Reckoning γ-Glutamyl-S-allylcysteine as a potential main protease (mpro) inhibitor of novel SARS-CoV-2 virus identified using docking and molecular dynamics simulation.

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2 or COVID-19), outbreak was first reported in December 2019 in the Wuhan, China. COVID-19 managed to spread worldwide and so far more than 9.1 million cases and more than 4.7 lakh death has been reported globally. Children, pregnant women, elderly population, immunocompromised patients, and patients with conditions like asthma, diabetes, etc. are highly vulnerable to COVID infection. Currently, there is no treatment available for COVID-19 infection. Traditional medicinal plants have provided bioactive molecules in the past that are efficiently used during conditions like cancer, malaria, microbial infections, immune-compromised states, etc. AYUSH India has recommended the use of Curcuma longa, Allium sativum, Ocimum tenuiflorum, and Withania somnifera for immune-boosting during SARS-CoV-2 infection. In the present study, we investigated the potential of 63-major bioactive molecules of these plants against SARS-CoV-2 main protease (Mpro) through docking studies and compared the results with known inhibitor 11a. Our results proposed cuscohygrine, γ-Glutamyl-S-allylcysteine, anahygrine, and S-allylcystein as the potent inhibitors against Mpro identified using molecular docking and molecular simulation dynamics. Interestingly, these molecules are from A. sativum, and W. somnifera, which are known for their antimicrobial and immunomodulatory potential. None of the proposed molecules have earlier been reported as antiviral molecules. Our results predict very strong potential of these four-molecules against SARS-CoV-2 Mpro, especially γ-glutamyl-S-allylcysteine, as all four form hydrogen bonding with Glu166 that is a crucial residue for the formation of the biologically active dimeric form of Mpro. Therefore, we strongly recommend further research on these biomolecules against SARS-CoV-2.

Introduction

Past two decades have witnessed the widespread health complications due to infections caused by Coronaviruses. These viruses are non-segmented positive sense RNA viruses belonging to the family of Coronaviridae and the subfamily of Orthocoronaviridae under superfamily of Nidovirales [1]. Initially, Coronaviruses were divided into the three groups (1, 2, and 3) based on their antigenic reactivity and variations in the genomic sequences. Recently, International Committee on Taxonomy of Viruses proposed new classification for the Coronaviruses in which they are divided into three genera, α-coronaviruses, β-coronaviruses and γ-coronaviruses, which corresponds to the group 1, 2 and 3 in the previous classification [2,3].

Coronaviruses are known to cause variety of infections in both, animals and humans, primarily targeting the respiratory system. In humans, coronaviruses are known to cause variety of minor infections such as common cold. HCoV-229E and HCoV-OC43 were the first two coronaviruses reported in humans [4], followed by HCoV-HKU1 and HCoV-NL63 [5]. All these four viruses rarely cause medical emergencies (unless associated with comorbid disease conditions) and are known to cause acute upper respiratory tract infections more frequently than lower respiratory tract infections [6,7]. Human coronaviruses (HCoV) are also capable of causing severe respiratory tract infection, of which zoonotic Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) and The Middle East Respiratory Syndrome Coronavirus (MERS-CoV) are best known [8]. SARS-CoV and MERS-CoV are known to cause serve health complications in the past one and a half decade and were designated as global epidemics due to heavy loss of life caused due to their infections globally [9,10]. These viruses may prove lethal in special category patients, particularly infants, elderly, pregnant females, immunocompromised patients, and patients having history of respiratory disease [5]. Currently, there is no specific treatment and drugs available for countering coronavirus infections and human immunity remains the primary defense against these viruses [11,12].

Recently, world witnessed yet again global pandemic due to the coronavirus which has still not reached its peak if seen in terms of global mortality. Novel Human Coronavirus Disease (COVID-19) outbreak was first reported in the December 2019 in the Wuhan city of the Hubei province in the central China [13,14]. Since then COVID-19 has managed to spread worldwide across 213 countries and territories with more than 9.1 million cases and more than 4.7 lakh death globally (as on 23rd June 2020). These figures are continuously increasing and are expected to go much higher in coming days, especially in countries like the United States of America, Brazil, Russia, Spain, United Kingdom and India, which have witnessed community spread of COVID-19. Due to the lack of effective treatment so far, the World Health Organization (WHO) has been actively involved in the development of diagnostic tools and issuing guidelines for treatment procedures, patient care, prevention and monitoring, besides, promoting research for developing an effective drug/vaccine for the treatment of COVID-19 [15]. The virus, believed to be originated in bats, has now been known to primarily infect human respiratory system through human-to-human transmission via respiratory droplets and direct contact [6,14].

Since, there is no drug or vaccine available for the treatment of COVID-19, the concept of drug repurposing with oral Hydroxychloroquine has shown some positive effects in some reports. Moreover, the hunt for identifying novel targets to prevent replication of the SARS-CoV-2 is ongoing. This has led to the identification of several nonstructural SARS-CoV-2 proteins, which includes main protease (Mpro), Nonstructural protein 10 (nsp10), nsp12, helicase (nsp13), N-terminal exoribonuclease and C-terminal guanine-N7 methyl transferase (nsp14), Uridylate-specific endoribonuclease (nsp15), 2′-O-methyltransferase (nsp16), and RNA-dependent RNA polymerase (RdRp). These SARS-CoV-2 proteins function to replicate and pack the viral genome [16,17]. Blocking the functioning of these proteins which adversely affect the replication and packing of the SARS-CoV-2 genome, and therefore, these proteins are now being considered as the key targets for drug development against SARS-CoV-2 [17].

