Virtual screening of plant-derived compounds against SARS-CoV-2 viral proteins using computational tools.

The new SARS-CoV-2, responsible for the COVID-19 pandemic, has been threatening public health worldwide for more than a year. The aim of this work was to evaluate compounds of natural origin, mainly from medicinal plants, as potential SARS-CoV-2 inhibitors through docking studies. The viral spike (S) glycoprotein and the main protease Mpro, involved in the recognition of virus by host cells and in viral replication, respectively, were the main molecular targets in this study. Molecular docking was performed using AutoDock, which allowed us to select the plant actives with the highest affinity towards the viral targets and to identify the interaction molecular sites with the SARS-CoV2 proteins. The best energy binding values for S protein were, in kcal/mol: -19.22 for glycyrrhizin, -17.84 for gitoxin, -12.05 for dicumarol, -10.75 for diosgenin, and -8.12 for delphinidin. For Mpro were, in kcal/mol: -9.36 for spirostan, -8.75 for N-(3-acetylglycyrrhetinoyl)-2-amino-propanol, -8.41 for α-amyrin, -8.35 for oleanane, -8.11 for taraxasterol, and -8.03 for glycyrrhetinic acid. In addition, the synthetic drugs umifenovir, chloroquine, and hydroxychloroquine were used as controls for S protein, while atazanavir and nelfinavir were used for Mpro. Key hydrogen bonds and hydrophobic interactions between natural compounds and the respective viral proteins were identified, allowing us to explain the great affinity obtained in those compounds with the lowest binding energies. These results suggest that these natural compounds could potentially be useful as drugs to be experimentally evaluated against COVID-19.

Introduction

New strains of viruses are constantly emerging and represent a challenge to public health. Among the recent emerging viruses were the Severe Acute Respiratory Syndrome coronavirus (SARS-CoV) in 2003 influenza A (H1N1) virus in 2009, Middle East Respiratory Syndrome coronavirus (MERS-CoV) in 2012, and Severe Acute Respiratory Syndrome coronavirus 2 (SARS-CoV-2), which causes the 2019-nCoV infection (COVID-19) that appeared at the end of 2019 in China (Gorbalenya et al., 2020). The disease started in Wuhan City, central China, mainly linked to workers from the Wuhan South China Wholesale Seafood Market, which sells, among other products, various types of live and dead exotic animals. After 14 months of the outbreak, the global acceleration in case incidence has slowed down, with around 3 million new cases reported; however, death rates continue to increase with over 88,000 new deaths reported last week (World Health Organization, 2020). As of February 9, there have been almost 105 million cases and over 2.3 million deaths reported globally since the beginning of the pandemic. These data demonstrate the urgent need to find quick solutions to the problem.

To date, very vaccines against SARS-CoV-2 have been developed showing good percentages of effectiveness in clinical trials; however, the world population has not been massively vaccinated and the period of immunity that these vaccines will generate is yet unknown (Detoc et al., 2020; Kaur & Gupta, 2020; Malik et al., 2020). On the other hand, it is important to mention that therapies with monoclonal antibodies proved to be effective as modulators of the immune response against SARS-CoV-2 (Jahanshahlu & Rezaei, 2020; Pecetta et al., 2020; Zhao, 2020). Some monoclonal antibodies such as infliximab are already in the stage of clinical trials in patients with COVID-19 (Clinical (Clinical Trials, 2020)). Also, many synthetic antiviral molecules, already existing in the market were evaluated at different levels, many of which reached the stage of clinical trials. However, none of them showed to be effective enough against SARS-CoV-2 (Clinical (Clinical Trials, 2020; Frediansyah et al., 2021; Whang et al., 2020)). In this context, the search for new therapeutic agents is imperative as another possible alternative for the treatment of COVID-19; indeed, great effort and important resources are nowadays being dedicated worldwide to achieve that goal. Antivirals are among the compounds used to prevent the development of viruses, useful in restricting viral infection. Of special interest are the secondary metabolites of plant origin since many of them not only act as powerful antivirals but also have the ability to stimulate the immune system, which is essential in a response to SARS-CoV-2 infection (Cinatl et al., 2003; Luengo-Fernandez et al., 2008; Rawat et al., 2020).

Compounds known for their effectiveness against other viruses are being investigated in search of a solution to the rapid and lethal advance of the new SARS-CoV-2. Examples of that are some protease inhibitors used for treatment of patients with Human Immunodeficiency Virus infection and acquired immune deficiency syndrome (HIV/AIDS) or neuraminidase inhibitors or nucleotide analogues like Remdesivir (GS-5734™) used for other pathogenic CoVs (Khaerunnisa et al., 2020). Another interesting molecule is chloroquine/hydroxychloroquine, an acidotropic amine form of quinine, extensively used against malaria and previously described as a potent inhibitor of most coronaviruses in vitro, including SARS-CoV-1 (Devaux et al., 2020; White et al., 2014). It is interesting to note that quinone is a metabolite produced by cinchona trees (Cinchona officinalis), native from Peru. Plants are great producers of secondary metabolites and as such, especially those from poorly studied regions, they captured the interest of many researchers trying to develop new medicines (Yang et al., 2018). Antiviral bioactivities have been described in numerous medicinal plants and associated with compounds like flavonoids, heterosides, terpenes, and triterpenes, organic acids, alkaloids, saponins, and quaternary ammonium salts, among others (Cecil et al., 2011; Khaerunnisa et al., 2020; Roy & Saraf, 2006; Wang et al., 2017).

