Molecular docking and dynamic simulations of some medicinal plants compounds against SARS-CoV-2: an in silico study

ABSTRACT

COVID-19 pandemic has poses urgent health challenge, and this project aims to identify potential inhibitors to combat this virus. We screened 198 bioactive compounds from five selected medicinal plants previously reported to be antiviral against SARS-CoV-2 protease and two co-receptors followed by molecular dynamics simulations. From the screened compounds, Astragalin demonstrated very strong molecular interactions with the molecular docking binding energies −8.5, −8.0, −7.6 kcal/mol for 6LU7, 6LZG, and 6VXX proteins of SARS-CoV-2, respectively. Hydrogen bonding interaction with the active site catalytic residue HIS-41 or CYS-145 of the main protease SARS-CoV-2 was observed. Binding free energies (ΔG_bind) from MM-GBSA after 50 ns MD simulations showed that Astragalin has the highest energy of −33.00 and −34.89 kcal/mol in complex with the main protease and spike glycoprotein of SARS-CoV-2, respectively. The study identifies Astragalin as a better inhibitor for the inactivation of COVID-19 and should be pursued as a potential drug candidate for this virus.

KEYWORDS:

• SARS-CoV-2
• ACE2 human
• protease
• molecular docking
• molecular dynamics
• medicinal plants

1. Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has recently emerged in December 2019 in Wuhan, China [1] causing millions of death worldwide. Common symptoms of SARS-CoV-2 infection are fever, dry cough, breathing difficulties (dyspnoea), diarrhoea, headache, pneumonia, and ventilator support. Disease onset may cause progressive respiratory failure owing to alveolar damage and even death [1,2]. The World Health Organization (WHO) declared it as pandemic on 11 March 2020 because of its fast-growing human-to-human transmission and reported cases in nearly every region [3,4]. So far, there is no effective cure for the treatment of the virus, thus the need to identify potent drugs/vaccines.

SARS-CoV-2 is an enveloped spherical or pleomorphic particle containing single-stranded (positive-sense) RNA paired with a capsid nucleoprotein containing a matrix protein. The envelope contains glycoprotein in the form of clubs with a genome of nearly 30,000 bases belonging to the genus betacoronavirus (β-CoV). The viral genome encodes four structural proteins, the envelope (E), the nucleocapsid (N), the membrane (M), and the spike (S) proteins [5,6]. The spike glycoprotein attaches the virion to the host cell membrane in two functional units, S1 and S2. S1 domain of SARS-CoV-2 spike glycoprotein potentially interacts with the human host cell targets such as (ACE2, cyclophilins, CD26, Ezrin, and other cell adhesion factors) a key immune-regulatory factor for hijacking and virulence [7,8]. SARS-CoV-2 attach to the host cells through Spike protein and human Angiotensin-Converting Enzyme 2 (hACE2) receptor interaction.

For the virus receptor (spike protein) to bind to its cellular ligand (ACE2), activation by transmembrane protease TMPRSS2 is needed [5]. One way to control the infection is to seek protease inhibitors to prevent the S-protein of SARS-CoV-2 cleavage to the S1 domain. Inhibiting the ACE2 receptor is likely to induce a conformational change of ACE2, therefore blocking the interaction between spike protein and ACE2 receptor.

Also, since SARS-CoV-2 is a single-stranded RNA virus, SARS-CoV-2 3-chymotrypsin-like protease (3CL^pro) also called main protease (M^pro) is a key CoV enzyme, which plays a critical role in facilitating viral replication and transcription, making it an attractive drug target for the virus [9]. Despite this, some protease inhibitors, and antiviral agents namely Chloroquine, Remdesivir and Hydroxychloroquine have been reported to be highly effective in the control of COVID-19 infection in in vitro under investigation as agents to treat COVID-19 [10–13].

Medicinal plants are important sources of structurally novel compounds that functions as a lead for the development of novel drugs [14]. Herbal medicine has the benefits of being stable at room temperature, non-injectable, and with no serious side effects [15]. The herbal preparation also can be used as a health drink, health food or dietary supplement, to prevent immune system impairment [15], which is needed for the recovery from SARS-CoV-2 infection [16].

