Finding potential inhibitors against RNA-dependent RNA polymerase (RdRp) of bovine ephemeral fever virus (BEFV): an in-silico study

Abstract The bovine ephemeral fever virus (BEFV) is an enzootic agent that affects millions of bovines and causes major economic losses. Though the virus is seasonally reported with a very high morbidity rate (80–100%) from African, Australian, and Asiatic continents, it remains a neglected pathogen in many of its endemic areas, with no proper therapeutic drugs or vaccines presently available for treatment. The RNA-dependent RNA polymerase (RdRp) catalyzes the viral RNA synthesis and is an appropriate candidate for antiviral drug developments. We utilized integrated computational tools to build the 3D model of BEFV-RdRp and then predicted its probable active binding sites. The virtual screening and optimization against these active sites, using several small-molecule inhibitors from a different category of Life Chemical database and FDA-approved drugs from the ZINC database, was performed. We found nine molecules that have docking scores varying between −6.84 to −10.43 kcal/mol. Furthermore, these complexes were analyzed for their conformational dynamics and thermodynamic stability using molecular dynamics simulations in conjunction with the molecular mechanics generalized Born surface area (MM-GBSA) scheme. The binding free energy calculations depict that the electrostatic interactions play a dominant role in the RdRp-inhibitor binding. The hot spot residues, such as Arg565, Asp631, Glu633, Asp740, and Glu707, were found to control the RdRp-inhibitor interaction. The ADMET analysis strongly suggests favorable pharmacokinetics of these compounds that may prove useful for treating the BEFV ailment. Overall, we anticipate that these findings would help explore and develop a wide range of anti-BEFV therapy. Communicated by Ramaswamy H. Sarma


Introduction
Bovine ephemeral fever virus (BEFV) is an evolving arbovirus and is the causative agent of bovine ephemeral fever disease. It spans endemically into Africa, Australia, and Asiatic continents across the globe (Walker, 2005). A seasonal arbovirus causes outbreaks during the monsoon onset that accompanies high humidity, conducive for vector (includes Culicoids and mosquitoes) propagation, and facilitates transmission and spread (Walker & Klement, 2015). The clinical symptoms include polyphasic fever, synovitis, orthostatic hypotension, muscle stiffness, lameness, respiratory problems, abortion, and predisposition towards secondary bacterial infections (Barigye et al., 2016). The BEFV causes severe economic loss, particularly (MacFarlane, 1955) high yielders. The losses projects into reduced milk yield (by 10-15%) (Davis et al., 1984), poor reproductive performance (St George, 1986), body weight loss in beef (Lee, 2019), and associated trade restrictions (He et al., 2016). With a morbidity rate approaching 80-100% and a mortality rate of 2-29% (Lee, 2019), the livestock industries have long sought effective preventive and mitigation strategies. In contrast, no effective antivirals or vaccines are available yet.
The BEFV belongs to Rhabdoviridae and is placed under the genus Ephemerovirus (Murphy et al., 1972;Wu & Lehane, 1999). The viral genome is a non-segmented, negative single-stranded RNA (-ssRNA) of $14.9 KB and encodes ten different structural and non-structural proteins. Among the structural proteins, nucleocapsid (N), phosphoprotein (P), glycoprotein (G), and RNA-dependent RNA polymerase (L) are indispensable for viral genome transcription and replication (Walker et al., 1991). Among all the above proteins, L protein is the largest and most conserved viral protein, and its RdRp domain possesses RNA catalyzing-enzymatic activity required for RNA synthesis (Cheng et al., 2011;Dhillon et al., 2000). Proof-of-concept studies suggested the RdRp as a validated antiviral target for numerous RNA viruses such as West Nile virus (Malet et al., 2007), ebolavirus (Mirza et al., 2019), hepatitis C (Varshney et al., 2012), Zika virus (Singh & Jana, 2017), and coronaviruses (Elfiky, 2020), etc.
The RdRp domain adopts a canonical right-handed core structure as the 'fingers-palm-thumb' sub-domain (Liang et al., 2015) that lie in close spatial proximity to each other. The highly conserved core encompasses the catalytic active binding site within the protein. The BEFV-RdRp shares a strikingly similar structure with highly conserved motifs within the core structure as is seen with other RNA viruses such as vesicular stomatitis virus (VSV) (Liang et al., 2015), bovine viral diarrhea virus (BVDV) (Tiong-Yip et al., 2014), respiratory syncytial virus (RSV) (Fearns & Deval, 2016), etc. In all the above viruses, the central palm region retains the most typical conserved motifs, namely pre-motif A, motif A, B, and C. In contrast, the last motif, D, lies within the thumb region (Dhillon et al., 2000). At the junction of the three subdomains, a deep cavity is formed that incorporates RdRp active site with GDNQ residues (present within the motif C) that is certainly involved in RNA polymerization, primer-template entry-exit channel, incoming NTPs, in addition to forming a large priming loop for regulating RNA transcription and elongation (Peersen, 2017). Thus, the large cavity in RdRp has been a significant target for developing antiviral drugs (Duan et al., 2017;Pattnaik et al., 2018).
Despite annual outbreak reports for more than a century, BEFV remains a poorly understood virus with no effective treatment options. The high morbidity rate with increased fatality rates among the susceptible bovines enhanced geographical expansions, and enormous impact on the dairy sector's economy is alarming. Though a few studies have been conducted in some endemic areas to develop vaccines (Walker & Klement, 2015), not much effort has been made to investigate anti-viral BEFV-targeted interventions (Tseng et al., 2020). Symptomatic-based medication involving analgesic drugs, non-steroidal anti-inflammatory drugs (like ketoprofen, flunixin meglumine, and phenylbutazone), certain corticosteroids, and herbal medicines are administered to the diseased animal (DiK et al., 2019;St George, 1986). Often, calcium borogluconate injections are administered to treat hypocalcemia signs (e.g. ruminal stasis, paresis, loss of reflexes). Therefore, there is a deemed necessity to discover some new compounds for the treatment and control of BEFV in all of its endemic areas. Hence, to the best of our knowledge, this is the first attempt to develop BEFV-RdRp targeted antiviral drugs with minimal side effects to the host.
In the present study, an extensive in-silico based antiviral drug discovery approach was performed to investigate the potential RdRps inhibitors against BEFV. Due to the unavailability of the crystal structure of BEFV-RdRp, we predicted the model by using the comprehensive computational approach that included structure prediction, homology modeling, and sequence database searches. Further, the determination of its catalytic domains and active site (GDNQ) based on the consensus sequences and template superimposition structural analysis was performed. Consequently, we performed the structure-based virtual screening against the predicted catalytic residues of the core domain to explore the potential drug candidates for anti-BEFV therapy, followed by druggability and toxicity analysis. Additionally, molecular dynamics simulation and binding free energy calculations provided a cost-effective and time-efficient strategy for ascertaining potential hits. Thus, our study represents the first attempt to build the 3D structural model of BRFV-RdRp, determine the binding site, and search for the antiviral drugs that could target RdRp of BEFV. It provides an effective therapeutic strategy against this overlooked virus to resolve one of the emerging hindrances in the nation's economy.