Development of an effective treatment will require months, which will result in a lot more mortality globally and will have disastrous economic effect. Therefore, the search for the effective molecules for the management of COVID-19 should be directed toward traditional medicinal plants [18,19]. The Ministry of Ayurveda, Yoga & Naturopathy, Unani, Siddha and Homoeopathy (AYUSH), India has taken an initiative to promote the use of medicinal plants used in traditional Indian system of medicines for the management of infections for the management of COVID-19. These plants include Curcuma longa (Turmeric), Allium sativum (Garlic), Ocimum tenuiflorum (Tulsi), and Withania somnifera (Ashwagandha). All these traditional medicinal plants are known to have variety of medicinal values which includes antiviral, anti-inflammatory, and immune-boosting potential [20-23]. Currently, here is lack of effective curative and preventive therapeutic measures against COVID-19 crises and the immune system of the infected person is the most effective remedy to fight against this virus during non-sever early stage and severe later stage of the viral infection. It is now evident that during the early non-severe stage of COVID infection immune system prevent the spread of virus and during sever later stages, strong immunity is known to restrict the harmful effects of cytokine storm syndrome. All the traditional medicinal plants recommended by the AYUSH to be used during COVID pandemic are reported in Ayurveda as Rasayana or immunomodulatory [24]. Considering the potential of these plants to curb COVID infection based on the traditionally accumulated knowledge, AYUSH India issued an advisory to use these plants to control the health crisis arising due to COVID infection [25]. However, the experimental evidence for their effectiveness against the COVID-19 is not yet established.

In our present study we selected these plants as per the recommendations of AYUSH, raditional medicinal values of these plants and for their oral efficacy (suggesting their use orally can be intended). We further identified the bioactive molecules present in these plants and subjected them to molecular docking study to identify molecules with potential anti-coronavirus effect by investigating their interaction energy and type of interactions in concerned enzyme of SARS-CoV-2. For this study, we selected the Mpro of SARS-CoV-2. The Mpro of SARS-CoV-2 is crucial for the viral gene replication and gene expression via proteolytic processing of replicase polyproteins. Targeting Mpro can result in the inhibition of viral replication and therefore is an attractive target for drug design against SARS-CoV-2 [26-29].

Material and methods

Datasets

The x-ray crystal structure of the Mpro of SARS-CoV-2 was selected from the PDB database. We selected the crystal structure of Mpro (PDB Id: 6M0K) [26] with resolution 1.50 Å. The standard inhibitor for this binding cleft is identified as the 11a ({N}-[(2 ∼ {S})-3-(3-fluorophenyl)-1-oxidanylidene-1-[[(2 ∼ {S})-1-oxidanylidene-3-[(3 ∼ {S})-2-oxidanylidenepyrrolidin-3-yl]propan-2-yl]amino]propan-2-yl]-1 ∼ {H}-indole-2-carboxamide). 11a has already been established as the inhibitor of Mpro and this was used as a positive control in our study. Mpro is a homo dimer of which only one chain is available in this crystal structure. Table S1 depicts the dataset of potential bioactive molecules of Curcuma longa, Allium sativum, Ocimum tenuiflorum, and Withania somnifera, which were used for molecular docking investigations.

Ligand-Preparation

The selected bioactive molecules were prepared using the protocol “prepare protein” of Accelrys Discovery studio package (BIOVIA, 2016). The 3-D conformers of these bioactive molecules were generated and saved as mol2 format for further analysis. DFT minimization protocols of Gaussian16 was used to optimize the ligand geometry [30,31]. The optimized molecules were then used for the docking analysis.

Molecular Docking

Molecular docking studies were performed using the SeeSAR version 9.2. Ligand binding site and the coordinates for the docking of the test compounds was kept same as for the positive control of the Mpro, 11a. The binding affinities of the selected bioactive molecule with the Mpro of SARS-CoV-2 were calculated by the Hydrogen dehydration (HYDE) scoring function of the software [32,33]. The estimated binding affinity in SeeSAR ranges from mM < μM < nM < pM. In our study, the selection of the best scoring models was dependent on various parameters, such as Ligand Efficiency (LE), estimated binding affinity, and Torsion [34]. Moreover, the scoring function of HYDE depends on atom type-specific hydration and desolvation terms, which were conservatively aligned by Octanol to water partition coefficients (Kow) of small molecules. The water molecules which were residing inside the ligand were expelled out and water molecules which were around the ligand were removed. This removed the chances of H-bond interactions between water molecules and protein or ligands, which may have resulted in a disfavored enthalpic contribution. Further, this formation of new H-bonds between the ligands and the target protein resulted in the compensation of the loos of energy as a result of water molecule removal. Unfavorable energy due to hydrophobic interactions between ligands and target protein prompt the breakage of H-bonds. Moreover, the energy of the ligand-protein complex was increased due to the removal of water molecules from the hydrophobic pockets of the complex and this phenomenon known as the hydrophobic impact [32,33]. An important drug-like parameter, the Lipophillic Ligand Efficiency (LLE), was calculated by merging the in-vitro binding strength of a ligand and lipophilicity [35]. The LE indicate the potency per atom and it is numerically calculated as the quotient of ΔG and the number of non-hydrogen atoms of the compound [36]. After completion of the molecular docking analysis, all the complexes were analyzed on basis of their binding affinity, LE, and Torsion values and the top most molecules were selected.