Molecular docking is a useful tool in structural molecular biology to evaluate the possibility of interactions and to predict binding modes of ligand-proteins from knowledge of the three-dimensional structure. The search is conducted in a way that all the candidate dockings are ranked through the Gibbs free energy (negative values indicate spontaneous chemical reactions; the greater the magnitude, the higher is the affinity of the ligand-macromolecule interaction), allowing us to make hypotheses or decisions about how, or which, ligands inhibit the target. This tool can play a key role in computer-assisted drug design or in the selection of the best candidates from a library of compounds with a certain activity (Morris & Lim-Wilby, 2008). Analyzing these virtual protein-ligand interactions becomes critical when rushing in the search of metabolites that can cure or ease symptoms of diseases when time and resources for extensive experimental work are scarce. Considering that in the last months many SARS-CoV-2 macromolecules have been crystallized, more effective simulations may be performed cross-out if model proteins from other closely related viruses have to be used (Jin et al., 2020; Wrapp et al., 2020; Zhang et al., 2020).

In the initial stage of the SARS-CoV-2 replication cycle, the attachment of the virion to the host cell is initiated by interactions between the S protein (spike) and human angiotensin-converting enzyme (ACE2 receptor) (Hoffmann et al., 2020). The S protein-receptor interaction is the primary determinant for a coronavirus to infect a host species and governs the tissue tropism of the virus. Structurally, the coronavirus spicule is a trimer arranged asymmetrically. Each monomeric unit has three segments: an extensive ectodomain, a single-pass transmembrane segment, and a small intracellular tail (Walls et al., 2016). The ectodomain contains two subunits; the S1 subunit allows binding to the ACE2 receptor while the S2 allows the fusion of viral and host membranes. Two major domains in S1, N-terminal domain (S1-NTD) and C-terminal domain (S1-CTD), have been identified. The S1-CTD, also known as the receptor-binding domain (RBD), is responsible for recognizing protein receptors ACE2 (Liu et al., 2015).

The next step in the infection cycle is the translation of the replicase gene from the viral genome (positive polarity RNA). Replicase gene encodes two overlapping polyproteins, pp1a and pp1ab, required for viral replication and transcription (Baranov et al., 2005). Polyproteins pp1a and pp1ab contain the non-structural proteins (nsps) 1–11 and 1–16, respectively. These polyproteins are subsequently cleaved into the individual nsps (Ziebuhr et al., 2000). SARS-CoV-2 encodes two proteases that cleave the polyproteins; the papain-like proteases, PLpro (encoded by nsp3) and chymotrypsin-like cysteine protease known as 3C-like protease (3CLpro). The 3CLpro enzyme, also called main protease (Mpro), hydrolyses the polyprotein at no less than 11 conserved sites between nsp4 and nsp16, starting with the autolytic cleavage of this enzyme itself (Hegyi & Ziebuhr, 2002).

The aims of this work were to evaluate by molecular docking selected compounds obtained from natural sources for antiviral activity against SARS-CoV-2. Synthetic drugs that are currently being evaluated for the treatment of COVID-19 were used as controls. Two viral proteins with different and essential functions were tested: the SARS-CoV-2 main protease (Mpro) and the S protein (spike) in its prefusion conformation.

Materials and methods

Viral proteins

SARS-CoV-2 Mpro (PDB ID: 6LU7) and S protein (PDB ID: 6VSB) were used in this study. Three-dimensional structure are shown in Table 1 .

Table 1

SARS-CoV-2 main protease (PDB: 6LU7) and S protein (PDB: 6VSB) structures. Active sites for main protease and interaction domains with the ACE2 receptor for S protein.

PDB ID	Macromolecule	Target-ligand binding sites
6LU7	Image 1	Active sitesTHR24aTHR26aPHE140aASN142aGLY143aCYS145aHIS163aHIS164aGLU166aHIS172a
6VSBCHAIN A (blue)CHAIN B (red)CHAIN C (grey)	Image 2	Receptor binding domain (RBD)bresidues340 to 521 (light blue) of chain B

Active site residues

Recognition sites.

Antiviral compounds

Secondary metabolites with different structures were analyzed to select the most promising to assess a possible relationship between structure and activity. Many of them already had proven antiviral activity against other viruses (Alonso & Desmarchelier, 2015; Cecil et al., 2011; Enkhtaivan et al., 2018; Shin et al., 2015).

