Discovery of Small Molecules from Echinacea angustifolia Targeting RNA-Dependent RNA Polymerase of Japanese Encephalitis Virus

The Japanese encephalitis virus (JEV), a mosquito-borne flavivirus that causes viral encephalitis leading to neural damage, is a major threat in most Asian countries. The RNA-dependent RNA polymerase (RdRp) present in the viral genome is the key component for genome replication, making it an attractive target for antiviral drug development. In this study, the natural products from Echinacea angustifolia were retrieved for structure-based virtual screening against JEV–RdRp. The top six compounds (Echinacoside, Echinacin, Rutin, Cynaroside, Quercetagetin 7-glucoside, and Kaempferol-3-glucoside) were obtained based on the highest negative docking score, ADMET (absorption, distribution, metabolism, excretion, and toxicity), and molecular interaction. The computational analysis of these selected compounds against the co-crystallized ligands, i.e., ATP and GTP, were performed. Further, 100 ns molecular dynamic simulation and post-free binding energy calculation of all the selected compounds complexed with JEV–RdRP were performed to check the stability of the complexes. The obtained results showed considerable stability and intermolecular interaction with native ligand-binding site residues of JEV–RdRp. Hence, selected natural compounds are admissible inhibitors of JEV–RdRp protein and can be considered for future antiviral drug development studies.