Materials and methods
The methods involved in the study include homology modeling of protein, prediction of active sites, and virtual screening of libraries of compounds from the ZINC and Life Chemical (LC) databases. The schematic representation of the workflow for selecting lead molecules against BEFV-RdRp is shown in Figure 1.

Protein homology modeling
In the absence of any crystal structure of BEFV-RdRp, structure elucidation was done by homology modeling. The complete sequence of BEFV-RdRp (protein ID: QOU09210) was retrieved from the NCBI database, and the sequence from 48 to 887 amino acids containing the catalytic domain was modeled. The BLASTp (Altschul et al., 1990) search was performed against the PDB database (Berman et al., 2000) to identify its structural homologs. BEFV-RdRp model was selected using a combinatorial approach of motif search, fold recognition, and structure-based sequence alignment coordinated by the protein structure prediction by various built-in algorithms such as I-TASSER (Iterative Threading Assembly Refinement) (Zhang, 2008) and SWISS-MODEL (Guex & Peitsch, 1997). Each software predicted the same template of VSV-RdRp against BEFV-RdRp, suggesting that the best available template was chosen for the generation of the 3D protein model of the target protein. Superimposition of the BEFV-RdRp model and template was done for structural comparison using TM-align (Zhang & Skolnick, 2005) and FATCAT (Ye & Godzik, 2004) online web tools. TM-align evaluates the optimal structural similarity based on TM-score (score > 0.5 indicates two structures likely to have the same fold in CATH/SCOP). FATCAT allows flexible protein structure comparison between two structures by evaluating the p-value (pvalue < 0.005 indicates significant similar pairs). Overall, the model was built based on various attributes, such as similar fold, comparable secondary structure, RMSD between the model and template, sequence identity, and C a -RMSD of backbone atoms.

Model refinement and validation
The best-fitted model was evaluated in verify3D (Eisenberg et al., 1997) and PROCHECK (Laskowski et al., 1993) in the RAMPAGE server based on the stereochemical properties and its geometry. The model's quality was evaluated by calculating the dihedral angles (/, w) of each residue in the energetically allowed and disallowed regions by constructing the Ramachandran plot. Finally, the energy minimization was conducted for optimizing the model using the 3D-refine webserver (Bhattacharya et al., 2016). It provides the initial and refined structures and the statistics related to the quality scores of refined models, similarities, and deviations of the refined structures with respect to the input structure, RMSD scores, and the Molprobity scores. Lastly, the structural characterization of the predicted model was accomplished utilizing the ProSA-webserver (Wiederstein & Sippl, 2007) that is configured to find native structure consistency. Finally, the 3D-refine minimized structure was used to study the virtual screening and MD simulations.

Prediction of binding site
Identification and characterization of active sites of the target protein are of prime importance in drug discovery. Identifying active sites could be possible by using various analytical tools based on the available known structures of proteins. Initial identification was made by performing the multiple sequence alignment (MSA) of closely related members in the Clustal Omega, searching for their conservation pattern among the RdRp core sequence. Further, analytical tools such as COACH-metaserver (Yang et al., 2013) and 3DLigSite (Wass et al., 2010) were utilized to predict the binding sites after superimposing the closest template (PDB ID: 5A22) to the generated BEFV-RdRp model.

Ligand preparation
In the present study, we used 43,174 compounds for the computational screening from the ZINC15 and Life Chemical databases against the target RdRp of BEFV. From the ZINC15 database, FDA-approved 1615 drugs were used. In the case of the Life Chemical database (LC), four different categories were used, namely, LC_antiviral (10,158), LC_Polymerase (18,279 ligands), LC_Building blocks (12,209 ligands), and LC_Natural compounds (913 ligands). Broadly we have categorized all-LC ligands as Antiviral (AV) and ligands from the ZINC database as FDA in the subsequent analysis.
Ligands were imported in the Schrodinger Glide (Friesner et al., 2006) and subjected to the Ligprep module  to generate all the ligands conformer along with their tautomeric combination under pH 7.0 6 2.0. The ligands were minimized and optimized using the OPLS3 (Harder et al., 2016) force field after the addition of hydrogen atoms.

