key: cord-0840130-qlz2q3d4
authors: Chikhale, Rupesh; Sinha, Saurabh K.; Wanjari, Manish; Gurav, Nilambari S.; Ayyanar, Muniappan; Prasad, Satyendra; Khanal, Pukar; Dey, Yadu Nandan; Patil, Rajesh B.; Gurav, Shailendra S.
title: Computational assessment of saikosaponins as adjuvant treatment for COVID-19: molecular docking, dynamics, and network pharmacology analysis
date: 2021-01-25
journal: Mol Divers
DOI: 10.1007/s11030-021-10183-w
sha: d52f8f6014882346e5c83cbbc65c707c6ebf99bc
doc_id: 840130
cord_uid: qlz2q3d4

Saikosaponins are major biologically active triterpenoids, usually as glucosides, isolated from Traditional Chinese Medicines (TCM) such as Bupleurum spp., Heteromorpha spp., and Scrophularia scorodonia with their antiviral and immunomodulatory potential. This investigation presents molecular docking, molecular dynamics simulation, and free energy calculation studies of saikosaponins as adjuvant therapy in the treatment for COVID19. Molecular docking studies for 23 saikosaponins on the crystal structures of the extracellular domains of human lnterleukin-6 receptor (IL6), human Janus Kinase-3 (JAK3), and dehydrogenase domain of Cylindrospermum stagnale NADPH–oxidase 5 (NOX5) were performed, and selected protein–ligand complexes were subjected to 100 ns molecular dynamics simulations. The molecular dynamics trajectories were subjected to free energy calculation by the MM-GBSA method. Molecular docking and molecular dynamics simulation studies revealed that IL6 in complex with Saikosaponin_U and Saikosaponin_V, JAK3 in complex with Saikosaponin_B4 and Saikosaponin_I, and NOX5 in complex with Saikosaponin_BK1 and Saikosaponin_C have good docking and molecular dynamics profiles. However, the Janus Kinase-3 is the best interacting partner for the saikosaponin compounds. The network pharmacology analysis suggests saikosaponins interact with the proteins CAT Gene CAT (Catalase) and Checkpoint kinase 1 (CHEK1); both of these enzymes play a major role in cell homeostasis and DNA damage during infection, suggesting a possible improvement in immune response toward COVID-19. [Image: see text] SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11030-021-10183-w.

Since the outbreak of COVID-19 pandemic, constant attempts are being made worldwide to develop a vaccine and utilize repurposing approaches to find out an effective cure and preventive measure [1] . Drug development using antiviral and immunomodulatory botanicals or traditional medicines is no exception to this. In these directions, two types of targets are mainly exploited; first is SARS-CoV-2 proteins to prevent the entry of virus in host cells to reduce viral load, and second is to modulate human proteins to control the replications of virus and/or to counter the immunopathogenicity [2] . A plethora of evidence suggests that viral proteins like spike proteases have been extensively studied within a short period [3] . With continuous research, a lot has been known about the COVID-19 right from SARS-CoV-2 entry in human beings till the severity of pathology it imposes and mechanism involved therein [4] .

In the human transmission of COVID-19, the angiotensin-converting enzyme 2 (ACE2) transmembrane protein provides the pathway for SARS-CoV-2 to infect host epithelial cells. ACE2 is highly expressed on lung alveolar epithelial cells; hence, the lungs are particularly vulnerable to injury following SARS-CoV-2 infection [4] . During attachment, SARS-CoV-2 binds to ACE2 via the surface glycoprotein, which is cleaved into S1 and S2. S1 subunit contains a receptor-binding domain and gives easy access to COVID-19 for direct attachment to the peptidase domain of ACE2 [2, 4] . ACE2, being a membrane receptor, needs to be glycosylated to internalize and active.

SARS-CoV-2 is a single-stranded RNA virus and its intracellular uptake is mediated by endosomes [1, 2, 4] . After incorporation into endosomes within endothelial cells, a viral protease enzyme 3CLpro of SARS-CoV is entangled in apoptosis of cell causing activation of endosomal NADPH oxidase (a pentameric complex of gp91phox, p67phox, p47phox, p22phox, and Rac) via TLR7. This catalyzes the production of superoxide and hydroxyl radicals increasing the generation of ROS which further leads to activation of NF-κB for inflammation. The viral protease enzyme also activates a NF-kB factor-dependent reporter gene which is linked with ROS generation [5] .

On the other hand, 3a protease of SARS-CoV is linked with mitochondrial cell death pathways. The mechanism of this protease depends on the activation of p38 mitogenactivated protein kinase (p38MAPK) which is a key protein kinase to regulate the transcriptional activity of NF-κB and expression of proinflammatory molecules [5] . Thus, SARS-CoV causes the phosphorylation and activation of p38MAPK which then activates NF-κB, a proinflammatory transcription factor. Together with activated p300 (a proinflammatory histone acetyltransferase), the NF-κB subsequently leads to transcription and increased production of different proinflammatory cytokines and chemokines interleukin-2 (IL-2), IL-7, IL-10, granulocyte colony-stimulating factor (G-CSF), tumor necrosis factor α (TNF-α), as well as IL-6-induced genes and TNF-α pathway-related genes. Cytokines are important mediators of the mediator of inflammation. Persistently, high levels of these cytokines, especially IL-6, as observed in severely ill patients, are called cytokine storm which is responsible for the state of hyperinflammation and ultimate lung damage. Janus Kinases are important protein kinases that are known to play an important role in IL6 and other cytokine signaling and hyperimmune activation exacerbating lung damage important features of the acute respiratory distress syndrome.

