key: cord-0924684-8apy731y
authors: Mosquera-Yuqui, Francisco; Lopez-Guerra, Nicolas; Moncayo-Palacio, Eduardo A.
title: Targeting the 3CLpro and RdRp of SARS-CoV-2 with phytochemicals from medicinal plants of the Andean Region: molecular docking and molecular dynamics simulations
date: 2020-10-21
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1835716
sha: edee9cb52f6a5e3d12fca1b67bb044d6ee52f1e8
doc_id: 924684
cord_uid: 8apy731y

Given the highly contagious nature of SARS-CoV-2, it has resulted in an unprecedented number of COVID-19 infected and dead people worldwide. Since there is currently no vaccine available in the market, the identification of potential drugs is urgently needed to control the pandemic. In this study, 92 phytochemicals from medicinal plants growing in the Andean region were screened against SARS-CoV-2 3 C-like protease (3CLpro) and RNA-dependent RNA polymerase (RdRp) in their active sites through molecular docking. The cutoff values were set from the lowest docking scores of the FDA-approved drugs that are being used to treat COVID-19 patients (remdesivir, lopinavir, and ritonavir). Compounds with docking scores that were lower than cutoff values were validated by molecular dynamics simulation with GROMACS, using root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and intermolecular hydrogen bonds (H-bonds). Furthermore, binding free energies were estimated using the MM-PBSA method, and ADMET profiles of potential inhibitors were assessed. Computational analyses revealed that the interaction with hesperidin (theoretical binding energies, ΔG(bind) = −15.18 kcal/mol to 3CLpro and ΔG(bind) = −9.46 kcal/mol to RdRp) remained stable in both enzymes, unveiling its remarkable potential as a possible multitarget antiviral agent to treat COVID-19. Importantly, lupinifolin with an estimated binding affinity to 3CLpro higher than hesperidin (ΔG(bind) = −20.93 kcal/mol) is also a potential inhibitor of the 3CLpro. These two compounds displayed suitable pharmacological and structural properties to be drug candidates, demonstrating to be worthy of further research. Communicated by Ramaswamy H. Sarma

Coronaviruses (CoVs) are a group of enveloped, large nonsegmented positive-sense RNA viruses able to cause enteric, respiratory, and central nervous diseases in animals, including humans (McIntosh, 1974; Weiss & Navas-Martin, 2005) . Whereas most of human CoVs cause mild respiratory infections, the outbreaks of severe acute respiratory syndrome (SARS), Middle East respiratory syndrome (MERS), and coronavirus disease 2019 which emerged in 2003, 2012, and 2019, respectively, have demonstrated how devastating and life-threatening they can be.

COVID-19 is caused by a new type of transmissible pathogenic SARS-CoV (SARS-CoV-2) -formerly 2019-nCoV -which belongs to the zoonotic CoVs of the genus Betacoronavirus in the family Coronaviridae. Symptoms include fever, dry cough, fatigue, shortness of breath, and in some cases gastrointestinal distress (Guo et al., 2020) .

The SARS-CoV-2 genome size of about 30 kb encodes multiple structural and non-structural proteins (Ahmed et al., 2020) . The spike (S), envelope (E), membrane (M), and nucleocapsid (N) proteins are considered essential to produce structurally complete viral particles (Ahmed et al., 2020; Mangar et al., 2020) . The 3 C-like protease (3CLpro), also known as main protease, is one of the two proteases that SARS-CoVs utilize for replication and infection processes (Pillaiyar et al., 2016) . Unlike the catalytic tryad found in serine proteases and other cysteine proteases, a catalytic dyad is present in the active site of SARS-CoVs main protease (Anand et al., 2003; Zhang et al., 2020) . The RNA-dependent RNA polymerase (RdRp), also named nsp12, is the most highly conserved protein among RNA viruses and its role is to catalyze the viral genome replication and transcription (J acome et al., 2020) . RdRp activity is dependent on magnesium ions and requires the non-structural proteins nsp7 and nsp8 for complete activity (Kirchdoerfer & Ward, 2019) . The active site of the SARS-CoV-2 RdRp is formed by conserved polymerase motifs (A-G), where the motif A and motif C contain the divalent-cation-binding amino acid D618, and the catalytic residues 759SDD761, respectively .

Taken together, these enzymes are attractive drug targets to treat these human CoVs.