Protein preparation and receptor grid generation
Next, the optimization and refinement of the BEFV-RdRp model were performed using the standard protocol in the protein preparation wizard Sastry et al., 2013;Velazquez-Campoy et al., 2002) of Schr€ odinger software, with the help of Maestro platform. The structure was refined and optimized at pH 7.0 using PROPKA (Olsson et al., 2011). Finally, under the default settings and the presence of the OPLS3 (Harder et al., 2016) force field, the target protein was minimized and ready for docking. The docking grid was prepared using Glide around the predicted binding site residues, preserving the conserved structural motifs (centering G 739 D 740 N 741 Q 742 residues) to reduce the search space ligand optimization keeping a 12 Å cubic space around it.

Virtual screening
Virtual screening (VS) evaluates different ligand categories searching for potential compounds that could be used to combat a drug target. The sequential molecular docking of the selected ligands from the ZINC: 1,615 ligands and Life Chemical databases: 41,559 ligands from its different categories were done against BEFV-RdRp. The molecular docking and screening of potential lead molecules were done utilizing the virtual screening workflow under the Glide module (Friesner et al., , 2006Halgren et al., 2004;Roux et al., 1996) of the Schr€ odinger suite. During VS, the filtering criteria of QikProp (Briggs et al., 1996) and the Lipinski's rule of five were taken into consideration. The VS accompanies sequential high throughput virtual screening (HTVS), standard precession (SP), and extra precession (XP) to achieve a potential set of lead molecules with high accuracy. A flexible docking method was chosen while performing the molecular docking in the virtual screening workflow. For each library, the top 50% of docked complexes from HTVS were used in the SP docking method, followed by the top 20% of SP docking to the XP docking procedure, keeping the default parameters for the rest. Finally, based on the Glide score, negative binding free energy and the interacting residues involved at the active site were retained after XP docking.

Pharmacokinetic and toxicological profiling
Computational validation of the virtually predicted drug candidates is extensively performed to generate potential lead substances with efficient drug likeliness and pharmacokinetics. To accomplish this, the initial screening of filtered compounds was done by utilizing the Lipinski's rule of 5 (Lipinski, 2004) and QikProp (Briggs et al., 1996) module of the Schrodinger suite with set thresholds during the virtual screening. Additionally, for the selected six molecules from the LC database in-silco tools such as ADMET-SAR (Yang et al., 2019) and pkCSM (Pires et al., 2015) were utilized to predict the essential physiochemical, pharmacodynamics, and toxicity properties of the screened compounds. The various properties include CNS activity, blood-brain penetration, Madin-darby Canine Kidney permeability, %human oral absorption, and certain undesired effects, including mutagenicity, tumorigenicity, and toxicity.

Molecular dynamics simulations
Molecular dynamics (MD) simulation is the most widely applied method for studying the binding of protein-ligand at the atomic level. In this study, we used MD simulation to consider the dynamic behavior of our top hit compounds. The MD simulations of all systems were performed by using the pmemd.cuda module of the AMBER18 (Case, 2018). Post MD analysis of the trajectories was performed using the Cpptraj (Roe & Cheatham 2013) module. AMBER ff14SB forcefield (Maier et al., 2015) and the generalized Amber force field (GAFF2) (Wang et al., 2004) were used in the simulations for assigning parameters to proteins and small molecules, respectively. The Leap module (Salomon-Ferrer et al., 2013) of AMBER was used to add all missing hydrogen atoms. Counterions were appropriately added to neutralize the whole system. The solvation of complexes was performed in an octahedron periodic box with the TIP3P water model (Jorgensen et al., 1983). The buffer distance was kept around 10 Å along each side. The Langevin thermostat (Pastor et al., 1988) and the Berendsen Barostat (Berendsen et al., 1984) was used for controlling the temperature at 300 K and pressure at 1 bar, respectively. The SHAKE algorithm (Kr€ autler et al., 2001) was used to constrain the highly fluctuating hydrogen bonds. Two femtoseconds of time-step were used in the simulations. The long-range electrostatic interactions were considered by using the Particle mesh Ewald (PME) (Darden et al., 1993) scheme, and the nonbonded cut off was taken as 10 Å. Lennard-Jones potential was used to treat the van der Waals interactions.
Two steps of minimization were performed using the steepest and conjugate algorithm for 500 steps each. Initial restraint minimization was done with a harmonic force constant of 2.0 kcal mol À1 Å À2 . Followed the minimization, systems were gradually heated from 0 K to 300 K in the NVT ensemble. Subsequently, 50 ps MD simulations with a restraint force constant of 2.0 kcal mol À1 Å À2 were performed at a constant pressure of 1 atm, and 300 K. Complexes were then equilibrated for 1 ns. Finally, production simulations of 100 ns were performed at the NPT ensemble for each receptor-inhibitor complex. The cartesian coordinates were saved for every 10 ps, resulting in 10,000 configurations for each simulation.