Molecular dynamics (MD) simulation

The best two ligands predicted from the docking study were used for the molecular dynamics simulation of Mpro in the presence of 11a using GROMACS 2019 software [37-39]. All the test docked complexes were compared with the 11a-Mpro docked complex (taken as control). SwissParam was used to generate the topology of the ligands. This provided us with the parameters and topology for test molecules compatible with CHARMM all atoms force field, for use with CHARMM and GROMACS [40]. The topology of the protein was created using GROMACS utilities using CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins) with the water model set to TIP 3-point. The Ligand-protein complex were defined with unit cell box under periodic boundary conditions using 1.0 nm distance from the protein to the box faces with triclinic shape and was filled with water [41]. Further, Na+ and Cl- counter ion neutralization was performed along with the energy minimization to equilibrate the system under NVT (constant particle number, constant volume and constant temperature) for 25 ns at 300 K. After completion of the NVT run, the system was proceeded with NPT (constant particle number, constant volume and constant temperature) simulation and the molecular dynamics run was performed for 25 ns. LINCS or Linear Constraint Solver algorithm was used to constrain all the covalent bonds [42] and the PME or Particle Mesh Ewald method was used to treat the electrostatic interactions. The cutoff radii for the van der Waals and Coulomb interactions were set at 14.0 Å and 10.0 Å respectively. Trajectories were recorded after completion of NPT and NVT simulations and were analyzed for root-mean-square fluctuation (RMSF), root-mean-square deviation (RMSD), number of Hydrogen bonds formed between the ligand and Mpro protein, and radius of gyration (Rg) by using the ’gmx rmsf’, ‘gmx rms’, ‘gmx hbond’ and gmx gyrate’of GROMACS utilities [37-39]. The stabilkity of the Ligand–protein complex was determined by the dynamics of ligand-protein hydrogen bonds with respect to time. Graphs were prepared using the XMgrace tool.

Binding free energy calculations MM/GBSA

The single trajectory approach was used for the binding free energy calculation using molecular mechanics generalized Born surface area (MM/GBSA) [43-46]. MM/GBSA in PRIME module of Maestro 11.4 was used to calculate the thermodynamic data of Coulomb energy (ΔGCoulomb), total free energy change (ΔGBind), Hydrogen-bonding correction (ΔGHbond), Pi-Pi packing correction (ΔGPacking), Lipophilic energy (ΔGLipo), and Van der Waals energy (ΔGvdW) for the trajectories obtained on performing MD simulations. MM/GBSA calculations were performed with the OPLS_2005 force field that use the VSGB 2.0 solvation model.

The free energy values were calculated using the following Equations (1) and (2): and where Δ is the conformation entropic contribution, and ΔG is the molecular mechanics' interaction energy (electrostatic + van der Waals interaction) between protein and ligand. ΔG and ΔG depict the polar solvation energy and the nonpolar solvation energy, respectively.

In silico pharmacokinetic properties using ADMET analysis

In our study, we used the pkCSM - pharmacokinetics server [47] for predicting the absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties of the top hits (Chloroquine, Cuscohygrine, γ-Glutamyl-S-allylcysteine, Anahygrine and S-allylcystein). This server predicted physiochemical as well as pharmacological properties. Simplified Molecule Input Line Entry Specification (SMILES) of the selected molecules were retrieved from PubChem, followed by uploading them to pkCSM - pharmacokinetics server. The server computed in-vivo absorption parameters like water solubility in the buffer system (SK atomic types, mg/L), Human intestinal absorption (HIA, %), in-vivo Caco2 cell permeability (Human colorectal carcinoma), in-vivo skin permeability (logKp, cm/hour), and in-vivo P-glycoprotein inhibition. We determined the metabolic parameters by using in-vivo Cytochrome P450 2C9 inhibition, in-vivo Cytochrome P450 2C19 inhibition, in-vivo Cytochrome P450 2D6 inhibition, in-vivo Cytochrome P450 3A4 inhibition, in-vivo Cytochrome P450 2D6 substrate, and in-vivo Cytochrome P450 3A4 substrate. For the distribution properties we included tests like, Blood-Brain Barrier (BBB) penetration, Central Nervous System (CNS) permeability and Lipinski’s Rule (Rule of Five). To access the toxicity of compounds, a range of vital endpoints such as, ames test, acute algae toxicity, 2 years carcinogenicity bioassay in rat, 2 years carcinogenicity bioassay in mouse, in-vivo Ames test result in TA100 strain (Metabolic activation by rat liver homogenate) were computed. Moreover, many drugs are often withdrawn at clinical trial stages due to their poorer renal clearance, which makes excretion a very important parameter. Therefore, in this study we also included total renal clearance and renal OCT2 Substrate to identify the excretion efficacy of the molecules under study.