Potential COVID-19 treatments include ivermectin, a broadspectrum anti-parasitic agent (Caly et al., 2020) ; type I interferons, a group of cytokines with a broad spectrum activity against RNA viruses (Mantlo et al., 2020) ; favipiravir, a viral RdRp inhibitor ; and convalescent plasma. The protease inhibitors lopinavir and ritonavir, used to treat HIV infections, and recommended by the National Health Commission (NHC) of China to treat COVID-19, did not significantly accelerate clinical improvement, reduce mortality, or decrease throat viral RNA detection in a randomized, controlled, open-label, not blinded trial involving hospitalized patients with severe COVID-19 . However, the work had several limitations and larger studies with greater variety of patients are needed to evaluate the efficacy of this treatment. Remdesivir, a nucleoside analog used to block viral RdRps including filoviruses, paramyxoviruses, pneumoviruses, as well as animal and human CoVs, was used to treat severely ill patients with COVID-19 in a randomized, double-blind, placebo-controlled, multicentre trial, showing a numerically faster time to clinical improvement but not statistically significant . Using insect cells, it was found to be a direct-acting antiviral inhibiting with the same potency and mechanism of action the RdRp of SARS-CoV-2 and its ortholog in SARS-CoV (Gordon et al., 2020) . Remdesivir was approved on May 1, 2020, by the United States Food and Drug Administration (FDA) to treat hospitalized COVID-19 patients. More recently, in a large trial including hospitalized patients critically ill with COVID-19, the steroid dexamethasone reduced deaths by one-third in patients receiving invasive mechanical ventilation and by one-fifth in patients receiving oxygen but were not on ventilators. Importantly, this steroid had not effect on patients with less severe cases those not receiving respiratory support (Horby et al., 2020) .

The high diversity and complex molecular structures of natural compounds make them an abundant biological source for drug discovery, even more, considering the limited number of plant species that have been explored for pharmaceutical purposes (Saklani & Kutty, 2014) . Flavonoids are naturally occurring compounds that have been proposed as attractive antiviral agents given their broad mechanisms of action including transcription and translation blocking, as well as their general safety and non-cytotoxicity to human cells (Lalani & Poh, 2020; Paduch & Kandefer-Szerszen, 2014) . Herbacetin, rhoifolin and pectolinarin were demonstrated to inhibit the protease activity of the SARS-CoV 3CLpro in vitro (Jo et al., 2020) . We focused our research to the Andean region considering that it is one of the most biodiverse regions in the world due to the different environmental niches resulting from its great elevational and latitudinal diversity gradient (Anthelme et al., 2014; Mutke et al., 2014; Pennington et al., 2010) . The present work was aimed to perform molecular docking, molecular dynamics simulations, binding free energy estimations, and to predict ADMET profiles of natural compounds identified in medicinal plants growing in the Andean region as potential inhibitors of the 3CLpro and RdRp of SARS-CoV-2.

Two SARS-CoV-2 proteins were selected as targets based on the key role they played in the replication/transcription machinery. The 3 D structures of the 3CLpro (PDB ID: 6LU7) and RdRp (PDB ID: 6M71) were downloaded from the Protein Data Bank. Water molecules and inhibitor were removed from 6LU7 using Discovery Studio Visualizer (2016). Similarly, co-factors nsp7 and nsp8 were removed from 6M71. Missing atoms for residues F70, K73, R74, E83, K98, F101, F102, I114, R365, and D824, in the RdRp file, were modeled using the Swiss-Pdb Viewer v4.1.0 (Guex & Peitsch, 1997) . AutoDock Tools v1.5.6 (Sanner, 1999) was used to add polar hydrogens to the structures and convert them in pdbqt format.