Sixteen plant-derived compounds were assessed for the interaction with Mpro: α-amyrin (CID 73170), berberine (CID 2353), geranial (CID 638011), glycyrrhetinic acid (CID 10114), lapachenole (CID 174859), lapachol (CID 3884), limonene (CID 22311), maytansine (CID 5281828), meliacine (Andrei et al., 1990; Hassan et al., 2015), menthol (CID 1254), rutin (CID 5280805), oleanane (CID 9548717), sambunigrin (CID 91434), sauroxine (CID 5462447), spirostan (CID 6857439), taraxasterol (CID 115250), and theophylline (CID 2153). The semisynthetic compound analyzed was the N-(3-acetylglycyrrhetinoyl)-2-amino-propanol (Zígolo et al., 2018). The synthetic nelfinavir (CID 64143) and atazanavir (CID 148192), were used as controls for the interaction with Mpro.

Regarding the interaction with the S protein (PDB: 6VSB), seven natural compounds: aloin (CID 9866696), diosgenin (CID 99474), amygdalin (CID 34751), gitoxin (CID 91540), glycyrrhizin (CID 14982), delphinidin (CID 68245), and dicumarol (CID 54676038), and also three synthetic compounds: chloroquine (CID 2719), hydroxychloroquine (CID 3652), and umifenovir (CID 131411), were assessed.

Docking analysis

Docking calculations were carried out with Autodock 4.2 program, which uses Lamarckian genetic algorithm (LGA) for calculations (Morris et al., 2009). Number of genetic algorithm (GA) runs was set to 200 for each case analyzed and binding affinities were determined and reported in kcal/mol.

The Autodock 4.2 program was applied considering all rotatable bonds for ligands and the whole protein as a rigid structure. For the location and extent of the 3D area the search space was defined by specifying a center, the number of points at each dimension and the points between spaces to focus the search space in the enzyme active site or putative binding region.

For protease, the grid box center coordinates corresponded to the Cα atom of the histidine residue 163: x (−19.055), y (16.637), and z (64.100). Dimensions of the grid box were 90 × 90 × 90 and point spacing was 0.492 Å. For the S protein, the grid box center coordinates corresponded to the oxygen atom of the alanine 419 residue for chain B: x (232.020), y (245.760), and z (265.659). Dimensions of the grid box were 124 × 124 × 124 and point spacing was 0.869 Å.

Prior to docking Mpro, the inhibitor N3 and water molecules were removed from the protein structure. Polar hydrogen atoms were added, and Gasteiger atom charges were assigned to protein atoms. Other parameters were set to default values.

The 200 conformers found for each compound with the Autodock program were grouped in clusters that were ordered according to a ranking, which was determined by the stability of the enzyme-ligand complexes within each cluster. Clustering of the 200 conformers was done according to the similarity with the conformation adopted inside the enzyme. The criterion used to evaluate such similarity was the residue mean quadratic square root deviation (RMSD). The residues were obtained through the difference between the atom coordinates of a given conformer respect to the cluster to which the most stable conformer belongs. Each cluster grouped conformers with RMSD lower or equal to 2.0 Å.

For visualization of protein-ligand complexes the software Visual Molecular Dynamics 1.9.1 (VMD) was used (Theoretical and Computational Biophysics Group, University of Illinois).

Data processing using data mining and statistics

The frequency distribution (number in cluster) of plant-derived compounds were first plotted vs. Gibbs free energy of conformers for a specific target or protein. Then, the best compounds for each protein were selected considering the following criteria for each compound: i) clusters with highest frequencies, ii) clusters with lowest Gibbs free energies, iii) conformational dispersion associated with the number of clusters (the higher the number of clusters was, there was more conformational dispersion for the interaction studied), iv) intra-cluster dispersion given by the variation around the mean binding energy, and v) intra-cluster dispersion given by the percentile of energy values close to the minimum binding energy (more negative free energy). Statistical software R (https://cran.r-project.org/) with the Rstudio IDE (https://rstudio.com/) was used for processing.

Different clusters were ordered according to the energy of the most stable conformer within each group: the ranking was ordered according to the increasing stability of the lowest energy conformers in each cluster.

Results

Selection of conformers

From the 200 docking runs for all the evaluated compounds, both from natural origin and the synthetic controls, the binding energies for each conformer were obtained which, in turn, were grouped into clusters. The binding energy for each conformer was plotted as a function of the frequencies (number of conformers in each cluster) for the Mpro protein (Fig. 1, A and B; S1A) and for S glycoprotein (Fig. 2; Fig. S2).

Fig. 1

A) Cluster analysis of conformations. Frequency (number of conformers in clusters) of each compound after 200 runs considering the binding energy of conformations corresponding to the docking of α-amyrin, berberine, geranial, glycyrrhetinic acid, lapachenole, lapachol, limonene, meitansine, meliacine, menthol, N-(3-acetylglycyrrhetinoyl)-2-amino-propanol, nelfinavir, oleanane, sambunigrin, sauroxine, spirostan, taraxasterol, and theophylline with SARS-CoV-2 main protease (Mpro). B) Frequency (number of conformers in clusters) and binding energy of the most promising compounds regarding affinity with SARS-CoV-2 main protease (Mpro): α-amyrin, N-(3-acetylglycyrrhetinoyl)-2-amino-propanol, oleanane, spirostan, and nelfinavir and atazanavir as controls.

Fig. 2