The plant Ageratum conyzoides has been specified for the treatment of “craw-craw” better known as rashes, pneumonia by rubbing it on the patient's chest, diarrhoea, wound healing and skin diseases, infectious and mental diseases as well as headaches and dyspnoea [17]. Andrographis paniculata (Burm. F.) Nees herb is well explored therapeutically and effectively used as an immunostimulant, for asthma, influenza, respiratory tract infections, anti-pathogenic bacteria, anti-HIV, antimalarial, anti-snake venom, fever among others [18,19]. Andrographolide showed good inhibition of CHIKV infection, and it reduces the virus production by approximately 3log10 with EC₅₀ of 77 µM without cytotoxicity [20,21].

Azadirachta indica (neem) has proven to be multifunctional and correctly stated as “Tree for Solving Global Problems”. Neem has been extensively used conventionally as a remedy for many diseases. Studies have identified active compounds from different part of neem plant as antimalarial, anti-inflammatory, anti-hyperglycaemic, anti-mutagenic, anti-ulcer, anti-carcinogenic, anti-oxidant, anti-viral, and immune-modulator drugs [22–24].

Momordica charantia (Bitter gourd or Bitter melon) with similar pharmacological importance of A. indica has been widely used traditionally as unorthodox medicines for numerous sicknesses such as antitumor, anti-diabetic, anthelmintic, antimalarial, and antiviral agents [25]. Studies have confirmed that the protein of M. charantia has an effective and wide antiviral activity against different subtypes of influenza A virus [26] and Dengue virus [27,28].

Phyllanthus amarus is a widely distributed tropical medicinal plant, which belongs to the family Euphorbiaceae. The compounds extract of P. amarus show a wide spectrum of therapeutic activities such as antibacterial, antiplasmodial, antiviral, antimalarial, anti-inflammatory, anticancer, anti-microbial, anti-diabetic, anti-oxidant, hypolipidemic, nephroprotective, hepatoprotective, and diuretic properties [25,29–32]. Methanol extract of Phyllanthus species was tested for anti-HIV-1 reverse transcriptase (RT) mechanism using the HIV-RT assay by inhibition of the HIV-1 RT enzyme based on their IC50 values [28]. Some suggest that P. amarus may be helpful in the cure of hepatitis B virus infection [29].

Molecular docking can predict the ligand–protein interactions by recognizing the apt active sites in the protein, getting the best geometry of the ligand–receptor complex, and computing the interaction energy for various ligands to design additional effective ligands.

Owing to the numerous activities reported for these class of compounds, we pursue a virtual screening against receptors of SARS-CoV-2 to determine the inhibitory activity of the bioactive constituents of Ageratum conyzoides, Andrographis paniculata (Burm. F.) Nees, Azadirachta indica, Momordica charantia, and Phyllanthus amarus, which have been used for folk medicine for various diseases.

To discover lead compounds for clinical use, the aim of this article is to indicate the drug candidate from medicinal plants that can act as an antiviral, protease inhibiting agents, and RNA interfering agents for SARS-CoV-2. We initiated virtual screening of compounds, molecular docking, and molecular dynamics to single out new drug leads that target the main protease (MPro) and Spike glycoproteins (S) of SARS-CoV-2.

2. Materials and methods

2.1. Ligand preparation

Two hundred and eight (208) compounds (i.e. 198 extracted compounds from 5 selected plants and 10 extra compounds) were obtained from PubChem [33], in SDF format. Some of these ligands are essentially bioactive compounds in medicinal plants and information on their bioactivities were retrieved from Dr. Duke's Phytochemical and Ethnobotanical Database [34]. The ligands were opened in Spartan ’14 version 1.1.4 and converted to a 3D structure. The ligands were energetically minimized and geometrically optimized using the Semi-Empirical method with a basis set of PM3. Drug-like properties were calculated using Lipinski's rule of five [35] and Veber's standards [36] for oral bioavailability of drug candidates. The considered Lipinski's filters mass (≤500 g/mol), clogP (≤5), Hydrogen bond acceptors (≤5), hydrogen bond donors (≤10), Molar refractivity (), and polar surface area (≤140 Å${}^{\wedge }$2) [35,37]. Adherence with Lipinski's rule of five was calculated using SWISSADME prediction [38].

