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

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 non-segmented 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 (COVID-19) 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ácome 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 (Gao et al., 2020). Taken together, these enzymes are attractive drug targets to treat these human CoVs.

Potential COVID-19 treatments include ivermectin, a broad-spectrum 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 (Cai et al., 2020); 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 (Cao et al., 2020). 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 (Wang et al., 2020). 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.

Methodology

Preparation of SARS-CoV-2 protein targets

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.

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 docking

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

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) (Vanommeslaeghe et al., 2012; Vanommeslaeghe & MacKerell, 2012). 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 (ΔGbind) of a protein-ligand complex in solvent is calculated as the free energy difference between the complex (Gcomplex) and the summation of the free energy of the protein (Gprotein) and ligand (Gligand):

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 () 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.

Results and discussion

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 (). The low affinities of chloroquine and hydroxychloroquine to both enzymes () 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)

Table 1.

Theoretical MM-PBSA free energies (kcal/mol) of the stable complexes.

Complex	Van der Waal energy	Electrostatic energy	Polar solvation energy	SASA energy	Binding energy
3CLpro-hesperidin	−38.79 (2.99)	−14.16 (3.23)	42.23 (5.99)	−4.46 (0.31)	−15.18 (4.82)
3CLpro-lupinifolin	−43.44 (3.26)	−2.78 (2.76)	30.31 (4.40)	−5.02 (0.28)	−20.93 (3.45)
RdRp-hesperidin	−31.96 (3.94)	−19.92 (-8.40)	47.20 (11.58)	−4.77 (0.54)	−9.46 (7.23)

Standard deviations are in parentheses.

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 (). 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-α-rhamnoside) showed significant affinities to both viral targets (). In addition, the quercetin glycosides guaijaverin (quercetin-3-O-α-arabinoside) and hyperin (quercetin-3-O-β-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

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.

Figure 1.

Backbone root mean square deviation (RMSD) of the 3CLpro (A) and RdRp (B) of SARS-CoV-2 in complex with phytochemicals, during 200 ns of the molecular dynamics simulation period. Backbones with mean RMSD lower than free protein are plotted.

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.

Figure 2.

Backbone root mean square deviation (RMSD) of phytochemicals from complexes with 3CLpro (A) and RdRp (B), during 200 ns of the molecular dynamics simulation period.

Figure 3.

Snapshot of the superimposed structures of 3CLpro in complex with hesperidin (A), lupinifolin (B), and rutin (C). Structures were obtained from the trajectory file in the interval of 20 ns.

Figure 4.

Snapshot of the superimposed structures of RdRp-hesperidin complex. Structures were obtained from the trajectory file in the interval of 20 ns.

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 β-sheet content remained stable (Figure S3, Figure S4).

Figure 5.

Root mean square fluctuation (RMSF) of backbone atoms of 3CLpro (A) and RdRp (B) as free proteins and forming complexes with phytochemicals, during the simulation period.

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

Figure 6.

Backbone radius of gyration of 3CLpro (A) and RdRp (B) as free proteins and forming complexes with phytochemicals, during 200 ns of the molecular dynamics simulation period.

Hydrogen bonds

It is well known that H-bonds are responsible for the secondary and tertiary structural protein motifs. Formation of H-bonds 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).

Figure 7.

Total number of hydrogen bonds interactions between phytochemicals (hesperidin and lupinifolin) and proteins 3CLpro (A) and RdRp (B).

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 (ΔGbind) 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 = 2 was set for the MM-PBSA calculations of these complexes. On the other hand, a value of = 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).

Figure 8.

Molecular interactions of hesperidin (blue) and lupinifolin (green) in complex with the 3CLpro of SARS-CoV-2 taken from the molecular docking (A,D) and during the simulation at 100 ns (B,E) and 200 ns (C,F).

Figure 9.

Molecular interactions of hesperidin in complex with the RdRp of SARS-CoV-2 taken from molecular docking (A) and during the simulation at 100 ns (B) and 200 ns (C).

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 ΔGbind = −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).

Figure 10.

MM-PBSA binding free energy contribution per residue of 3CLpro in complex with hesperidin and lupinifolin.

Figure 11.

MM-PBSA binding free energy contribution per residue of RdRp in complex with hesperidin.

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.

Table 2.

Predicted ADMET profile of hesperidin and lupinifolin.

Descriptor	Hesperidin	Lupinifolin
Molecular weight (g/mol)	610.56	406.47
Number of rotatable bonds	7	3
Number of H-bond donors	8	2
Number of H-bond acceptors	15	5
Log P	−0.72	4.38
Gastrointestinal absorption	Low	High
BBB permeability	No	No
Substrate of P-gp	Yes	No
CYP1A2 inhibitor	No	No
CYP2C19 inhibitor	No	Yes
CYP2C9 inhibitor	No	Yes
CYP2D6 inhibitor	No	No
CYP3A4 inhibitor	No	Yes
Bioavailability score	0.17	0.55
Log CLtot (mL/min/kg)	0.211	0.31
Hepatotoxicity	No	No
Ames toxicity	No	No
hERG I channel inhibition	No	No
hERG II channel inhibition	Yes	No

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-go-go-related gene

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 anti-inflammatory 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 (Gao et al., 2020; 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, ΔGbind = −15.18 kcal/mol to 3CLpro and ΔGbind = −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 (ΔGbind = −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.

Click here for additional data file.