Cluster analysis of conformations. Frequency (number of conformers in clusters) of each compound after 200 runs considering the binding energy of conformations corresponding to the docking of aloin, diosgenin, amygdalin, gitoxin, chloroquine, glycyrrhizin, delphitinidin, dicumarol, hydroxichloroquine and umifenovir with SARS-CoV-2 S protein.

The interaction of all the natural compounds evaluated with each viral protein showed, in general, a low number of clusters although they were highly frequent, suggesting that there was a high probability for them to happen. Conversely, the synthetic antiviral nelfinavir and atazanavir used as controls for the Mpro protein and umifenovir, chloroquine, and hydroxychloroquine for the spike protein, showed broad conformational dispersion or high number of clusters (Figs. 1B and 2, S1A).

Each compound was represented by the most populated clusters, which in some cases were also those with the lowest binding energies. Then, statistical criteria were applied to select the best clusters when more than one was considered for one compound, and the binding energies associated with their interactions with the viral proteins. Three different statistics were calculated: a) the standard deviation of the Gibbs free energy of all the cluster's conformers (SD binding energy), b) the variation coefficient of the Gibbs free energy of all the cluster's conformers (CV), and c) the intra-cluster dispersion given by the 1% percentile of the Gibbs free energy close to the lowest binding energy (Lowest 1% perc). These metrics were applied in order to assign to each compound the Gibbs free energy corresponding to the interaction with the viral proteins based on the following criteria: i) the Mean binding energy was assigned when CV was lower than or equal to 10% and the value of the Lowest 1% perc was not higher than 60%; ii) the Lowest binding energy was assigned when the Lowest 1% perc reached at least 60%, despite the CV value was lower than or equal to 10% (Table 2 for Mpro and Table 3 for S protein).

Table 2

Compounds evaluated by docking with Mpro. Ranking; frequency; lowest and mean and standard deviation (SD) binding energies (ΔG); variation coefficient (CV); Lowest 1% ΔG and Lowest 1% perc, for selected clusters. Synthetic antivirals used as controls are in bold. Lowest 1% perc is an intra-cluster dispersion given by the 1% percentile of the Gibbs free energy close to the lowest binding energy.

Compound	Ranking	Frequency (Num in clus)	Lowest ΔG (kcal/mol)	Mean ΔG (kcal/mol)	SD ΔG (kcal/mol)	CV (%)	Lowest 1% ΔG (kcal/mol)	Lowest 1% perc
α-Amyrin	1	125	−8.41	−8.40	0.05	0.65	−8.33	96
Atazanavir	2	13	−10.56	−9.47	0.91	9.60	−10.45	15
Berberine	1	66	−7.35	−7.29	0.04	0 0.49	−7.28	80
Geranial	1	108	−5.60	−5.20	0.19	3.75	−5.54	5
Glycyrrhetinic acid	1	46	−8.15	−8.03	0.07	0.91	−8.07	43
Lapachenole	1	79	−7.60	−7.45	0.09	1.27	−7.52	18
Lapachenole	3	77	−6.37	−6.35	0.03	0.50	−6.31	96
Lapachol	3	97	−6.91	−6.71	0.12	1.83	−6.84	12
Lapachol	1	61	−8.17	−7.65	0.43	5.60	−8.09	8
Limonene	1	141	−6.31	−6.21	0.11	1.84	−6.25	55
Meitansine	7	18	−6.23	−5.46	0.40	7.24	−6.17	6
Meliacine	5	39	−7.36	−7.05	0.10	1.27	−7.29	3
Menthol	4	64	−5.38	−5.19	0.13	2.55	−5.33	20
Menthol	1	47	−6.55	−6.33	0.35	5.54	−6.48	66
N-(3-Acetylglycyrrhetinoyl)-2-amino-propanol	19	14	−7.58	−7.38	0.22	2.97	−7.50	43
N-(3-Acetylglycyrrhetinoyl)-2-amino-propanol	3	10	−8.37	−7.68	0.54	7.06	−8.29	10
N-(3-Acetylglycyrrhetinoyl)-2-amino-propanol	1	7	−8.75	−8.23	0.43	5.27	−8.66	29
Nelfinavir	1	6	−9.66	−8.25	0.85	10.26	−9.56	17
Nelfinavir	4	5	−7.83	−7.18	0.41	5.66	−7.75	20
Oleanane	2	59	−8.35	−8.34	0.01	0.10	−8.27	100
Rutin	6	2	−4.65	−4.30	0.49	11.51	−4.60	50
Sambunigrin	6	20	−5.35	−4.59	0.48	10.54	−5.30	10
Sauroxine	6	65	−6.06	−6.05	0.01	0.20	−5.99	98
Spirostan	4	71	−8.50	−8.50	0.01	0.04	−8.42	100
Spirostan	1	59	−9.36	−9.35	0.01	0.08	−9.27	100
Taraxasterol	6	72	−8.21	−8.11	0.04	0.59	−8.13	38
Theophylline	5	56	−4.49	−4.48	0.01	0.11	−4.45	100
Theophylline	1	32	−5.04	−5.02	0.02	0.32	−4.98	97

Table 3