Looking into these pathways of COVID-19 pathologies involving host proteins, the drugs/molecules can be developed into prospective drug candidates for COVID-19 treatment through possible influence on mechanisms, viz. suppression of the glycosylation and membrane attachment of ACE2; inhibition of the phosphorylation of p38 to inhibit the production of different pro-inflammatory cytokines; prevention of the translocation of gp91phox to the membrane to inhibit the activation of NADPH oxidase and reduce oxidative stress and inflammation; and inhibition of JAKdependent cytokine signaling to prevent hyperinflammation ( Fig. 1) .

Currently, chemoinformatics and molecular modeling method approaches are becoming vital components in the drug discovery due to reduction in cost and time. Computational methods signify a competent approach to efficiently alter large and diverse compound libraries to potential candidates for drug development. In this sense, there are many compounds suggested by computational methods that could be evaluated quickly with in vitro techniques [6, 7] .

Saikosaponins are major biologically active triterpene glucosides, mostly occurring in Traditional Chinese Medicines (TCM) like Bupleurum spp., Heteromorpha spp., and Scrophularia scorodonia. The traditional uses of these TCM plants and saikosaponins are well known and reported in Chinese Pharmacopeia and Japanese Pharmacopeia. Saikosaponins possess antiviral, anti-hepatitis, antihepatoma, anti-inflammatory, anti-depressant, anti-neurodegenerative, antitumor, nephron-protective, and antibacterial effects, both in vitro and in vivo [8] [9] [10] . Saikosaponins (a, b, c, and d) reported to have sedative and analgesic, antiviral, anticancer, anti-inflammatory, and antibacterial properties and protect the liver and kidney [11] [12] [13] .

Saikosaponins demonstrated immunomodulatory potential with enhanced production of immunosuppressive mediators (Th1/Th2 cells, IL-4, and IL-10) and suppression of proimmune mediators (IFN-c and TNF-a in splenocytes) [12] . As compared to lamivudine, saikosaponins revealed better HBV-DNA replication and HBsAg secretion, without inhibiting the cell proliferation [14] . Moreover, Saikosaponins (A, B2, C, and D) also inhibited the early stage of viral replication including absorption and penetration of the coronaviruses [1] . In support of established in vitro anticoronaviral potential, our recent in silico studies of 23 saikosaponins on NSP15 endoribonuclease (PDB ID: 6W01) and prefusion 2019-nCoV spike glycoprotein (PDB ID: 6VSB) proteins of 2019-nCoV2 reported the plausible role of Saikosaponin_U and Saikosaponin_V in the treatment for COVID-19 [15] . Recently, in silico studies have reported repurposing of existing drugs and plausible role of natural lead in the treatment for COVID-19 [16] [17] [18] .

In view of the association of cytokine storm in COVID-19 and the major role of inflammatory responses in the pathogenesis of SARS-CoV-2, targeting the coronavirus alone with antiviral therapy might not be sufficient to reverse highly pathogenic infections. The current prophylactic measures are insufficient, and as of now, there are no vaccines or antiviral drugs declared to be officially useful against the infection.

The present investigation demonstrates molecular docking, simulation, and MM-PBSA studies of TCM as adjuvant therapy in the treatment for COVID-19. Three different proteins targets involved in inflammation and oxidative stress were selected for molecular docking and simulation studies. 

The AMBER18 software package was used for MD simulations, and ligands were parameterized with ANTECHAM-BER. The protein-ligand complexes prepared in xLeap were subjected to 100 ns MD simulations on Nvidia V100-SXM2-16 GB Graphic Processing Unit using the PMEMD. CUDA module. Further, the 100 ns trajectories were subjected to MM-GBSA analysis using Amber18 and Amber18 tools on all the 10,000 frames (Supplementary File).

