Introduction

Coronavirus disease-19 (COVID-19) has become pandemic in humans by targeting the severe acute respiratory syndrome-coronavirus-2 (SARS-CoV-2) entry and replication in the viral host. Numerous studies have reported the potential druggable targets to control the SARS-CoV-2 infection (Gupta et al. 2021; Wu et al. 2020; Kushwaha et al. 2021a; Maurya et al. 2020). SARS-CoV-2 interacts with host cellular membrane ACE2 receptors (angiotensin converting enzyme-2) and TMPRSS2 co-receptors (trans-membrane serine protease-2) with its spike glycoprotein and initiate host-pathogen interaction. One of the Spike-protein’s co-receptor viz. DPP4/CD26 (dipeptidyl peptidase-4) involved in COVID-19 pathophysiology has received much attention. Solerte et al. (2020) determined that CD26 inhibitors may helpful to reduce the viral entry, replication, cytokine storm and inflammation in COVID-19 patients. Recently, genetics-based study from our group revealed the coding and regulatory variations of TMPRSS2 and CD26 in prolonging SARS-CoV-2 infection in humans (Senapati et al. 2020). The spike glycoprotein S1 and S2 domains involved in the protein-protein interaction and membrane fusion establish viral entry into the cell. Binding of cellular receptor to the S1 domain induce conformational change in the S2 domain which create cleavage site for TMPRSS2, a human serine protease. TMPRSS2-induced cleavage generates S-fragments and releases them into the cellular supernatant. After processing in the early and late endosome, viral-genome is being replicated with the support of RNA dependent RNA polymerase-RdRp and transcribed into viral mRNA. Following translation, synthesized viral proteins undergoes maturation and generation of viral particle. Targeting entry of virus particle into the cell by inhibiting spike and/or host membrane receptor/co-receptors is a better strategy to develop SARS-CoV-2 anti-viral drug discovery. Moreover, targeting TMPRSS2 and RdRp proteins may decrease the SARS-CoV-2 mediated virulence and host pathogenicity.

Type II transmembrane serine proteases is an important protein family involved in the pathophysiology of cancer and viral disease. TMPRSS2 is a member of hepsin/transmembrane protease, serine (TMPRSS) subfamily. TMPRSS2 is known to prime S2 domain of SARS-CoV-2 which is a key step in viral infection and virulence. Various in vitro, and in vivo studies showed TMPRSS2 inhibition mediated decrease in coronavirus (SARS-CoV, SARS-CoV-2 and MERS-CoV) induced immunopathology, virulence and airway spread of the disease (Iwata-Yoshikawa et al. 2019; Hoffmann et al. 2020) reported that a clinically proven TMPRSS2 inhibitor has potential to inhibit viral infection in SARS-CoV-2 infected lung cells. Cluster of differentiation 26 (CD26), also known as dipeptidyl peptidase-4 (DPP4) is a human cell surface glycoprotein having dipeptidyl peptidase activity at extracellular domain. Besides, it is also present in soluble form in body fluids and possesses significant serine exopeptidase activity (Casrouge et al. 2018). DPP4 inhibitors (DPP4i) are clinically approved antidiabetic drug for type 2 diabetes mellitus. Pharmacologically active CD26 inhibitors are known to decrease pro-inflammatory cytokines and lower the incidence risk of autoimmune diseases (Kim et al. 2015). Acute respiratory failure and deregulated inflammation is the major cause related to SARS-CoV-2 death. Reports suggest that CD26 inhibitor sitagliptin could be utilized as a therapeutic candidate for acute respiratory distress syndrome (ARDS). Clinically approved CD26 inhibitors including vildagliptin and saxagliptin either act as pseudo-substrates or substrate competitive inhibitors that binds to the same catalytic site of CD26 receptor (Berger et al. 2017). Importantly, the active site of the receptor is different from the viral-human (Spike glycoprotein-CD26) protein-protein interaction interface (Lu et al. 2013; Vankadari and Wilce 2020), where CD26 interacts with S1 domain of SARS-CoV-2 spike glycoprotein (Vankadari and Wilce 2020). SARS-CoV-2 like other coronaviruses utilizes CD26 other than ACE2 as a functional receptor to enter into the host cell.

Natural products have been used overtime for the treatment of various diseases (Kushwaha et al. 2019, 2020a, b, 2021b; Kumar et al. 2019; Waseem et al. 2017). Phytochemicals present in medicinal plants are known to possess antiviral potential by targeting different proteins of Ebola, Herpes Simplex Virus 1 and 2, Human Immunodeficiency Virus, Adenovirus, Respiratory Syncytial virus, including SARS-CoV-2 virus (Li et al. 2005; Naveed et al. 2018; Kushwaha et al. 2021c). Withania somnifera (Solanaceae) is a well-known Indian medicinal plant. Phytochemicals present in W. somnifera possess DPP4 inhibitory and immunomodulatory potential (Kempegowda et al. 2018). Besides, the plant is well documented for its antiviral potential against several viruses (such as H1N1 influenza and HIV virus). Research demonstrates that W. somnifera phytochemicals inhibit viral protein, alters inflammation and immune system (Cai et al. 2015; Jain et al. 2018; Ganguly et al. 2018). Previously, our team identified W. somnifera phytochemicals having potential to bind SARS-CoV-2 main protease. In the present study, for the first time, we studied the SARS-CoV-2 (Spike-protein, RdRp) and human host protein (CD26 and TMPRSS2) binding efficacy of W. somnifera phytoconstituents. These phytochemicals were docked against test proteins to identify their binding potential in comparison to respective standard inhibitors. Further, the MM-GBSA and molecular dynamics (MD) simulation was performed to identify the binding free energy, stability and protein conformation of the lead phytochemical-protein complexes. Moreover, the druglikeness and physiochemical properties of the lead compounds were also studied.