Results and discussion

Docking and interaction studies

The search for the bioactive molecules for the management of COVID-19 infection is ongoing throughout the world. The traditional drug screening process is highly time consuming with huge financial constraints, and has low degree of success rate. The availability and utilization of docking techniques has significantly reduced the screening time and have delivered good results in recent past [48,49]. Bioactive molecules from various plants show extensive diversity in their chemical and structural properties, which opens up a vast area of research for exploration of their effects on the target sites of various proteins involved in several ailments, including COVID-19. For example, a recent study demonstrated that the bioactive molecules from the Tea plant, viz. Theasinensin-D, Oolonghomobisflavan-A, and Theaflavin-3-O-gallate, possess better binding affinity toward Mpro site of SARS-CoV-2, when compared to atazanavir, darunavir and lopinavir [48], making these molecules as point of interest for COVID research.

Molecular docking techniques are used to predict the binding sites on the target proteins (or receptors) and then explore the binding affinity of the ligand (drug under consideration) with the target protein [19,50]. Docking studies can simultaneously investigate thousands of molecules and rank them according to their binding affinities by taking into consideration multiple factors involved in ligand-drug interaction [49,51]. Further, high resolution crystallographic structure of RdRp-RNA complex has unlocked vast possibilities in the screening of antiviral drugs targeting a specific protein of SARS-CoV-2 [48,52]. The present study was aimed to explore the potential of bioactive molecules from the medicinal plants, which were recommended to be beneficial against COVID-19 by the AYUSH India, for their potential against SARS-CoV-2 Mpro. We identified the bioactive molecules (Table S1) of C. longa, A. sativum, O. tenuiflorum, and W. somnifera and screened for their potential against SARS-CoV-2 Mpro.

A total of 63 bioactive molecules along with Chloroquine were identified and compared with internal standard 11a for their binding affinities against the active site of SARS-CoV-2 Mpro protein (Table S2). Reports suggest that chloroquine is having beneficial effect during SARS-CoV-2 infection and is having a good potential to target SARS-CoV-2 Mpro protein [53-56]. The docking results of best four molecules, compared to Chloroquine, are depicted in Table 1 in terms of estimated affinity (Ki), ligand efficiency (LE), lipophilic ligand efficiency (LLE), and Torsion values. Cuscohygrine, γ-Glutamyl-S-allylcysteine, anahygrine, and S-allylcystein were selected on the basis of their estimated binding affinity against Mpro of SARS-CoV-2 (Table S2 and Table 1), and their chemical structures are depicted in Table 2.

Table 1.

Interaction of the natural bioactive molecules with main protease (Mpro) of SARS-CoV-2.

Ligand Name	LLE	LE	Range of Estimated Affinity (Ki)	Torsion
Chloroquine (standard drug)			553.5 μM − 54.9 mM
Cuscohygrine			1.20 μM − 120.1 μM
γ-Glutamyl-S-allylcysteine			18.7 μM − 1.86 mM
Anahygrine			14.3 μM − 1.42 mM
S-Allylcystein			313.6 μM − 31.5 mM

Table 2

The chemical structures and docking results of the top compounds.

S. No.	Name	Structure	Range of Estimated Affinity (Ki)
1.	Chloroquine		553.4 µM − 54.9 mM
2.	Cuscohygrine		1.20 μM − 120.1 μM
3.	γ-Glutamyl-S-allylcysteine		18.7 μM − 1.86 mM
4.	Anahygrine		14.3 μM − 1.42 mM
5.	S-allylcystein		313.6 μM − 31.5 mM

SARS-CoV-2 Mpro is having a total of four binding sites, S1′, S1, S2, and S4, and the 11a molecule is known to inhibit the activity of the SARS-CoV-2 Mpro by interacting at the S1′ site of this protein [26]. Docking analysis for the test molecules was performed on the coordinates corresponding to S1′ site of Mpro. We observed that 11a molecule is capable of interacting with the multiple amino acids at the S1′ site of Mpro, which include hydrogen bonding with Gly143, Cys145, and Glu166; Pi-Pi stacking with Leu141, Met165, and Pro168; Pi-sigma interaction with His41, and Met49; and Halogen interactions with Arg188, and Gln189. These results are depicted in the Figure 1. Interaction with the Cys145 is regarded as the crucial amino acid for inhibiting the activity of Mpro, as demonstrated in our results for the 11a-Mpro interaction [26]. It is worth mentioning that the 11a is a large molecule and was developed using computer aided drug design. It can inhibit the SARS-CoV-2 in in-vitro conditions, but may pose a threat to human life if administered. This is a major limitation of 11a which has led to divert the search for potential safer inhibitors toward the molecules of natural origin and/or molecules which are already being used clinically (drug repurposing). Chloroquine is one of the example of drug repurposing. It was used during SARS-CoV-2 infection and has shown some positive results to suppress the progression of viral infection. The docking pose of chloroquine within the active site of SARS-CoV-2 Mpro are depicted in Figure 2. Chloroquine binds to the active pocket of Mpro with an estimated affinity between 553.5 μM–54.9 mM. Both the LLE and LE parameters of molecular docking were very low for chloroquine in the binding site of Mpro. The binding was supported by favorable Torsion values. Moreover, Chloroquine formed a hydrogen bond with residues Asn142, His164 and Pi/Pi-alkyl interactions with residue His41 while formed Pi-Sulfur bond with Met165 of the Mpro of SARS-CoV-2.

