{'fraction_non_alphanumeric': 0.03434290229128259, 'fraction_numerical': 0.016697392678430338, 'mean_word_length': 5.085096153846154, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reckoning γ-Glutamyl-S-allylcysteine as a potential main protease (mpro) inhibitor of novel SARS-CoV-2 virus identified using docking and molecular dynamics simulation

Abstract 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, a-coronaviruses, b-coronaviruses and c-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  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 23 rd 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 (M pro ), Nonstructural protein 10 (nsp10), nsp12, helicase (nsp13), N-terminal exoribonuclease and C-terminal guanine-N7 methyl transferase (nsp14), Uridylate-specific endoribonuclease (nsp15), 2 0 -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, antiinflammatory, and immune-boosting potential [20][21][22][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 M pro of SARS-CoV-2. The M pro of SARS-CoV-2 is crucial for the viral gene replication and gene expression via proteolytic processing of replicase polyproteins. Targeting M pro can result in the inhibition of viral replication and therefore is an attractive target for drug design against SARS-CoV-2 [26][27][28][29].

Material and methods
Datasets The x-ray crystal structure of the M pro of SARS-CoV-2 was selected from the PDB database. We selected the crystal structure of M pro (PDB Id: 6M0K) [26] with resolution 1.50 Å. The standard inhibitor for this binding cleft is identified as the 11a (fNg- 11a has already been established as the inhibitor of M pro and this was used as a positive control in our study. M pro 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 M pro , 11a. The binding affinities of the selected bioactive molecule with the M pro 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 < lM < 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 typespecific 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 Hbond 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 DG 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 M pro in the presence of 11a using GROMACS 2019 software [37][38][39]. All the test docked complexes were compared with the 11a-M pro 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 Clcounter 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 M pro protein, and radius of gyration (Rg) by using the 'gmx rmsf', 'gmx rms', 'gmx hbond' and gmx gyrate'of GROMACS utilities [37][38][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][44][45][46]. MM/GBSA in PRIME module of Maestro 11.4 was used to calculate the thermodynamic data of Coulomb energy (DGCoulomb), total free energy change (DGBind), Hydrogen-bonding correction (DGHbond), Pi-Pi packing correction (DGPacking), Lipophilic energy (DGLipo), and Van der Waals energy (DGvdW) 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 D TDS is the conformation entropic contribution, and DG MM is the molecular mechanics' interaction energy (electrostatic þ van der Waals interaction) between protein and ligand. DG PB and DG SA 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, c-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-Ogallate, possess better binding affinity toward M pro 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 M pro . 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 M pro .
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 M pro 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 M pro protein [53][54][55][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, c-Glutamyl-S-allylcysteine, anahygrine, and S-allylcystein were selected on the basis of their estimated binding affinity against M pro of SARS-CoV-2 (Table S2 and Table 1), and their chemical structures are depicted in Table 2.
SARS-CoV-2 M pro is having a total of four binding sites, S1 0 , S1, S2, and S4, and the 11a molecule is known to inhibit the activity of the SARS-CoV-2 M pro by interacting at the S1 0 site of this protein [26]. Docking analysis for the test molecules was performed on the coordinates corresponding to S1 0 site of M pro . We observed that 11a molecule is capable of interacting with the multiple amino acids at the S1 0 site of M pro , 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 M pro , as demonstrated in our results for the 11a-M pro 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. Based on the estimated affinity of Chloroquine (553.5 lM À 54.9 mM), we selected four molecules, viz. cuscohygrine, c-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 M pro are depicted in Figure 2. Cuscohygrine, c-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein binds to the active pocket of M pro with an estimated affinity between 1.20 lM À 120.1 lM, 18.7 lM À 1.86 mM, 14.3 lM-1.42 mM, and 313.6 lM-31.5 mM respectively. The LLE and LE parameters of molecular docking for all molecules were higher than chloroquine. LLE and LE parameters for c-glutamyl-S-allylcysteine were observed to be higher than the rest of the molecules in the binding site of M pro . 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 c-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 M pro . c-glutamyl-S-allylcysteine showed the most optimum interaction with M pro 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).

Molecular dynamics simulation
We conducted molecular dynamics simulations for different protein-ligand complexes, which includes positive control 11a-M pro , S-Allylcysteine-M pro , c-glutamyl-S-allylcysteine-M pro individually and their results were compared. S-Allylcysteine and c-glutamyl-Sallylcysteine were selected for molecular dynamics simulation based on the docking interaction observed with the M pro . Our results demonstrated the good overall stability of the M pro proteins in presence of the inhibitors being tested, S-Allylcysteine and c-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 ligandprotein 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 M pro protein backbone against the simulation time scale (0-25) ns) during its interaction with the internal ligand 11a, S-Allylcysteine and c-glutamyl-S-allylcysteine (Figure 4).
The RMSD values of 11a-M pro trajectories was least and could be sought to be the most stable interaction which ranged below  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-M pro (as a positive control) and (ii) S-Allylcysteine-M pro and (iii) c-glutamyl-Sallylcysteine-M pro 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.
In the present study, we observed almost identical RMSF of the M pro 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 M pro was determined throughout the simulation at 300 K for the ligand system ( Figure 6). The 11a-M pro complex exhibited maximum of five H-bonds formation during the course of simulation. While for S-Allylcysteine-M pro , 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 c-glutamyl-S-allylcysteine, where this ligand could consistently make four to five H-bonds with M pro 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 M pro , but also make appropriate hydrogen bonds with the protein. Further, c-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 M pro as made by 11a. Therefore the ratio of hydrogen bond to surface area is much greater for c-glutamyl-Sallylcysteine compared to 11a, suggesting very strong interaction with M pro in the S1' cleft. c-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.
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-M pro , S-Allylcysteine-M pro , and c-glutamyl-S-allylcysteine-M pro . Further, Figure 7(b) suggests the deviation in the temperature from the constant 300 K, where it is evident that the 11a-M pro system has maximum deviation, this can be due to the large size of the ligand (11a) attached to the protein (M pro ). Since 11a is relatively larger in size and has more possible interactions with M pro , 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, c-glutamyl-Sallylcysteine can be considered as the best contender to serve as inhibitor of M pro since this molecule showed maximum H-bond formation with M pro along with stability corroborated by RMSD and RMSF values of a protein.  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 (DGbind) 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 DGBind 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 DGvdW. These findings suggest that these compounds tends to stay in the vicinity of the interacting amino amides of M pro . 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 c-Glutamyl-S-allylcysteine is at par with native ligand, 11a (wellestablished inhibitor).

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 c-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 c-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, c-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   Table 3 along with other essential ADMET properties.
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. M pro 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][59][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 M pro , as indicated by the findings of our docking studies.
c-glutamyl-S-allylcysteine, and S-allylcystein are the major bioactive molecules of A. sativum (garlic). Although c-glutamyl-Sallylcysteine is not explored much for its pharmacological actions, S-allylcystein is known to possess anticancer, antidiabetic, nephroprotective, and neuroprotective potential [62][63][64]. Our findings suggest that both molecules are having good activity against SARS-CoV-2 M pro 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, c-glutamyl-Sallylcysteine, anahygrine, and S-allylcystein derived from Ashwagandha and Garlic binds to the active pocket of the M pro of SARS-CoV-2. The binding was stabilized by hydrogen bonds, Pi/Pialkyl 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 M pro . Further, Pi/Pi-alkyl interactions with Cys145 also appears to be important for the interaction with M pro of SARS-CoV-2 as this interaction was depicted by all molecules except cuscohygrine. Moreover, these results were confirmed through molecular dynamics simulations, where c-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 M pro 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, c-glutamyl-S-allylcysteine, anahygrine, and S-allylcystein demonstrated potential against M pro 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 M pro , suggesting that these molecules can target M pro very efficiently. Further, molecular dynamics simulation studies confirmed our results and predicted c-glutamyl-S-allylcysteine as the best molecule which could act on the M pro 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, c-glutamyl-Sallylcysteine, anahygrine, and S-allylcystein in these plants may provide protection against SARS-CoV-2 infection by targeting M pro . The proposed molecules, especially c-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.