Initially, the Isomeric SMILES of individual bioactives were retrieved from PubChem (https ://pubch em.ncbi.nlm. nih.gov/) or ChemSpider (www.chems pider .com) wherever applicable and their protein-based targets were identified at the minimum pharmacological activity > 0.5 using DIGEP-Pred [19] . The compounds with negative Mol Charge and molecular weight 1250 were excluded from the study as the inputs of these compounds are not recognized by the cheminformatic tool. The list of predicted targets was then queried using STRING [20] Version 11.0, and the protein-protein interaction was enriched concerning Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database.

The human IL-6R is a 468 amino acid protein which comprises a signal peptide (19 residues), an extracellular region (339 residues), a transmembrane domain (28 residues), and a short cytoplasmic domain (82 residues). The N-terminal domain D1 (residues 1-93) is distinctive of the Ig subfamily, and domains D2 (residues 94-194) and D3 (residues 195-299) are similar to cytokine binding domains [21, 22] . SARS-CoV-2 induces fatal acute respiratory disease, caused by combative inflammation initiated by viral replication [23] . Viral import and cell infection causes enormous autophagy of the epithelial and endothelial cells, vascular leakage and initiates the immune response of the host, which produces the ebullient proinflammatory cytokines and chemokines [24, 25] . SARS-CoV-2 patient's clinical results indicate that IL-6, IL-10, and tumor necrosis factor (TNF)-α increase during sickness and decrease after recovery [26] [27] [28] . Based on the docking score (Table 1 , Supplementary File), the best five saikosaponins and appropriate reference drug with respect to IL-6 and other proteins have been discussed.

The docking results of extracellular domains of human interleukin-6 receptor alpha chain revealed that Saikosaponin_V binds mainly into the domains D1 and D2 key residues. Saikosaponin_U, K, B1, and O showed a similar interaction (Fig. 2b , c, d, e) with docking scores of − 6.987, − 6.683, − 6.136, and − 6.072 kcal/mol respectively. All the saikosaponins showed better binding modes and affinity to the binding pocket than the standard 11-ketoβ-boswellic acid, which was found to interact with the binding pocket assembled by polar, Ser 109, Asn 226, Ser 152, Hie 223, Ser 156, Ser 224 charged Glu 114, Lys 154, Lys 105 hydrophobic Val 112, Phe 103 amino acids with docking score − 3.723 kcal/mol. The residues Lys 154, Lys 105 exhibited H-bonding with carboxyl group, and Glu 114 showed H-bonding with hydroxyl group of 11-keto-βboswellic acid (Fig. 2f ).

Janus Kinase 3 (JAK3) is a cytosolic tyrosine kinase associated with cytokine receptors and has a pivotal role in several chronic inflammations [29, 30] . Docking studies are found helpful to examine the plausible role of saikosaponins in inhibiting the JAK3. The JAK family comprises JAK1, JAK2, JAK3, and Tyk2 with over 1000 amino acids and between 120 and 140 kDa of molecular weight. There are four domains in each JAK protein: an N-terminal domain, an SH2 domain, a pseudokinase domain, and a tyrosine kinase domain [31] . The ATP-binding pocket in the tyrosine kinase domain from the crystal structure of JAK3 was chosen in docking studies. Saikosaponin_B4 was found to be positioned into the binding pocket assembled by polar Ser . All the saikosaponins showed lower binding affinity to the binding pocket than the cocrystallized ligand peficitinib, which was found to interact with catalytic site assembled by polar Asn 954 . 3f) . Superoxide, a reactive oxygen species, is known to cause cellular damages through oxidative stress [32] . The NADPH oxidases (NOXs) are the only enzymes that are devoted exclusively to the production of reactive oxygen species [33] . NOXs are membrane proteins and consist of six transmembrane domains (TM) and a cytosolic C-terminal dehydrogenase (DH) domain. DH includes the FAD and NADPH binding sites, while TM binds with two hemes [34] . For the development of superoxide and reactive oxygen molecules, DH transfers the electron from NADPH to oxygen bound to the moieties of the heme. Among these, Pro 694 was the most significant for function changes and reduction in NOX5 actions [35] . The docking study was performed on the crystal structure of the dehydrogenase domain of NOX5. In NOX5 binding pocket, Saikosaponin_C was found to be the best fit ligand, constructed by polar Thr 541, Hie 476, (Fig. 4a) .

The hydroxyl group of oxane ring of Saikosaponin_BK1 exhibited H-bonding with Thr 462, Hie 459, Hie 476, Arg 478, and Glu 691 with a docking score of − 7.813 kcal/mol (Fig. 4b) . The hydroxyl groups substituted on oxane rings of Saikosaponin_N exhibited H-bonding with Hie 459, Pro 460, and Thr 462; while hydroxymethyl groups on oxane ring and the aglycone showed H-bonding with Arg 478 and Lys 690 respectively with docking score of − 7.743 kcal/mol (Fig. 4c ). Saikosaponin_R and Saikosaponin_K showed a little bit low binding affinity as compared to the rest of saikosaponins with a docking score of − 7.646 and − 7.278 kcal/mol ,respectively (Fig. 4d, e) . In this study, we used 2-acetylphenothiazine as a standard that showed good binding affinity and was stabilized by π-π stacking interaction with Phe 461 and Trp 695 and exhibited docking score of − 8.338 kcal/mol (Fig. 4f) . Saikosaponin_C was the only ligand that showed better binding affinity than the standard, but other saikosaponins also exhibited comparable binding energy with the standard 2-acetylphenothiazine.