Figure 1.

The interaction of 11a in the binding pocket of SARS-CoV-2 Mpro shown in two-dimensional, and in three-dimensional space.

Figure 2.

Superimposition of all the ligands at the S1’ site of Mpro.

Based on the estimated affinity of Chloroquine (553.5 μM − 54.9 mM), we selected four molecules, viz. cuscohygrine, γ-glutamyl-S-allylcysteine, anahygrine, S-allylcystein, based on their binding affinities as depicted in Table 2. The docking poses of these molecules within the active site of SARS-CoV-2 Mpro are depicted in Figure 2. Cuscohygrine, γ-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein binds to the active pocket of Mpro with an estimated affinity between 1.20 μM − 120.1 μM, 18.7 μM − 1.86 mM, 14.3 μM- 1.42 mM, and 313.6 μM- 31.5 mM respectively. The LLE and LE parameters of molecular docking for all molecules were higher than chloroquine. LLE and LE parameters for γ-glutamyl-S-allylcysteine were observed to be higher than the rest of the molecules in the binding site of Mpro. The binding for all the natural molecules was also supported by favorable Torsion values. Amino acid interaction study of Anahygrine showed that it poorly interacts with the protein. It makes five bonds of which it forms only one H-bond. Cuscohygrine formed a hydrogen bond with residues Gly143 and Asn142, Pi Sulfur interaction Glu166. S-Allylcysteine and γ-glutamyl-S-allylcysteine showed the most promising results. S-Allylcysteine formed five H-bonds, namely with Asn142, Ser144, Gly143, Cys145 and His163; formed pi-alkyl interaction with Leu27. Here more importantly, this compound formed a hydrogen bond with Cys145, making this compound an impressive candidate as an inhibitor of Mpro. γ-glutamyl-S-allylcysteine showed the most optimum interaction with Mpro as it could make multiple interactions. It made in total six hydrogen bonds, involving five amino acids. These amino acids are His164, Glu166, Gln189, Arg188 and Thr190. Gln129 made unfavorable donor interaction, His163 formed Pi-alkyl interaction and Met165 formed carbon-hydrogen bond (Figure 3).

Figure 3.

Interactions made by the Ligands under investigation with Mpro.

Molecular dynamics simulation

We conducted molecular dynamics simulations for different protein–ligand complexes, which includes positive control 11a-Mpro, S-Allylcysteine-Mpro, γ-glutamyl-S-allylcysteine-Mpro individually and their results were compared. S-Allylcysteine and γ-glutamyl-S-allylcysteine were selected for molecular dynamics simulation based on the docking interaction observed with the Mpro. Our results demonstrated the good overall stability of the Mpro proteins in presence of the inhibitors being tested, S-Allylcysteine and γ-glutamyl-S-allylcysteine, in comparison 11a. The stability of the model system was checked using the RMSD, RMSF, hydrogen bonds formed, and by evaluating the total energy of the ligand-protein complex along the changes in pressure and temperature during the course of simulations.

RMSD is very crucial parameter to investigate the equilibrium of MD trajectories. RMSD of the protein backbone atoms is used to check the stability of ligand-protein complex during the simulations with respect to the time function. We calculated the RMSD values of the Mpro protein backbone against the simulation time scale (0–25) ns) during its interaction with the internal ligand 11a, S-Allylcysteine and γ-glutamyl-S-allylcysteine (Figure 4).

Figure 4.

Representation of ligand RMSD of Mpro backbone during interaction with internal standard 11a, S-Allylcysteine, and γ-glutamyl-S-allylcysteine during their interaction with Mpro of COVID-19 derived from NVT Simulation at 300 K.

The RMSD values of 11a-Mpro trajectories was least and could be sought to be the most stable interaction which ranged below 0.15 nm during simulation for the backbone of protein. While, S-Allylcysteine-Mpro interactions showed the RMSD to fluctuate for the first 5 ns of the simulation thereafter achieving stabilization and in this case the values never surpassed 0.6 nm. γ-glutamyl-S-allylcysteine-Mpro system showed much stabilized interactions as its RMSD values were much lower than that of S-Allylcysteine-Mpro complex and values under this case never surpassed 0.35 nm. However, despite there being differences in the ranges of RMSD fluctuations of all three systems, they all could be considered stable as none of the system showed RMSD values above 1 nm (Figure 4).

With respect to the average molecular dynamics simulation conformation, the RMSF reveals the means of portraying flexibility differences amongst various residues under investigation. The RMSF of the backbone atoms of each residue of (i) 11a-Mpro (as a positive control) and (ii) S-Allylcysteine-Mpro and (iii) γ-glutamyl-S-allylcysteine-Mpro individually was calculated to reveal the flexibility of the backbone structure in presence of both the ligands. This study indicates the behavior of protein in the presence of ligands under investigation. In this study, high RMSF value suggests greater flexibility and low RMSF value suggest the limited movements during simulation in relation to its average position. The RMSF of the residues are shown in Figure 5.

Figure 5.