Ligands used for molecular docking included 92 phytochemicals from 20 medicinal plants growing in the Andean region in South America, as well as some FDA-approved drugs with potential anti-SARS-CoV-2 (darunavir, nelfinavir, paritaprevir, saquinavir, setrobuvir, simeprevir, and sofosbuvir), the aminoquinolines chloroquine and hydroxychloroquine, and cinnamaldehyde as negative controls due to its poor affinity to these SARS-CoV-2 viral targets (Elfiky, 2020) . Natural compounds and synthetic drugs were downloaded from PubChem and DrugBank, respectively. Ligands were converted to PDB format using Open Babel (O' Boyle et al., 2011) and later prepared to pdbqt format using AutoDock Tools v1.5.6 (Sanner, 1999).

The active pockets on the target proteins were identified using the CASTp (Computed Atlas of Surface Topography of proteins) server (Tian et al., 2018) .

Docking analysis was performed using the AutoDock Vina software (Trott & Olson, 2010) . Grid boxes of 24 Å x 26 Å x 30 Å centered at x, y, z ¼ À10, 12, 69 and 36 Å x 34 Å x 19 Å centered at x, y, z ¼ 113, 113, 130 were set for 3CLpro and RdRp, respectively. A grid spacing of 1 Å and default exhaustiveness were used. After docking, the structures were examined in Discovery Studio Visualizer (2016).

The resulting protein-ligand complex structures from molecular docking were used to perform molecular dynamics (MD) simulations in GROMACS 2018 (Van Der Spoel et al., 2005) with MPI support at the Oklahoma State University HPC system. Ligand and protein parameters were generated using CHARMM general force-field (Vanommeslaeghe et al., 2010) and CHARMM36 force-field (Huang & MacKerell, 2013) , respectively; the topologies for both 3CLpro and RdRp were generated using the pdb2gmx utility included in GROMACS, while the CGenFF online program was used for the ligands (https://cgenff.umaryland.edu) . In order to prepare the system for MD simulations, all protein-ligand and protein with no ligand systems were centered in a rhombic dodecahedron box with a box-system distance of 1.0 nm and solvated with three-point transferable intermolecular potential (TIP3P) water (Bjelkmar et al., 2010) .

Charge neutralization was carried out for the 3CLpro system by adding four sodium ions, while six sodium ions were added for the RdRp system. Then, the systems were relaxed through 50000 steps of the steepest descent algorithm for energy minimization calculations at a tolerance value of 1000 kJ/(mol.nm). This was followed by the equilibration with position restraint on the protein and ligand molecules for 0.1 ns using NVT and NPT ensembles, wherein the systems were heated to 300 K using Berendsen thermostat (Berendsen et al., 1984) with a coupling time of 0.1 ps, and the pressure was maintained with a coupling to a reference pressure of 1 bar. For energy minimization, NVT, and NPT relaxation simulation, a smooth force-switch 1.2 nm cutoff was used in short-range interactions, and long-range electrostatics were evaluated using the PME (Particle-Mesh-Ewald) (Darden et al., 1993) ; additionally, hydrogen-bonds were restrained with the LINCS algorithm (Hess et al., 1997) .

Final MD simulations of 200 ns were performed without restraints using an integration time-step of 2 fs, and trajectory snapshots were captured at every 1 ps. Root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and number of intermolecular hydrogen bonds (H-bonds) were chosen as parameters to analyze MD trajectories using GROMACS utilities and plotted with Xmgrace. To further confirm the stability of protein-ligand complexes, changes in secondary structures during MD simulations were assessed using the Dictionary of Secondary Structure of Proteins (DSSP) algorithm.

The Molecular Mechanics Poisson-Boltzmann surface area (MM-PBSA) method (Kollman et al., 2000) , implemented in the g_mmpbsa tool (Kumari et al., 2014) of GROMACS, was used to estimate the relative binding free energy of each protein-ligand complex. Briefly, the binding free energy (DG bind ) of a protein-ligand complex in solvent is calculated as the free energy difference between the complex (G complex ) and the summation of the free energy of the protein (G protein ) and ligand (G ligand ):

The analyses were performed for 500 snapshots collected consecutively at an interval of 40 ps from the last 20 ns of MD simulations. The solute dielectric constant (e in ) was assigned based on the protein-ligand interfaces by using the scheme: 1 for non-polar residues, 2 for polar not charged residues, and 4 for polar charged (positively and negatively) residues (Ravindranathan et al., 2011; Wang et al., 2019) .

A total of 19 molecular descriptors were calculated to predict the absorption, distribution, metabolism, elimination, and toxicity (ADMET) profile of the selected potential inhibitors of the 3Clpro and RdRp of SARS-CoV-2, using the Web servers SwissADME (Daina et al., 2017) and pkCSM (Pires et al., 2015) which have been extensively validated with experimental data. The Lipinski's rule of five considering key physicochemical parameters (molecular weight, lipophilicity, polar surface area, hydrogen bonding, and charge) was used to evaluate the drug likeness of the potential inhibitors.

Docking scores of FDA-approved protease and RdRp inhibitors which are used to treat COVID-19 patients were calculated in order to set cutoff values. Remdesivir docking scores (-8.2 kcal/mol to 3CLpro and À7.5 kcal/mol to RdRp) were lower than or equal to those found for lopinavir (-7.9 kcal/ mol to 3CLpro and À6.9 kcal/mol to RdRp) and ritonavir (-8.2 kcal/mol to 3CLpro and À 7.2 kcal/mol to RdRp). Then, À8.2 kcal/mol to 3CLpro and À7.5 kcal/mol to RdRp were set as cutoff values.

The docking scores to RdRp of all tested FDA-approved drugs which are being proposed as candidates against SARS-CoV-2, excepting sofosbuvir, were lower than remdesivir. Moreover, paritaprevir, setrobuvir, simeprevir, and nelfinavir also showed higher affinity to 3CLpro than ritonavir and remdesivir (Supplementary Table 1 ). The low affinities of chloroquine and hydroxychloroquine to both enzymes (Supplementary Table 1 ) indicate that they do not interact effectively with these viral targets. In fact, the mode of action of these aminoquinolines is not related to RdRp or 3CLpro inhibition. CQ impairs the autophagosome fusion with lysosomes (Mauthe et al., 2018) , and interferes with the glycosylation of viral or host proteins. More specifically, CQ inhibits the glycosylation of HIV viral particles or SARS-CoV human receptor angiotensin-converting enzyme 2 (ACE2) (Savarino et al., 2004; Vincent et al., 2005) From the analysis of the 92 phytochemicals, 10 exhibited docking scores lower than the cutoff for 3CLpro (corosolic acid, friedelin, hesperidin, guaijaverin, hyperin, lupinifolin, mangiferin, pinocembrin-7-O-rutinoside, quercitrin, and rutin) while 12 had docking scores lower than the cutoff for RdRp (Supplementary Table 1 ). These compounds showing such affinities were selected to perform MD simulations in order to evaluate the stability of protein-ligand complexes.

Interestingly, unlike their aglycones, hesperidin (hesperetin-7-O-rutinoside), rutin (quercetin-3-O-rutinoside), and quercitrin (quercetin-3-O-a-rhamnoside) showed significant affinities to both viral targets (Supplementary Table 1 ). In addition, the quercetin glycosides guaijaverin (quercetin-3-O-a-arabinoside) and hyperin (quercetin-3-O-b-galactoside) revealed significant affinities to 3CLpro. By comparing the protein-ligand interactions, the glycoside moiety was involved in molecular interactions.

Molecular Dynamics (MD) is a computational approach used to analyze the dynamic behaviour of complex systems where atoms and molecules interact as a function of time. In general, computational methods improve the efficiency in drug discovery, Importantly, unlike general molecular docking methods, MD simulations consider the flexibility of the targets and combined with binding energy calculations a more accurate prediction of potential inhibitors can be developed (Liu et al., 2018) . Structural parameters including RMSD, RMSF, Rg, and number of intermolecular H-bonds were used to evaluate the stability, dynamic behaviour, and compactness of protein-drug complexes.

Root mean square deviation (RMSD) of the protein backbone was used to analyze the stability of 3CLpro and RdRp in complex with the selected phytochemicals, where a value lower than that calculated to the free protein indicates a stable complex (Liu et al., 2017) . The mean RMSD values for free 3CLpro and 3CLpro in complex with hesperidin, hyperin, lupinifolin, and rutin were 2.37 Å, 2.13 Å, 2.11 Å, 2.04 Å, and 2.00 Å, respectively. Similarly, the mean RMSD for free RdRp and RdRp in complex with hesperidin, pinocembrin-7O-rutinoside, speciophylline, and vitexin were 2.72 Å, 2.44 Å, 2.69 Å, 2.38 Å, and 2.67 Å, respectively. Then, these phytochemicals did not disturb the structural stability of these enzymes (Figure 1 ). Backbones with mean RMSD higher than free protein are plotted in Figure S1 .

To ensure the binding stability of these compounds in the active site, ligand positional RMSD was calculated. In complex with 3CLpro, high fluctuations were found for hyperin, corosolic acid, friedelin, guaijaverin, mangiferin, pinocembrin-7-O-rutinoside, and quercitrin (Figure 2(A) , Figure S2 ). In contrast, hesperidin and lupinifolin were stable throughout the simulation (Figure 2(A) ). Since rutin fluctuated from 0.2 nm to 2.3 nm at the end of the simulation, ten snapshots were downloaded in the interval of 20 ns during the entire simulation wherein it was observed that rutin had moved considerably from the active site, while hesperidin and lupinifolin remained firmly bound (Figure 3) . It suggests the inability of rutin to inhibit efficiently the 3CLpro of SARS-CoV-2. Among the ligands that formed stable complexes with RdRp, only hesperidin showed stable binding in the active site ( Figure  2 (B), Figure 4 ). Following analyses were done for hesperidin and lupinifolin which did not disturb the protein stability and remained firmly bound to the active site.

In order to calculate the residual mobility in the absence and presence of the phytochemicals, RMSF analysis was calculated and plotted against the residue number and was supported with the DSSP analysis. Part of the noticeable fluctuations occurring to backbones of 3CLpro complexes around position 50 corresponds to residues belonging to a loop surrounding the active site ( Figure 5(A) ). It has been demonstrated that loops can be involved in ligand binding and conformational changes can often take place in these sites (Lee et al., 2012; Teague, 2003) . In contrast, residues interacting with stable ligands including the catalytic dyad H41 and C145 for 3CLpro, as well as the catalytic residues 759SDD761 and the divalent-cation-binding residue D618 for RdRp were among the most stable residues. Missing residues between K50-Y69, F102-P112, and L895-M906 in the crystal structure of RdRp were evidenced by ambiguous fluctuations in those positions ( Figure 5(B) ). In all complexes, there were no significant structural changes during the throughout simulation. The helical and b-sheet content remained stable ( Figure S3 , Figure S4 ).

Radius of gyration (Rg) is a measure of the compactness of a protein which allows to understand its folding properties (Lobanov et al., 2008) . Small Rg values indicate a tight packing whereas high Rg values show a floppy packing. A relative constant Rg value through time indicates that the ligands hold the folding behavior of the protein whereas abrupt fluctuations of the Rg values denote protein folding instability (Khan et al., 2020) .

The 3Clpro in complex with hesperidin (mean Rg ¼ 2.24 nm) and lupinifolin (mean Rg ¼ 2.23 nm) presented more compactness compared with the free protein (mean Rg ¼ 2.27 nm). On the other hand, the RdRp in complex with hesperidin (Rg ¼ 3.04 nm) was slightly less compact than the free RdRp (Rg ¼ 3.02 nm). In all these complexes, the Rg kept constant with no abrupt fluctuations through the time, indicating that hesperidin, rutin, and lupinifolin maintain the folding behavior of 3Cl pro while hesperidin maintains RdRp folding ( Figure 6 ).

It is well known that H-bonds are responsible for the secondary and tertiary structural protein motifs. Formation of Hbonds between a ligand and a protein motif explains the binding affinity of a drug towards a protein target in molecular dynamics simulations; so, the more number of H-bonds the stronger interactions (Menendez et al., 2016) .

For 3CLpro, a mean number of hydrogen bonds of 5-6, and 1, were found for hesperidin and lupinifolin, respectively (Figure 7(A) ). The stability of lupinifolin can be explained by the six hydrophobic interactions revealed from the docking analysis. In addition, the contribution of van der Waals interactions to the binding energy of lupinifolin was the lowest among the complexes (Table 1) . Even though pinocembrin 7-O-rutinoside formed 4 hydrogen bonds, it was discarded because of the backbone instability obtained in the RMSD analysis. On the other hand, for RdRp in complex with hesperidin, an average of 2-3 hydrogen bonds was found (Figure 7(B) ). The formation of these hydrogen bonds explains the high affinities of the ligands to 3Clpro and RdRp, as well as the stability of these complexes over time.

In general, small non-covalent or reversible covalent inhibitors, such as phytochemicals, display several advantages regarding side effects and toxicity compared with covalent inhibitors (Pillaiyar et al., 2016) .

By rescoring docked complexes, the MM-PBSA method has demonstrated to be a valuable tool in drug discovery to remove possible false positive compounds obtained by docking analysis performed with standard methods (Kumari et al., 2014; Rastelli et al., 2009) . Binding energies (DG bind ) of the stable complexes with ligands firmly bound to the active site of 3CLpro and RdRp were estimated according to the MM-PBSA method (Table 1) . Molecular interactions of hesperidin and lupinifolin in complex with 3CLpro taken from molecular docking and during the simulation at 100 ns and 200 ns included mostly polar not charged and non-polar residues (Figure 8) , then a value of e in ¼ 2 was set for the MM-PBSA calculations of these complexes. On the other hand, a value of e in ¼ 4 was used for the RdRp-hesperidin complex because the highly charged protein-ligand interface including aspartic acid, glutamic acid, arginine, and lysine residues (Figure 9 ).

Results indicated that lupinifolin showed more affinity to 3CLpro compared to hesperidin (Table 1) . However, hesperidin also formed a stable complex with RdRp showing a DG bind ¼ À9.46 kcal/mol. In these complexes, the contribution of van der Waals interactions to the binding energy was lower than electrostatic and non-polar solvation energies. To understand protein-ligand associations, the total binding energies were decomposed into the contribution made by each residue. For 3CLpro complexes, H41, C145, and M165 were found to strongly interact with both hesperidin and lupinifolin, (Figure 10 ). For RdRp-hesperidin complex, residues S561 and T565 surrounding the active site as well as D703 were part of those contributing negatively to the binding energy (Figure 11 ).

In silico ADMET analysis is used to predict the pharmacokinetic profile of compounds, before experimental procedures (Jayaraj et al., 2020) . The bioavailability of any compound is highly influenced by its physicochemical properties of the compound. According to Lipinski's rule, a poor absorption and a low permeation are more likely when the molecule meets the following properties: the molecular weight of the compound is greater than 500 g/mol, there are more than 10 H-bond donors and more than 5 H-bond acceptors, and the logarithm of octanol-water partition coefficient (Log P) is greater than 5 (Benet et al., 2016) . Whereas the values for lupinifolin lie between the ranges, hesperidin violated molecular weight and number of H-bond acceptors ( Table 2) . The low lipophilicity of hesperidin (Log P ¼ À0.72) and poor water solubility of lupinifolin (Log p ¼ 4.38) implying a low oral bioavailability can be overcome by using oral drug delivery systems or intravenous administration.

The analysis of distribution took into consideration both the permeability across the blood-brain barrier (BBB) and whether the compounds could be substrates of the P-glycoprotein (P-gp) (Prachayasittikul & Prachayasittikul, 2016) . Since the inability of hesperidin and lupinifolin to permeate across the BBB (Table 2) , the central nervous system is protected from their action. Importantly, the distribution of hesperidin could be limited by the P-gp-mediated efflux.

Metabolism of drugs is carried out by the cytochrome P450 (CYP) superfamily. Among the analyzed CYPs (Table 2) , whereas hesperidin might not inhibit any of these enzymes, lupinifolin could inhibit CYP2C19, CYP2C9, and CYP3A4, but not neither CYP1A2 nor CYP2D6. Because of these inhibitions, lupinifolin may affect the metabolism and clearance of the potential drugs resulting into bioaccumulation (Terao & Mukai, 2014) , so optimal drug doses must be defined.

Excretion analysis considered the total clearance (CLtot) because it influences the half-life and bioavailability of drugs. Moreover, CLtot provides a framework for the initial dose for first in human studies (Berellini et al., 2012) . The CLtot of lupinifolin was higher than hesperidin (Table 2) . Neither hesperidin nor lupinifolin was found to be carcinogenic or to present hepatotoxicity. However, hesperidin revealed weak inhibition of the human ether-a-go-go-related gene II (hERG II) potassium channel.

ADMET: Absorption, distribution, metabolism, excretion, toxicity; Log P: logarithm of octanol-water partition coefficient; BBB: Blood-brain barrier; P-gp: P-glycoprotein; Log CLtot: logarithm of total clearance; hERG: human ether-a-gogo-related gene

From the combination of molecular docking to score binding poses, molecular dynamics to evaluate the interaction and stability of protein-ligand complexes, and MM-PBSA to estimate binding energies, among 92 phytochemicals, hesperidin was found to be a promising multitarget antiviral against SAS-CoV-2. More specifically, this flavonoid rutinoside showed high affinities and stability in complex with the 3CLpro and RdRp in their active sites. Structurally, hesperidin is constituted of the flavanone hesperitin (aglycone) and the disaccharide rutinoside (rhamnose linked to glucose), and it is mainly found in citrus fruits. Pharmacological effects of hesperidin include antimicrobial, antiviral, antihyperlipidemic, cardioprotective, antihypertensive, antidiabetic, and antiinflammatory activities (Man et al., 2019; Zanwar et al., 2014) . By performing cell-free and cell-based assays, hesperitin was found to significantly inhibit the cleavage processing of SARS-CoV 3CLpro while being less toxic to Vero cells when compared to other phenolic compounds (Lin et al., 2005) .

Considering that SARS-CoV-2 3CLpro is 96% identical with its SARS-CoV ortholog Zhang et al., 2020) , the estimated binding energies, and the stability with both 3CLpro and RdRp revealed by molecular dynamics simulation, hesperidin becomes a remarkable candidate to act as a multitarget antiviral, inhibiting two essential enzymes of the causative agent of COVID-19.

Lupinifolin was found to be a potential inhibitor of the main protease of SARS-CoV-2 because of its affinity and stability in complex with 3CLpro. This prenylated flavonoid has been identified in several medicinal plants belonging mostly to the family Fabaceae. Other sources include plants in the family Rutaceae and Apocynaceae (Ganapaty et al., 2006; Joycharat et al., 2016; Lin et al., 1991; Mahidol et al., 1997; Smalberger et al., 1974; Soonthornchareonnon et al., 2004) . The wide range of pharmacological effects of lupinifolin includes antimicrobial activities against multidrug-resistant enterococci and Mycobacterium tuberculosis (Sianglum et al., 2019; Soonthornchareonnon et al., 2004) , cytotoxicity against breast cancer, human small-cell lung (NCI-H187) and oral human epidermoid carcinoma (KB), with limited or not determined cytotoxicity against Vero cells (Soonthornchareonnon et al., 2004; Sutthivaiyakit et al., 2009) , and antiviral activity against herpes simplex virus type 1 (HSV-1) (Soonthornchareonnon et al., 2004) . Lupinifolin was found to inhibit bacterial growth via cell membrane disruption, increasing its permeability while decreasing salt tolerance (Sianglum et al., 2019; Yusook et al., 2017) . On the other hand, its antiviral mechanism has not been studied yet.

Protein-interacting residues taken from molecular docking and during the simulation at 100 ns and 200 ns were annotated in order to obtain deeper insight into the interaction pattern. In the 3CLpro complexes, residues forming the catalytic dyad, C145 and H41, were found to interact with hesperidin and lupinifolin (Figure 8 ). In the stable RdRp complex, residues of the catalytic site (759SDD761) as well as the divalent-cation-binding residue D618 interacted with hesperidin ( Figure 9 ).

It is worth mentioning that the oral bioavailability of flavonoids is restricted by the extensive first-pass metabolism in the intestine and liver (Gonzales, 2017) . For instance, flavonoid rutinosides such as hesperidin are metabolized via methylation, demethylation, glucuronidation, sulfation, and sulfoglucuronidation in the intestine and liver, and hydrolyzed by microbial enzymes in the large intestine (Boyle et al., 2000; Nectoux et al., 2019) . Moreover, the transport of prenylated flavonoids to blood circulation is lowered by prenyl groups and a long-term dietary intake can derive in tissue accumulation (Terao & Mukai, 2014) . Hence, the route of administration and delivery are crucial to consider. Moreover, in vitro and in vivo assays are required to evaluate the efficacy of these natural compounds.

Using a virtual screening approach, the present study aimed to elucidate natural compounds as potential inhibitors of the SARS-CoV-2 3CLpro and RdRp which are enzymes that play a key role in the virus replication cycle. After performing molecular docking to estimate the binding affinities and molecular dynamics simulation to evaluate the stability of the protein-ligand complexes, hesperidin (estimated binding affinities, DG bind ¼ À15.18 kcal/mol to 3CLpro and DG bind ¼ À9.46 kcal/mol to RdRp) was identified as a lead candidate because of its potential to efficiently inhibit both 3CLpro and RdRp. Furthermore, lupinifolin exhibited potential to inhibit the main protease showing a higher affinity than hesperidin (DG bind ¼ À20.93 kcal/mol) while displaying suitable pharmacokinetics and physicochemical properties. Based on the safety and low toxicity to normal cells reported in the literature, the predicted ADMET profile, and its affinity and stability when forming complexes with these SARS-CoV-2 enzymes, hesperidin becomes an attractive multitarget drug candidate to treat COVID-19. However, further in vivo/in vitro assays are required to evaluate the proposed therapeutic use.

Preliminary identification of potential vaccine targets for the COVID-19 Coronavirus (SARS-CoV-2) based on SARS-CoV Immunological studies

Coronavirus main proteinase (3CLpro) structure: Basis for design of anti-SARS drugs

Biodiversity patterns and continental insularity in the tropical high andes

BDDCS, the rule of 5 and drugability

In silico prediction of total human plasma clearance

Molecular dynamics with coupling to an external bath

Implementation of the CHARMM force field in GROMACS: Analysis of protein stability effects from correction maps, virtual interaction sites, and water models

Bioavailability and efficiency of rutin as an antioxidant: A human supplementation study

Experimental treatment with favipiravir for COVID-19: An open-label control study

The FDA-approved drug ivermectin inhibits the replication of SARS-CoV-2 in vitro

A trial of lopinavir-ritonavir in adults hospitalized with severe Covid-19

SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Particle mesh Ewald: An NÁlog(N) method for Ewald sums in large systems

Ribavirin, Remdesivir, Sofosbuvir, Galidesivir, and Tenofovir against SARS-CoV-2 RNA dependent RNA polymerase (RdRp): A molecular docking study

Anti-inflammatory activity of Derris scandens

Structure of the RNAdependent RNA polymerase from COVID-19 virus

In Vitro bioavailability and cellular bioactivity studies of flavonoids and flavonoid-rich plant extracts: Questions, considerations and future perspectives

Remdesivir is a direct-acting antiviral that inhibits RNA-dependent RNA polymerase from severe acute respiratory syndrome coronavirus 2 with high potency

SWISS-MODEL and the Swiss-PdbViewer: An environment for comparative protein modeling

The origin, transmission and clinical therapies on coronavirus disease 2019 (COVID-19) outbreak -an update on the status

LINCS: A linear constraint solver for molecular simulations

Effect of dexamethasone in hospitalized patients with COVID-19: Preliminary report

CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data

Sofosbuvir as a potential alternative to treat the SARS-CoV-2 epidemic

Structural insights on Vitamin D receptor and screening of new potent agonist molecules: Structure and ligand-based approach

Inhibition of SARS-CoV 3CL protease by flavonoids

Chemical constituents and biological activities of Albizia myriophylla wood

Targeting SARS-CoV-2: A systematic drug repurposing approach to identify promising inhibitors against 3C-like proteinase and 2 0 -O-ribose methyltransferase

Structure of the SARS-CoV nsp12 polymerase bound to nsp7 and nsp8 co-factors

Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models

g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

Flavonoids as antiviral agents for Enterovirus A71 (EV-A71)

Conformational sampling of flexible ligand-binding protein loops

Anti-SARS coronavirus 3C-like protease effects of Isatis indigotica root and plant-derived phenolic compounds

Three new flavonoids, 3 0 -methoxylupinifolin, laxifolin, and isolaxifolin from the roots of Derris laxiflora Benth

Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations

Molecular dynamics simulations and novel drug discovery

Radius of gyration as an indicator of protein structure compactness

Prenylated flavanones from Derris reticulata

Benefits of hesperidin for cutaneous functions

Comparative analysis based on the spike glycoproteins of SARS-CoV2 isolated from COVID 19 patients of different countries

Antiviral activities of Type I interferons to SARS-CoV-2 infection

Chloroquine inhibits autophagic flux by decreasing autophagosome-lysosome fusion

Coronaviruses: A comparative review

Hydrogen bond dynamic propensity studies for protein binding and drug design

Diversity patterns of selected Andean plant groups correspond to topography and habitat dynamics, not orogeny

Absorption and metabolic behavior of hesperidin (Rutinosylated Hesperetin) after single oral administration to sprague-dawley rats

Open Babel: An open chemical toolbox

Antitumor and antiviral activity of pentacyclic triterpenes. Mini-Reviews in Organic Chemistry

Contrasting plant diversification histories within the Andean biodiversity hotspot

An overview of severe acute respiratory syndrome-coronavirus (SARS-CoV) 3CL protease inhibitors: Peptidomimetics and small molecule Chemotherapy

pkCSM: Predicting small-molecule pharmacokinetic and toxicity properties using graphbased Signatures

P-glycoprotein transporter in drug development

Binding estimation after refinement, a new automated procedure for the refinement and rescoring of docked ligands in virtual screening

Improving MM-GB/SA Scoring through the Application of the Variable Dielectric Model

Herbal medicine as inducers of apoptosis in cancer treatment

Python: A programming language for software integration and development

Anti-HIV effects of chloroquine: Inhibition of viral particle glycosylation and synergism with protease inhibitors

Mechanism of action and biofilm inhibitory activity of lupinifolin against multidrug-resistant enterococcal clinical isolates

Flavonoids from Tephrosia-VII: The constitution and absolute configuration of lupinifolin and lupinifolinol, two flavanones from Tephrosia lupinifolia Burch (DC)

Lupinifolin, a bioactive flavanone from Myriopteron extensum (Wight) K. Schum. stem

Cytotoxic and antimycobacterial prenylated flavonoids from the roots of Eriosema chinense

Implications of protein flexibility for drug discovery

Prenylation modulates the bioavailability and bioaccumulation of dietary flavonoids

CASTp 3.0: Computed atlas of surface topography of proteins

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

GROMACS: Fast, flexible, and free

Automation of the CHARMM General Force Field (CGenFF) I: Bond perception and atom typing

CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

Automation of the CHARMM General Force Field (CGenFF) II: Assignment of bonded parameters and partial atomic charges

Chloroquine is a potent inhibitor of SARS coronavirus infection and spread

Assessing the performance of the MM/PBSA and MM/GBSA methods. 10. Impacts of enhanced sampling and variable dielectric model on protein-protein Interactions

Remdesivir in adults with severe COVID-19: A randomised, double-blind, placebo-controlled, multicentre trial. The Lancet

Coronavirus pathogenesis and the emerging pathogen severe acute respiratory syndrome coronavirus

Lupinifolin from Derris reticulata possesses bactericidal activity on Staphylococcus aureus by disrupting bacterial cell membrane

Chapter 76 -cardiovascular effects of hesperidin: A Flavanone glycoside

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved a-ketoamide inhibitors

The authors thank Dr. Francisco Javier Flores for revising and improving our manuscript. Portions of the computing for this project was performed at the OSU High-Performance Computing Center at Oklahoma State University.

The authors declare no conflict of interest.