Compounds evaluated by docking with the spike (S) glycoprotein. Ranking; frequency; lowest and mean and standard deviation (SD) binding energies (ΔG); variation coefficient (CV); Lowest 1% ΔG and Lowest 1% perc, for selected clusters. Lowest 1% perc is an intra-cluster dispersion given by the 1% percentile of the Gibbs free energy close to the lowest binding energy.

Compound	Ranking	Frequency (Num in clus)	Lowest ΔG (kcal/mol)	Mean ΔG (kcal/mol)	SD ΔG (kcal/mol)	CV(%)	Lowest 1% ΔG (kcal/mol)	Lowest 1% perc
Aloin	2	5	−6.28	−5.94	0.29	4.86	−6.22	20
Amygdalin	1	2	−7.37	−6.52	12.02	18.44	−7.30	50
Chloroquine	30	8	−5.96	−5.41	0.42	7.81	−5.90	12
Delphinidin	1	16	−8.12	−7.84	0.21	2.67	−8.04	25
Dicumarol	1	41	−12.05	−11.29	0.50	4.43	−11.93	7
Diosgenin	1	58	−10.75	−10.64	0.27	2.51	−10.64	72
Gitoxin	1	32	−17.84	−17.33	0.50	2.87	−17.66	22
Glycyrrhizin	1	64	−19.22	−19.00	0.42	2.22	−19.03	81
Hydroxy-chloroquine	1	3	−5.30	−4.64	0.79	17.08	−5.25	33
Umifenovir	22	4	−6.15	−5.36	0.85	15.83	−6.09	25
Umifenovir	1	3	−7.47	−6.95	0.47	6.83	−7.40	33

Synthetic antivirals used as controls are in bold.

When two or more clusters were initially selected for the same compound (for having similar frequencies and binding energies), the following criteria were applied: i) if the clusters compared had Lowest 1% perc value equal to or greater than 60%, the one with the highest Lowest 1% perc value was chosen, and ii) when the clusters compared had a Lowest 1% perc value lower than 60%, the CVs were compared and that with the lowest CV value was chosen.

Following the statistical criteria for the docking results, the binding energies for the interactions of each compound with the corresponding protein were selected (Table 4 ).

Table 4

Gibbs free energy (ΔG, in kcal/mol) calculated for the interaction between SARS-CoV-2 different molecular targets Mpro and S protein and selected active compounds.

Compound	ΔG (kcal/mol) Mpro	Compound	ΔG (kcal/mol) S protein
Berberine	−7.35	Aloin	−5.94
α-amyrin	−8.41	Amygdalin	−7.37
Menthol	−6.55	Delphinidin	−8.12
Limonene	−6.21	Dicumarol	−12.05
N-(3-Acetylglycyrrhetinoyl)-2-amino-propanol	−8.75	Diosgenin	−10.75
Glycyrrhetinic acid	−8.03	Gitoxin	−17.84
Lapachol	−6.71	Glycyrrhizin	−19.22
Taraxasterol	−8.11	Umifenovir	−7.47
Spirostan	−9.36	Chloroquine	−7.25
Oleanane	−8.35	Hydroxychloroquine	−5.30
Lapachenole	−7.45
Maytansine	−5.46
Meliacine	−7.05
Rutin	−4.03
Geranial	−5.20
Sambunigrin	−4.59
Sauroxine	−6.06
Theophylline	−4.49
Atazanavir	−9.47
Nelfinavir	−8.25

Synthetic antivirals used as controls are in bold.

When Mpro protein interacted with natural compounds, similar values than those observed with the synthetic antivirals atazanavir and nelfinavir, were found. However, for synthetics control compounds greater dispersion and lower frequencies were observed (Fig. 1A and B; S1A). The strongest affinities for spirostan and for N-(3-acetylglycyrrhetinoyl)-2-amino-propanol, with ΔG of −9.36 and − 8.75 kcal/mol, were obtained, respectively. Other compounds showed interesting affinities for Mpro, for example: α-amyrin, oleanane, taraxasterol, and glycyrrhetinic acid (Table 4).

On the other hand, the plant-derived natural compounds evaluated for S protein, yielded very high affinity values compared to controls, except for aloin (−5.94 kcal/mol). The highest affinities were observed for glycyrrhizin, gitoxin, dicumarol, diosgenin, and delphinidin, which also showed low dispersion and high frequencies (Fig. 2; Fig. S2). The most interesting results were observed for glycyrrhizin and gitoxin (−19.22 and −17.84 kcal/mol, respectively), which presented more than twice the affinity for the viral S protein compared to the synthetic antivirals umifenovir, chloroquine, and hydroxychloroquine (−7.47, −7.25, and −5.30 kcal/mol, respectively) (Table 4).

Analysis of the molecular interactions of candidate compounds with the proteins evaluated