Molecular dynamics simulations (MDS) were performed to explore the dynamic nature and stability of the protein-ligand complex in water. It helps to understand the interaction between ligand and protein amino acid residues at the atomistic level. In the present work, the integrated workflow of molecular docking, molecular dynamics simulations and free energy calculations was employed to understand the properties of various small molecule inhibitors with the target human proteins [17, 36, 37] . The MDS studies on protein-ligand complexes for 100 ns allow exploring a wide space of conformational optimization and binding affinity. Herein, 100 ns MDS was performed on each of the complexes: human interleukin-6 bound to Saikosaponin_V and Saikosaponin_U; Janus Kinase 3 (JAK3) bound to Saikosaponin_B4 and Saikosaponin_I; and dehydrogenase domain of NOX5 bound to Saikosaponin_C and Saikosaponin_BK1. The selection of these compounds was based on their binding affinity to target receptors compared to the other phytochemicals screened in the present study. Preprocessing of the complex was performed in three stages with two steps of minimization of the waters and completes complex, simulated annealing, NPT, and NVT equilibration of 5 ns each followed by 100 ns of production stage [18] . Trajectory analysis of these complexes gives information on their binding modes, the formation of hydrogen bond, pi-pi interactions and van der Waals interactions. RMSD, The MDS trajectories of human interleukin-6 with Saikosaponin_U and Saikosaponin_V were analyzed, and the protein RMSD, ligand RMSD, and per amino acid (aa) fluctuations, RMSF were recorded (Fig. 5a, b,  c) . The RMSD analyses of human interleukin-6 bound to Saikosaponin_U showed an increase in RMSD from the beginning of the MDS with initial equilibration at around 4 Å for the first 25 ns after which it had gradual rise to 7.5 Å between 25 and 50 ns and later frequent changes increasing to 12.5 Å. This could be attributed to the relative motion of the loop region between aa25-aa60 during MDS; this is supported by the RMSF for the aa in this region of the receptor. The RMSF of aa50 is about 7.5 Å in the complex of Saikosaponin_U [ Fig. 5b (black) ]. These RMSD fluctuations are not observed in the second complex with Saikosaponin_V but the RMSF has similar trends [ Fig. 5b (red and black) ]. The Lig-RMSD of Saikosaponin_U remains between 5 and 10 Å for about 75 ns and then sharply rises to 60 Å [ Fig. 5c (black) ]. To investigate this sharp increase in Lig-RMSD, visual inspection of the MDS trajectory was performed (Fig. 5d , e, f, g). Saikosaponin_U is a large molecule, and its binding to the shallow cavity of the receptor allows it to move into a conformationally more stable position and hence even a small change in the ligand conformation is reflected with high RMSD in the plot (Fig. 5c, black) . The hydrogen bonds initially present between the residues Leu100 and Ser101 of human interleukin-6 with Saikosaponin_U (Saikosaponin_U-O6-Leu100-O; bond length = 2.69 Å. Saikosaponin_U-O8-Leu101-OG; bond length = 3.17 Å) breaks and new hydrogen bond with aa Glu114 is formed (Saikosaponin_U-O23-Glu114-OE1; Bond length = 2.76 Å) (Fig. 5d, e, f, g) . This making and breaking of bond happens during the last 25 ns of the MDS which is reflected in the Lig-RMSD. The ligand superposition also shows some structural conformational changes (Fig. 5h, i, j, k) . During the MDS, it forms a new hydrogen bond with the Arg5 while the Glu144 and Pro7 are conserved (Saikosaponin_V-O14-Glu10-N; bond length = 2.74 Å and Saikosaponin_V-O19-Pro7-O; bond length = 2.79 Å) (Fig. 5h, i, j, k) . The ligand RMSD is also nominal with fluctuations between 0 and 5 Å for about 90 ns and an increase to 10 Å toward the end of MDS. The ligand superposition does not show a huge change in conformations, and the surface analysis also confirms the stability of the complex. The MDS studies on the protein-ligand complex for the human interleukin-6 bound to Saikosaponin_U and Saikosaponin_V suggest that these complexes are quite stable.

The MDS trajectories of Janus Kinase 3 (JAK3) bound to Saikosaponin_B4 and Saikosaponin_I were analyzed, and the protein RMSD, ligand RMSD, and per amino acid (aa) fluctuations, RMSF were recorded (Fig. 6a, b, c) . The protein RMSD for Saikosaponin_B4 rises to about 2 Å during the first 25 ns and then equilibrated at around 2 Å till 85 ns after which it fluctuates slightly toward the end of MDS. The Saikosaponin_I bound to JAK3 protein equilibrates around 30 ns between 2.5 and 3.0 Å with a slight fluctuation between 70 to 75 ns. The long-equilibrated states in RMSD suggest that these two saikosaponins can form stable complexes with the JAK3 (Fig. 6a) . This finding also reflects in the RMSF and Ligand RMSD. The aa RMSF for both the complexes has a similar pattern with the fluctuations ranging between 0.5 and 2.5 Å; however, the Saikosaponin_B4 bound to residues aa226, aa227 has slightly higher RMSF compared to Saikosaponin_I bound to JAK3 (Fig. 6b) . The Interaction analysis of the human interleukin-6 (PDB: 1N26) bound to Saikosaponin_V during the molecular dynamics simulation; h equilibrated structure of Saikosaponin_V (green) bound to the human interleukin-6 before MDS production phase; i Saikosaponin_V (purple) bound to the human interleukin-6 post-MDS production phase; j receptor surface analysis; k ligand analysis ligand RMSD of Saikosaponin_B4 gradually increased to 5 Å till 50 ns of the MDS after which it equilibrated at 3-4 Å for the rest of the MDS. Saikosaponin_I RMSD increases gradually to 7 Å during the first 25 ns; afterward, it equilibrated in the same conformation for the rest of the MDS (Fig. 6c) . These observations support the stability of these complexes. The protein-ligand interactions for JAK3 bound to Saikosaponin_B4, and Saikosaponin_I during the MDS are shown in Fig. 6d , e, f, g.

