Sars-cov-2 host entry and replication inhibitors from Indian ginseng: an in-silico approach.

COVID-19 has ravaged the world and is the greatest of pandemics in modern human history, in the absence of treatment or vaccine, the mortality and morbidity rates are very high. The present investigation identifies potential leads from the plant Withania somnifera (Indian ginseng), a well-known antiviral, immunomodulatory, anti-inflammatory and a potent antioxidant plant, using molecular docking and dynamics studies. Two different protein targets of SARS-CoV-2 namely NSP15 endoribonuclease and receptor binding domain of prefusion spike protein from SARS-CoV-2 were targeted. Molecular docking studies suggested Withanoside X and Quercetin glucoside from W. somnifera have favorable interactions at the binding site of selected proteins, that is, 6W01 and 6M0J. The top-ranked phytochemicals from docking studies, subjected to 100 ns molecular dynamics (MD) suggested Withanoside X with the highest binding free energy (ΔGbind = -89.42 kcal/mol) as the most promising inhibitor. During MD studies, the molecule optimizes its conformation for better fitting with the receptor active site justifying the high binding affinity. Based on proven therapeutic, that is, immunomodulatory, antioxidant and anti-inflammatory roles and plausible potential against n-CoV-2 proteins, Indian ginseng could be one of the alternatives as an antiviral agent in the treatment of COVID 19. Communicated by Ramaswamy H. Sarma.

Introduction

The whole world is currently facing the health crisis in the form of novel corona virus (2019-nCoV) outbreak and forced every nation to face challenges such as testing, quarantining and treating people affected by a coronavirus. Among the reported family of coronaviruses, 2019-nCoV is a novel strain not identified in humans previously (Ji et al., 2020). The initial outburst of 2019-nCoV in Wuhan spread briskly and greatly and affected other parts of China. The spread was so rapid that it took the shape of an epidemic in no time with coverage in several other countries of all the five continents Asia, Europe, Australia, Africa and the Americas (Wu et al., 2020). The range of 2019-nCoV epidemic is considerably wider than the severe acute respiratory syndrome (SARS) epidemic. World Health Organization (WHO) has declared 2019-nCoV as pandemic on 11 March 2020, and has also announced COVID-19 as a synonym for this new coronavirus disease. The widespread outbreak of a SARS-CoV-2 is declared as ‘Public Health Emergency of International Concern’ (PHEIC), affecting around 212 countries and almost half a million deaths worldwide. In this context, WHO has released the COVID-19 advice for the public and the same is available on the WHO website (WHO, 2020).

The coronaviruses, highly enveloped single-stranded RNA play a critical role in initial RNA synthesis of the infectious cycle, template for replication and transcription and also act as a substrate for packaging into the progeny virus. In all types of coronaviruses, two-thirds of the genome encodes a replicase polyprotein, pp1ab, that consists of two overlapping open reading frames (ORFs), that is, ORF1a and ORF1b and are processed by viral proteases to undergo cleavage into 16 different nonstructural proteins (NSPs), involved in transcription and replication (Boopathi et al., 2020; Cotten et al., 2013; Gupta et al., 2020). Another NSP that helps 2019-nCoV to get entry into the host cells involves a densely glycosylated, homotrimeric class I fusion spike (S) protein (Sarma et al., 2020). This S protein present in a metastable prefusion conformation undergoes structural rearrangement for the viral membrane to get fuses with the membrane of the host cell. The process is accelerated as a result of the binding of S1 subunit to the host-cell receptor and transition of the S2 subunit to a highly stable postfusion conformation (Hasan et al., 2020). Further, recent studies have reported that spike (S) glycoprotein possesses a human angiotensin-converting enzyme 2 (ACE2) binding site and has a 10 to 20-fold higher binding affinity towards 2019-nCoV S as compared with SARS-CoVs (Pillaiyar et al., 2020; Sinha et al., 2020).

The most effective methods to combat viral infection (specific vaccines and antiviral drugs) are being explored and developed and the world is waiting for the most effective treatment to combat COVID-19. Recent in-silico studies have reported the plausible role of new molecules and repurposing of existing drugs in the treatment of COVID-19 (Arya & Dwivedi, 2020; Beura & Prabhakar, 2020; A. Kumar et al., 2020; Mahanta et al., 2020). However, it may take months or years in the development of such effective treatments and thus exploration of prompt treatment options is also crucial. During the 2003 SARS outbreak, the effectiveness of herbal treatments was demonstrated . Therefore, it is the need of the hour to find the remedy for minimizing the morbidity and mortality due to COVID-19, based on complementary traditional medicines, since to date there is no official treatment available to treat COVID-19.