Representation of Molecular Dynamics (a) RMSF values of Mpro backbone during its interaction with 11a, (b) S-Allylcysteine (c) γ-glutamyl-S-allylcysteine during simulation.

In the present study, we observed almost identical RMSF of the Mpro protein backbone in presence of all three molecules studied, suggesting the similar stability of protein. Moreover, RMSF values are higher when the binding is poor, which was definitely not the case in our studies (Figure 5). The inter-molecular hydrogen bonding between the ligand and protein is important for the stabilization on ligand-protein complex. The stability of the hydrogen bonds formed by all three molecules under investigation with Mpro was determined throughout the simulation at 300 K for the ligand system (Figure 6). The 11a-Mpro complex exhibited maximum of five H-bonds formation during the course of simulation. While for S-Allylcysteine-Mpro, there were consistently maximum of three hydrogen bonds being formed throughout the simulation study. But the results of docking were validated here in case of γ-glutamyl-S-allylcysteine, where this ligand could consistently make four to five H-bonds with Mpro throughout the course of simulation. Overall scenario suggests both natural ligands not only fits very well in to the binding cleft of 11a on S1’ site of Mpro, but also make appropriate hydrogen bonds with the protein. Further, γ-glutamyl-S-allylcysteine is a much smaller molecule in comparison to 11a and is observed to make almost same number of H-bonds with the Mpro as made by 11a. Therefore the ratio of hydrogen bond to surface area is much greater for γ-glutamyl-S-allylcysteine compared to 11a, suggesting very strong interaction with Mpro in the S1’ cleft. γ-glutamyl-S-allylcysteine, therefore, by making enough hydrogen bonds to settle and block the space in the cleft might render the protein ineffective and this could be the key rationale behind its mode of action.

Figure 6.

Representation of hydrogen bond formation of (a) 11a, (b) S-Allylcysteine (c) γ-glutamyl-S-allylcysteine, with Mpro during simulation.

Molecular dynamics simulations performed at NVT reflects the stability of the ligand-protein complex in terms of the deviations from the pressure and temperature. In these studies, higher energy of the ligand-protein complex is undesirable and reflects the instability of the system.

From Figure 7(a) it is evident that all the three systems have almost identical total energy suggesting the similar stability of all the three systems i.e. 11a-Mpro, S-Allylcysteine-Mpro, and γ-glutamyl-S-allylcysteine-Mpro. Further, Figure 7(b) suggests the deviation in the temperature from the constant 300 K, where it is evident that the 11a-Mpro system has maximum deviation, this can be due to the large size of the ligand (11a) attached to the protein (Mpro). Since 11a is relatively larger in size and has more possible interactions with Mpro, this could have resulted in the increased deviation in the overall temperature during simulations. Figure 7(c) suggests identical deviation in the pressure changes for all the three protein-ligand systems, suggesting their identical behavior in simulation. The results of the molecular dynamics simulations suggest that the protein behave in a similar way in presence of the all three ligand tested. Moreover, γ-glutamyl-S-allylcysteine can be considered as the best contender to serve as inhibitor of Mpro since this molecule showed maximum H-bond formation with Mpro along with stability corroborated by RMSD and RMSF values of a protein.

Figure 7.

Representation of (a) total energy changes (b) Changes in the temperature of the system (c) changes in the pressure for the systems 11a-Mpro, S-Allylcysteine-Mpro (c) γ-glutamyl-S-allylcysteine-Mpro .

MM/Gbsa binding free energy calculations

MD trajectory analysis of both the trajectories with (i) Mpro-11a and (ii) Mpro- γ-Glutamyl-S-allylcysteine were performed to understand the thermodynamic stability of the complexes which allows us to understand the spontaneity of the ligand receptor binding (Table 4). Binding free energy change calculation gives an estimate of that how strongly a ligand interact with the amino acids of the target protein. The energy released (ΔGbind) due to bond formation, or rather interaction of the ligand with protein is in the form of binding energy. This energy determines the stability of the protein–ligand complex under investigation. In our study, we observed that both ligands have ΔGBind in the negative range and shares the identical values. Moreover, both the compounds, tend to have similar energies corresponding to van der Waals interactions represented as ΔGvdW. These findings suggest that these compounds tends to stay in the vicinity of the interacting amino amides of Mpro. Further, both the compounds show negative values for Coulomb energy, suggests that these ligands while interacting with RBD has poor potential energy, suggesting better stability as the ligands do not have enough potential energy to get destabilized. In addition to the total energy, the contributions to the total energy from different components such as Lipophilic energy, Hydrogen-bonding correction, Pi-Pi packing correction and Van der Waals energy is depicted in Table 4. MM/GBSA analysis suggests that spontaneity of interaction forming capabilities of γ-Glutamyl-S-allylcysteine is at par with native ligand, 11a (well-established inhibitor).

Table 4.

MM/GBSA profiles of 11a, and γ-Glutamyl-S-allylcysteine, while interacting with Mpro.

Ligand	ΔGBind(Kcal/mol)	ΔGCoulomb(Kcal/mol)	ΔGHbond(Kcal/mol)	ΔGLipo(Kcal/mol)	ΔGPacking(Kcal/mol)	ΔGvdW(Kcal/mol)
11a	−61.834	−29.63	−2.59	−15.44	−2.25	−34.22
γ-Glutamyl-S-allylcysteine	−72.45	−28.71	−3.32	−20.53	−3.29	−37.67