Key interactions of possible candidate drugs with Mpro were analyzed: α-amyrin, spirostan, taraxasterol, oleanane, lapachol, glycyrrhetinic acid, and N-(3-acetylglycyrrhetinoyl)-2-amino-propanol (Table 5 ). All the natural and semisynthetic active compounds bound to the enzyme active site, as well as the controls nelfinavir and atazanavir. Only interactions of nelfinavir with Mpro with lower binding energy were shown (Table 5). The semi-synthetic compound N-(3-acetylglycyrrhetinoyl) -2-amino-propanol presented five hydrogen bond interactions, four with the residues of the active site and one with the tyrosine residue 54. It also showed hydrophobic interactions with six residues of the protease in the vicinity of the catalytic site. All these interactions increased the binding affinity with respect to the control compound, which presented only two hydrogen bond interactions with the active site residues and hydrophobic interactions with other four protease residues. Glycyrrhetinic acid, a precursor of the above-mentioned semisynthetic compound, had four hydrogen bond interactions with catalytic residues and hydrophobic interactions with leucine residue 4.

Table 5

Three-dimensional (3D) molecular interaction of candidate compounds with Mpro.

Active site residues are represented with balls and sticks: THR24 (blue), THR26 (violet), PHE140 (pink), ASN142 (brown), GLY143 (white), CYS145 (yellow), HIS163 (orange), HIS164 (red), GLU166 (grey), and HIS172 (green).

It is important to note that five of the six candidate compounds (Table 5) have in common pentacyclic triterpene structures: taraxasterol, N-(3-acetylglycyrrhetinoyl) -2-amino-propanol, glycyrrhetinic acid, α-amyrin, and oleanane. All of them but oleanane presented hydrogen bond interactions; however, its high affinity, comparable with the compounds that showed hydrogen bond interactions highlighted the importance of hydrophobic interactions that could help stabilize oleanane binding with the main protease.

Spirostan interacted with the Mpro mainly by hydrophobic interactions, presenting only one hydrogen bond interaction with the serine 46 catalytic residue. Similarly, taraxasterol only had one hydrogen bond interaction.

In addition to glycyrrhizin and gitoxin, the interaction of dicumarol with the spike protein showed one of the most favorable Gibbs free energies (−12.05 kcal/mol) (Table 6 ). Glycyrrhizin, with the highest affinity for the glycoprotein, interacted with a great number of residues from A and C chains. There was a total of four hydrogen bonds with residues of the A chain (PHE855, VAL976 and ASP979), and another two with the B chain (ASP571 and ASP574). Hydrophobic interactions occurred with the remaining of the 14 residues.

Table 6

Three-dimensional (3D) molecular interaction of candidate compounds with S protein.

The S protein residues that interacted by hydrogen bonds are represented with lines and with balls and sticks those that interacted hydrophobically with the compounds.

Discussion

Plants have a long evolutionary history of developing resistance against viruses. This adaptation, together with the ability to produce secondary metabolites, has generated a wide range of possible sources for antiviral drugs. Herbal infusions and essential oils from plants have been used ancestrally in different civilizations. Their low toxicity and minimal adverse effects at the doses used have long been known, representing an important advantage as compared with new synthetic drugs, which need to be evaluated for safety through strict clinical trials (Abbas Zaidi et al., 2016; Casais de Corne et al., 1997; Horden, 2016). Other advantages of natural vs. synthetic (pure) compounds are that: i) they are mixtures of active compounds named phytocomplexes, formed by isomers and analogue compounds, with the capacity for synergistic activity, and ii) the mixture usually contains compounds that help mitigate the toxicity of the main active drugs. Also, microorganisms are less likely to develop resistance to a phytocomplex than to a unique compound in a synthetic drug (Alonso, 2008; Bruneton, 2001; Camponovo & Bandoni, 1955). To simplify the analyses in this study, only the main active compounds were used since the chemical structures are known.

Liu and coworkers (2015) successfully crystallized the SARS-CoV-2 main Mpro protease, which is a potential drug target for the inhibition of SARS-CoV-2 as it plays an essential role in viral replication and transcription (Liu et al., 2015). From the docking of Mpro with natural compounds, the most favorable binding energies were obtained for compounds with steroid and pentacyclic triterpene structures. However, it is important to note that compounds with different structures extracted from plants were evaluated. Thereby alkaloids (berberine, theophylline, sauroxin), terpenoids (limonene, geranial), phenols (lapachol, lapachenol), alcohols (menthol), steroids (spirostan), flavonoids (rutin), glycosides (sambunigrin), pentacyclic triterpenes (amirin, taraxasterol, oleanaane, glycyrrhetinic acid, N-(3-acetylglycyrrhetinoyl)-2-amino-propanol) and macrolides (maytansine), were assayed.

Berberine is an isoquinoline alkaloid extracted from plants of the Berberis family (Coptis chinensis Berberis buxifolia, Berberis vulgaris, Argemone mexicana, etc.). It has antimalarial, antibacterial, and antiviral effects, among other biological activities (Wang et al., 2017). Berberine can also inhibit the replication of influenza A in vitro in several different cell types and with two influenza A (H1N1) strains (Cecil et al., 2011). As in this study, other researchers found that berberine displayed binding energy of −7.91 kcal/mol and had the ability to interact with the protease through the amino acid ASP187 producing different conformational changes in the target enzyme (Agrawal et al., 2020).