In the initial conformation of Saikosaponin_B4 bound to JAK3 before MDS, a hydrogen bond between the 4-hydroxymethyl group of tetradecahydropicen-3-yl ring and the Pro172 of the active site was observed (Saikosapo-nin_B4-O5-Pro172-O; bond length = 3.13 Å) (Fig. 6d) . The first oxane ring of Saikosaponin_B4 forms a hydrogen bond with the Arg139 (Saikosaponin_B4-O7-Asp153-OD1; bond length = 2.63 Å) and the second oxane ring hydrogen bonds with Asp153 (Saikosaponin_B4-O7-Asp153-OD1; bond length = 2.90 Å). The post-MDS analysis shows retention of the hydrogen bond with Pro172 (Saikosaponin_ B4-O7-Asp153-OD1; bond length = 2.90 Å), but interactions between the oxane rings of Saikosaponin_B4 and the receptor are lost (Fig. 6e) . In this process, the ligand stabilizes and gains an energetically lower conformation (Fig. 6f, g) . However, the positional conformation is retained by the hydrogen bond with Pro172 and the channel shape of the active site was shown by the white surface in Fig. 5f .

Before performing the MDS ligand-receptor interaction analysis, it shows a hydrogen bond interaction between the hydroxy group of oxane ring in Saikosaponin_I with the Lys158 of JAK3 (Saikosaponin_I-O11-Lys158-NZ; bond length = 2.75 Å) (Fig. 6h) . The post-MDS analysis suggests the breaking of this bond. Figure 6c shows a sharp rise in the RMSD around 20 ns which suggests a sudden conformational change in the ligand. This could have caused the breaking of the old hydrogen bond and formation of a new hydrogen bond with Glu57 (Saikosaponin_I-O16-Glu57-OE1; bond length = 2.77 Å) (Fig. 6i) . The surface analysis of the superposed MDS trajectory structures shows the shift from the initial structure (Fig. 6j, k) . These observations suggest the stabilization of Saikosaponin_I occur during the MDS. It achieves energetically and conformationally stable binding position during the MDS and further stabilizes the protein-ligand complex which is evident from equilibrated protein and ligand RMSDs (Fig. 6a, c) .

The MDS trajectories of NADPH-Oxidase 5 (NOX5) bound to Saikosaponin_BK1 and Saikosaponin_C were analyzed, and the protein RMSD, ligand RMSD and per amino acid (aa) fluctuations, RMSF were recorded (Fig. 7a,  b, c) . The Saikosaponin_BK1 bound to NOX5 protein shows convergence after 20 ns of MDS, and the receptor protein is an oxidase domain with more than 260 amino acid residues (Fig. 7a) . Initially, the RMSD increases to around 2.25 Å and then 1.75 to 2.5 Å, which is highly acceptable for a protein of this size. It suggests that the binding of Saikosaponin_BK1 stabilizes the complex. The NOX5 bound to Saikosaponin_C shows fluctuations between 1.5 and 2.5 Å throughout the MDS. The RMSF of NOX5 bound to both the ligands was analyzed (Fig. 7b) . The RMSF for NOX5 is similar for both the complexes with a high fluctuation for residues between Gly95-Val100 of about 3.25 Å. The ligand RMSD of Saiko-saponin_BK1 and Saikosaponin_C shows interesting information (Fig. 7c) . The Saikosaponin_BK1 RMSD increases for 1 to 2 Å between 15 and 20 ns of MDS and stabilizes to 2 Å for the rest of the simulation. This suggests a stable complex and strong binding between the protein and ligand. On the other hand, RMSD of Saikosaponin_C fluctuates between 1 and 4 Å, but it has a pattern in fluctuations suggesting a flapping motion at an interval of every 20 ns. It suggests that the NOX5-Saikosaponin_C complex is stable but the ligand achieves various conformations within the active site of the protein which allows for some sort of flexibility to Saikosaponin_C. To understand the ligand interactions, visual assessment of the trajectories was performed.