Materials and methods

Withania somnifera phytochemical retrieval and preparation

.

Phytochemicals present in Withania somnifera were retrieved from previous published literature Kushwaha et al. 2021c.) The 3D or 2D structures of standard compounds against respective targeted proteins were retrieved from the NCBI PubChem in .sdf format. The .sdf structure files were converted into .mol to get the respective 3D conformation of the compound with the assistance of Open babel molecule format converter (O’Boyle et al. 2011). Energy of the test ligands were minimized by applying mmff94 force field and conjugate gradients optimization algorithm using PyRx-Python prescription 0.8 for 200 steps (Dallakyan 2008).

Protein retrieval and preparation

Crystal structure of target proteins viz., CD26 (PDB ID: 4PNZ), Spike glycoprotein (PDB ID: 6VSB), and RdRp (PDB ID: 6M71) were obtained from Protein Data Bank (https://www.rcsb.org/) with a resolution of 1.90 Å, 3.46 Å, and 2.90 Å (Berman et al. 2000). TMPRSS2 modeled protein (PDB ID: 5CEL was used as template) was adapted from our previous publication (Kushwaha et al. 2021c). Three dimensional structure of the targeted proteins were prepared for molecular docking using UCSF Chimera (Pettersen et al. 2004). These proteins were cleaned and optimized by removing the ligands and other heteroatoms including water. Following this step, the energy minimization of the proteins was performed using steepest descent method having 100 steps (step size 0.02 Å), and a conjugate gradient method with 10 steps (step size 0.02 Å) using UCSF Chimera.

Molecular docking studies

For docking analysis, receptors and ligands were loaded in Auto Dock Tools 1.5.6 (ADT) (Trott and Olson 2010; Singh et al. 2021; El Khatabi et al. 2021; Verma et al. 2020). Following the merging of non-polar hydrogens and torsions applied to the ligands by rotating all rotatable bonds, the Gestgeiger partial charges were assigned. Docking calculations were performed for the entire protein model. Auto Dock tools were used to assign polar hydrogen atoms, Kollman charges, and solvation parameters. For the exploration of the space of active binding with different efficacy, the Lamarckian Genetic Algorithm option of Auto Dock 4.2 was used. The grid boxes were prepared for the entire binding site (Spike glycoprotein) or the active site (CD26, RdRp and TMPRSS2) of the target proteins to provide sufficient space for the ligand’s translational and rotational walk. A maximum number of 27,000 GA operations were generated on a single population of 150 individuals for each of the 30 independent runs. UCSF Chimera (v.1.10.2) was used for visualization and LigPlot+ (v.1.4.5) used for the analysis of interface between receptors and ligands (Laskowski et al. 2011).

MM-GBSA analysis

The binding free energies of the targeted protein(s) complexed with the lead phytochemicals were calculated by using Prime MM-GBSA (Molecular Mechanics/Generalized Born Model and Solvent Accessibility) calculations. The calculation includes the VSGB solvent model, OPLS_2005 force field, and rotamer search algorithms. The docked PDB file was used to Prime MM-GBSA simulation and subsequent computation of the binding free energy.

Molecular dynamics simulation

All-atom MD simulations were performed for TMPRSS2 and RdRp protein in unbound and lead phytochemical/standard inhibitor bound complex. Simulations were performed using GROMOS 54A7 force field in the GROMACS 5.1.1 package (Van Der Spoel et al. 2005; Ezebuo et al. 2019; Kushwaha et al. 2021c). Ligand topologies and parameters generated using PRODRG server and were combined with the protein topology to make simulation complexes (Schüttelkopf and van Aalten 2004; Singh et al. 2020a, b). All systems were solvated using the SPC water model in a cubic box. Energy minimization was performed to get rid of steric clashes by using 1500 steps of the steepest descent approach. Following energy minimization, the temperature of all systems was increased from 0 to 300 K in equilibrium phase of 100 ps at a constant volume and stable pressure of 1 bar. The final MD run was set to 10,000 ps for all six systems. Trajectory analysis was performed using various GROMACS utilities such as gmx rms, gmx rmsf, gmx gyrate, gmx sasa, gmx hbond and gmx do_dssp for the analysis of rmsd, rmsf, radius of gyration, solvent accessible surface area, hydrogen bonds and secondary structure respectively.

Pharmacological and physiochemical property prediction

Swiss ADME (Swiss Institute of Bioinformatics, Switzerland), a widely used web tool was utilized for the prediction of physicochemical properties, pharmacokinetics, drug likeness and water solubility by bioavailability radar generation. Bioavailability radar confers six physicochemical properties i.e. flexibility, lipophilicity, saturation, size, polarity and solubility. The Egan BOILED-Egg method of SwissADME tool was applied to determine the absorption of the inhibitors in the gastrointestinal tract and brain (Daina et al. 2017). BOILED-Egg (Brain or Intestinal Estimated permeation predictive model), also called Egan egg, was used to predict a threshold and a clear graphical representation of how far a molecular structure is from the ideal one for good absorption. In 2D graphical representation, the molecules present in yolk sac will passively permeate through the blood–brain barrier (BBB), whereas the molecules present in the white region are expected to be passively absorbed by the gastrointestinal (GI) tract.

Result

In the present study, a total of 93 Withania sominfera compounds were searched from the literature and docked against SARS-CoV-2 (spike glycoprotein and RdRp) and human host proteins (CD26 and TMPRSS2). Binding energy of the lead compounds (cut-off value ≤− 4.5 kcal/mole) against targeted test proteins is shown in Table 1. Two dimensional structure of the lead compounds against all targeted protein with cut-off value  ≤ binding energy of standard inhibitors are depicted in Fig. 1.

Fig. 1
figure 1

Structure of lead phytochemicals present in W. somnifera plant identified in the present study against targeted SARS-CoV-2 and Human host targeted proteins. (Duke 1992; Kumar and Patnaik 2016; Singh et al. 2017)

Full size image

In the present study, galactitol, chlorogenic acid, quercetin-3-rutinoside-7-glucoside, isochlorogenic acid, rutin, quercetin, and caffeoyl-quinic acid present in W. somnifera showed potential binding (− 10.1, − 9.92, − 9.66, − 8.51, − 8.38, − 8.25, and − 8.21 kcal/mole respectively) against SARS-CoV-2 Spike glycoprotein S1 domain. Standard inhibitor VE607 showed binding energy − 7.181 kcal/mole. Docking poses all lead (binding energy cutoff value ≤ standard inhibitor) and standard compound at S1 domain is shown in Fig. 2 A. Surface structure of top two lead phytochemicals and standard inhibitor at S1 domain is depicted in Fig. 2B−D. Binding energies and type of interaction of standard inhibitor and lead W. somnifera phytochemicals with Spike glycoprotein S1 domain is provided in Tables 1 and 2. VE607, galactitol, chlorogenic acid, quercetin-3-rutinoside-7-glucoside, isochlorogenic acid, rutin, quercetin, and caffeoyl-quinic acid interact with the active site of the protein using 3, 4, 3, 3, 4, 7, 2 and 4 hydrogen bonds respectively. Notably VE607 and galactitol showed hydrogen bonding with similar amino acid residues and both ligands interact with Asn422 and Leu492 via hydrophobic interaction. Rutin interacted with minimum (01 residue), while Chlorogenic acid showed maximum (11 residue) hydrophobic interaction with the S1 domain of SARS-CoV-2 spike glycoprotein.

Table 1 Binding energy of W. somnifera phytochemicals against SARS-CoV-2 and human host targeted proteins with cutoff  ≤ − 4.5 kcal/mole
Full size table
Fig. 2
figure 2

Withania somnifera lead phytochemical docking pose and interaction with SARS-CoV2 Spike glycoprotein (S-protein) protein. A Docking pose of standard inhibitor and lead phytochemicals at S1 domain of S-protein. B Surface structure of S-protein interacted with lead compound galactitol. C Surface structure of S-protein interacted with lead compound chlorogenic acid. D Surface structure of S-protein interacted with standard inhibitor VE607. Green and orange color represents the amino acid involved in hydrogen bonding and hydrophobic interaction respectively. The molecular docking was performed as per methodology discussed in material and method section

Full size image

Quercetin 3-rutinoside-7-glucoside, rutin, caffeoylquinic acid, sitoindoside X, and dulcitol present in W. somnifera plant showed potential binding (binding energy − 11.69, − 9.49, − 7.91, − 7.15, and − 6.61 kcal/mole) against human CD26 active site. Standard inhibitor sitagliptin showed binding energy of − 6.6 kcal/mole. Docking pose of the lead phytochemicals (binding energy cutoff value ≤ standard inhibitor) and the standard compound, at the human CD26 active site is depicted in Fig. 3A. Surface structure of the top two lead phytochemicals and standard inhibitor at human CD26 active site is directed in Fig. 3B−D. Binding energies and type of interaction of standard inhibitor and lead W. somnifera phytochemicals with CD26 active site is shown in Tables 1 and 2. Sitagliptin, quercetin 3-rutinoside-7-glucoside, rutin, caffeoylquinic acid, sitoindoside X, and dulcitol interact with the active site of the protein using 5, 8, 6, 8, 3 and 6 hydrogen bonds respectively. These six phytochemicals against human CD26 showed binding with mostly similar type of amino acid residues (Table 2). Furthermore, W. somnifera compounds were docked on protein-protein interaction interface of human CD26 to find the CD26-S1 domain protein-protein interaction inhibition potential of the compounds. Result showed tight binding of dulcitol, quercetin-3-rutinoside-7-glucoside, and rutin phytochemicals on CD26-Spike protein (S1 domain) protein-protein interaction interface with − 8.19, − 7.13, and − 7.08 kcal/mole binding energy. Dulcitol, quercetin-3-rutinoside-7-glucoside, and rutin formed hydrogen bond with Ser239 (two bond), Tyr241(two bond); Ala707, Asp192, Asp739, Gln123, Lys122, Lys250; and Arg253, Asp192, Glu237, Leu235, Pro249, Thr231, Thr251 amino acid residues at CD26-Spike protein-protein interaction interface.

Fig. 3
figure 3

Withania somnifera lead phytochemical docking pose and interaction with human host CD26 protein active site. A Docking pose of standard inhibitor and lead phytochemicals at CD26 active site. B Surface structure of CD26 active site with lead compound Quercetin-3-rutinoside-7-glucoside. C Surface structure of CD26 active site with lead compound rutin. D Surface structure of S-protein interacted with standard inhibitor sitagliptin. Green and orange color represents the amino acid involved in hydrogen bonding and hydrophobic interaction

Full size image

In the present study, rutin, isochlorogenic acid, caffeoyl-quinic acid, 4-2-3-didehydrosomnifericin, and chlorogenic acid present in W. somnifera plant showed potential binding (− 7.71, − 6.68, − 5.89, − 5.58, and − 5.05 kcal/mole) against human TMPRSS2 active site. Standard inhibitor camostat mesylate showed binding energy − 5.146 kcal/mole. Docking poses for all lead (binding energy cutoff value ≤ standard inhibitor) and standard compound at human TMPRSS2 active site is depicted in Fig. 4A. Surface structure of top two lead phytochemicals and standard inhibitor at human TMPRSS2 active site is depicted in Fig. 4B−D. Binding energies and type of interaction of standard inhibitor and lead W. somnifera phytochemicals with TMPRSS2 is provided in Tables 1 and 2. Camostat mesylate, rutin, isochlorogenic acid, caffeoyl-quinic acid, 4-2-3-didehydrosomnifericin, and chlorogenic acid interacted with the active site of the protein using 2, 4, 5, 3, 4 and 2 hydrogen bonds. Camostat mesylate and top five W. somnifera compounds (rutin, isochlorogenic acid, caffeoyl-quinic acid, 4-2-3-didehydrosomnifericin, and chlorogenic acid) against human TMPRSS2 showed binding with three (Ile279, Phe231, Thr324) common amino acid residues (Table 2). Moreover, isochlorogenic acid and caffeoyl-quinic acid showed maximum binding similarity with standard inhibitor against host TMPRSS2 among top five lead molecules.

Fig. 4
figure 4

Withania somnifera lead phytochemical docking pose and interaction with human host TMPRSS2 protein. A Docking pose of standard inhibitor and lead phytochemicals at TMPRSS2 protein active site. B Surface structure of TMPRSS2 protein with lead compound rutin. C Surface structure of TMPRSS2 protein with lead compound Isochlorogenic acid B. D Surface structure of TMPRSS2 interacted with standard inhibitor camostst mesylate. Green and orange color represents the amino acid involved in hydrogen bonding and hydrophobic interaction respectively

Full size image

In the present study, quercetin-3-rutinoside-7-glucoside, rutin, and chlorogenic acid present in W. somnifera plant showed potential binding (− 11.61, − 9.07, and − 8.05 kcal/mole) against SARS-CoV-2 RdRp active site. Standard inhibitor remdesivir showed binding energy − 7.56 kcal/mole. Docking poses all lead (binding energy cutoff value ≤ standard inhibitor) and standard compound at SARS-CoV-2 RdRp active site which is depicted in Fig. 5A. Surface structure of top two lead phytochemicals and standard inhibitor at SARS-CoV-2 RdRp active site is provided in Fig. 5B−D. Binding energies and type of interaction of standard inhibitor and lead W. somnifera phytochemicals with RdRp is shown in Tables 1 and 2. Remdesivir, quercetin-3-rutinoside-7-glucoside, rutin, and chlorogenic acid interacted with the active site of the protein using 5, 9, 6 and 6 hydrogen bonds. Remdesivir and quercetin-3-rutinoside-7-glucoside interacted with 6 similar amino acid residue (Arg553, Tyr619, Asp623, Asp618, Asp760 and, Lys621) of SARS-CoV-2 RdRp (Table 2).

Fig. 5
figure 5

Withania somnifera lead phytochemical docking pose and interaction with SARS-CoV-2 RdRp active site. A Docking pose of standard inhibitor and lead phytochemicals at SARS-CoV-2 RdRp active site. B Surface structure of RdRp active site with lead compound Quercetin-3-rutinoside-7-glucoside. C Surface structure of RdRp active site with lead compound rutin. D Surface structure of RdRp active site interacted with standard inhibitor remdesivir. Green and orange color represents the amino acid involved in hydrogen bonding and hydrophobic interaction respectively

Full size image
Table 2 Amino acid residues and type of interaction involved in binding of standard ligands and W. somnifera lead phytochemicals with target SARS-CoV-2 and human host proteins
Full size table
Fig. 6
figure 6

Binding energy of W. somnifera lead phytochemical groups against the targeted SARS-CoV-2 (Spike-glycoprotein and RdRp) and Human host (CD26 and TMPRSS2) test proteins. The lead (more potent binding in comparison to standard inhibitors of the targeted proteins) phytochemicals were groped as A’-Other phytochemicals; B’-Flavonoids/Polyphenoilcs; C’-Saponin; D’-Alkaloids; E’-Steroidal/Benzenoid lactones. Other phytochemicals included pentose phosphate, glucinolates, amino acids, fatty acids and cyclitol. The analysis was done on GraphPad Prism software using data provided in Table 1

Full size image
Fig. 7
figure 7

MM-GBSA analysis of the lead W. somnifera phytochemical complexed with the SARS-CoV-2 (Spike-glycoprotein and RdRp) and Human host (CD26 and TMPRSS2) test proteins. The graph represents the contribution of different types of interaction energy and the algebraic sum of the total energy of the ligand-protein complex. The total binding free energy involves different types of interactions (coulomb, covalent, lipo (lipophilic), Hbond (hydrogen bond), SolvGB (generalized born electrostatic solvation energy), vdW (Van deer Waals), and packaging interactions). A Binding free energies of CD26 and Quercetin-3-rutinoside-7-glucoside complex. B Binding free energies of TMPRSS2 and Rutin complex. C Binding free energies of RdRp and Quercetin-3-rutinoside-7-glucoside complex. D Binding free energies of Spike-protein and Galactitol complex

Full size image

On the basis of docking results, we selected proteins that are involved in the viral entry (TMPRSS2) and virulence (RdRP) to perform MD simulation study. The simulations were carried out on unbound state and lead phytochemical bound state of the respective proteins to elucidate the dynamic behavior of both proteins. Simultaneously, the protein bounded with their biologically proven inhibitors also ran to compare the results. MD simulation results for RdRP-Quercetin-3-rutinoside-7-glucoside (RQ) complex, TMPRSS2-Rutin (TR) complex and their respective unbound proteins and complex with standard inhibitors are shown in Figs. 8 and 9. Favipiravir and camostat mesylate were used as biologically approved RdRP and TMPRSS2 protein inhibitors in MD simulation study. The MD trajectory of both proteins was stable under most parts of 10 ns simulation period. Consequently, these trajectories were used for further analysis. RMSD analysis revealed that RQ and TR complexes did not shown any significant deviation from their respective unbound proteins (Figs. 8A and 9A). RMSD of both complexes were comparable to their respective biological known inhibitors complex. To understand the effect of inhibitor binding on the flexible areas of the protein structure, root-mean-square fluctuation (RMSF) was analyzed. Little fluctuation was observed in the RQ and TR complexes in comparison to unbound proteins (Figs. 8B and 9B). Binding of standard inhibitors to their respective proteins showed similar pattern of fluctuation during simulation period.

Fig. 8
figure 8

Simulation results of W. somnifera lead phytochemical targeting RdRP protein. Plot of molecular dynamic simulation trajectories of COVID-19 RdRP protein and protein-ligand complexes during 10 ns simulation. A The root mean square deviation (RMSD) of solvated SARS-CoV2 RdRP protein and RdRP -favipiravir and RdRP-Quercetin-3-rutinoside-7-glucoside complex during 10 ns molecular dynamics simulation. B The root mean square fluctuation (RMSF) values of solvated SARS-CoV2 RdRP protein and RdRP-favipiravir and RdRP-Quercetin-3-rutinoside-7-glucoside complex plotted against residue numbers. C Plot of solvent accessible surface area (SASA) during 10 ns molecular dynamics simulation of SARS-CoV2 RdRP protein and RdRP-favipiravir and RdRP-Quercetin-3-rutinoside-7-glucoside complex. D Plot of radius of gyration (Rg) during 10 ns molecular dynamics simulation of SARS-CoV2 RdRP protein and RdRP-favipiravir and RdRP-Quercetin-3-rutinoside-7-glucoside complex. E Plot of number of hydrogen bond in the SARS-CoV2 RdRP protein in unbound state and bound with favipiravir and Quercetin-3-rutinoside-7-glucoside. F Plot of number of hydrogen bonds formed between RdRP and favipiravir along with RdRP and Quercetin-3-rutinoside-7-glucoside. Unbound protein parameters are depicted in black color. Parameters for RdRP-favipirvavir complex and RdRP-Quercetin-3-rutinoside-7-glucoside complex are represented in red and green color respectively

Full size image

Solvent accessible surface area analysis revealed that there were no significant difference in the SASA value found for RQ and TR complexes, unbound proteins and their complex with respective standard inhibitors (Figs. 8C and 9C). The SASA value for RQ and RdRP- Favipiravir complexes as well as the unbound protein progressively decreased during the 10ns simulation period (Fig. 8C). Radius of gyration for RQ and TR complexes are compared with their respective unbound protein and standard inhibitor bound protein in Figs. 8D and 9D. TQ complex showed lesser Rg value during the 10 ns simulation period in comparison to unbound and camostat mesylate protein complex (Fig. 9D). Hydrogen bond analysis was performed for RQ and TQ complexes during simulation period and results are shown in Figs. 8E, F, 9E, F. Binding of Quercetin-3-rutinoside-7-glucoside to RdRP protein resulted in increase in the number of hydrogen bonds in the RdRP protein (Fig. 8E). In an event of TMPRSS2, binding of identified lead inhibitor had no overall effect on intra-molecular hydrogen bonding network of TMPRSS2 (Fig. 9E). In both cases (RQ and TQ complexes), number of hydrogen bonds between protein and respective identified lead inhibitor were significantly higher than their standard inhibitors (Fig. 8F and 9F).

Fig. 9
figure 9

Simulation results of W. somnifera lead phytochemical targeting TMPRSS2 protein. Plot of molecular dynamic simulation trajectories of Human TMPRSS2 protein and protein-ligand complexes during 10 ns simulation. A The root mean square deviation (RMSD) of solvated Human TMPRSS2 protein and TMPRSS2-Camostat mesylate and TMPRSS2-rutin complex during 10 ns molecular dynamics simulation. B The root mean square fluctuation (RMSF) values of solvated Human TMPRSS2 protein and TMPRSS2-Camostat mesylate and TMPRSS2- rutin complex plotted against residue numbers. C Plot of solvent accessible surface area (SASA) during 10 ns molecular dynamics simulation of Human TMPRSS2 protein and TMPRSS2- Camostat mesylate and TMPRSS2- rutin complex. D Plot of radius of gyration (Rg) during 10 ns molecular dynamics simulation of Human TMPRSS2 protein and TMPRSS2-Camostat mesylate and TMPRSS2- rutin complex. E Plot of number of hydrogen bond in the Human TMPRSS2 protein in unbound state and bound with Camostat mesylate and rutin. F Plot of number of hydrogen bonds formed between TMPRSS2 and Camostat mesylate along with TMPRSS2 and rutin. Unbound protein parameters are depicted in black color. Parameters for TMPRSS2- Camostat mesylate complex and TMPRSS2- rutin are represented in red and green color respectively

Full size image

We further performed secondary structure analysis for both RQ and TR complex, unbound RdRP and TMPRSS2 proteins and standard inhibitor bound complexes during the 10 ns simulation (Supplementary Fig. 1). Secondary structure analysis was conducted to assign the secondary structure content of the proteins as a function of time. The various secondary structures of protein such as, helix, sheet, bridge, and turn were divided for particular residues for every step in 10 ns simulation. The number of amino acid residues participating in formation of a particular secondary structure can be seen in secondary structure plots of both RdRP and TMPRSS2 protein simulations in unbound and lead/standard inhibitor bound complex (Supplementary Fig. 1). The average number of residues participating in the formation of secondary structure in RdRP-favipiravir complex and RdRP-Quercetin-3-rutinoside-7-glucoside complex seems to increase slightly in comparison to unbound RdRP protein. In case of TMPRSS2, there was a slight decrease in the number of residues involve in the secondary structure formation in TMPRSS2- Camostat mesylate due to the increase in the percentage of bend (Supplementary Fig. 1). However in TMPRSS2- rutin complex, there is a slight increase in the number of residues involved in the formation B-sheet along with a slight increase in the number of residues forming the structured region.

Drug-likeness and toxicity properties of all the lead phytochemicals (binding energy ≤ respective standard inhibitor) against targeted proteins were performed and results are shown in Fig. 10. Compound itoindoside IX ; caffeoyl quinic acid; dulcitol; sitoindoside X, chlorogenic acid and galactitol were predicted as not absorbed and not brain penetrant as they were located outside the egg (shown as white circle) whereas compound 2,3-didehydrosomnifericin and quercetin were predicted as well-absorbed (inside the egg) (Fig. 10). Moreover, these compounds (2, 3-didehydrosomnifericin and quercetin) showed blood brain barrier inaccessibility (out of egg yolk). It is noteworthy that compound 2, 3-didehydrosomnifericin and quercetin showed higher probability (white circle area of the egg) of passive absorption by the gastrointestinal tract. Anahygrine, cuscohygrine, (-)-anaferine molecules showed potential to cross blood brain barrier (yellow colour circle as egg yolk). Notably, the compound (-)-Anaferine and 2, 3-Didehydrosomnifericin were predicted to be actively effluxed by P-glycoprotein (blue dot) while compound Anahygrine, cuscohygrine, and quercetin were predicted as non-substrate of P-gp (red dot). Moreover, among the not absorbed compounds sitoindoside IX and sitoindoside X were predicted actively effluxed by P-gp while Caffeoyl quinic acid, dulcitol, 2, 3-didehydrosomnifericin and galactitol showed P-gp non-substrate. Bioavailability radar analysis showed that galactitol has all the drug likeness properties within optimal range (pink colour web part) (Fig. 10C). In contrast, quercetin-3-rutinoside-7-glucoside and rutin possessed slightly higher size and polarity (white portion of the egg) (Fig. 10B, D). Further, we predicted the physiochemical properties and toxicological profile of the lead phytochemicals and result are summarized in Table 3.

Table 3 Physiochemical properties and toxicological parameters of lead W. somnifera phytochemicals against SARS-CoV-2 and human host target protein
Full size table
Fig. 10
figure 10

Boiled egg diagram and bioavailability radar map of W. somnifera compounds. A Boiled egg diagram of compound M1, M2, M3, M6, M7, M8, M9, M11, M12, M13 and M14. Bioavailability radar map of B Rutin, C Galactitol and D Quercetin-3-rutinoside-7-glucoside depicting the LIPO (lipophilicity), SIZE (molecular weight), POLAR (polarity), INSOLU (insolubility) INSATU (insaturation) and FLEX (rotatable bond flexibility) parameters. M1 Anahygrine, M2 Cuscohygrine, M3 (-)-Anaferine, M6 Sitoindoside IX, M7 Caffeoyl quinic acid, M8 Dulcitol, M9 Sitoindoside X, M11 2,3-Didehydrosomnifericin, M12 Chlorogenic acid, M13 Galactitol, M14 Quercetin

Full size image

Discussion

In the present study, we found that lead W. somnifera phytochemicals interacted with two of the critical amino acids (Gln493 and Ser494) by forming hydrogen bonding (Table 2). Beside these, the phytochemicals showed hydrophobic interaction with other amino acids (Tyr495, Asn422, and Pro491) involved in ACE-2 binding (Gordon et al. 2020; Shang et al. 2020) reported that residue 455, 482–486, 493, 494, and 501 are critical SARS-CoV-2 S1 domain amino acids involved in interaction with human ACE-2 protein. Lead phytochemicals also showed hydrophobic interaction (Leu492, Pro491, Asp442, Asn422, Ser494, Tyr495, and Gln493) with SARS-CoV-2 receptor binding domain (RBD) and ACE-2 protein-protein interaction interface residues. The results indicate that W. somnifera phytochemicals have potential to disrupt spike glycoprotein-ACE-2 protein-protein interaction by binding at SARS-CoV-2 RBD. Thus, the lead molecules might inhibit the viral entry into the cell.

The results demonstrate that quercetin 3-rutinoside-7-glucoside and other lead compounds (except sitoindoside X) exhibited greater hydrogen bond (˃ 5) interaction with the CD26 protein in comparison to standard inhibitor. Although, quercetin 3-rutinoside-7-glucoside and caffeoylquinic acid formed same number of hydrogen bond with the target active site of the protein, but former showed comparatively tight binding which might be due to the involvement of more hydrophobic interaction in the binding (Table 2). Different binding efficacy of various lead phytochemicals towards CD26-SARS-CoV-2 S-protein-protein interface and CD26 active site further corroborate the previous reporting that these two binding sites on CD26 are different (Lu et al. 2013; Vankadari and Wilce 2020). Several in vitro studies reported CD26 or DPP4 activity inhibition potential in W. somnifera extracts which corroborate with our present findings that phytochemical present in the plant possess CD26 inhibitors (Singh et al. 2020a, b). CD26/DPP4 possesses dipeptidyl peptidase activity and protein interface for viral protein attachment. Inhibition of peptidase activity by small molecules is known to decrease inflammation and associated disease occurrence. Moreover, disruption of viral-host protein-protein interaction by inhibitors has suggested a treatment strategy against SARS-CoV-2 (Fig. 11). Binding of W. somnifera phytochemicals to active/protein-protein interaction interface sites propose it as potential candidate inhibiting SARS-CoV-2 entry into the cell.

Fig. 11
figure 11

Lead phytochemicals present in W. somnifera and their SARS-CoV-2 and human host protein targets. The red color arrow in front of dulcitol, rutin, and QRG indicates their potential binding with the target SARS-CoV-2 and human host proteins involve in viral entry and replication. The detailed function of each protein is discussed in introduction part. Blue color text indicates the target proteins. QRG-quercetin-3-rutinoside-7-glucoside

Full size image

Human TMPRSS2 is an important protease involved in the priming of SARS-CoV-2 S protein that makes it a probable drug target against viral virulence. Similar to standard inhibitor, isochlorogenic acid B and caffeoyl-quinic acid showed hydrophobic interaction with Arg277 residue of TMPRSS2 protein (Table 2). Tight binding of W. sominfera phytochemicals with critical amino acid residues, in comparison to standard inhibitor(s) indicate the SARS-CoV-2 virulence inhibitory potential of the plant. Present study corroborate with the finding of Hoffmann et al. (2020) which reported the TMPRSS2 inhibition mediated decreased SARAS-CoV-2 viral virulence by using clinically proven protease inhibitor. Gordon et al. (2020) reported that remdesivir, a nucleotide analogue possess RNA synthesis inhibition potential in SARS-CoV-2 by delaying the RNA chain termination mechanism. Several clinical trials were conducted to study the efficacy of RdRp inhibitors in SARS-CoV-2 patients. Binding of W. somnifera lead molecules at RdRp active site showed interaction with similar amino acids (Arg553, Asp623, Tyr619, Lys621, Asp618, and Asp760) which indicated their SARS-CoV-2 RdRp inhibition potential.

The present study showed no significant deviation in RMSD values after binding of the lead molecules to RdRP and TMPRSS2 proteins indicating the stable protein ligand complex formation. RMSF results indicate that lead inhibitor binding provide more flexibility to both proteins under the 10 ns simulation period investigation which means that the bound confirmations of both proteins were more stable during simulation. SASA analysis explores the surface of proteins that can be assessed by the solvent molecules. RQ and TR complex showed favorable SASA value in comparison to unbound and standard inhibitor bound protein indicating the stable binding of lead molecule to the respective protein. Radius of gyration provides valuable information about protein conformation and its stability during the simulation period. Relatively stable Rg value during and at the end of the simulation indicate favorable binding of ligand to respective proteins without affecting its stability. Hydrogen bonds are an important type of interaction and play significant role in providing the stability to the molecular complexes. Increased number of hydrogen bonds in the RdRP protein after lead molecule binding indicates that binding of inhibitor has made structure more stable in comparison to unbound RdRP protein and RdRP–favipiravir complex (Fig. 8E). Analysis of hydrogen bonding network between protein and inhibitor can provide clear indication of stability of binding between protein and inhibitor. Our results indicate that more number of hydrogen bonds were formed by the lead inhibitors in comparison to standard inhibitors, which might be responsible for their more negative docking scores. Over, these MD simulation analysis (RMSD, RMSF, SASA, Rg, HBond and Secondary structure) and free energy calculations showed that the binding of W. somnifera lead phytochemicals to their respective targeted proteins stabilized the protein structure.

Three phytochemicals viz. quercetin-3-rutinoside-7-glucoside, galactitol and rutin possess anti-SARS-CoV-2 potential by targeting viral and host proteins involved in disease infection and virulence (Fig. 7). Rutin, a plant flavonoid, contain flavonol (quercetin) and a disaccharide (rutinose). Recently, in silico SARS-CoV-2 main protease inhibitory potential of rutin has been reported (Das et al. 2021). Galactitol is a sugar alcohol found in plants. Recently, in silico antiviral activity of galactitol against Ebola virus has been reported (Nagrajan et al. 2019). Quercetin 3-rutinoside-7-glucoside is a flavonoid glycoside identified in W. somnifera aerial part (Mundkinajeddu et al. 2014). Antiviral potential of flavonoid glycoside has been reported against human coronaviruses (Schwarz et al. 2014). Furthermore, chlorogenic acid, isochlorogenic acid B, caffeoyl-quinic acid, sitoindosite X, and 4-2-3-didehydrosomnifericin phytochemicals showed better binding efficacy than respective standard inhibitors of the targeted proteins (Tables 1 and 2). Published literature showed potent antiviral efficacy of these phytochemicals against herpes simplex virus 1 and 2; human immunodeficiency virus, adenovirus, respiratory syncytial virus, and Ebola virus (Li et al. 2005; Naveed et al. 2018). Sitoindosite X is a glycowithanolides identified in W. somnifera known to possess anti-stress, and memory boosting properties in vivo. In pharmacological parameters, galactitol showed better potential (such as low molecular weight) in comparison to rutin and quercetin-3-rutinoside-7-glucoside. Toxicological profile prediction of lead molecules is an important step in in silico drug discovery. The result showed that all the three lead compounds are not able to pass blood brain barrier. Galactitol was predicted as substrate for P-Glycoprotein which indicates its possible efflux by the Pgp expressing cells. The low absorption through the gastrointestinal tract of the lead compounds reflects the general limitation in therapeutic natural compounds. Over all better pharmacological, ADMET and toxicological profile (Fig. 10; Table 3) of the lead molecules further substantiate their potential as anti SARAS-CoV-2 drug agents.

Conclusions

The COVID-19 pandemic has challenged the scientific community for development of effective molecules and/or herbal formulation for their use as antiviral remedy against SARS-CoV-2. Computer aided drug designing approaches such as virtual screening, molecular docking, and molecular dynamics simulation could provide information on such lead compounds in isolated form and/or present in herbal formulations. Our results uncovers that phytochemicals present in W. somnifera have potential to bind at active/protein-protein interaction interface of the SARA-CoV-2 (Spike glycoprotein, and RdRP) and human host (CD26 and TMPRSS2) target proteins involved in the viral entry and pathogenesis. The MD simulation investigation established that the protein-lead compound complex possess stable conformation and lower protein-ligand interaction energy. MM-GBSA analysis also approves the docking and simulation study of protein-lead compound complex. Over, these studies show that quercetin-3-rutinoside-7-glucoside is a potent inhibitor of human CD26 and SARS-CoV-2 RdRp while galactitol and rutin showed highest binding efficacy against viral spike glycoprotein and human TMPRSS2 protein. We conclude that W. somnifera compounds possess inhibitory potential against viral attachment, host immune hijacking, and virulence. We further recommend that the in vitro and in vivo studies should be performed to establish the anti SARS-CoV-2 potential of the W. somnifera compounds.