Introduction
Japanese encephalitis (JE) is the most commonly diagnosed epidemic encephalitis-a vector-borne acute central nervous system infection, caused by the Japanese encephalitis virus (JEV) that belongs to the genus Flavivirus of the family Flaviviridae and comprises five genotypes (GI-GV) [1][2][3]. JEV is transmitted to humans by the bite of infected mosquitoes (Culex spp.) and causes severe neurological manifestations [4,5]. JEV is an endemic in Life 2022, 12, 952 3 of 17 database (https://www.rcsb.org/) [54]. Additionally, a total of 53 natural compounds, reported in Echinacea angustifolia, were collected from the literature, and their 3D atomic coordinates were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih. gov/) [55] to form the ligand library for further computational analysis.

Virtual Screening, Redocking Analysis, and ADME
For the structure-based virtual screening, the JEV-RdRp domain was processed and minimized using the Dock Prep tool in UCSF Chimera [56], with default parameters. Here, the receptor molecule was prepared by removing the native ligand, water molecules, heteroatoms, and non-polar hydrogen atoms, along with the addition of Gasteiger charges and polar hydrogen atoms, followed by structure minimization, using 100 and 50 steps of steepest descent and conjugate gradient methods, respectively. Furthermore, for the ligand library preparation, i.e., the collected natural compounds of Echinacea angustifolia and structure-based virtual screening, PyRx0.8 (Virtual Screening Tool) [57] was used, under default parameters. Herein, each ligand was initially added for Gasteiger partial charge, and energy was minimized under MMFF94 force field using default steepest descent and conjugate gradient with 1000 steps, followed by the addition of polar hydrogens, as reported earlier [58]. The prepared ligands were docked in the binding pocket of JEV-RdRp covering essential residues (Lys 463 , Lys 471 , Arg 474 , Asp 541 , Ser 604 , Asp 668 , Arg 734 , Arg 742 , Ser 799 , Try 800 , and Ser 801 ) in the docking grid 39.3513 Å × 29.2100 Å × 28.2521 Å and center at −37.5873 Å, −5.1747 Å, −35.1844 Å along the X, Y, and Z axes.
Subsequently, the top 10 ranked ligands based on binding scores were selected for redocking analysis by comparison to the ATP and GTP as reference molecules using the Chimera-AutoDock Vina Plugin setup [56,59]. Herein, the processed protein structure was docked with the selected compounds in the same docking pocket as considered for virtual screening (center coordinates: −37.5873 Å, −5.1747 Å, −35.1844 Å; grid size: 39.3513 Å × 29.2100 Å × 28.2521 Å) under default parameters. Then, at least 10 docked poses were generated for each ligand, and the docked conformations with the highest negative docking scores and least root-mean-square deviation (RMSD) values (by default 0 in AutoDock Vina) were extracted for binding pose and intermolecular interaction analysis, under default parameters, in Maestro tool of Schrödinger suite 2018-4 [60]. All 3D and 2D figures were generated in Maestro tool of Schrödinger suite2018-4 [60]. SwissADME online server (http://www.swissadme.ch) was used to understand the properties related to the pharmacokinetics and drug likeliness of the selected compounds such as absorption, distribution, metabolism, and excretion (ADME) [61]. The bioavailability is another important parameter that makes any drug a promising therapeutic, and therefore, several parameters were observed during the ADME analysis such as blood-brain barrier (BBB), permeability, and cytochrome inhibition activities.

Molecular Dynamics Simulation
The molecular dynamics (MD) simulations for the selected best poses of the top six protein-ligand complexes were performed using with academic Maestro-Desmond v5.6 suite [62,63]. Initially, the docked complex was placed in the center of the orthorhombic box (10 Å × 10 Å × 10 Å) and amended with an explicit TIP4Pwater model using the system builder module. Thereafter, the complete system was neutralized using the counter sodium and chloride ions, while being aced at a distance of 20 Å around the docked ligand in the binding pocket of protein, and then minimized under default parameters using the minimization tool. Next, the whole system was simulated using an NPT ensemble, maintained by the Nose-Hoover thermostat and the Martyna-Tobias-Klein barostat method [64], with the temperature set at 300 K and pressure set at 1.013 bar under default parameters. For the simulations, the cutoff radius in Coulomb interactions was set at 9.0 Å, and the longrange electrostatic interactions were computed via the particle mesh Ewald method [65]. Each complex was simulated for 100 ns under optimized potentials for liquid simulations (OPLS)-2005 force-field parameters, and a total of 5000 frames were saved for analysis. The generated JEV-RdRp trajectories with selected ligands were analyzed for statistical parameters, including root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) analysis, and intermolecular interactions (protein-ligand contact mapping) as a function of 100 ns using the Simulation Interaction Diagram (SID) tool implemented in free academic Desmond module with Maestro-Schrödinger suite 2018-4 interface [62,63].

Essential Dynamics
The essential dynamics, defined in terms of principal component analysis (PCA), to study the dynamic motion of the protein was assessed in the presence of docked ligands from the respective MD simulation trajectories using the Bio3d package [66]. For this purpose, a total of 5000 frames of the Cα atoms from the respective MD simulation trajectories were superimposed on the initial pose to reduce the root-mean-square variations between similar residues by utilizing the fit.xyz function. These superimposed poses were then analyzed to generate respective plots of PCA components from each simulation trajectory using the pca.xyz function of the plot (pc) in the Bio3d package [66].

Prime MM/GBSA Binding Free Energy Calculations
To calculate the binding free energy and ligand strain energy for the docked poses, the endpoint MMGBSA calculation was applied to the complete 100 ns MD simulation trajectories using the thermal_mmgbsa.py python script in the Prime MMGBSA module of the Schrödinger suite (Schrödinger Release 2020-4: Prime, Schrödinger, LLC, New York, NY, USA, 2020), where each snapshot was treated for the removal of solvent and ions, as well as split into individual protein and ligand conformation for the free energy calculation. The net free binding energy (∆G) was calculated using the following equation: where ∆G Bind denotes the binding free energy; ∆G complex indicates the binding free energy of the complex; ∆G Receptor and ∆G Ligand exhibit the energy for receptor and ligand, respectively. Finally, the computed binding free energy, along with energy dissociation components, was provided as average with the standard deviation for each simulation trajectory.

Structure-Based Virtual Screening and ADME
The natural products reported in the E. angustifolia were collected from the reported literature, and the respective 3D structures were collected retrieved from the PubChem database. The collected compounds were used as ligand libraries for the structure-based virtual screening against the nucleotide GTP-binding pocket in the crystal structure of JEV-RdRp (PDB ID: 4HDG). The top conformations of the screened ligands exhibited binding affinity ranging from −11.1 to −1.9 kcal/mol, targeting the binding pocket of the selected JEV-RdRp protein (Table S1). Thus, the top docked poses of the first six compounds-namely, Echinacoside, Echinacin, Rutin, Cynaroside, Quercetagetin 7glucoside, and Kaempferol-3-glucoside-with significant binding scores (>−9 kcal/mol) were considered for the ADME/Tox and redocking analysis to find drug-likeness properties and the most ideal docking conformation, respectively, in the selected binding pocket of the viral protein.
Drug likeness and pharmacological characters are considered to be important factors for understanding the medicinal application of the compounds. Therefore, the ADME analysis of the top six compounds, i.e., Echinacoside, Echinacin, Rutin, Cynaroside, Quercetagetin 7-glucoside, and Kaempferol-3-glucoside, were predicted using SwissADME, resulting in the properties related to the pharmacokinetics and toxicity (Table S2). All of the selected compounds were found to be non-inhibitors of Cytochrome P450 2D6 (CYP2D6), as the inhibition of CYP2D6 leads to drug-drug interactions. Additionally, each Life 2022, 12, 952 5 of 17 of these natural compounds was found to be impermeable to the Blood-brain barrier (BBB) ( Table S2). Based on these findings, as well as other properties such as pharmacokinetics and drug-likeness, the selected compounds have considerable medicinal properties.

Redocking and Molecular Contact Analysis
From the virtual screening, six potential natural compounds-Echinacoside, Echinacin, Rutin, Cynaroside, Quercetagetin 7-glucoside, and Kaempferol-3-glucoside-were selected for redocking in the selected binding pocket of JEV-RdRp against ATP and GTP as reference compounds using AutoDock Vina. Afterward, the docked poses with the highest negative docking energy values corresponding to zero RMSD values for each natural compound were considered for further computational analysis ( Figure 1). In this study, among the selected natural compounds, Echinacoside docked with JEV-RdRp showed maximum docking energy (−11.1 kcal/mol), while JEV-RdRp-Kaempferol-3-glucoside docked complex was noted for lowest docking scores (10 kcal/mol) by comparison to reference compounds ATP (8.6 kcal/mol) and GTP (−9.0 kcal/mol). Additionally, a substantial number of intermolecular interactions, including hydrogen bond formation, π-π stacking, π-cation, hydrophobic, polar, negative, positive, glycine, and salt bridge interactions were also noted in the docked complexes by comparison to the reference docked complexes (Table 1, Figure S1). Therefore, the calculated docking scores and intermolecular contact profiling between the docked natural compounds and JEV-RdRp indicates the stability of the respective docked complexes.

Molecular Dynamics Simulation Analysis
A molecular dynamics simulation is used to predict the stability of protein-ligand complexes with respect to simulation intervals. In this study, the molecular dynamics simulation was analyzed in terms of the last pose from the simulation trajectory, protein and protein-fit-ligand root-mean-square deviation (RMSD), protein and protein-fit-ligand root-mean-square fluctuation (RMSF), and protein-ligand contact mapping as a function of 100 ns interval. Figure 2 shows the last poses of the MD trajectory were extracted and compared with the initially docked poses to monitor the relative occupation of the docked ligands in the protein structure. Notably, all of the docked natural compounds showed relatively substantial residence in the selective pocket of the JEV-RdRp as a function of 100 ns MD simulation interval, except for an acceptable deviation in the ligand conformations, which was noted against reference compounds ( Figure S2). These results suggested the docked natural compounds as substantial inhibitors by comparison to the reference compounds, i.e., ATP and GTP. also noted in the terminal regions and residues forming direct contact or adjacent residues to the interacting residues with the docked ligands. Moreover, the calculated protein-fit-ligand RMSD values showed considerable deviations (<4.5 Å) (Figure 3) throughout the 100 ns simulation interval against reference compounds, i.e., ATP (<4.0 Å) ( Figure 3g) and GTP (<4.0 Å) (Figure 3h). Notably, among the selected natural compounds, Rutin, Cynaroside, and Kaempferol-3-glucoside showed the most substantial stability and equilibrium in protein-fit-ligand values as a function of 100 ns against reference compounds, i.e., ATP and GTP. Moreover, the calculated RMSF values also showed variation within 2 Å ( Figure S4) during the 100 ns, supporting the observed stability of the docked ligands with the JEV-RdRp during MD simulation.  Moreover, the protein-ligand contacts were also extracted from the 100 ns MD simulation trajectories for all of the docked complexes. Specifically, each type of intermolecular interaction, hydrogen bonding, hydrophobic interactions, water bridge formation, and ionic interaction was extracted and plotted as a fraction of the total molecular contact formed during the 100 ns interval for all the simulated complexes.
In the case of the Echinacoside-JEV-RdRp docked complex, it exhibited hydrogen bond formation for more than 100% of simulation time in Ser 666 and Asp 669 residues, while Moreover, the protein-ligand contacts were also extracted from the 100 ns MD simulation trajectories for all of the docked complexes. Specifically, each type of intermolecular interaction, hydrogen bonding, hydrophobic interactions, water bridge formation, and ionic interaction was extracted and plotted as a fraction of the total molecular contact formed during the 100 ns interval for all the simulated complexes. In the case of the Echinacoside-JEV-RdRp docked complex, it exhibited hydrogen bond formation for more than 100% of simulation time in Ser 666 and Asp 669 residues, while Trp 800 and Ile 802 residues showed 50% and 35% of total simulation time for the hydrophobic interaction with the docked ligand; Gly 472 , Cys 714 , and Arg 474 residues participated in water bridges formation, in which Arg 474 residue showed more than 100% of water bridge formation during the simulation time in addition to a significant H-bond. Additionally, Arg 734 and Arg 742 residues formed an ionic bond for more than 20% simulation time (Figure 4a). Likewise, the Echinacin-JEV-RdRp complex exhibited hydrogen bond interaction in Asp 541 and Lys 463 residue for 90% and 70% of total simulation time; Trp 800 and Arg 474 residues were also noted to form hydrophobic interaction for 100% of simulation time. Gly 605 and Arg 460 residues showed more than 80% of the total interaction fraction in water bridge formation during the simulation period (Figure 4b). Additionally, protein-ligand contact analysis of Rutin-JEV-RdRp showed a high contribution of Asp 668, Asp 669 , and Ser 604 in hydrogen bond formation (100%) during the total simulation interval. Tyr 610 residue participated in hydrophobic interaction for more than 50% of the simulation time. A water bridge was formed by Gly 412 , Trp 477 , and Arg 484 for more than 50% of simulation time and by Gly 412 for 100% of simulation time (Figure 4c). Furthermore, the Cynaroside-JEV-RdRp complex exhibited a hydrogen bond with Glu 510 and Asp 668 residues for 100% of the simulation time. Tyr 610 also showed the hydrophobic interaction with the docked ligand at a 100% of simulation interval. Asp 541 and Asp 669 residues were involved in water bridge formation for more than 50%, along with the hydrogen bond (Figure 4d). In comparison, in the Quercetagetin 7-glucoside-JEV-RdRp complex, Glu 510 residue showed more than 100% of the total interaction fraction in hydrogen bond formation during the simulation interval, and Gly 605 and Asp 668 residues participated in hydrogen bond formation for more than 50% of simulation time; Asp 668 also participated in water bridge formation (50%) during the simulation interval. Ile 802 residue was noted for hydrophobic interaction with the docked ligand for more than 50% of the simulation time (Figure 4e). In the protein contact mapping analysis of the Kaempferol-3-glucoside-JEV-RdRp complex, Asp 541 exhibited more than 100% of hydrogen bond formation during the 100 ns simulation time, whereas Tyr 610 , Trp 800 , and Ile 802 residues showed 50% of hydrophobic interaction during the total simulation period. Notably, Asp 668 residue showed more than 50% of the water bridge interaction along with the hydrogen bond (Figure 4f). In contrast, the protein-ligand mapping of the JEV-RdRp with its native ligand ATP exhibited more hydrogen bond formation for 100% of simulation time with Arg 460 , Lys 463 , Lys 471 , Arg 474 , and Arg 734 residues, along with water bridge formation (50% of the interaction fraction). Ile 802 and Lys 471 residues participated in hydrophobic (50%) and ionic interaction (30%) during the 100 ns simulation period (Figure 4g). When the protein-ligand mapping of the JEV-RdRp with its native ligand GTP was analyzed, Arg 460 , Arg 474 , and Arg 742 participated in hydrogen bond formation for more than 100%, in addition to the water bridge interaction of the total simulation time. Additionally, Lys 463 and Arg 734 exhibited hydrophobic and ionic bonds for more than 90% and more than 50% of the total interaction fraction (Figure 4h). The solvent accessible surface area (SASA) of the selected natural compound and the reference molecules were also studied to track the protein's flexibility, stability, and folding in the presence and absence of ligands ( Figure S5) [43].
Collectively, the analysis of the MD simulations supported the selected natural compounds as inhibitors of JEV-RdRp against the reference compounds, i.e., ATP and GTP. Hence, based on the statistical analysis and intermolecular interaction profiling of the docked natural compounds in the selected pocket of JEV-RdRp, the potential ligands can be placed in the order of Rutin, Cynaroside, Kaempferol-3-glucoside, Echinacoside, Quercetagetin 7-glucoside, and Echinacin, exhibiting the most stability with the JEV-RdRp protein. Collectively, the analysis of the MD simulations supported the selected natural compounds as inhibitors of JEV-RdRp against the reference compounds, i.e., ATP and GTP. Hence, based on the statistical analysis and intermolecular interaction profiling of the docked natural compounds in the selected pocket of JEV-RdRp, the potential ligands can be placed in the order of Rutin, Cynaroside, Kaempferol-3-glucoside, Echinacoside, Quercetagetin 7-glucoside, and Echinacin, exhibiting the most stability with the JEV-RdRp protein.

Principal Component Analysis
Studying the essential dynamics based on principal component analysis provides a method to collect the domain dynamics and displacement of atoms in the protein, required for biological function. Figure 5 shows the percentage of variance (%) for the cal-

Principal Component Analysis
Studying the essential dynamics based on principal component analysis provides a method to collect the domain dynamics and displacement of atoms in the protein, required for biological function. Figure 5 shows the percentage of variance (%) for the calculated mean square position variations in the covariance matrix as a function of extracted 20 eigenmodes from the respective MD simulation trajectories of the docked complexes. Notably, all of the JEV-RdRp structures docked with natural compounds relatively showed similar drops in eigen fraction values by comparison to the protein structures docked with reference compounds GTP but not ATP ( Figure 5 and Figure S6). However, no significant changes in the eigen fraction were observed from 6 to 20 eigenvalues. These results indicated the significant conformational changes in the viral protein structure docked with natural compounds to attain the most stable complex formation against reference compounds. docked with reference compounds GTP but not ATP (Figures 5 and S6). However, no significant changes in the eigen fraction were observed from 6 to 20 eigenvalues. These results indicated the significant conformational changes in the viral protein structure docked with natural compounds to attain the most stable complex formation against reference compounds.   Figure S6). These observations suggested that the docked natural compounds have the potential to induce conformational fluctuations in JEV-RdRp to distribute its biological function for the replication of the JE virus.