The visual inspection of the NOX5-Saikosaponin_ BK1 complex before performing MDS revealed the presence of a hydrogen bond acceptor donor pair between the backbone of NOX5 and Saikosaponin_BK1 (Saikosapo-nin_BK1-O16-Arg68-N; bond length = 2.98 Å) (Fig. 7d ). This complex also formed a weaker side-chain donor pair (Saikosaponin_BK1-O15-Thr52-OG1; bond length = 3.17 Å). Post-MDS analysis shows the hydrogen bond observed before the simulation breaks and is replaced by the other one (Saikosaponin_BK1-O17-Arg68-N; bond length = 2.22 e Saikosaponin-BK1 (purple) bound to the NOX5 post-MDS production phase; f receptor surface analysis; g ligand analysis. Interaction analysis of the NADPH-Oxidase 5 (NOX5) (PDB: 5O0X) bound to Saikosaponin-C during the molecular dynamics simulation; h equilibrated structure of Saikosaponin-C (green) bound to the NOX5 before MDS production phase; i Saikosaponin-C (purple) bound to the NOX5 post-MDS production phase; j receptor surface analysis; k ligand analysis Å). The Saikosaponin_BK1 also forms other interactions; Saikosaponin_BK1-O14-Thr52-N; bond length = 3.12 Å and Saikosaponin_BK1-O17-Thr52-OG1; bond length = 1.77 Å (Fig. 7e) . The surface analysis of the complex before and after MDS suggests a good affinity of ligands to the active site (Fig. 7f, g) and hence low RMSD as discussed previously.

The visual inspection of the NOX5-Saikosaponin_C complex before performing MDS revealed the presence of a hydrogen bond acceptor donor pair between the backbone of NOX5 and Saikosaponin_BK1 (Saikosaponin_C-O7-Pro50-O; bond length = 2.76 Å) (Fig. 7h) (Fig. 7i) . The newly formed interactions aid the complex to remain stable and within the active site (Fig. 7j, k) , albeit higher fluctuation in the ligand RMSD observed during the MDS analysis.

