Targeting the 3CLpro and RdRp of SARS-CoV-2 with phytochemicals from medicinal plants of the Andean Region: molecular docking and molecular dynamics simulations

Abstract 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, ΔGbind = −15.18 kcal/mol to 3CLpro and ΔGbind = −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 (ΔGbind = −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


Introduction
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.

Library of ligands
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).

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

Molecular dynamics
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.

Binding free energy estimation
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).

Admet prediction
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.

Molecular docking
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.

Molecular dynamics
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
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.

Root mean square fluctuation and DSSP
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
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).

Hydrogen bonds
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).

Binding energy estimation
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).

Admet prediction
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.

Potential inhibitors of the 3CLpro and RdRp of SARS-CoV-2
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.

Conclusions
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.