2.2. Receptor preparation

The crystal structure of SARS-CoV-2, COVID-19 3CL^pro/M^pro (PDB ID: 6LU7) [9] and spike glycoproteins (PDB ID: 6VXX, 6LZG) [39,40] were retrieved from the RCSB Protein Data Bank [41] in PDB format. All water molecules and heteroatoms were removed from the crystal structures using Discovery Studio software v. 2020 and then saved in PDB format for further analysis. For the spike glycoprotein (6VXX), alpha – chain (A chain) was used for the docking study. This is because there are multiple alignments of the amino acid sequences, this means, they have identical residue.

2.2.1. Determination of active sites

The amino-acids in the active site of the receptor was determined as it was reported in the study of the proteins [9,39,40]. The determination of the amino acids in the active site was used to analyse the docking evaluation results. Meanwhile, a three–dimensional grid box map was made at maximum on M^Pro (6LU7) and on the Spike glycoproteins (6LZG and 6VXX), to allow the ligand to be docked to all parts of the receptor.

2.3. Molecular docking

The prepared receptors (PDB ID: 6LU7, 6LZG, 6VXX) were treated, i.e. removal of water molecules, ligand, and heteroatoms from the receptor, with the use of the Discovery Studio programme. AutoDock version 4.2.6 was used for protein optimization and the addition of the polar hydrogen group. Virtual screening using automated docking involves the preparation of the receptor (comprises of assigning Gasteiger charges and Kollman charges to all atoms and assigning AD4 forms to protein structure atoms) prior to docking [42], ligands, grid box, and the algorithm. With PyRx and AutoDock Vina [43], the correct grid box size was evaluated for the possible binding site of the ligands, and the receptor grid centre was placed on the receptor's active site residue. AutoDock 4.2.6, Pymol, and Discovery Studio were used for the docking analyses. Analysis of the docked complexes was done to obtain compounds with better binding energy (BE) and ligand's inhibitory constant (Ki) with the target protein. Analysis of the structural complex of the ligand–protein interaction was done using Discovery Studio 2020 programme.

2.4. Molecular dynamics simulations

The top three categories of compounds based on their binding affinities in complexation with the main protease and the two spike glycoproteins of SARS-CoV-2 from the docking studies were subjected to molecular dynamics (MD) simulations. A similar MD protocol in our previous studies [44] was applied in this work. The simulations were conducted on the GPU version of the AMBER 14 software package. The molecular geometries of the ligands and proteins were reproduced using the general AMBER force field (GAFF) and AMBER force field 99SB, respectively. The ligand–protein complexes were neutralized with the required counter ions. Thereafter, the complexes were solvated with TIP3P water box of 12 Å. The solvated parameters (“top and crd)” for MD simulations were saved. The complexes were minimized using steepest descent and conjugated gradient minimizations for 2500 and 10,000 steps respectively to eliminate atomic clashes. Stepwise heating of the complexes was applied until 300 K and relaxation at the equilibration stage before the actual production. A 50 ns MD simulations were performed; trajectories analysis was extracted on the CPPTRAJ module of AMBER 14.

3. Results and discussion

3.1. Docking score

The binding affinities of all the compounds tested against the main protease and the two spike proteins are presented in Table 1. Based on the docking score, it is obvious that Astragalin, Isoquercitrin, Momordicin, Friedelin, Nelfinavir, Nimbaflavone, Momordicoside-I, Kaempferol, Gedunin, and Meldenin interacted well with the main protease and Spike glycoprotein of SARS-CoV-2 than Remdesivir and Chloroquine.

Table 1. Docking score result for SARS–CoV-2 receptors binding with the potential inhibitors.