Withania somnifera (Solanaceae) (L.) Dunal, popularly known as ‘Ashwagandha’ and ‘Indian Ginseng’ is a prime medicinal plant of Ayurvedic and indigenous medicines from India and has been used as an herbal tonic and healthy food to treat various kinds of diseases and human ailments (Gurav & Gurav, 2014). W. somnifera contains alkaloids, flavonoids and steroidal lactones namely withanine, somniferine, somnine, somniferinine, withananine, pseudowithanine tropane, pseudo-tropine, choline, anaferine, anahydrine, isopelletierine, withaferin A and B, 27-deoxywithaferin, dihydrowithaferin A, 17-hydroxywithaferin A, withanolide A, B, C, D, E, G, J, L, M, N, O, P, Q, R and S, withanone, withanosine II, III, IV, V, VI, X and XI, Quercetin, and Quercetin-3-O-galactosyl-rhamnosyl-glucoside (QGRG; Abraham et al., 1968; Elsakka et al., 1990; Gupta & Rana, 2007; Krison & Glotter, 1980; Matsuda et al., 2001). Ashwagandha is widely claimed to have hepatoprotective, anxiolytic, antidepressant, nootropic, antimicrobial, anti-inflammatory, antioxidant, anti-stress, anticonvulsant, cardio-protective, antitumor, anti-genotoxic, anti-Parkinson and immunomodulatory properties (Dar et al., 2015; Gurav et al., 2020).

Ashwagandha holds an important place among the Ayurvedic Rasayana herbs (a preparation that works as a health tonic to children, a medicine to middle-aged persons and rejuvenator to the elderly). More importantly, the plant has been reported for its potent antiviral activity against different kinds of viruses such as H1N1 influenza, herpes simplex type-1, hepatitis, coxsackie virus, bursal disease virus, HIV, their infections and replications (Akram et al., 2018; Cai et al., 2015; Hattori et al., 1995; Kambizi et al., 2007; Mishra et al., 2013; Mukhtar et al., 2008; Pant et al., 2012). Further, studies have reported that drugs having early adaptogenic, immunomodulatory, anti-inflammatory and antioxidant potential along with its basic antiviral potential could prove quite effective against COVID 19 (Jayawardena et al., 2020; Sinha et al., 2020; Thevarajan et al., 2020; Tillu et al., 2020). A recent study also suggests that Ayurveda Rasayana such as Ashwagandha, that is, Indian Ginseng, can be a potential candidate for the management of COVID-19 and as a better and safer alternative to disease-modifying drugs such as HCQ (hydroxychloroquine; Patwardhan et al., 2020).

For complete viral replication 2019-nCoV requires four major structural proteins namely the spike (S) protein, nucleocapsid (N) protein, membrane (M) protein and the envelope (E) protein (Choudhury, 2020; Dewald & Burtram, 2019; Mittal et. al., 2020). Coronavirus spike (S) protein, a glycoprotein, uses ACE2 for entry of coronavirus into host cells (Walls et al., 2020). Here, the receptor-binding domain (RBD) of spike protein is involved in anchoring the protein on ACE2. Nucleocapsid (N) protein, an endoribonuclease, on the other hand, interacts with the viral genomic RNA and involved in viral replication and other mechanisms (Deng & Bak, 2018).