In particular, the structural family of pentacyclic triterpenes, with high affinity values for viral protease, could be suggested as possible candidates for the treatment of COVID-19. Taraxasterol is a pentacyclic triterpene only found in a few plants like arnica, chicory, and dandelion, but with remarkably interesting biological activity (Krajcovicova et al., 2004; Zígolo et al., 2018). Further, it would be also important to evaluate more compounds of plant origin with common steroid structures due to the high value found in this work for spirostan. Other compounds in this family with steroidal structures and potent antiviral activity have already been reported against poliovirus, measles virus, Influenza virus, herpes simplex virus, human Immunodeficiency virus, and other viral diseases (Kaliyaperumal et al., 2020).

Hydrophobic interactions played a particularly important role for the Mpro protein. In the case of oleanane, which showed a good binding energy value, the only type of interactions that presented were those of hydrophobic nature with the protease residues. Furthermore, for all the other compounds evaluated, despite forming hydrogen bonds there were also a significant number of hydrophobic interactions with protease residues. Khaerunnisa et al. (Khaerunnisa et al., 2020) also showed the importance of hydrophobic interactions for a series of compounds of plant origin in the interaction with the viral protease Mpro4. It is known that spirostan is the most abundant steroidal saponin, being the triterpene and steroid glycosides of important biological and medicinal applications (Pérez-Labrada et al., 2012). Several of these molecules have exhibited potent antiviral activities, making them promising candidates for drug development according to the results of this work and other studies (Cecil et al., 2011; Enkhtaivan et al., 2018; Shin et al., 2015).

In the molecular docking for Mpro protein of SARS-Cov-2, the semi-synthetic compound N-(3-acetylglycyrrhetinoyl)-2-amino-propanol performed better than glycyrrhizic acid. The former, obtained by biocatalytic reaction of aminolysis from glycyrrhetinic acid (derived from the plant Glycyrrhiza glabra L), already showed powerful antiviral activity against the herpes virus types I and II, including strains resistant to acyclovir (De Clercq, 2000; Zígolo et al., 2018). In addition, Zígolo et al. (Zígolo et al., 2018) showed low toxicity of the N-(3-acetylglycyrrhetinoyl)-2-amino-propanol in Vero cell cultures, making this compound a promising candidate for COVID-19 treatment27. In a molecular docking study of natural products that are commonly used in Indian and other cuisines against SARS-CoV-2 main protease (PDB: 6LU7 and PDB: 6Y2E), it was shown that several of these compounds had high binding affinity against Mpro, and they were comparable with a known anti-HIV drug, saquinavir. In our study, we report plant-derived active compounds with much higher binding affinities than other natural compounds reported as possible antivirals for SARS-CoV-2. Glycyrrhizin (binding energy of −8.40 kcal/mol) was one of the best potential candidates assessed. It is remarkable that a diammonium salt of glycyrrhizin was recommended for the cure of COVID-19 in China. The binding affinity of the positive control, saquinavir, was −7.90 kcal/mol (Sampangi-Ramaiah et al., 2020).

In a recent study, several compounds, mainly flavonoids (kaempferol, quercetin, luteolin-7-glucoside, demethoxycurcumin, naringenin, apigenin-7-glucoside, oleuropein, curcumin, catechin, epicatechingallate, zingerol, gingerol, and allicin) and the synthetic compounds nelfinavir and lopinavir were analyzed as potential inhibitors of the SARS-CoV-2 Mpro (PDB: 6LU7) by docking studies (Khaerunnisa et al., 2020). In our hands, the binding energy found (−8.25 kcal/mol) for the interaction of nelfinavir with Mpro protein was similar to that reported by Khaerunnisa et al. (Khaerunnisa et al., 2020) (−8.37 kcal/mol). Interestingly, despite the similarity in the ΔG obtained in both studies, differences in the hydrogen bonds formed were observed (Khaerunnisa et al., 2020). While residues Glu166, Gln189, and Gln192 were used according to Khaerunnisa et al. (Khaerunnisa et al., 2020), we only found Glu166 and His172 in the interaction between nelfinavir and Mpro protein by docking (Khaerunnisa et al., 2020). Differences observed would be consequences of the parameters considered by different authors, for example, the rotatable links in the coupling. Furthermore, this difference may also be due to the selection of different interaction conformations made by the authors. In this work, rigorous statistical criteria to make the best decision in choosing the correct cluster (detailed in Section 3.1), were applied. In this sense, nelfinavir presented a broad conformational dispersion (number of clusters) and similar frequencies in the clusters with better binding energies.

All of the compounds evaluated against S glycoprotein showed very favorable binding energies, higher than those for the synthetic antiviral compounds, umifenovir (−7.47 kcal/mol), chloroquine (−7.25 kcal/mol), and hydroxychloroquine (−5.30 kcal/mol), directed to that viral target. The most favorable binding energies were obtained for glycyrrhizin (−19.22 kcal/mol), gitoxin (−17.84 kcal/mol), dicumarol (−12.05 kcal/mol), diosgenin (−10.75 kcal/mol), and delphinidin (−8.12 kcal/mol). These results are very promising compared to those obtained in recent works, where the coupling of natural and synthetic compounds with viral glycoprotein S, were analyzed. Thus, the best binding energies obtained for S protein were: −11.55 kcal/mol to coenzyme A, −11.089 kcal/mol to flavin adenine dinucleotide, and −9.36 kcal/mol to tiludronate (Hall & Ji, 2020).