Ligand	Binding affinity (kcal/mol) 6LU7	Binding affinity (kcal/mol) 6ZLG	Binding affinity (kcal/mol) 6VXX
5-hydroxy-7,8,2′-trimethoxyflavone	−7.5	−8.4	−7.4
Andrographolide	−7.1	−7.8	−7.0
Ascorbigen	−7.4	−8.3	−7.5
Astragalin	−8.5	−8.8	−7.6
Azadiradione	−7.6	−8.8	−7.9
Beta-amyrin	−7.4	−9.7	−7.8
Chloroquine	−5.9	−7.6	−5.3
Cucurbitacin-B	−7.5	−9.0	−6.6
Cucurbitacin-K	−7.6	−8.6	−6.9
Diosgenin	−7.3	−9.2	−8.3
Epiazadiradione	−7.5	−8.2	−7.6
Epoxyazadiradione	−7.4	−9.1	−7.7
Friedelin	−8.2	−9.0	−8.5
Gallocatechin	−7.2	−9.2	−6.8
Gamma-cadinene	−6.1	−8.6	−7.4
Gedunin	−7.7	−8.4	−7.4
Hydroxychloroquine	−6.1	−7.2	−5.7
Isoquercitrin	−8.5	−8.1	−6.7
Ivermectin	−5.6	−9.9	−8.5
Kaempferol	−7.8	−7.7	−7.1
Lectin	−7.4	−9.7	−6.8
Lopinavir	−7.1	−8.5	−6.4
Meldenin	−7.7	−7.9	−7.5
Meliantriol	−7.4	−8.3	−8.2
Momordicin	−8.4	−8.6	−7.7
Momordicoside-I	−7.9	−8.9	−8.0
Multiflorenol	−7.4	−9.2	−8.2
Myricetin	−7.3	−8.5	−7.2
Nelfinavir	−8.2	−9.2	−6.5
Neoandrographolide	−7.6	−8.7	−7.1
Nimbaflavone	−8.1s	−8.3	−7.8
Nimbinin	−7.5	−8.8	−7.6
Nimbolin-B	−7.6	−9.0	−7.2
Nimolicinol	−7.2	−9.1	−8.3
Nimolinone	−7.6	−8.7	−7.8
Oseltamivir	−6.1	−5.9	−6.0
Phyllanthusin E	−7.2	−7.7	−6.8
Quercetin	−7.3	−9.5	−7.1
Quercitrin	−7.7	−8.6	−7.2
Remdesivir	−7.7	−8.6	−6.9
Taraxerol	−7.4	−9.0	−8.0

Thirty-three active compounds in medicinal plants used as an antioxidant, antiviral, and antimalarial, which passed the Lipinski thumb rule alongside seven other compounds were considered for this study. In total, we docked two hundred and eight (208) compounds. Eight out of the forty-one compounds showed a good binding affinity (ΔG ≤ −7.8) with the Main protease (6LU7); these compounds are Astragalin, Isoquercitrin, Momordicin, Friedelin, Nelfinavir, Nimbaflavone, Momordicoside-I, and Kaempferol.

It was also found that thirty-four compounds have a strong interaction (ΔG ≤ −8.0) with the spike glycoprotein (6LZG); these are Ivermectin, Lectin, Beta-amyrin, Quercetin, Diosgenin, Multiflorenol, Nelfinavir, Gallocatechin, Nimolicinol, Epoxyazadiradione, Friedelin, Taraxerol, Cucurbitacin-B, Nimbolin-B, Momordicoside-I, Astragalin, Azadiradione, Nimbinin, Neoandrographolide, Nimolinone, Momordicin, Remdesivir, Cucurbitacin-K, Quercitrin, Gamma-cadinene, Myricetin, Lopinavir, Gedunin, 5-hydroxy-7,8,2′-trimethoxyflavone, Nimbaflavone, Meliantriol, Ascorbigen, Epiazadiradione, and Isoquercitrin. Besides, seventeen compounds were found to have a strong interaction (ΔG ≤ −7.6) with the spike glycoprotein (6VXX); these are Ivermectin, Friedelin, Diosgenin, Nimolicinol, Multiflorenol, Meliantriol, Taraxerol, Momordicoside-I, Azadiradione, Beta-amyrin, Nimolinone, Nimbaflavone, Epoxyazadiradione, Momordicin, Astragalin, Nimbinin, and Epiazadiradione. The key residue(s) for each ligand interaction with the SARS-CoV-2 receptors through hydrogen bond interaction and their respective distances were presented in Table S2. Hydrogen bonding interactions are essential for proteins, as they provide the organization for distinct folding and the selectivity that supports molecular recognition in the ligand–protein interface [45].