Binding free energy using the MM/GBSA scheme
The binding affinity of the different lead molecules with the BEFV-RdRp was carried out using the molecular mechanics Generalized-Born surface area (MM/GBSA) (Jonniya et al., 2021a(Jonniya et al., , 2021bKar et al., 2007Kar et al., , 2011Kar & Knecht, 2012a, 2012cKollman et al., 2000;Thurakkal et al., 2020;. We used the single trajectory approach to calculate binding free energy (DG bind ), which involved the simulations of only the proteininhibitor complex. The details of the MM/GBSA scheme was previously discussed in our previous papers for various protein-ligand complexes (Jonniya et al., 2019;Roy et al., 2020Roy et al., , 2020aSk et al., 2020aSk et al., , 2020bSk et al., , 2020cSk et al., , 2021.The same protocol has been followed here, and the binding free energy is calculated using the following Equation: DG bind is the summation of the enthalpy and entropy terms, which further represents the summation of the internal energy (DE internal ), solvation free energy (DG solv ), and the configurational entropy (-TDS). Furthermore, DE internal represents the summation of the DE covalent (bond, dihedral, and angle energies), electrostatic energy (DE elec ), and the van der Waals interaction energy (DE vdW ). In contrast, the solvation free energy (DG solv ) includes the polar solvation energy (DG pol ) and the non-polar solvation energy (DG np ). The polar part of the solvation energy DG pol was calculated using the Generalized-Born model (GB) (Gohlke et al., 2003), while the non-polar part was computed using the SASA. Overall, 6000 structural frames from the last 60 ns were used for the calculation of the MM/GBSA. In addition, the per-residue of decomposition of energy was also calculated to reveal the contributions from each residue by the molecular mechanics generalized-Born surface area (MM-GBSA) scheme (Gohlke et al., 2003).

Protein modeling and refinement
The BEFV-L protein comprises of RdRp domain and the capping domain, as shown in Figure 2(A). Herein, by considering it as the potential target, the BEFV-RdRp domain was modeled. The PDB-BLAST analysis of the RdRp protein sequence revealed a 47.95% sequence identity from the VSV-RdRp (PDB id: 5A22) with 99% query coverage. Due to low % identity, further accuracy of the predicted 3D structure was analyzed using different servers. The structure of the BEFV-RdRp was modeled in the I-TASEER and SWISS-MODEL. The most favorable template reported for each case was from the VSV-RdRp (PDB: 5A22). The quality of obtained structures from different servers was checked and assessed in their Ramachandran plots, as shown in Supplementary Information Figure S1(A). The best model was selected from the SWISS-MODEL. Further, refinement was done in the 3DRefine web server, followed by the energy minimization, Figure S1(B).
Additionally, the quality of the model was checked based on the MolProbity score, which was obtained as 2.2; lower the score indicates a more physically realistic model, and Z score which was obtained as À10.78 as given by the ProSA server and ProSa designated energy plot which revealed a high-quality structure Figure S1(C) and S1(D). The structural alignment by superimposing the energy minimized refined structure of BEFV-RdRp onto the VSV-RdRp revealed an RMSD of 0.32 Å, signifying their comparable structural similarity shown in Figure 2(B). The residues depicted in the binding pocket sites of BEFV-RdRp are the functionally important conserved residues, shown in Figure 2(C). The core structure comprises three characteristic sub-domains as a finger (343-551 and 594-664 aa), palm (551-594 and 664-772 aa), and thumb (773-845 aa) as per the sequence similarity of VSV (Liang et al., 2015). The conserved catalytic active site comprising 739-GDNQ-742 resides positioned within the palm domain (see Figure 2D). The sequence alignment of BEFV-RdRp with other closely related members of the genus Ephemerovirus and Vesiculovirus revealed a high level of sequence conservation as reflected by its conserved motifs the structural and sequence levels, as shown in Figure 3(A).
All the conserved motifs, typically A to D motifs, the active site GDNQ motif (lying in the motif C of the palm domain), and their localization within the modeled protein are shown in Figure 3(B).

Virtual screening
The different compounds, 41,559 from the libraries of Life Chemical and 1615 FDA drugs, were docked against the BEFV-RdRp protein utilizing the virtual screening (VS) workflow in the Schrodinger suite. During the VS, filtering criteria from QikProp and the Lipinski's rule of five were taken into consideration. The obtained screened molecules with a cutoff docking score of À6.0 kcal/mol were retained after XP docking and shown in Table S1. The docking score from the antiviral group (AV) varies from À6.0 to À10.43 kcal/mol. However, in the FDA group, it ranges from À6.47 to À8.58 kcal/mol.
The top lead molecules based on the Glide score from each category were finally chosen for further analysis, which comprises six top lead molecules from the antiviral (AV) category; AV1-AV6, in which AV1 was from LC_Antiviral, AV3 from LC_Polymerase, and remaining AV2, AV4, AV5, and AV6 were from the LC_building blocks and the top three lead molecules from the FDA category. The docked complex of each selected ligand inside the binding pocket is given in Figure 4. Therefore, a total of nine molecules were selected for subsequent analysis, as shown in Table 1. The histogram drawn based on the docking score for all selected compounds is shown in Figure S2 for an easy interpretation at a glance.

ADMET properties
Selected ligands were subjected to in-silico drug-likeness prediction and further ADMET analysis for studying their pharmacokinetics properties and listed in Table S2. The physiochemical parameters were evaluated based on the Lipinski's rule of five, which in all compounds were strictly followed with no violation (except AV4 that showed one violation). Also, the molecular weight that for all compounds showed <500 g/mol, anticipating rapid transport, diffusion, and absorption. The topological polar surface area (TPSA) analysis showed values range from 63.32 to 167.10 Å for all compounds, as a lower TPSA value may provide better CNS penetration (Blake, 2000). The water solubility of the compounds showed all to be soluble (except AV1, with poor solubility). Furthermore, the predicted key absorption parameters such as human colon adenocarcinoma-2 cell line (CaCo2) permeability and human intestinal absorption (HIA), demonstrating high permeability of AV1, AV2, AV5 relative to AV3, AV4, and AV6 and could be easily absorbed via the intestine (except for the AV3).
A balance between hydrophobicity and hydrophilicity for good absorption and distribution within the body is another requirement (Duffy & Jorgensen, 2000). The values of log Po/ w and logS provide the estimated octanol/water partition coefficient and the drug's aqueous phase solubility, respectively. All the active ligands exhibit within the recommended range of logPo/w as 5 and logS as À6.5-0.5 implying that they are neither too hydrophilic nor too hydrophobic (Table S2).
The drug distribution is regulated by many parameters such as concentration in plasma and its binding ability to plasma proteins, etc. The volume of distribution at steady-state (VDss) indicates that the expected theoretical dose for uniform distribution in the plasma was low for all six compounds. Additionally, as determined by the fraction in the unbound state, the degree of diffusion through the plasma membrane increases in the order AV1 < AV2 < AV4 < AV3 < AV5 < AV6.
Furthermore, ligands' interaction with cytochromes P450 (CYP) has a key role in drug metabolism. This superfamily consists of five major isoforms (CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4) (Di, 2014). Each isozyme is a key player in the drug elimination through metabolic biotransformation as 50 to 90% of these predicted drugs may act as a substrate of these isoenzymes (Blake, 2000). Hence, except AV1, all compounds did not inhibit any of the CYP450 isoforms, indicating the expected drug metabolism capability.
The blood-brain barrier (BBB) permeation was measured using the BOILED-Egg model (Daina & Zoete, 2016), which showed few BBB permeate in the studied compounds (such as AV3 and AV5). The CNS permeability highlights very minute penetration in the brain. However, the medium degree of lipophilicity of the drugs indicated no detrimental effect on the exposure to the nervous system. For a predicted ligand to attain potent CNS activity, permeability glycoprotein (P-gp), an active efflux membrane protein, may play a pivotal role in the active efflux of the ligand from the brain back into the blood if it acts as a substrate for it.
However, such ligands must be co-administered with a P-gp inhibitor (Ghose et al., 2012). Among the six lead hits, only two compounds, AV1 and AV3, are predicted as P-gp substrates and require P-gp inhibitor co-administration for efficient drug activity. However, AV1 and AV2 themselves were  predicted to be P-gp I and II inhibitors, respectively. Thus, all six lead hits may regulate the physiological functions of Pglycoprotein.
On evaluating most of the toxicity dimension values, almost all compounds were predicted to overrule each of the toxicity modules. Although AV1 was projected as the substrate of the renal organic cation transporter-2 (Renal OCT2), all the remaining compounds may be cleared by other available pathways, like air, bile, and sweat. Additionally, none of the compounds were associated with skin sensitization. While most of the compounds did not show any hepatotoxicity, AV3 was proposed to cause the minor effect. Strikingly, no unwanted effects such as mutagenicity, tumorigenicity, and irritability were demonstrated by any of the compounds studied. Overall, the in-silico ADME report revealed that all compounds except for a few that showed minor violations, all others fulfill the desired properties with high bioavailability and likeliness to be the lead compound.

Structural stability and flexibility of the system simulations
To assess the stability of docked complexes and their interacting residues, we conducted the molecular dynamics simulations of these complexes (Chandra et al., 2021;Parvez et al., 2020;Ribaudo et al., 2020;. The top selected nine compounds: three FDA and six antiviral compounds were subjected to 100 ns of MD simulations to monitor the dynamic, structural, and energetic features, which is further compared with its apo (without ligand). The time evolution of root means square deviation (RMSD) of the backbone atoms relative to the respective energy minimized structures are shown in Figure 5(A). Average RMSD values for the last 50 ns simulations are reported in Table 2. It is evident from Figure 5(A), the RMSD values for all protein backbone atoms for each system initially increase up to 20 ns; after that, all nine systems are reached a well equilibrium  The data are reported as average ± standard error of the mean (SEM).
state throughout the entire last 80 ns production simulations. As seen from RMSD values compared to the rest of the antiviral compounds. These observations show that the MD simulation convergency is reliable and can be used for further dynamic analyses and energetic calculations. Furthermore, the time evolution of the RMSD values of the backbone atoms of residues within 5 Å around the ligand in the binding pocket relative to the respective initial conformation is shown in Figure 5(B). The last 50 ns trajectory average values are listed in Table 2. It is clear from Figure  5(B) that RMSD fluctuations are minimal throughout the entire 100 ns simulations for each system. AV4 shows the lowest average changes at the binding pocket region, and AV2 and AV3 bound complexes show higher RMSD values of more than 2 Å. The RMSD values of RdRp/FDA complexes range between 1.47 Å to 1.73 Å. This implies that the binding pocket conformation of AV4 and FDA3 is slightly different from the other inhibitors in some regions in the pocket. Furthermore, to illustrate receptor backbone C a atoms flexibility, the root means square fluctuations (RMSFs) in nine complexes were monitored ( Figure 5C). According to Figure  5(C), nine complexes have more or less similar fluctuation tendencies. The notable high dynamic fluctuations are observed at both C-terminal and N-terminals and the nonligand binding sites (around Leu151, Glu229, Leu498, Ala529, Ser752, and Phe851) and loops regions. Moreover, the lower fluctuations are observed at active sites, and the average RMSF of these residues is less than 1 Å, such as Leu397, Val555, Glu564, Ile631, Asn634, Asp636, Gln707, Gly710, Leu 714, and Asn741. This suggests that after binding the FDA and small antiviral molecules, the active site of RdRp still rigid and interacts strongly with each inhibitor.
Finally, we calculate the heavy atoms RMSD based 2D potential of the mean force (PMF) profile of each small molecules (both FDA and AV), see in Figure 5(D). It is evident from Figure 5(D) that FDA3 and AV4 have a single free minimum found at $0.5 and $0.8 Å. The movement of each small molecule is dissimilar inside the binding site, as seen in Figure 5(D). In AV1, AV3, and AV6, the free energy minimum is shifted to a higher RMSD value (>1 Å) compared to other small molecules. Similarly, FDA1, FDA2, AV2, and AV5 show high probable free energy minimum at RMSD values of <1 Å. These results agree with the binding pocket RMSD analysis. Moreover, we also monitored the center of mass (CoM) distances between ligand and RdRp active site residues, and the time evolution of these distances is shown in Figure S3. The CoM distance analysis clarifies that the distance between FDA-approved small molecules and binding pocket residues lies within $5 Å. On the other hand, in antiviral inhibitors, only AV1 and AV4 molecules CoM distance lie within $5 Å and the rest of the molecules (AV2, AV5, and AV6) within $10 Å. In AV3, we observed a different scenario; up to 25 ns distance is increasing, then it reached a stable equilibrium at around 15 Å for the rest 75 ns simulation. This analysis also suggests that FDA3 and AV4 are likely bound to RdRp more tightly than the rest of the inhibitors.
Besides, we computed the RMSD of the finger, palm, and thumb region of RdRp, and shown in Figure S4. Mainly the binding site residues lie in the palm and finger regions. As depicted in Figure S4, each complex's palm region exhibited more stability than the finger and thumb regions. For RdRp/ FDA complex, the trend of higher structural deviations was seen for the finger region, being highest in the RdRp/FDA1 compared to others. In contrast, all of their palm regions remain stable throughout the simulations. Among the RdRp/ AV complexes, again, the finger region showed high deviations in all except for RdRp-AV4, where the thumb region exhibited the most deviations. However, the palm region that contains most of the binding site residues exhibited the lowest variation in the RdRp/AV4 while the RdRp/AV3 complex showed the largest. Overall, it showed that the palm Standard errors of the mean are provided in the parenthesis. region's binding site region stabilizes upon binding of the ligands, especially in RdRp/FDA complexes and RdRp/ AV4 complex. We also predicted the radius of gyration (R g ) and solvent accessible surface area (SASA) of RdRp of all nine complexes. The temporal evolutions of these parameters are shown in Figure S5(A). The last 50 ns average values of these parameters are recorded in Table 2. According to Table 2, the average values of R g range from 29.37 to 30.61 Å, suggesting that all the complexes are relatively compact and not much changed after inhibitor binding. Similarly, the average SASA values vary between 39,689 Å 2 to 42,413 Å 2 (see Figure S5B).

Binding free energy
To give insight into the binding affinity of the best-screened molecules, three FDA drugs and six antiviral molecules (AV1-AV6) against the RdRp were analyzed using the MM-GBSA scheme. Among all, except AV2, AV3 and AV5 displayed a comparatively high binding affinity against RdRp, as shown in Table 3 and Figure S6. The different components of the binding free energy, such as electrostatic (DE elec ), van der Waals (DE vdW ), and solvation free energy (DG pol and DG np ), were estimated and depicted in Figure S6. It is evident from Figure S6 that DE elec , DE vdW, and DG np favor the binding of the ligand to RdRp of BEFV, while DG pol opposes the complex formation. Among various energetic components, the intermolecular electrostatic interaction energy contributes most favorably to the complexation of all cases except AV1 and AV3. In the case of the AV1-RdRp complex, DE elec and DE vdW are comparable; however, the AV3/RdRp complex exhibited an opposing behavior of DE elec as given by the positive energy; instead, the complex formation is driven by the van der Waals interactions (see Table 3).
It is evident from Table 3, the values of DE elec interaction vary from À128.67 to 166.72 kcal/mol in FDA drug complexes, while for AV complexes, it ranges between 0.73 and À111.39 kcal/mol. In contrast, DE vdW interaction varies from À19.42 to À49.16 kcal/mol in all the complexes. Moreover, the total energetically flavoring components (DE elec þ DE vdW ) for RdRp/FDA complexes were high compared to RdRp/AV complexes. Furthermore, it is to be noted that in the RdRp/ AV1, RdRp/AV2, and RdRp/AV3 complexes, the DE elec is overcompensated by the desolvation polar energy (DG pol ) as given by 3.1 to 22.92 kcal/mol. However, for complexes RdRp/AV4, RdRp/AV5, RdRp/AV6 and RdRp/FDA1, RdRp/ FDA2, RdRp/FDA3, the total electrostatic interaction energy (DE elec þ DG pol ) favours the complex formation which ranges from À1.63 to À17.28 kcal/mol and the RdRp/FDA3 complex exhibited the most significant total electrostatic interactions. Moreover, the total hydrophobic interactions (DE vdW þ DG np ) for RdRp/AV1, RdRp/AV3, and RdRp/AV4 were more energetically favorable than RdRp/FDA complexes. Among them, RdRp/AV1 exhibited the highest total hydrophobic interactions, even better than FDA drugs. However, it is overcompensated by the total electrostatic interactions (DE elec þ DG pol ¼ 22.92 kcal/mol), due to which it binds less favorable than FDA drugs against RdRp. Similarly, the RdRp/AV3 complex exhibited favorable total hydrophobic interactions (DE vdW þ DG np ), being À39.0 kcal/mol compared to FDA drugs. However, it also gets overcompensated by the total electrostatic interactions (DE elec þ DG pol ¼ 11.10 kcal/mol) that disfavor the complex formation.
The estimated binding free energies (DG) for RdRp/FDA were in the order of FDA3 (À43.55 kcal/mol) > FDA2 (À37.53 kcal/mol) > FDA1 (À34.29 kcal/mol), and for the complexes RdRp/AV it follows: AV4 (À40.11 kcal/mol) > AV1 (À31.97 kcal/mol) > AV6 (À29.39 kcal/mol) > AV3 (À27.8 kcal/mol) > AV5 (À24.65 kcal/mol) > AV2 (À19.53 kcal/mol). Overall, it suggests that the screened FDA drugs bind strongly to the BEFV-RdRp compared to the AV molecules. Among them, FDA3 could be a potent inhibitor to BEFV/RdRp as both the DE elec and DE vdW contribute more favorably than other complexes. The second most promising molecule which binds to RdRp of BEFV was AV4. Its total molecular interactions (DE ele þ DE vdW ¼ À124.88 kcal/mol) favors the binding, including the total electrostatic interaction and total hydrophobic interactions among the other AV complexes. Hence, the screened FDA3/FDA2/FDA1 and AV4 could be exploited as potent inhibitors against the BEFV RdRp. Our study correlated well with the disease ailments and the corresponding resultant drugs screened after the virtual screening. The FDA drugs showed relevance against the associated diseased condition, thus signifying the usefulness of these drugs that may prove potent to treat the disease by targetting the BEFV-RdRP.

Per-residue energy contributions to binding
Furthermore, quantitively each residual contribution in the binding of ligands to RdRp of BEFV was analyzed by performing ligand-residue pair decomposition using the MMGBSA scheme. It reveals the hot spot residues from the binding pocket, which imparts knowledge regarding the mutational study of BEFV-RdRp while designing any inhibitors. Different components of each residual pair with FDA drugs and AV molecules were listed in Table S3, including the residues' backbone and side-chain contributions. Figure  6 highlighted those residues having more than À1.0 kcal/mol of energy towards the binding contributions.
Further, the comparative analysis of the highest contributing residues from the BEFV-RdRp complexes was depicted in Figure 7. It suggests that key contributing residues such as Glu633, Glu563, Glu707, and Asn636 from RdRp/FDA1 and RdRp/AV4 favor the binding to its high electrostatic contributions. Also, the residue Asp740 contributed its high electrostatic contributions as shown in RdRp/FDA2, RdRp/FDA1, and RdRp/AV4 complexes reveal that in the conserved GDNQ motif, Asp740 play a significant role in the binding of BEFV-RdRp.
In addition, we complement the above findings by providing information on H-bond and hydrophobic residues from the final conformation of the production simulation by plotting a 2D Ligplot diagram, as shown in Figure 8. The hydrogen bonds are shown with green dotted lines, and the red spike represents the hydrophobic residues. The RdRp/FDA3 and RdRp/AV4 complex showed the highest contributions of hydrophobic residues, including the most stable H-bonds of Asn636, Glu633, Asn637, Glu707, and Trp635. The significant hydrophobic interacting residues in RdRp/FDA3 comprise Glu563,Leu398,Val402,Leu553,Gly705,Leu706,Phe566,and in RdRp/AV4 includes Leu398,Phe566,Lys634,Arg565,and Asp740. Although the RdRp/AV1 complex comprises the most hydrophobic interactions, which provides for Leu398, Lys634, Glu707, Glu633, Lys804, Gln738, Gly803, Asp740, Arg565, and Gly564, however, it showed only the Asp631 Hbond interactions, suggesting that the hydrophobic rather than electrostatic interactions govern the primary interaction compared to other complexes. Overall, it indicates that the RdRp/FDA3 complex has a higher binding affinity than other complexes due to stable H-bond and hydrophobic interactions.

Discussion
BEFV is a neglected livestock pathogen that has devastating economic consequences for the dairy sector. Its enhanced pathogenicity and transmission possess an alarming situation to the livestock industry worldwide (Yanase et al., 2020). The threat of the currently emerging BEFV necessitates safe and effective antivirals therapeutic strategies, in outbreak situations, in priority. Presently, none of the FDA-approved medications is prescribed to treat or prevent BEFV. Neither other compound has yet been reported explicitly for the BEFV treatment. Thus, it is critical to explore compounds that are efficient and cost-effective.
Computational approaches provide a rapid and scientifically valid framework for establishing highly selective inhibitors against critical viral proteins and assisting in developing antiviral drugs (Maggiora & Shanmugasundaram, 2011). Across the diverse RNA viruses, the RdRp, with its high degree of structural conservation, is one of the critical proteins in the BEFV life cycle and has been suggested as a promising antiviral target for drug development (Tchesnokov et al., 2019;Tian et al., 2021). Current BEFV treatment is generalized, including analgesics and antiinflammatory drugs as the first choice of action, often followed by corticosteroids in severe cases, leading to the immunosuppression that may increase the chances of undesirable secondary infections.
In the absence of any BEFV-RdRp crystal structure, the present study employed a comprehensive computational approach to predict the model of BEFV-RdRp and followed by its binding site prediction. The VSV-RdRp (PDB: 5A22) served as the most suitable template to generate a 3D model of BEFV-RdRp for homology modeling, further refined and minimized. In addition, the RdRp sequence and structural alignment showed that motifs (A-D) with druggable binding sites and an active site with GDNQ residues within the catalytic center are highly conserved among Rhabdoviridae members. Additionally, computer-aided drug design may assist in finding new therapeutics due to its rapid and accurate screening capability from large libraries of compounds. Due to the lack of any previous drug discovery reports in BEFV, there were no compounds to take is as reference compounds. Therefore, we adopted the virtual screening workflow to identify potential inhibitors against BEFV-RdRp from the diverse, extensive compound library. In total, we took 1,615 FDA-approved drugs from the ZINC database and 141,559 compounds from the different categories of Life Chemical (LC) databases, such as LC_antiviral (10,158), LC_Polymerase (18,279 ligands), LC_Building blocks (12,209 ligands), and LC_Natural compounds (913 ligands). In the screening module of Schrodinger, we considered the combinatorial approaches such as HTVS, SP, and XP along with the initial filtering criteria from the QikProp module and the Lipinski's rule of five. Based on the initial filters and XPdocked score of binding energies, we selected nine compounds; three from the ZINC database, named as FDA-FDA3, and six from the LC database, named as AV1-AV6. Further, extensive pharmacological and bioavailability analyses of those selected compounds were performed. Such in-silico drug-likeness prediction and further ADMET profiling present various opportunities that help accelerate the new target discovery and predict their bioavailability as potent drugs. Moreover, to understand and endorse the stability of the docked complex, molecular dynamic simulations in conjunction with binding free energy calculations were conducted (Panda et al., 2020;Wang, 2013).
Our study presented the first effort on the structural modeling of BEFV-RdRp. Further, molecular docking and MD simulation were employed to anticipate potent anti-BEFV drug candidates targeting this enzyme. Moreover, the search for typical drug surface hotspots was investigated to assess their affinity and interaction pattern.
The structural parameters were assessed by molecular stability analysis of the protein-ligand complexes. The RMSD profile of the protein specifies the overall integrity of each complex during 100 ns of MD simulations. Also, analysis of RMSD within 5 Å around the ligand in the binding pocket relative to the respective initial conformation indicates minimal RMSD fluctuations, ranges between 1-2 Å. The lowest value was obtained for AV4 and highest in cases of AV2 and AV3, while all FDA complexes ranged between this range. Further, flexibility analysis (RMSF) revealed lower fluctuations at the active site upon ligand binding. Subsequent explorations were done by MM-GBSA to evaluate the binding free energy and affinities between the protein and ligands (Kumar et al., 2020;Tortorici et al., 2019). It contributed to different components of binding free energy. It evidenced the electrostatic interactions (DE elec ) as the highest contributors in most of the lead compounds, indicating their spatial stability at the active site of RdRp, except for AV3 having positive DE elec energy. The most favorable molecular activity was shown by FDA3 (DG bind ¼ À43.55 kcal/mol) and AV4 (DG bind ¼ À40.11), proposing these as the most potent candidates against RdRp in BEFV therapy.
Interestingly, the drug surface hotspot analysis revealed that almost all compounds have a typical pattern of molecular binding sites, implying these as RdRp hotspots that prefer direct interaction with inhibitors. The top contributing residues include Arg565, Asp631, Glu633, Asp740, and Glu707, although Lys634, Asn636, Glu792, and Tyr632 residues also depicted significant contribution. Also, the aspartate residues (Asp740) within the conserved active site and another (Asp636) were divulged as highly interacted residue in almost all the compounds. Furthermore, H-bond findings indicated that the H-bonds were stable during the simulation and would most likely hold a significant role. However, hydrophobic residues were also significantly complemented in the complex stabilization.
Lastly, the biological properties of all screened FDA compounds also sufficed their established medical implications as drugs. For example, Droxidopa (FDA1) used to treat neurogenic orthostatic hypotension and extreme fatigue (Kaufmann, 2017). Epinephrin (FDA2) works by relaxing the muscles in the airways and tightening the blood vessels during respiratory and heart-related issues and anaphylaxis (Dretchen et al., 2020), while L-norepinephrine (FDA3) is used in vasodilatory shock states, including septic shock and neurogenic shock (Beloeil et al., 2005). Since all diseased conditions are evident during BEF disease manifestations, it corresponds well with possible BEFV therapeutics.
Overall, the combinatorial docking and molecular dynamic simulations suggest three FDA drugs (FDA3, FDA2, and FDA1) and one Antiviral compound (AV4) as lead compounds against BEFV-RdRp. Hence, it revealed that FDA drugs exhibited better efficacy than the AV category against BEFV-RdRp. Our study correlated well with the disease ailments and the corresponding resultant drugs.

Conclusion
The present study marked the first effort to identify possible inhibitors of BEFV-RdRp for targeted antiviral treatment by utilizing an integrated computational approach. It demonstrated the RdRp domain as an effective anti-BEFV therapeutic agent to inhibit BEFV propagation within the host with these biologically active compounds. The in-depth structural elucidation of interacting residues and the dynamic conformational exploration of the identified compounds offered a way to design selective inhibitors against BEFV. All the inhibitors, majorly FDA3, AV4, FDA2, and FDA1, occupied the binding site that consists of critical catalytic residues concluding their inhibitory potential for RdRp enzymatic activity. However, this in-silico study can be easily coupled with in vitro and in vivo experiments to confirm their effectiveness against BEFV.

Disclosure statement
The authors declare no conflict of interest.

Data availability statement
The data that support the findings of this study are available from the corresponding author (PK) upon reasonable request.

Authors contributions
PK and DN conceived and supervised the project. SP and NAJ conducted molecular dynamic simulations, performed data analysis, and wrote the manuscript. MFS performed data analysis and took part in manuscript writing. PK and DN edited the manuscript. All authors have approved the final version of the manuscript.

Ethics statement
There is no human or animal experiment in this study.