Recent in-silico studies have reported the plausible antiviral role of Withanone-N in COVID-19 by inhibiting the functional activity of SARS-CoV-2 protease Mpro (V. Kumar et al., 2020). Being Ayurvedic Rasayana, recently Indian Government (Ministry of AYUSH) along with Council of Scientific and Industrial Research (CSIR) and Indian Council of Medical Research (ICMR) has approved Ashwagandha for clinical trials against SARS-CoV-2 (https://www.nhp.gov.in). Thus, keeping the above view into consideration, the present investigation was undertaken to determine the efficacy of 48 different phytoconstituents from W. somnifera in comparison to reference drugs (hydroxychloroquine, lopinavir and remdesivir) against two different protein targets, that is, NSP15 endoribonuclease and prefusion spike RBD from SARS-CoV-2 using in-silico docking simulation and molecular dynamics (MD) study.

Material and methods

Molecular docking

Protein preparation

For the docking simulation, the crystal structures of NSP15 and spike protein were used. The crystal structure of NSP15 endoribonuclease of SARS CoV-2 (PDB ID: 6W01) with 1.9 Å resolution was retrieved from protein data bank (https://www.rcsb.org). This protein has bound citrate in its crystal structure. The crystal structure of SAR-CoV-2 spike protein with RBD bound with the ACE 2 (PDB ID: 6M0J) with 2.45 Å resolution was used in docking simulation. Overall protein inaccuracies were handled through the protein preparation wizard of Schrödinger maestro 2018-1 MM share version and missing side chains were modeled by the prime module. While preparing the protein for docking, hydrogen atoms were added, water molecules and other nonstandard residues were removed and partial charges were assigned by using the OPLS-2005 force field. Further, the protonation states of residues were assigned through PROPKA and the protein structure was subjected to restrained minimization with 0.3 Å root mean square deviation (RMSD). The sitemap module was employed to analyze the prospective binding sites. These binding sites were used in generating the grid box large enough to accommodate the structures of phytochemicals. For this, the dimension of the grid box was chosen as 20 × 20 × 20 Å3 and the docking simulation was carried out. The results of the docking simulations were analyzed from the analysis of the docking score, the types of interactions at the binding site residues.

Ligand preparation

The 2 D structures of a total of 48 bioactive molecules from W. somnifera were downloaded from the PubChem compound database (https://pubchem.ncbi.nlm.nih.gov) in SDF format. The 3 D structures and generation of lowest energy conformations of these phytochemicals were carried out with the Ligprep module of the maestro. The OPLS-2005 force field was employed while generating the stable conformations of the phytochemicals and the most stable conformations were used in docking studies.

Protein–ligand docking

Schrödinger Glide module was used during docking simulation of each bioactive molecule of W. somnifera into the binding site of respective optimized protein structures. The binding pose with the lowest docking score was retained and the docking results were analyzed using glide XP visualize.

MD simulation and molecular mechanics-generalized born solvent accessibility (MM-GBSA) analysis

System preparation

The top three binding energy scorer ligand to each protein, that is, QGRG, Withanoside X and Ashwagandhanolide for spike RBD and QGRG, Dihydrowithaferin A and Withanolide N for NSP15 Endoribonuclease were selected for further MD study. All the MD simulations were done on AMBER 18 software package (Lee et al., 2018). ANTECHAMBR was used for ligand preparation and to determine the charges on the ligand and further GAF force field was used for parametrization (Wang et al., 2001). Complexes of protein and ligand were prepared with the help of xleap. The SARS-CoV-2 spike protein RBD and NSP 15 endoribonuclease were solvated separately in truncated octahedron of TIP3P (Price & Brooks, 2004) box giving a total of 24,515 and 20,364 water molecules, respectively. A sufficient number of counter ions Na+ and Cl− were added to neutralize the simulation system and 0.1 M of ionic strength was achieved. To parameterize the amino acids and to model the proteins FF14SB force field was used (Maier et al., 2015).

Unbiased MD simulation

Simulations were performed for 100 ns of time step on Nvidia V100-SXM2-16GB Graphic Processing Unit using the PMEMD.CUDA module (Peramo, 2016). Simulations were run at 1 atm constant pressure using Monte Carlo barostat and 300 K constant temperature by using Longevin thermostat with a collision frequency of 2 ps−1 and the volume exchange was attempted for every 100 fs. An integration step of 2 fs was also used for simulation the hydrogen atoms involving bonds were constrained by using the SHAKE algorithm (Andersen, 1983). Long-range electrostatic interactions were computed by using Particle Mesh Ewald method while for short-range interaction a cutoff of 8 Å was used (Essmann et al., 1995). Equilibration consisted of rounds of NVT and NPT equilibration for 10 ns in total. CPPTRAJ (Roe & Cheatham, 2013) was used to analyze the interactions over full trajectory after taking configuration at every 4 ps. Dynamic behavior of the entire simulated systems was investigated in detail through various analyzing parameters, such as RMSD, Root Mean Square Fluctuations (RMSF), protein-ligand contacts, the radius of gyration and protein secondary structure elements (SSE). All the ligand RMSDs were graphically analyzed to see the stability of the ligand to the protein. RMSD, RMSF and MM-GBSA binding free energy (Rastelli et al., 2010) were determined after analyzing the trajectories.

MM-GBSA analysis

The MM-GBSA was performed on Amber18 and Amber18 tools. After simulation of the protein–ligand complexes, all the trajectories of 100 ns covering all the 10,000 frames were used for MM-GBSA analysis. The MM-GBSA based binding free energy (ΔGbind) calculations were performed on the molecular dynamics simulations (MDS) trajectories. The major energy components, such as H-bond interaction energy (ΔGbind_H-bond), Coulomb or electrostatics interaction energy (ΔGbind_Coul), covalent interaction energy (ΔGbind_Cov), lipophilic interaction energy (ΔGbind_Lipo), electrostatic solvation free energy (ΔGbind_Solv) and van der Waals interaction energy (ΔGbind_vdW) altogether contribute to the calculation of MM-GBSA-based relative binding affinity.

Results and discussion

Molecular docking studies

SARS-CoV-2 uses a homotrimeric class I fusion spike glycoprotein to make entry into the host cells ACE-2 receptor. This binding interaction is accelerated by binding of the S1 subunit to the host-cell receptor and by a transition of the S2 subunit to a highly stable postfusion conformation. Besides, the RBD of S1 subunit is composed of five twisted β sheets β1, β2, β3, β4 and β7 which are antiparallel to each other. An extended inclusion was found between β4 and β7 containing some α loops called receptor binding motif (RBM). This RBM contains most of the residues which are critical for connection between n-COVID-19 and ACE-2. A recent study has confirmed that among all the residues Arg 319 to Phe 541 only 17 residues Lys 417, Gly 446, Tyr 449, Tyr 453, Leu 455, Phe 456, Ala 475, Phe 486, Asn 487, Tyr 489, Gln 493, Gly 496, Gln 498, Thr 500, Asn 501, Gly 502 and Tyr 505 are crucial for binding with ACE-2 make contact with the 20 residues of N-terminal peptidase domain of ACE-2. Among these 17 residues Gln 493, Asn 501, Tyr 449, Tyr 489 and Tyr 505 are strongly connected with the help of H-bonding and Lys 417 by salt bridge interaction. NSP15 endoribonuclease bears a catalytic C-terminal domain, which constitutes a set of trimers, that is, hexamers and responsible for cutting the double-stranded (ds) RNA substrates with specificity via the Mn + 2-dependent endoribonuclease activity. Additionally, each monomeric unit is comprised of ∼345 amino acids which are folded into three domains namely N-terminal, middle and nidoviral RNA uridylate-specific endoribonuclease (NendoU) C-terminal catalytic domain, where the NendoU C-terminal catalytic domain comprises of two β-sheets which are antiparallel. These β-sheets contain six key amino acids namely His 235, His 250, Lys 290, Thr 341, Tyr 343 and Ser 294. Among them, His 235, His 250 and Lys 290 constitutes the catalytic triad, His 235 serves as a general acid, His 250 acts a base, while Ser 294 together with Tyr 343 have been found to govern U specificity. The binding site residues involved in the interaction with citrate are key residues. To check the potency and insights possible mechanism of 48 herbal-based ligands, the docking study was performed on two proteins.

The docking results of 2019-nCoV spike glycoprotein revealed that QGRG was properly positioned into the binding pocket surrounded by polar amino acid residues (Gln 493, Ser 494, Gln 498 and Asn 501), charged residues (Glu 406, Arg 403 and Lys 417), hydrophobic residues (Tyr 453, Tyr 495, Tyr 505, Phe 497, Leu 455 and Ile 418) and neutral amino acid residues (Gly 496 and Gly 502) with binding energy −9.246 kcal/mol. The key interactions at the binding site of spike glycoprotein are shown in Figure 1. The docking scores of top scorer W. somnifera bio-actives are given in Table 1. Whereas the docking scores of other phytochemicals are given in supplementary material Table S1.

Figure 1.

Binding-interaction analysis of (a) QGRG; (b) ashwagandhanolide; (c) Withanoside IV; (d) Withanoside X; and (e) withanolide III with SARS-CoV-2 spike RBD (PDB ID: 6M0J).

Table 1.

List of bio-actives with binding interaction parameter, that is, binding energy with the PDB: 6W01 and PDB: 6M0J of SARS-CoV-2.

Sr. No.	Pubchem CID	Name	Binding affinity kcal/mol
			6W01	6M0J
1	132519418	Quercetin-3-O-galactosyl-rhamnosyl-glucoside	−6.70	−9.25
2	15411208	Dihydrowithaferin A	−5.98	−2.82
3	23266147	Withanolide N	−5.97	−0.57
4	101168807	Withanoside X	−5.94	−7.07
5	443143	Anaferine	−5.82	−2.55
6	5280343	Quercetin	−5.79	−4.41
7	10700345	Withanoside V	−5.59	−4.66
8	442877	Withasomnine	−5.46	−2.57
9	92987	Pelletierine_Isopelletierine	−5.40	−2.50
10	8424	Tropine	−5.22	−2.54
11	16099532	Ashwagandhanolide	−4.08	−6.50
12	71312551	Withanoside IV	−4.79	−6.12
13	101168810	Withanoside III	−4.75	−5.76
14	10952344	Withanoside XI	−4.55	−5.65
15	3652	Hydroxychloroquine	−5.23	−3.57
16	92727	Lopinavir	−5.17	−4.22
17	121304016	Remdesivir	−5.02	−4.65

The structures of all W. somnifera phytochemicals investigated in the present investigation are given in supplementary material Figure S1.

Hydroxyl group of oxane ring (2H-tetrahydropyran ring) of QGRG exhibited H-bonding with Glu 406, Gln 493, Gly 496 and Gln 498 and hydroxyl group of terminal phenyl ring showed H-bonding with Gly 502. The phenyl ring of QGRG exhibited π–π interactions with Tyr 505. These types of interactions for all top ligands are tabulated in Table 2. Withanoside X, a glucoside of pubesenolide, showed a similar type of interactions with polar, charged and hydrophobic residues as shown in Figure 1.

Table 2.

Binding interaction of different bioactive herbal ligands with the active site of SARS-CoV-2 spike RBD (PDB ID: 6M0J).

		Interaction (PDB-6M0J)
			No. of H-Bond
S. No.	Ligand	H-Bond	HBA	HBD	π–π Stacking	Salt bridge
1	QGRG	Glu 406, Gln 493, Gly 496, Gln 498, Gly 502	1	5	Tyr 505	-
2	Withanoside X	Gln 409, Glu 406, Glu 484, Tyr 453	1	4	-	-
3	Ashwagandanolide	Asp 427, Arg 466, Pro 463, Glu 516	2	4	–	–
4	Withanoside IV	Gln 498, Ser 494, Gln 493, Arg 403, Gly 496, Tyr 505	2	4	–	–
5	Withanoside III	Gln 409, Glu 406, Lys 417, Arg 403	3	2	–	–
6	Hydroxychloroquine	Gln 409, Asn 501, Glu 406	1	2	Tyr 505	Glu 406
7	Lopinavir	Gln 493, Arg 403, Tyr 453	2	1	Tyr 505	–
8	Ramdisvir	Arg 403, Glu 406, Gln 493	1	3	Tyr 505	–

H-Bond, hydrogen bonding; HBD, hydrogen bond donner; HBA, hydrogen bond acceptor in respect to residues; QGRG, Quercetin-3-O-galactosyl-rhamnosyl-glucoside.

The binding energy of −7.07 kcal/mol is slightly higher than the binding free energy of QGRG, owing to the type of hydrogen bond interactions and hydrophobic interactions. Ashwagandhanolide, a dimeric thiowithanolide constitutes many hydrophobic and hydrophilic parts in its structure. It also shows many key interactions at the binding site with various amino acid residues.

The docking results of NSP15 Endoribonuclease revealed that QGRG has the most favorable interactions at the binding pocket surrounded by polar residues Gln 245, Thr 341 and His 338; charged residues Glu 340, Asp 240, Glu 234, Lys 335, Hip 250, Lys 290 and Hip 235 and hydrophobic residues Val 339, Trp 333, Met 219, Ala 232 and Gly 230, Gly 247 and Gly 248. The interactions QGRG and other ligands produce at the binding site are shown in Figure 2.

Figure 2.

Binding-interaction analysis of (a) QGRG; (b) Withanoside X; (c) anaferine; (d) dihydrowithaferin A; and (e) withanolide N with NSP15 endoribonuclease (PDB ID: 6W01).

The docking scores are given in Table 2 and the key interactions are given in Table 3. The higher docking score and corresponding binding free energy estimate of QGRG may be due to H-bonding of the hydroxyl group of oxane ring of QGRG with Glu 340, Gln 245, and the hydroxyl group of terminal phenyl ring showed H-bonding with Val 339. Dihydrowithaferin A, a non-glycosidic with anolide, was found to bind at the binding pocket similarly with a slightly lower docking score than QGRG. The higher docking score may be due to the propensity of hydrogen bond formation with protonated His 250 and His 235 residues and Thr 341 and Lys 290. In the case of Withanolide N, another non-glycosidic withanolide, the hydrogen bonds are formed with the same residues. Recent in-silico studies and our previous reports on NSP15 Endoribonuclease (6W01) and RBD of spike protein (6M0J) from SARS-CoV-2 evidenced the involvement of same key residues in the binding pocket (Agata & Piotr; Jonathan & Kevin, 2020; Sinha et al., 2020).

Table 3.

Binding interaction of different bioactive herbal ligands with the active site of SARS-CoV-2 Nsp15 endoribonuclease (PDB ID: 6W01).

		Interaction (PDB-6W01)
S. No.	Ligand	H-bond	No. of H-bond		π–π Stacking	Salt bridge
			HBA	HBD
1	QGRG	Glu 340, Gln 245, Val 339	2	3	–	–
2	Dihydrowithaferin A	Hip 250, Hip 235, Thr 341, Lys 290	4	–	–	–
3	Withanolide N	Hip 250, Hip 235, Thr 341, Lys 290	4	–	–	–
4	Withanoside X	Lys 290, Hip 250, Asp 220	2	2	–	–
5	Anaferine	Glu 261,Gly 239, Arg 258	1	2	–	Glu 261
6	Hydroxychloroquine	Glu 261, Gly 239	–	2	–	–
7	Lopinavir	Lys 290, Hip 235	2	–	His 243	–
8	Remdesivir	Hip 250, Hip 235, Thr 341	3	–	–	–

H-Bond, hydrogen bonding; HBD, hydrogen bond donner; HBA, hydrogen bond acceptor in respect to residues; QGRG, Quercetin-3-O-galactosyl-rhamnosyl-glucoside.

The docking results of 2019-nCoV spike glycoprotein with reference drugs revealed higher docking scores than the W. somnifera phytochemicals. This may be in part due to a lesser number of hydrogen bonds formed at the binding site. Specifically, Hydroxychloroquine, a 4-aminoquinoline derivative was found to bind at the binding pocket and forms bi-furcated H-bonding with Gln 409 and Glu 406. The quaternary nitrogen of hydroxychloroquine formed a salt bridge with Glu 406 and both aromatic and heteroaromatic ring of quinoline ring exhibited π-π interactions with Tyr 505 (Figure 3).

Figure 3.

Binding-interaction analysis of (a) hydroxychloroquine; (b) lopinavir; (c) remdesivir IV with NSP15 endoribonuclease (PDB ID: 6W01) and (d) hydroxychloroquine; (e) lopinavir; and (f) remdesivir with SARS-CoV-2 spike RBD (PDB ID: 6M0J).

The carboxyl group of amide in Lopinavir exhibited H-bonding with Gln 493, Arg 403 and nitrogen of amide showed H-bonding with Tyr 453. The phenyl ring of lopinavir exhibited π–π interactions with Tyr 505 and another phenyl ring showed Pi-cation interaction with Lys 417. In the case of remdesivir, the carbonitrile group showed H-bonding with Arg 403 and the hydroxyl group of tetrahydrofuran ring exhibited H-bonding with Glu 406. The nitrogen of phosphoramide depicted interactions with Gln 493. The phenyl ring of remdesivir exhibited π–π interactions with Tyr 505.

The docking results of NSP15 Endoribonuclease with reference drugs revealed higher docking scores than the W. somnifera phytochemicals. Particularly, hydroxychloroquine was found to bind at the binding pocket surrounded by charged residues Hip 235, Glu 261, Lys 257, Asp 240, Glu 234, Arg 258 hydrophobic residues Phe 241, Ala 218, Met 219, Tyr 238 and Gly 239 and the resultant binding energy −5.227 kcal/mol. Lopinavir and remdesivir were found to interact in near analogues way at the binding pocket but with consequent higher binding energies.

These results suggest that the QGRG is the best fit ligand to each selected protein among all the tested herbal ligands. QGRG has a dihydroxy chromenone ring, substituted with dihydroxy phenyl and connected with 2 trihydroxy oxane rings by an ether linkage. In the case of protein 6M0J, the binding free energy in terms of docking scores was found to be much lesser than all the reference drugs. The oxane ring of QGRG which is connected with other linkage was found to be properly buried into the outer concave surface of RBM and stabilized by 5 H-binding interaction, and most of the crucial amino acids of RBM was also found in the binding pocket of QGRG docked complex. A total of 10 residues out of 17 were found in the binding pocket, in which Gln 493, Gly 502 and Tyr 505 exhibited strong interaction with ligand. Withanoside X, Ashwagandhanolide, Withanoside IV and Withanoside III contain a dodecahydro-cyclo-pentaphenanthrene ring. Binding energy differs depending on the number and distribution of oxane rings and hydrophilic sites in these phytochemicals. All these ligands showed lesser binding free energy and more favorable binding interactions compared to the reference drugs. In the case of NSP15 protein, QGRG is the best withanolide in terms of its least binding free energy. Dihydrowithaferin A, Withanolide N and Withanoside X ranked 2nd, 3rd and 4th, respectively, but showed almost similar binding energy. Hydroxymethyl substituted oxenone ring of Withanolide N only showed 4 H-bonding interaction with the key amino acid of the binding site.

From the docking studies on these key 2019-nCoV viral proteins, the withanolide glycosides, containing D-glucose or D-galactose were found producing very favorable binding at the binding site due to various key hydrogen bond and hydrophobic interactions. Though docking studies give reasonable predictions of binding modes of ligands and estimates of binding free energy, the effects of biological environments such as aqueous medium, temperature, pressure and ionic composition of the system are not accounted for docking results. In such situations, MD studies can give accurate predictions of binding modes and estimates of binding free energies.

MD simulation

In biomolecular modeling studies, molecular docking combined with MDS of biological macromolecules is a successful and well-established method towards drug discovery and drug delivery process. It is a widely employed process for understanding structural stability and molecular interactions profile between protein and ligands or peptides. Therefore, in the present study 100 ns of MD simulations were performed for all selected best six docked complexes to observe how the interaction pattern of the binding site of NSP15 endoribonuclease and spike RBD from SARS CoV-2 adapts to the docked bioactive. The simulations trajectories provide the behavioral parameters, such as binding orientation and strength, and each bioactive-protein induced fit effects of individual particle motions as a function of time for each docked complex. Dynamic behavior of the entire simulated systems was investigated in detail through various analyzing parameters, such as RMSD, RMSF, protein–ligand contacts, the radius of gyration and protein SSE. All the ligand RMSDs were graphically analyzed to see the stability of the ligand to the protein. The protein–RMSD, protein–RMSF, ligand–RMSD are given in Figures 4–6, and protein and ligand RMSD values for selected complexes are given in Table 4 and binding free energy components for the protein-ligand complexes calculated by MM-GBSA analysis are given in Table 5.

Table 4.

Protein and ligand RMSD values of six selected protein-ligand complexes.

Bioactive compounds	RMSD (Å)			Ligand RMSD (Å)
	Protein backbone atoms			Ligands atoms
	Maxa	Minb	Avgc	Maxa	Minb	Avgc
Ashwagandhanolide	2.65	0.63	1.77	4.62	0.39	2.14
Dihydrowithferin A	3.82	0.64	2.36	2.36	0.16	0.91
QGRG	2.55	0.63	1.85	8.74	0.72	3.62
Quercetin	3.82	0.61	2.24	2.88	0.21	1.54
Withanolide N	2.86	0.62	2.10	2.18	0.23	0.91
Withanoside X	3.55	0.63	2.75	16.46	0.66	12.17

QGRG, Quercetin-3-O-galactosyl-rhamnosyl-glucoside.

aMaximum.

bMinimum.

cAverage.

Table 5.

Binding free energy components for the protein ligand complexes calculated by MM-GBSA analysis.

Compounds	MM-GBSAa
	ΔEVDW	ΔEELE	ΔGGB	ΔGSurf	ΔGgas	ΔGSol	ΔGbind
Ashwagandhanolide	−48.56 (4.87)	−35.56 (8.54)	52.69 (6.83)	−6.48 (0.47)	−88.30 (8.40)	46.20 (6.83)	−42.09 (5.02)
Dihydrowithferine A	−30.93 (2.67)	−21.80 (4.33)	29.01 (3.14)	−3.94 (0.23)	−54.46 (4.44)	25.07 (3.13)	−29.38 (2.81)
QGRG	−28.14 (3.63)	−64.14 (16.19)	71.17 (11.64)	−5.96 (0.53)	−114.62 (17.09)	65.20 (11.20)	−49.41 (6.88)
Quercetin	−15.45 (3.49)	−52.42 (16.33)	56.93 (12.22)	−3.80 (0.41)	−89.67 (16.66)	53.13 (11.94)	−36.54 (5.68)
Withanolide N	−34.16 (2.82)	−36.47 (7.66)	45.20 (6.94)	−4.45 (0.29)	−71.29 (8.37)	40.75 (6.80)	−30.53 (3.38)
Withanoside X	−19.38 (6.29)	−16.91 (10.39)	30.29 (10.05)	−2.90 (0.98)	−116.81 (19.89)	27.38 (9.50)	−89.42 (16.09)

QGRG, Quercetin-3-O-galactosyl-rhamnosyl-glucoside; ΔEVDW, van der Waals contribution from MM; ΔEELE, electrostatic energy as calculated by the MM force field; ΔGGB, the electrostatic contribution to the solvation free energy calculated by GB; ΔGSurf, solvent-accessible surface area; ΔGSol, solvation free energy; ΔGgas, gas phase interaction energy; ΔGbind, binding free energy.

aAll energies are in kcal/mol with standard deviation in parenthesis.

The analysis was performed on the molecular interactions before and after the MDS production step. Before the MD studies, dihydrowithaferin A shows to form two hydrogen bonds between the O1 atom of oxypentacyclo ring with the backbone nitrogen of Gly291 and Gly249 (Figure 4(d)).

Figure 4.

SARS-CoV-2 spike RBD (PDB:6M0J) interactions: (a) protein–RMSD; (b) protein per residue RMSF during simulation; (c) ligand–RMSD during MD simulation; (d) pre- and post-MDS interaction of Ashwagandhanolide with the spike protein; (e) pre- and post-MDS interaction of QGRG with the spike protein; and (f) pre- and post-MDS interaction of Withanoside X with the spike protein.

Towards the end of the simulation, the compound optimizes the conformation with an average ligand RMSD of 0.91 Å (Figure 4(a), Table 4) to form a hydrogen bond between the hydroxyl group of the oxypentacyclo ring and residues Thr342 and His251. The ligand RMSD shows a stable transition initially till 30 ns after the RMSD fluctuated by 2 Å between 30 and 65 ns after which the RMSD was stable for the rest of the simulation. Quercetin shows various intermolecular interactions with the endoribonuclease residues owing to the presence of sugar moieties. However, the most significant interaction is observed between Lys 346 and the O17 atom chromenone ring of Quercetin (Figure 4(e)). This hydrogen bond is a key interaction between the ligand and receptor which is evident from the bond length; at the beginning, the bond length was 2.02 Å and fluctuates to 2.39 by the end of MDS simulation. The average ligand RMSD of 1.54 Å also supports this finding (Figure 4(a), Table 4). The compound Withanolide N, after minimization and equilibration steps of MDS, show a hydrogen bond with the residue Ser295 of the endoribonuclease and O3 of the phenanthrene ring structure. The bond length was 2.19 Å at the beginning of MD production and shortens to 2.04 Å, the ligand also forms a new hydrogen bond (2.43 Å) between Gly249 and O4 of pyran ring. The ligand RMSD was stable throughout the MD production phase with low RMSD fluctuation at an average of 0.91 Å. This indicates that Withanolide N forms a very stable complex with an endonuclease. The protein RMSD and RMSF for these complexes do not show a huge fluctuation individually and compared to each other. Figure 4(b) and (c) shows protein RMSD stays low throughout the MDS between 0.5 and 2.75 indicating the overall stability of dihydrowithaferin A, Quercetin and Withanolide N bound to the SARS-CoV-2 NSP15 endoribonuclease complexes.

The phytochemicals Ashwagandhanolide, QGRG and Withanoside X in complex with the SARS-CoV-2 spike protein (PDB: 6M0J) were selected for the MDS based on molecular docking study results. The spike protein-Ashwagandhanolide complex shows three main hydrogen bond interactions at the pre-production stage: Tyr495 (O—O13; 1.86 Å), Lys417 (N—O6; 1.88 Å) and Arg403 (NH1—O4; 2.05 Å).

The hydrogen bonds Lys417 (N—O6; 1.98 Å) and Arg403 (NH1—O4; 1.83 Å) were preserved during the MDS, the Tyr 495 interaction broke down a weaker but new hydrogen bond was formed at the end of the MDS, Asn501 (N—O13; 2.98 Å; Figure 5(d)). The ligand RMSD remains stable for the first 30 ns but then fluctuates between 1.5 and 4.5 Å with an average of 2.14 Å. The QGRG and spike protein complex shows several intermolecular interactions but it fails to present any strong hydrogen bonds with the backbone of the receptor. The sugar moieties of QGRG forms hydrogen bonds with Ser494, Gly496, Lln493, Arg403, Gln498, Asn501 and Gln493 of the spike protein, these interactions are found to stabilize this complex throughout the MDS production step. The ligand RMSD remains stable up to 50 ns with little fluctuations (Figure 5(c)), however towards the end of the simulation the fluctuations are quite high, and this is also reflected in the protein RMSD and RMSF (Figure 5(a) and (b)). The Withanoside X complex with spike protein is interesting (Figures 5(f) and 6).

Figure 5.

SARS CoV-2 NSP15 endoribonuclease (PDB:6W01) protein interactions: (a) protein–RMSD; (b) protein per residue RMSF during simulation; (c) ligand–RMSD during MD simulation; (d) pre- and post-MDS interaction of dihydrowithaferin A with the spike protein; (e) pre- and post-MDS interaction of Quercetin with the spike protein; and (f) pre- and post-MDS interaction of withanolide N with the spike protein.

Figure 6.

MDS analysis for the Withanoside X bound to spike protein, the Withanoside X changes conformation to obtain higher conformational stability and in-process changes the conformation showing high binding energy.

After the equilibration step, the Withanoside X complex was analyzed, it shows a hydrogen bond between the side chain of Lys417 and O5 of the ligand. During the MDS this ligand undergoes conformational changes and thus optimizes its binding o the receptor. Figure 6 shows the transition of Withanoside X during the MDS and corresponding ligand RMSD. Initially, the RMSD increases to 15 Å during the first 40 ns of simulation, later it stabilizes and stays around 15 Å for all of the 60 ns. The protein RMSD and RMSF also show the same trend during the simulation (Figure 5(a) and (b)).

MM-GBSA analysis

MM-GBSA analysis was performed on all of the six protein-ligand complexes to evaluate the affinity of ligands to the target protein receptors. The MM-GBSA based binding free energy (ΔGbind) calculations were performed on the MDS trajectories. The binding energies measured by this method are more efficient than the GlideScore values for the selection of protein–ligand complexes. The major energy components, such as H-bond interaction energy (ΔGbind_H-bond), coulomb or electrostatics interaction energy (ΔGbind_Coul), covalent interaction energy (ΔGbind_Cov), lipophilic interaction energy (ΔGbind_Lipo), electrostatic solvation free energy (ΔGbind_Solv) and van der Waals interaction energy (ΔGbind_vdW) altogether contribute to the calculation of MM-GBSA-based relative binding affinity. The binding energies and the contributing factors calculated for the MDS trajectories are mentioned in Table 5. Out of the six complexes studied, three complexes showed high binding free energies. The group of compounds that bound to the SARS-CoV-2 spike RBED, QGRG and Ashwgandhanolide showed high binding free energies. The QGRG showed ΔGbind of −49.41 kcal/mol with a high contribution from the ΔGgas (−114.62 kcal/mol) and ΔGGB (71.17 kcal/mol). Ashwagandhanolide had the second-highest binding free energy (ΔGbind of −42.09 kcal/mol). The ΔGSurf (−6.48 kcal/mol) for Ashwagandhanolide indicates its better surface accessibility compared to other molecules. In the second group of compounds binding the SARS-CoV-2 NSP15 endoribonuclease, Withanoside X shows a very good binding free energy (ΔGbind = −89.42 kcal/mol). The high binding affinity is also justified by the fact that during MDS the molecule optimizes its conformation for better fitting with the receptor active site.

Conclusion

The antiviral activity, adaptogenic and immunomodulatory potential of Ashwagandha is well documented in the literature. In present work, the possible inhibitory potential of phytochemicals from Ashwagandha was analyzed through in silico methods. The docking studies and the MD pointed out the possible lead-like properties to some phytoconstituents from Ashwagandha. Specifically, the phytochemicals/bio-actives such as QGRG, Withanoside X, Ashwagandanolide, Dihydrowithaferin A and Withanolide N hold promise in inhibiting the SARS-CoV-2 key viral proteins. The present study could be the starting point for the future ligands from natural sources in 2019-nCoV spike glycoprotein and NSP15 endoribonuclease.

Author contributions

SSG was involved in conceptualization, supervision, investigation; RBP, SKS and RVC were involved in methodology, software data curation and software validation; SKP and AS were involved in writing—original draft preparation; SKS, NSG and RSP contributed for formal analysis, visualization, illustration, reviewing and editing.