3.2. Post molecular dynamics simulations

The post MD analysis such as root mean square deviation fluctuation (RMSF), root mean square deviation (RMSD), and free binding energy were investigated. For clarity, the analysis in this section will be focused on the main protease while the results for the spike proteins will be provided in the supplementary information. The compounds simulated are Astragalin (Astra), Kaempferol (Kaem), and Nimbaflavone (Nimb) in complex with the main protease of the SARS-CoV-2.

3.3. Binding free energy calculations

The association of ligand–protein interaction is quantified using the MM/GBSA method [46,47]. The Generalized Born (GB) solvation model [48] is used to compute the electrostatic component of the solvation binding free energies. Estimation of the binding free energies are extracted from the equations below (1) $\begin{array}{rl}\mathrm{\Delta }{G}_{bind}& =G_{ligand\text{-}protein}-G_{protein}-G_{ligand}\end{array}$(1) (2) $\begin{array}{rl}\mathrm{\Delta }{G}_{bind}& =\mathrm{\Delta }{E}_{MM}+\mathrm{\Delta }{G}_{solv}-\mathrm{T}\mathrm{\Delta }\mathrm{S}\end{array}$(2) (3) $\begin{array}{rl}{E}_{MM}& =E_{int}+E_{ele}+E_{vdw}\end{array}$(3) (4) $\begin{array}{rl}\mathrm{\Delta }{G}_{solv}& =\mathrm{\Delta }{G}_{GB}+\mathrm{\Delta }{G}_{SA}\end{array}$(4) (5) $\begin{array}{rl}\mathrm{\Delta }{G}_{SA}& =\gamma \text{SASA}+\text{b}\end{array}$(5)

Where $G_{ligand\text{-}protein}$, $G_{protein}$, $G_{ligand}$ is the energy of the complex, protein and ligand respectively. The energy components such as $E_{MM}$ is the molecular mechanics energy of the system, $E_{int}$, electrostatic energy, $E_{ele}$, and van der Waals term, $E_{vdw}$. Also $\mathrm{\Delta }{G}_{solv}$, solvation energy which is divided into the polar ($\mathrm{\Delta }{G}_{GB}$) and non-polar ($\mathrm{\Delta }{G}_{SA}$) contributions. $\mathrm{\Delta }{G}_{SA}$, is the non-polar contribution to solvation-free energy that is linearly dependent on the solvent accessible surface area (SASA).

3.4. System stability and flexibility

RMSD value for the heavy backbone atoms was calculated for the 50 ns simulation for all the systems to validate their stabilities. The RMSD plots are represented in Figure 1 with average values of 1.85, 1.25, and 0.82 Å for Astra-6LU7, Kaem-6LU7, and Nimb-6LU7, respectively. The small average RMSD values of the complexes indicate that the simulations were already converged [49].

Figure 1. RMSD plots for the studied complexes.

It can be observed that the systems were equilibrated after 5 ns, although some deviations were seen during the simulation period. A stable protein structure is judged with a corresponding smaller deviation [50]. The residual flexibility of the studied proteins on binding with the ligands was evaluated using the RMSF parameter. The RMSF (Figure 2) of the α-Carbon backbone for the systems showed a similar evolution. However, the Kaem-6LU7 system displayed more flexible residues at residue number 152–155 in comparison to Astra-6LU7 and Nimb-6LU7. In general, different protein regions showed varying levels of flexibility, which impacted the binding interactions with the ligands as represented on Table 2.