ADMET analysis

All the ADMET properties of the top screened compounds along with GRL0617 is depicted in Table 3. For a compound to be classified as an oral drug, it is important to predict its mobility through the intestinal epithelial layers of cells that predicts the bioavailability. The theoretical model makes the use of Caco-2 permeability and its value higher than 0.90 means the compound has high permeability. Under present study all the compounds shows the Caco-2 permeability values in positive integer suggesting them to be absorbing through intestinal epithelial layers. The only exception is γ-Glutamyl-S-allylcysteine where the value is −0.517. Intestinal absorption (human) value is another parameter that calculates the absorption of the drug from human gut when administered orally. All the phytochemicals showed efficient intestinal absorption except γ-Glutamyl-S-allylcysteine. The next important parameter is skin permeability, where it was observed that all the compounds under study have values smaller than −2.5 log Kp, which means these compounds have poor permeability. All the compounds under study showed identical skin permeability. The ATP-binding cassette (ABC) transporter is important for transport of molecules through cell membrane and P-glycoprotein is its component needed for efficient transport, and its value ‘yes’ predicts the compound to pass cell membrane through ABC transporter. Here, γ-Glutamyl-S-allylcysteine and Chloroquineare are predicted to pass through the cell membrane via ABC transporters. Total diffusion of drug in total blood volume is determined by “volume of Distribution (VDss)” and its value below −0.15 logVDss suggest poor diffusion while values higher than 0.45 LogVDss suggests faster and higher equal distribution of drug in the total blood volume. The ability of compound to travel to the brain is given by the values of Blood-Brain barrier (BBB) permeability. When the logBB values are greater than 0.3, they can pass BBB. Values of compounds under study predicted that all the compounds may not be able to cross BBB and this attribute was only found to be present by Chloroquine. Further, where none of the compounds except Chloroquine is predicted to show CYP1A2, CYP2C19, CYP2C9, CYP2D6 and CYP3A4 inhibition as analysis obtained from metabolism prediction, while Chloroquine is predicted to inhibit only CYP2D6. Renal excretion rate of all the compounds under study were different and no compound except Chloroquine predicted to show AMES toxicity which important to predict as this test suggest drug’s property of mutagenicity. Potassium flux in heart is controlled by hERG I and II, improper flux of potassium by these transporters can cause QT syndrome where the Q and T peaks of heart electrocardiogram gets altered. The impact of screened compounds under present study on hERG I and II transporters are shown in Table 3 along with other essential ADMET properties.

Table 3.

ADMET Analysis of the Top Hits.

Property	Model Name	Predicted Value					Unit
		Chloroquine	Cuscohygrine	γ-Glutamyl-S- allylcysteine	Anahygrine	S-allylcystein
Absorption	Water solubility		−1.108	−2.891	−1.121	−2.888	Numeric (log mol/L)
Absorption	Caco2 permeability	1.624	1.364	−0.517	1.349	0.704	Numeric (log Papp in 10−6 cm/s)
Absorption	Intestinal absorption (human)	89.95	94.096	8.312	93.917	79.971	Numeric (% Absorbed)
Absorption	Skin Permeability	−2.679	−2.984	−2.735	−3.054	−2.736	Numeric (log Kp)
Absorption	P-glycoprotein substrate	Yes	No	Yes	No	No	Categorical (Yes/No)
Absorption	P-glycoprotein I inhibitor	No	No	No	No	No	Categorical (Yes/No)
Absorption	P-glycoprotein II inhibitor	No	No	No	No	No	Categorical (Yes/No)
Distribution	VDss (human)	1.332	0.979	−0.48	0.957	−0.561	Numeric (log L/kg)
Distribution	Fraction unbound (human)	0.191	0.757	0.452	0.755	0.444	Numeric (Fu)
Distribution	BBB permeability	0.349	0.236	−1.124	0.225	−0.277	Numeric (log BB)
Distribution	CNS permeability	−2.191	−3.226	−4.02	−3.237	−3.417	Numeric (log PS)
Metabolism	CYP2D6 substrate	Yes	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP3A4 substrate	Yes	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP1A2 inhibitor	No	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP2C19 inhibitor	No	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP2C9 inhibitor	No	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP2D6 inhibitor	Yes	No	No	No	No	Categorical (Yes/No)
Metabolism	CYP3A4 inhibitor	No	No	No	No	No	Categorical (Yes/No)
Excretion	Total Clearance	1.092	1.159	0.3	1.218	0.591	Numeric (log ml/min/kg)
Excretion	Renal OCT2 substrate	Yes	No	No	No	No	Categorical (Yes/No)
Toxicity	AMES toxicity	Yes	No	No	No	No	Categorical (Yes/No)
Toxicity	Max. tolerated dose (human)	−0.167	0.093	1.119	0.115	1.126	Numeric (log mg/kg/day)
Toxicity	hERG I inhibitor	No	No	No	No	No	Categorical (Yes/No)
Toxicity	hERG II inhibitor	Yes	No	No	No	No	Categorical (Yes/No)
Toxicity	Oral Rat Acute Toxicity (LD50)	2.85	2.385	2.438	2.407	2.02	Numeric (mol/kg)
Toxicity	Oral Rat Chronic Toxicity (LOAEL)	1.026	0.799	2.29	0.808	2.635	Numeric (log mg/kg_bw/day)
Toxicity	Hepatotoxicity	Yes	No	No	No	No	Categorical (Yes/No)
Toxicity	Skin Sensitization	No	Yes	No	Yes	No	Categorical (Yes/No)
Toxicity	T. pyriformis toxicity	1.558	0.291	0.285	0.227	0.166	Numeric (log ug/L)
Toxicity	Minnow toxicity	0.747	2.072	2.928	2.067	2.088	Numeric (log mM)