MM-GBSA analysis was performed on all of the six protein-ligand complexes, viz. human interleukin-6 receptor, Janus Kinase-3, and NADPH-Oxidase 5 (NOX5) bound to various saikosaponins (Table 1 and Supplementary File) to evaluate the affinity of ligands to the target protein receptors. The MM-GBSA-based binding free energy (∆G bind ) calculations were performed on the 100 ns long 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 (∆G bind _H-bond), Coulomb or electrostatics interaction energy (∆G bind_Coul ), covalent interaction energy (∆G bind_Cov ), lipophilic interaction energy (∆G bind_Lipo ), electrostatic solvation free energy (∆G bind_Solv ) and van der Waals interaction energy (∆G bind_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 1 and Supplementary File. The human interleukin-6 receptor-Saikosaponin_U complex had binding energy of − 145.65 kcal/mol with the largest contribution from the gas phase interaction energy (ΔG gas ) of − 207.22 kcal/mol. The human interleukin-6 receptor-Saikosaponin_V complex had binding energy of − 155.79 kcal/mol, which is higher than the binding energy of Saikosaponin_U. The Human JAK3-Saikosaponin_B4 complex had binding energy of − 66.01 kcal/mol; however the JAK3-Saikosaponin_I complex had a large biding energy of − 300.07 kcal/mol with the largest contribution from the gas phase interaction energy (ΔG gas ) of − 375.61 kcal/mol. The NADPH--oxidase 5 (NOX5)-Saikosaponin_BK1 complex had binding energy of − 176.90 kcal/mol; however, the NADPH-oxidase 5 (NOX5)-Saikosaponin_C complex had lower binding energy of − 114.61 kcal/mol compared to Saikosaponin_ BK1. These results indicate that the JAK3 would be the best interacting partner for the saikosaponin compounds.

During network pharmacology analysis, out of 23 bioactives, Saikosaponin_U and Saikosaponin_2 could not be predicted for any targets because of molecular weight > 1250 and negative mol charge, respectively. Rest of the compounds collectively modulated six genes, i.e., CAT, CD83, CHEK1, NPPB, SMN2, and TNFRSF1A; where all the proteins except CHEK1 were upregulated. Further, KEGG analysis exhibited the regulation of amyotrophic lateral sclerosis (hsa05014) and HTLV-I infection (hsa05166) pathways via common modulation of TNFRSF1A ( Table 2 ). The interaction of each compound with its targets and respective modulated proteins is presented in Fig. 8 .

As recorded in KEGG, hsa05014 pathway is networked with 12 different pathways (https ://www.genom e.jp/dbget -bin/www_bget?hsa05 014) in which autophagy, apoptosis, PI3K signaling, oxidative phosphorylation, ubiquitin-proteasome system, cytoskeletal regulation, and RNA processing may play an important concern in the COVID-19 management. In COVID-19 as there is the invasion of the viral infection, it is important to remove the unwanted cellular components which are affected by the coronavirus which could play an important role in its management. Since corona virus is reported to cell necrosis by triggering inflammatory responses and activating caspase-8 [38] , triggering of apoptosis via hsa05014 pathway could play a tricky game as its regulation may inhibit the excessive cell necrosis. Similarly, it has been a great concern with the upregulation of immunity against COVID-19 and multiple approaches have been made to identify the lead hits against COVID-19 via this approach [39, 40] . Since PI3K signaling, oxidative phosphorylation, ubiquitin-proteasome system, cytoskeletal regulation, and RNA processing directly concerned with immune system regulation [41] [42] [43] [44] [45] are networked with hsa05014, modulation of this pathway would add a beneficial effect in COVID-19 management by upregulating the immune system. Additionally, hsa05166 is concerned with multiple pathways that are associated with the pathways that are involved in the regulation of viral infection, i.e., PI3K signaling (nt06114), TGFB signaling (nt06118), JAK-STAT signaling (nt06119), calcium signaling (nt06120), TNF signaling (nt06123), MHC presentation (nt06129), cell cycle (nt06130), cytokine-cytokine receptor interaction (nt06150) and human T-cell leukemia virus 1 (nt06160) (https ://www.genom e.jp/dbget -bin/www_ bget?hsa05 166). Since hsa05166 was affected by the saikosaponins, they may affect the multiple signaling pathways concerned with coronavirus. Additionally, PI3K signaling (nt06214), TGFB signaling (nt06218), JAK-STAT signaling (nt06219), calcium signaling (nt06220), TNF signaling (nt06223), and MHC presentation (nt06229) are associated with regulation of immune system [46] [47] [48] [49] [50] [51] and are concerned with hsa05166, regulated by saikosaponins. Hence, saikosaponins could upregulate the immune system and add a beneficial effect to the immune compromised subjects.

Gene CAT codes catalase, an enzymatic anti-oxidant biomarker [52] , that scavenges free radicals in multiple infectious and non-infectious pathogenesis which in turn minimizes the cell necrosis and inflammation. Since inflammation is one of the major causes of progressive necrosis and cell damage, activation of this enzyme could contribute to minimizing the above conditions via scavenging free radicals. CD83 plays a prime role important role in boosting the immune system through activation of lymphocytes and scavenging of exogenic components [53, 54] . This helps the immune system to fight back with infection which is beneficial to immune-compromised subjects suffering from multiple infectious and noninfectious diseases. CHEK1 is involved in the cell damage [55] and its downregulation could minimize necrosis and attenuate the cell damage. NPPB plays a key role in cardiovascular homeostasis through natriuresis, diuresis, vasorelaxation, and inhibition of renin and aldosterone secretion (https ://www.prosp ecbio .com). Activation of this protein could minimize the fluid retention, thereby facilitating O 2 and CO 2 exchange in COVID-19-infected subjects. SMN2 plays an important role in the splicing of cellular pre-mRNAs [56] . Upregulation of this protein could provide a Fig. 8 Interaction of bioactives with their targets and regulated pathways more favorable condition for the synthesis of homeostatic proteins and boosting of the immune system. TNFRSF1A contributes to inducing noncytocidal effects including antiviral state and activation of the acid sphingomyelinase (https ://www.prosp ecbio .com).

CAT (Identifier: ENSP00000241052, CAT; organism: Homo sapiens) codes the protein catalase and STRING database links it with ten different genes, i.e., HAO1, SCP2, DAO, ACOX1, HSD17B4, SOD2, AKT1, SOD1, SOD3 and GSR. The catalase promotes growth of cells including T-cells, B-cells, myeloid leukemia cells, melanoma cells, mastocytoma cells, and normal and transformed fibroblast cells which are concerned with the immune system regulation [57] . Additionally, catalase also acts as an antioxidant [58] which would play an important role in scavenging hydroxyl free radicals produced during cell necrosis due to corona virus infection. Further, STRING networks CD83 (Identifier: ENSP00000368450, CD83; Organism: Homo sapiens) with CD86, TNF, CD40, CD80, IL10, CD1A, CD1B, CD1C, CD1D, and CD1E and encodes CD83 antigen which could play a significant role in antigen presentation or the cellular interactions followed by the activation of lymphocytes and T cell response which are associated to the symptoms of COVID-19 [59] . Likewise, CHEK1 and encodes serine/threonine-protein kinase Chk1; links with ATM, EXO1, CDC45, CDC6, CDC25A, CDC25C, RAD51, TOPBP1, ATR, and CLSPN and is required for checkpoint-mediated cell cycle arrest and is also involved in immune regulation [60] which could add a beneficial effect to the immune compromised subjects.

In the present work, the role of saikosaponins as adjuvants in the treatment for COVID-19 is analyzed through in silico methods. The molecular docking studies and the dynamic simulations projected the plausible potential of Saikosaponins_U and Saikosaponin_V as promising adjuvants in the COVID-19 treatment. The present study could be the starting point for the future ligands from natural sources in 2019-nCoV as adjuvant treatment.

Authors' contributions SSG and MMW contributed to conceptualization, supervision, investigation. RVC, SKS, and RBP conducted molecular docking and dynamics studies. PK, NSG, and YDD performed network pharmacology and analysis. SKS, SKP, and MA performed methodology, software data analysis. SSG and NSG contributed to writing-original draft preparation. SKP and MA were involved in formal analysis, visualization, reviewing and editing.

Spread of SARS-CoV-2

Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune crossreactivity with SARS-CoV

Consensus virtual screening of dark chemical matter and food chemicals uncover potential inhibitors of SARS-CoV-2 main protease

Organ-protective effect of angiotensin-converting enzyme 2 and its effect on the prognosis of COVID-19

Synergistic activation of NF-B by nontypeable Haemophilus influenzae and tumor necrosis factor

In silico tools to study molecular targets of neglected diseases: inhibition of TcSir2rp3, an epigenetic enzyme of Trypanosoma cruzi

Computational drug design methods-current and future perspectives

A comprehensive review and perspectives on pharmacology and toxicology of saikosaponins

Saikosaponins: a review of pharmacological effects

Antiviral effects of saikosaponins on human coronavirus 229E in vitro

Genus Bupleurum: a review of its phytochemistry, pharmacology and modes of action

Triterpenoid saponins from Bupleurum fruticosum

Saikosaponin A inhibits influenza A virus replication and lung immunopathology

Anti-SARS coronavirus 3C-like protease effects of Isatis indigotica root and plantderived phenolic compounds

Identification of bioactive compounds from Glycyrrhiza glabra as possible inhibitor of SARS-CoV-2 spike glycoprotein and non-structural protein-15: a pharmacoinformatics study

An in silico evaluation of different Saikosaponins for their potency against SARS-CoV-2 using NSP15 and fusion spike glycoprotein as targets

In silico investigation of phytochemicals from Asparagus racemosus as plausible antiviral agent in COVID-19

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

DIGEP-Pred: web service for in silico prediction of drug-induced gene expression profiles based on structural formula

STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets

Structural design and molecular evolution of a cytokine receptor superfamily

Structure of the extracellular domains of the human interleukin-6 receptor α-chain

Advancements in the battle against severe acute respiratory syndrome

Understanding SARS-CoV-2-mediated inflammatory responses: from mechanisms to potential therapeutic tools

Virology, epidemiology, pathogenesis, and control of COVID-19

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Induction of pro-inflammatory cytokines (IL-1 and IL-6) and lung inflammation by Coronavirus-19 (COVI-19 or SARS-CoV-2): anti-inflammatory strategies

COVID-19 illness in native and immunosuppressed states: a clinical-therapeutic staging proposal

The JAK-STAT pathway: impact on human disease and therapeutic intervention

Involvement of the janus kinase/signal transducer and activator of transcription signaling pathway in multiple sclerosis and the animal model of experimental autoimmune encephalomyelitis

The role of JAK/STAT signaling pathway and its inhibitors in diseases

Combating oxidative stress in vascular disease: NADPH oxidases as therapeutic targets

Intramembrane bis-heme motif for transmembrane electron transport conserved in a yeast iron reductase and the human NADPH oxidase

Crystal structures and atomic model of NADPH oxidase

Design, synthesis, pharmacological evaluation and molecular docking studies of substituted oxadiazolyl-2-oxoindolinylidene propane hydrazide derivatives

LEDGF/p75 IN interaction inhibitors: in silico studies of an old target with novel approach

SARS-CoV-2 triggers inflammatory responses and cell death through caspase-8 activation

Network pharmacology of AYUSH recommended immune-boosting medicinal plants against COVID-19

Anthraquinone derivatives as an immune booster and their therapeutic option against COVID-19

Mitochondria and microbiota dysfunction in COVID-19 pathogenesis

Immune pathology associated with altered actin cytoskeleton regulation

RNA regulation of the immune system

The PI3K/Akt/mTOR pathway in innate immune cells: emerging therapeutic applications

Immunologic aspects of protein degradation by the ubiquitin-proteasome system

The role of the PI3K signaling pathway in CD4+ T cell differentiation and function

Transforming growth factor-β signaling in immunity and cancer

Regulation of JAK-STAT signalling in the immune system

Calcium signaling in immune cells

Tumor necrosis factor superfamily in innate immunity and inflammation

Major histocompatibility complex (MHC) class I and MHC class II proteins: conformational plasticity in antigen presentation

Antioxidant activity of SOD and catalase conjugated with nanocrystalline ceria

CD83: Activation marker for antigen presenting cells and its therapeutic potential

The CD83 molecule-an important immune checkpoint. Front Immunol

Checkpoint kinase 1 in DNA damage response and cell cycle regulation

Targeting RNA-splicing for SMA treatment

Introduction to T and B lymphocytes

Role of catalase in oxidative stress-and age-associated degenerative diseases

T cell responses in patients with COVID-19

Modulation of human checkpoint kinase Chk1 by the regulatory β-subunit of protein kinase CK2

Conflict of interest The authors declare no conflict of interests.