Figure 2. RMSF plots for the α-Carbon backbone of the main protease for three systems.

Table 2. The calculated binding free energies and their relative components for the complexes using MM-GBSA method obtained from the last 1000 frames of 50 ns. The energy components are in kcal/mol.

Compound	ΔEvdw	ΔEele	ΔGgas	ΔGpolar	ΔGnon-polar	ΔGsolvation	ΔGbind
Astra-6LU7	$-44.56$	$-8.63$	$-53.19$	$25.15$	$-5.28$	$19.87$	$-33.32$
Nimb-6LU7	$-42.69$	$-4.24$	$-46.92$	$19.98$	$-4.92$	$15.06$	$-31.86$
Kaem-6LU7	$-32.23$	$-6.34$	$-38.57$	$17.46$	$-4.24$	$13.22$	$-25.35$

3.5. Binding free energy calculations

To evaluate the binding free energy of the ligand–protein complexes, we employed the MM-GBSA method [51,52] and they were generated from the last 1000 frames of 50 ns. An in silico prediction of ligand–protein binding free energy has been reported to be the main objective of structure-based drug design [53]. The binding free energies for the systems are presented in Table 5 with the accompanying energy components. According to the calculated results presented in Table 2, ΔG_bind for Astragalin, Nimbaflavone, and Kaempferol in complex with the main protease of SARS-Cov-2 were −33.32, −31.86, and −25.35 kcal/mol respectively. The results for the Astra – 6LU7 complex demonstrate larger binding free energies compared to Nimb – 6LU7 and Kaem – 6LU7 with a difference of 1 and 8 kcal/cal respectively. Binding free energies for the spike proteins (6LZG and 6VXX) are presented on Tables S3 and S4.

The relative energy contributions for the systems are consistent with their binding free energies with Astra – 6LU7 having the highest increases in electrostatic energy when binding in both gas, polar and non-polar phases and a favourable change in van der Waals energy when binding, the change in electrostatic energy compensates significantly. Important residual interactions were observed as presented in Figure 3 for the systems. Hydrogen bonding interactions were seen for GLU166, GLY143, SER144, ASN142, there is interaction with the catalytic residue CYS145 (Figure 3).

Figure 3. Ligand-protein interactions resulting from MM-GBSA calculations.

4. Conclusion

Using various in silico methods, we have identified some potential antiviral inhibitors, which can serve as a starting point to stop the replication of SARS-CoV-2. These inhibitors were derived from known medicinal plants. The use of virtual screening and molecular dynamics simulations were utilized for this study. Top lead compounds (n = 41) showed that many of the compounds interact with the active site residues of SARS-CoV-2 from the virtual screening protocol. Molecular dynamics simulations of the top 3 compounds were performed, Astragalin, Nimbaflavone, and Kaempferol in complex with the main protease of SARS-CoV-2. The results of the study revealed that Astragalin and Nimbaflavone in complex with the main protease of SARS-Cov-2 showed the best binding free energies ΔG_bind of −33.32 and −31.86 kcal/mol, respectively. Important residual interactions with the active region of the receptors were observed, such as Cys 145. It has been reported that Astragalin has antiviral activity and has underlying anti-H1N1 influenza virus activity [53]. Studies have demonstrated that Astragalin functions as an effective lead compound antagonizing oxidative stress induced by endotoxin, which leads to airway dysfunction and inflammation [54]. The study which utilized state-of-the-art virtual screening and molecular dynamics simulations are crucial techniques required for drug design and identification of new inhibitors.

Acknowledgement

The authors are grateful for the computational resources provided by CHPC (www.chpc.ac.za).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

Where, ΔEvdw van der Waals energy, ΔEele electrostatic energy, ΔG_gas gas-phase free energy, ΔG_polar polar free energy, ΔG_non-polar non-polar free energy, ΔG_solvation solvation free energy. ΔG_bind total binding free energy.