Overall discussion of the docking and molecular dynamics results

The recent outbreak of the COVID-19 or the SARS-CoV-2 has initiated a worldwide hunt for the therapeutic molecule effective against the viral infection [14,16,17]. One of the most effective strategy for therapeutics of SARS-CoV-2 is to target viral transcription and replication. A biomolecule which is capable of inhibiting the viral replication and transcription would be the ideal molecule to control spread of COVID-19. Mpro has been identified as the play vital role in the viral transcription and replication, and therefore, ant biomolecule capable of this protein could possibly prevent the viral transcription and replication [26].

Cuscohygrine and anahygrine are the two major biomolecules present in the W. somnifera. Cuscohygrine is primarily used as a biomarker to distinguish a practice of cocaine abuse from coca chewing [57] and has never been identified or recognized for its antiviral potential. Likewise, anahygrine is also not known for any major health benefits or for antiviral potential. W. somnifera is known for variety of health benefits on variety of ailments, one of which being a potent immune booster and immune stimulator [58-60]. Further, a traditional Indian formulation (Amukkara Choornam) based on W. somnifera has been reported to show beneficial effects against Chikungunya Virus [61], suggesting the antiviral potential of this plant. The broad health benefits of this plant are due to the presence of diverse chemical constituents such as cuscohygrine and anahygrine. Investigating the efficacy of cuscohygrine and anahygrine against the SARS-CoV-2 might yield beneficial effects via its action on SARS-CoV-2 Mpro, as indicated by the findings of our docking studies.

γ-glutamyl-S-allylcysteine, and S-allylcystein are the major bioactive molecules of A. sativum (garlic). Although γ-glutamyl-S-allylcysteine is not explored much for its pharmacological actions, S-allylcystein is known to possess anticancer, antidiabetic, nephroprotective, and neuroprotective potential [62-64]. Our findings suggest that both molecules are having good activity against SARS-CoV-2 Mpro and need further investigation. Moreover, A. sativum is also recognized for its variety of health benefits, which includes anti-microbial effects and immune booster potential [63,65].

In this study we demonstrated that cuscohygrine, γ-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein derived from Ashwagandha and Garlic binds to the active pocket of the Mpro of SARS-CoV-2. The binding was stabilized by hydrogen bonds, Pi/Pi-alkyl interactions, and favorable LE, LLE and Torsion values. Moreover, all four molecules formed a hydrogen bond with residue Glu166, a crucial residue for the formation of biologically active dimeric form of Mpro. Further, Pi/Pi-alkyl interactions with Cys145 also appears to be important for the interaction with Mpro of SARS-CoV-2 as this interaction was depicted by all molecules except cuscohygrine. Moreover, these results were confirmed through molecular dynamics simulations, where γ-glutamyl-S-allylcysteine demonstrated best results. Our findings suggest that in addition to the immunoboosting potential, these bioactive molecules could actively participate in combating the notorious coronavirus within the host body.

Conclusion

The present study was aimed at identifying the bioactive molecules for their potential against Mpro main protease of SARS-CoV-2 through molecular docking and molecular dynamics studies. In this study, we looked into the major bioactive molecules present in four medicinally important plants of India, viz. C. longa, A. sativum, O. tenuiflorum, and W. somnifera, which are recommended by the AYUSH, India for their beneficial effects against SARS-CoV-2 infection. Out of the 63 bioactive molecules, cuscohygrine, γ-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein demonstrated potential against Mpro main protease and their activity was observed to be better than the chloroquine. All these molecules demonstrated hydrogen bonding with the Glu166 residue, which is a crucial residue for the formation of biologically active dimeric form of Mpro, suggesting that these molecules can target Mpro very efficiently. Further, molecular dynamics simulation studies confirmed our results and predicted γ-glutamyl-S-allylcysteine as the best molecule which could act on the Mpro of SARS-CoV-2. Interestingly, all these molecules are not known for their potential against viral infections. Moreover, these molecules are present in good quantity in A. sativum and W. somnifera and both these plants are known to possess immunomodulatory, immune-booster and immune-stimulatory potential, which could be because of these molecules. Our results also validate the recommendations of AYUSH India for the use of A. sativum and W. somnifera for SARS-CoV-2 and suggest the presence of cuscohygrine, γ-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein in these plants may provide protection against SARS-CoV-2 infection by targeting Mpro. The proposed molecules, especially γ-glutamyl-S-allylcysteine, need to be investigated further for their potency against SARS-CoV-2. These molecules may find a clinical application, however, an intensive and rapid research is required.

Click here for additional data file.