In agreement with the low Gibbs free energy, numerous molecular interactions were observed for the glycyrrhizin: six hydrogen bonds and many hydrophobic interactions, all of them with residues from chains A and C from the glycoprotein. This is despite no interactions with the RBD region (recognition of ACE human receptors), chain B, were observed. Therefore, glycyrrhizin would be a good candidate to be evaluated experimentally to learn how it would affect the glycoprotein activity. Interestingly, a dietary supplement of Glycyrrhiza glabra L. (SoriaNatural®) is already in the market demonstrating the safety of consuming this compound since it has already been approved by the corresponding organizations.

Gitoxin is a secondary cardiac glycoside from Digitalis purpurea and Digitalis lanata. This compound and others from Digitalis sp. have demonstrated action on heart muscle (Haustein et al., 1970). The aglycone that is the active part of the compound has a cyclopentaneperhydrophenanthrene structure, which is common in steroids. This hydrophobic portion of the molecule established many hydrophobic interactions with the spike protein residues. The glycone portion of the gitoxin interacted with two N-acetyl-d-glucosamine residues from spike protein. It also formed hydrogen bonds. Many of these interactions of different nature were established in the RBD region of S glycoprotein. For this reason, we consider that it is the most promising compound with potential to inhibit spike protein.

Dicumarol is a derivative from cumarin and was successfully used as an inhibitor against influenza A (An et al., 2014). Moreover, it could be a particularly good candidate for inhibiting the S glycoprotein by blocking the RBD region, which contains human ACE2 receptors, as the compound interacted with several residues in that region. One important aspect to consider regarding dicumarol and other cumarins, is that they have already been approved in the pharmaceutical market for their effectiveness as: antidiabetic (Dwivedi et al., 2008), antiviral (Yu et al., 2003), anti-inflammatory (Kontogiorgis & Hadjipavlou-Litina, 2004), and anticancer and antileukemia activities (Stefanachi et al., 2011; Yang et al., 2011); thus, they are considered safe compounds for human oral consumption. Considering that dicumarol and coumarins have powerful anticoagulant effects, their use in certain patients must be under medical control. In fact, being derived from medicinal plants with ancestral use proven for oral consumption is an important advantage of some natural origin compounds. The medicinal use of phytocomplexes allows not only increasing the desired action of the main active compounds but also decreasing adverse effects than if any of the pure active compounds were used for the same treatment. Phytocomplexes are plant extracts or essential oils that contain more than one main active compounds and other substances that can help decreasing the adverse effects of the main compounds, contributing to the better absorption and solubility of the main compounds. This important advantage would facilitate the production and reduce the cost because fewer downstream operations would be needed as there is no need to obtain the high purity product (Cortés Gallardo et al., 2004; Duraffourd et al., 1987; Hostettmann & Marston, 1995; Pellecuer, 1993; Rawat et al., 2020; Sharapin et al., 2000).

Conclusions

COVID-19 is a global health problem. Despite all the efforts dedicated to stop its pandemic dissemination there is no conclusive treatment to date.

Based on the binding energies obtained and on the analysis of the frequency and dispersion of clusters, interaction conformations between the proposed natural compounds and Mpro and spike proteins are highly possible which means an effective and strong bond. This could translate into a high possibility that these active compounds are efficient inhibitors of viral proteins, which must be corroborated experimentally. In addition, the tested natural compounds performed better than the synthetic antivirals used as controls: nelfinavir and atazanavir, for protease Mpro and umifenovir, chloroquine and hydroxichloroquine for S protein. Thus, the compounds evaluated as inhibitors for the viral protease and glycoprotein S could be potential active agents for treatment of COVID-19, either alone or as adjunct therapies with other medications demonstrated have antiviral activity in vivo.

This study also contains valuable information to increase the knowledge of certain structural families of compounds with great affinities by Mpro protease and SARS-CoV-2 spike proteins. Furthermore, the present study provides molecular details that set the basis to further propose structural modifications of some compounds to make the interaction between them and proteins even more effective.

The following are the supplementary data related to this article.

Fig. S1

Cluster analysis of conformations. Frequency (number of conformers in clusters) of each compound after 200 runs considering the binding energy of conformations corresponding to the docking of α-amyrin, berberine, geranial, glycyrrhetinic acid, lapachenole, lapachol, limonene, meitansine, meliacine, menthol, N-(3-acetylglycyrrhetinoyl)-2-amino-propanol, nelfinavir, oleanane, sambunigrin, sauroxine, spirostan, taraxasterol, and theophylline with SARS-CoV-2 main protease (Mpro).

Fig. S2

CRediT authorship contribution statement

María Antonela Zígolo: Conceptualization, Methodology, Formal analysis, Writing – original draft. Matías Rivero Goytia: Formal analysis, Data curation. Hugo Ramiro Poma: Methodology, Writing – review & editing. Verónica Beatriz Rajal: Visualization, Writing – review & editing, Supervision. Verónica Patricia Irazusta: Visualization, Writing – review & editing, Supervision.

Declaration of competing interest

The authors declare no conflict of interests.