Binding Free Energy Calculation
All of the generated MD simulation trajectories created as a function of 100 ns interval were treated for free binding energy calculation using MM/GBSA method against reference complexes. Figure 6 shows high and low binding energy values for all the docked complexes with respect to time. Notably, Echinacoside, Echinacin, and Rutin were found to exhibit the highest binding energy values (>−100 kcal/mol) for some poses, with net binding free energy values of −80 ± 8.04, 81.67 ± 8.31, and 80.33 ± 5.54 kcal/mol, respectively. In comparison, the other three natural compounds showed mean binding free energy between 66 to 57 kcal/mol. Interestingly, the calculated binding free energy values were relatively higher than those of the reference complexes, i.e., ATP (69.62 ± 8.62 kcal/mol) and GTP (47.98 ± 11.41 kcal/mol) ( Figure 6). Additionally, net energy dissociation components-namely, ∆G Bind Coulomb , ∆G Bind Covalent , ∆G Bind Hbond , ∆G Bind Lipo , ∆G Bind Packing , ∆G BindSolvGB , and ∆G Bind vdW -were determined using the Prime MM/GBSA method and were also studied on the complete 100 ns MD trajectories of each complex. Interestingly, substantial contributions of ∆G Bind Coulomb , ∆G Bind Lipo , and ∆G Bind vdW interactions were noted for favorable energy, while ∆G Bind Covalent and ∆G BindSolvGB exhibited unfavorable energy to the net binding free energy for all of the protein-ligand complexes during the 100 ns MD simulation (Table 2). Altogether, the natural compounds were noted for higher binding free energy, except Cynaroside, in the selective pocket of JEV-RdRp against reference compounds as a function of 100 ns interval.

Conclusions
Japanese encephalitis virus-RNA-dependent RNA polymerase (JEV-RdRp) is a potential target for antiviral drugs, as they are responsible for viral genome replication. The purpose of this study was to identify potent antiviral compounds from E. angustifolia by inhibiting JEV-RdRp using in silico approach. Based on the substantial docking energy (>−10 kcal/Mol) and pharmacokinetics analysis, six compounds-Echinacoside, Echinacin, Rutin, Cynaroside, Quercetagetin 7-glucoside, and Kaempferol-3-glucoside-were selected. Based on the hydrogen, hydrophobic, and other molecular interaction analyses, the stability of the respective compounds with JEV-RdRp's binding site was studied. The final confirmation of the stability of the protein-ligand complex was made with the analysis