Nature's therapy for COVID-19: Targeting the vital non-structural proteins (NSP) from SARS-CoV-2 with phytochemicals from Indian medicinal plants.

Background: Containing COVID-19 is still a global challenge. It has affected the "normal" world by targeting its economy and health sector. The effect is shifting of focus of research from life threatening diseases like cancer. Thus, we need to develop a medical solution at the earliest. The purpose of this present work was to understand the efficacy of 22 rationally screened phytochemicals from Indian medicinal plants obtained from our previous work, following drug-likeness properties, against 6 non-structural-proteins (NSP) from SARS-CoV-2.
Methods: 100 ns molecular dynamics simulations were performed, and relative binding free energies were computed by MM/PBSA. Further, principal component analysis, dynamic cross correlation and hydrogen bond occupancy were analyzed to characterize protein-ligand interactions. Biological pathway enrichment analysis was also carried out to elucidate the therapeutic targets of the phytochemicals in comparison to SARS-CoV-2.
Results: The potential binding modes and favourable molecular interaction profile of 9 phytochemicals, majorly from Withania somnifera with lowest free binding energies, against the SARS-CoV-2 NSP targets were identified. It was understood that phytochemicals and 2 repurposed drugs with steroidal moieties in their chemical structures formed stable interactions with the NSPs. Additionally, human target pathway analysis for SARS-CoV-2 and phytochemicals showed that cytokine mediated pathway and phosphorylation pathways were with the most significant p-value. Conclusions: To summarize this work, we suggest a global approach of targeting multiple proteins of SARS-CoV-2 with phytochemicals as a natural alternative therapy for COVID-19. We also suggest that these phytochemicals need to be tested experimentally to confirm their efficacy.

Introduction

The pandemic COVID-19 (coronavirus disease 2019) has spread to more than 203 countries and territories. Initial reports suggested that 20% cases are severe but its fast transmission rates has severely impacted world economy and the health care sector causing global panic (Guan et al., 2020; Huang et al., 2020; Wang et al., 2020). SARS-CoV-2 has been related to gastroenteritis, heart failure and respiratory complications. It is also severe in immune compromised patients, in men- independent of age and elderly patients with co-morbidities (Deng and Peng, 2020; Guan et al., 2020; Huang et al., 2020; Jin et al., 2020a; Islam et al., 2020). Given the severity of the infection, the proposed therapies include targeting viral replication and blocking the virus attachment to human cell receptors by its structural proteins (Wu et al., 2020). Unfortunately, coronavirus vaccines are not yet on the market and large-scale manufacturing capacity for these vaccines does not exist yet (Amanat and Krammer, 2020). Only two candidate vaccines have reached stage 2/3 clinical trials and one in phase 2 (COVID-19 vaccine tracker, 11 June 2020). Considering economic impact on stock markets in recent times, funding for vaccine production seems to be great investment soon (Amanat and Krammer, 2020). In this context development of natural drug therapy can help to contain the spread of the virus.

Nature provides a vast library of chemicals yet to be explored as drugs for treatment of various viral ailments (Denaro et al., 2020). Research have been conducted in screening natural drugs against SARS-CoV-2 targets (Jin et al., 2020b; Kumar et al., 2020; Luo et al., 2020; Pang et al., 2020; Xu et al., 2020; Yang et al., 2020; Zhang et al., 2020; Islam et al., 2020). In this context, we have recently presented 66 virtually screened phytochemicals from 55 Indian medicinal plants against 8 SARS-CoV-2 targets (Parida et al., 2020a). We also presented work related to molecular dynamics simulation of phytochemicalsfor the spike glycoprotein and the main protease of SARS-CoV-2 (Parida et al., 2020b). Interestingly from our previous work and the vast literature that exist for drug development against SARS-CoV-2, it was realized other proteins of SARS-CoV-2 must be targeted to avoid multi drug resistance strains of the SARS-CoV-2 virus as it is prone to high rate of mutations (Pachetti et al., 2020). In this context 6 non-structural proteins (NSPs) as targets for drug development is being suggested. These NSPs may have conserved roles within the viral lifecycle of SARS-CoV-2 (Littler et al., 2020). The NSPs selected in this study were due to their indispensable nature in the viral life cycle. Other than the main protease NSP5 (3CLpro) vital targets of SARS-CoV-2 considered in this work were NSP3, NSP9, NSP10, NSP12, NSP15 and NSP16.

NSP3 is the papain-like proteinase- reported to induce protective immunity by blocking host immune response by inhibiting the expression of innate immunity genes (Fehr et al., 2016). Its conserved macrodomain is involved in binding of ADP-ribose in most of the corona viruses (Han et al., 2011; Cho et al., 2016; Lei et al., 2018) and bears 91.8% sequence similarity with SARS-CoV NSP3 (Yoshimoto, 2020). Another promising target used in this study was NSP9 (RNA-binding protein) as it is thought to mediate viral replication. NSP9 from SARS-CoV-2 bears 99.1% sequence similarity with that from SARS-CoV (Yoshimoto, 2020; Sutton et al., 2004; Frieman et al., 2012).NSP10 was chosen as another important drug target as it acts as a cofactor for the activation of the replicative enzyme in SARS-CoV (Joseph et al., 2006; Matthes et al., 2006; Bouvet et al., 2014).Its sequence similarity (99.3%) with SARS-CoV makes it a promising drug target (Yoshimoto, 2020). NSP12 or the RNA dependent RNA polymerase has been chosen as another important drug development candidate as it is involved in viral replication (Neogi et al., 2020), however, but no specific inhibitors have been validated until now. Therefore, in this work other than the main catalytic pocket of NSP12, NSP12-NSP8 and NSP12-NSP7 interfaces were also chosen as potential targets. NSP8 synthesizes primer- up to six nucleotides in length, for NSP12-RdRp RNA synthesis. Further, the NSP7 acts as a cofactor to enhance the RdRps enzyme activity of NSP12 (Wu et al., 2020). Another interesting SARS-CoV-2 target used in this work was NSP15 with 97.7% sequence similarity with SARS-CoV (Yoshimoto, 2020). This protein was reported to degrade polyuridine sequences and viral dsRNA to prevent host recognition of the virus to prevent the host immune sensing system from detecting the virus (Hackbart et al., 2020; Decroly et al., 2012). It bears 99% sequence similarity with SARS-CoV NSP16.

Briefly it can be said that the aim of the present work was to explore various NSPs from SARS-CoV-2 as promising targets for phytochemicals as drug candidates. In this study, we present the analysis of 100 ns molecular dynamics trajectories for 6 NSPs from SARS-CoV-2 and their phytochemical complexes. The binding free energy and the interaction dynamics of the phytochemicals were analysed by various computational approaches. Further, the therapeutic efficacy of these phytochemicals was also experimented by biological pathway enrichment analysis. In summary it can be said that targeting these NSPs with phytochemicals can serve to overcome the predicted multi-drug resistance as SARS-CoV-2 which is fast evolving due to its propensity to high mutations.

Methods

Molecular dynamics simulations and trajectory analyses of the phytochemicals against the vital NSPs of SARS-CoV-2

Phytochemicals were successfully screened from Indian medicinal plants by molecular docking against 10 SARS-CoV-2 protein target binding sites in our previous work (Parida et al., 2020a). 22 phytochemicals with good docking scores that followed drug-like properties against SARS-CoV-2 NSP3, NSP9, NSP10, NSP12, NSP16 and NSP15 were analysed by 100 ns molecular dynamics (MD) simulations (Table 1 ). The repurposed drugs that had the lowest binding score were also subjected to MD simulations as a control group. All the simulations were carried out with Gromacs 2020.2 software package (Lindahl et al., 2020) with AMBER99SB-ILDN force field and TIP3P water molecules. The force field parameters for the phytochemicals and the repurposed drugs were generated by ACPYPE (AnteChamberPYthon Parser interface) (da Silva andVranken, 2012). Charges were neutralized and system was minimized with 5000 steps by steepest descent algorithm. The systems were equilibrated by 1 ns position restraint simulations of 1000 kJ mol−1 nm−2 in the NVT and NPT ensembles. Equilibrated systems were used to simulate a 100 ns with no restraint production run. Electrostatic interactions were calculated with Particle-Mesh-Ewald summation (PME) (Darden et al., 1993). Post-MD analyses were performed, this includes root mean square deviation (RMSD), root mean square fluctuations (RMSF), the radius of gyration (Rg) and hydrogen bond occupancy. Molecular Mechanics - Poisson Boltzmann Surface Area (MM-PBSA) was applied on snapshots obtained from MD trajectory to estimate the binding free energy using the GROMACS tool g_mmpbsa (Baker et al., 2001; Kumari et al., 2014). The Principal Component Analysis (PCA) was performed using GROMACS v2020.2. Further, Dynamic Cross Correlation Matrix (DCCM) analysis performed using the Bio3D package (Grant et al., 2006).

Table 1

The phytochemicals and their tartgets from SARS-CoV-2.

*Phytochemicals/ Repurposed drugs	Complexes	Docking Score (Kcal/mol)	#2D structure	*Plant Source
NSP 10 (PDB 6W4H)
27-Hydroxywithanolide B	N10P1	−8.32	Image, table 1	Withania somnifera
Anaferine	N10P2	−6.22	Image, table 1	Withania somnifera
Baricitinib	N10-R3	−9.22	Image, table 1	Repurposed drug
NSP12 D1 (RNA Binding Site) (PDB 6M71)
12-Deoxywithastramonolide	N12d1P1	−7.58	Image, table 1	Withania somnifera
Methylprednisolone	N12d1R2	−6.54	Image, table 1	Repurposed drug
NSP12 D2 (NSP12-NSP7 Interface) (PDB 6M71)
Withastramonolide	N12d2P1	−9.72	Image, table 1	Withania somnifera
Withanolide B	N12d2P2	−9.2	Image, table 1	Withania somnifera
12-Deoxywithastramonolide	N12d2P3	−9.05	Image, table 1	Withania somnifera
Withanolide R	N12d2P4	−8.94	Image, table 1	Withania somnifera
Withaferin A	N12d2P5	−7.9	Image, table 1	Withania somnifera
Hydroxychloroquine	N12d2R6	−7.73	Image, table 1	Repurposed drug
NSP12 D3 (NSP-NSP8 Interface) (PDB 6M71)
Withaferin A	N12d3P1	−9.94	Image, table 1	Withania somnifera
Lopinavir	N12d3R2	−8.33	Image, table 1	Repurposed Drug
NSP16 (PDB 6W4H)
Solvanol	N16P1	−10.98	Image, table 1	Solanum nigrum
(-)-Anaferine	N16P2	−10.91	Image, table 1	Withania somnifera
Limonin	N16P3	−10.73	Image, table 1	Nigella sativa
Chloroquine	N16R4	−10.11	Image, table 1	Repurposed Drug
NSP9 (PDB 6W4B)
27-Hydroxywithanolide B	N9P1	−8.28	Image, table 1	Withania somnifera
12-Deoxywithastramonolide	N9P2	−8.19	Image, table 1	Withania somnifera
Azadiradionolide	NPP3	−7.3	Image, table 1	Azadirachta indica
Baricitinib	N9R4	−6	Image, table 1	Repurposed drug
NSP15 (PDB 6W01)
Somniferine	N15P1	−7.21	Image, table 1	Withania somnifera
Vindolinine	N15P2	−6.87	Image, table 1	Withania somnifera
Dexamethasone	N15R3	−7.08	Image, table 1	Repurposed drug
NSP3 (PDB 6W02)
2,3-Dehydrosomnifericin	N3P1	−12.3	Image, table 1	Withania somnifera
Withanolide B	N3P2	−11.44	Image, table 1	Withania somnifera
24,25-dihydrowithanolide D	N3P3	−10.24	Image, table 1	Withania somnifera
27-Deoxy-14-hydroxywithaferin A	N3P4	−9.49	Image, table 1	Withania somnifera
Baricitinib	N3R5	−10.38	Image, table 1	Repurposed drug

The phytochemicals (docking scores) and interacting residues within the binding pocket were identified from a previous work to possess drug likeness properties (Parida et al., 2020a).

The 2D images of the phytochemicals were obtained from PubChem Database (https://pubchem.ncbi.nlm.nih.gov/).

High Performance Computing resources of Amazon EC2 G4 instances were employed for carrying out computational work on Amazon Web Services cloud computing platform (https://aws.amazon.com) with the help of RONIN interface (https://ronin.cloud/) as an award granted by the COVID-19 High Performance Computing (HPC) Consortium for carrying out this work (https://covid19-hpc-consortium.org/projects).

Biological pathway enrichment analysis of human protein targets identified for the phytochemicals

SwissTargetPrediction server was employed to identify top 15 human protein targets of the phytochemicals (Gfeller et al., 2014). Redundancy removal resulted into 93 target proteins for the phytochemicals. Further, 172 protein targets for human coronavirus were obtained by literature survey (Chakrabarty et al., 2020; Zhou et al., 2020). Enrichr and REVIGO webserver were used for gene ontology (GO) enrichment analysis (Supek et al., 2011; Kuleshov et al., 2016). The generated networks were visualized in Cytoscape v3.8 (Shannon et al., 2003). Enrichr was also employed to calculate the p- values using the Fisher exact test for pathway and database enrichment for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and VirusMint databaseichr. VirusMint is database of interactions between viral and human proteins (Chatr-aryamontri et al., 2009).

Results

Molecular dynamics simulation trajectory analysis of the NSP complexes

In this work, 6 NSP and 8 potential drug target sites docked with 22 phytochemicals and 10 repurposed drugs obtained from our previous work, were explored by 100 ns molecular dynamics simulations (Parida et al., 2020a). Throughout this paper the NSP complexes with the phytochemicals and the repurposed drugs have been abbreviated as presented in Table 2 . The binding sites have been illustrated in Fig. 1 . The binding site of NSP10 and NSP16 were the NSP10-NSP16 binding interface (PDB ID: 6w4h). The binding site in NSP10 (known to stimulate methyl transferase activity of NSP 16) comprised of residues Asn4293-Tyr4349 between helix 1 and helix 3. In NSP 16 it involved residues Lys 6836-Met6839 and Val6842- Asp6904 which lies in between helix 3 and β-sheet A edge β-strand 4. The inhibitor binding site in NSP9 (involved in viral replication, PDB ID: 6w4b) was the dimerization interface consisting of residue Arg100- Thr110 (helix 1). In NSP 12 (PDB ID: 6m71) three sites were chosen: D1 which is the catalytic pocket; D2: NSP12-NSP7 binding interface and D3: NSP12-NSP8 binding interface. In case of NSP15 (PDB ID: 6w01), an uridylate specific endoribonuclease, active site pocket with conserved residues His235, His250, Lys290, Thr341, Tyr343, and Ser294 was the site of inhibitor docking. In NSP3 (PDB ID: 6w02) which is the ATP ribose phosphatase of SARS-CoV-2, the adenosine-5-diphohoribose binding site was the inhibitor docked site. The site was in between β-strand 3 of β-sheet B- helix 2 and strand 3 of β-sheet A and helix 7. The details of the virtual screening of these 8 target sites can be found in our previous work (Parida et al., 2020a).

Table 2

Relative lowest free binding energy calculated by MM/PBSA for each NSP.

Inhibitor	van der Waal energy (kJ/Mol)	Electrostatic energy (kJ/Mol)	Polar solvation energy (kJ/Mol)	SASA energy (kJ/Mol)	Relative Binding Free energy (kJ/Mol)	#Residue contribution Binding Free energy (kJ/Mol)
NSP 10
N10P1	−184.657+/- 20.543	−5.743 +/- 4.589	97.672 +/−19.470	−20.243 +/- 1.783	−112.971 +/−16.527	Leu4345, Ile4308, Val4369
N10R3	−177.662 +/−18.014	−4.835 +/- 5.206	75.330 +/−17.752	−17.317 +/−1.421	−124.484 +/−25.220	Tyr4329, Leu4345, His4336
NSP12D1
N12d1P1	−124.475 +/- 21.162	−5.636 +/- 4.267	51.656 +/−16.145	−13.259 +/−1.881	−91.714 +/−16.660	Phe440,Phe415,Leu437
N12d1R2	−96.233 +/- 13.862	−3.993 +/−2.454	30.061 +/- 8.921	−10.827 +/−1.178	−80.992 +/ 13.123	Phe415,Phe440,Leu437
NSP12D2
N12d2P4	−142.880 +/−23.178	−6.687 +/−10.524	72.704 +/−36.840	−16.064 +/- 2.008	−92.926 +/−31.586	Lys798,Arg553,Tyr455
N12d2R6	−167.847 +/−10.302	−11.017 +/−2.508	108.884 +/−10.887	−15.784 +/- 0.797	−85.765 +/- 13.760	Trp800,Glu811,Trp617
N12D3
N12d3P1	−169.341 +/−18.584	−2.053 +/−3.186	65.330 +/−10.075	−17.814 +/- 1.169	−123.877 +/−15.553	Val398,Val330,Tyr273
N12d3R2	−325.095 +/- 24.003	−5.013 +/−2.809	138.228 +/- 11.958	−30.849 +/−1.728	−222.729 +/−21.472	Phe396,Tyr273,Val675
NSP16
N16P1	−152.891 +/−23.291	−1.058 +/−4.995	48.993 +/−10.961	−15.571 +/−1.899	−120.528 +/−23.072	Tyr6950, Leu6898,Phe6947
NSP9
N9P2	−125.731 +/−15.961	−2.591 +/- 2.944	37.451 +/−11.636	−13.235 +/- 1.491	−104.106 +/−12.125	Leu113,Phe76,Leu104
N9R4	−119.176 +/−13.521	−1.765 +/−2.925	30.670 +/- 8.861	−12.280 +/- 1.230	−102.552 +/−15.112	Leu104,Leu113,Leu107
NSP15
N15P1	−128.354 +/−16.899	−6.799 +/- 8.384	69.397 +/−21.569	−13.507 +/- 1.749	−79.263 +/−15.485	Glu340,Glu234,His235
NSP3
N3P1	−204.165 +/−16.108	−7.942 +/−3.650	107.923 +/−13.714	−21.812 +/- 0.982	−125.995 +/−13.066	Phe132,Ile131,Val49
N3P4	−183.720 +/- 23.432	−12.307 +/−4.755	89.847 +/−18.679	−18.589 +/−1.892	−124.769 +/−16.441	Ile131,Phe132,Ala50

For abbreviations please refer Table 1. The 500-step bootstrap +/- standard error have been presented.

The top 3 contributing residues are arranged in descending order of the contribution to the binding energy.

Fig. 1

The cartoon representation of the NSP from SARS-CoV-2 showing the binding sites of the inhibitors. A) NSP10, B) NSP16, C) NSP9, D) NSP12 D1, E) NSP12 D2, F) NSP12 D3, G) NSP15, H) NSP3. The amino acid residues are represented by one letter code.

Results of RMSD, RMSF and Rg were plotted and have been illustrated in Fig. 2 . In the case of NSP10 (Fig. 2A), slightly lower average RMSD value was obtained for N10P1 (0.49) and N10R3 (0.48) in comparison with NSP10 apo-protein (0.51 nm). The average RMSF fluctuations of residues were slightly lower for N10P1 and N10R3 (~0.2 nm) in comparison to the apo-protein (0.25 nm). It can also be observed from Fig. 2A that N10P1 showed comparatively decreased fluctuation (∆RMSF) in the residues interacting with the phytochemicals w.r.t. the apo-protein.

Fig. 2

Comparative trajectory analysis of the 8 NSP binding sites with the phytochemicals and the repurposed drugs along with their apo-proteins. For abbreviations please refer main text and Table 1. A) NSP10, B) NSP12d1, C) NSP12d2, D NSP12d3, D) NSP16, E) NSP15, F) NSP9, G) NSP3. Purple brackets in ∆RMSF plots represent residue atom numbers interacting with the phytochemicals and the repurposed drugs (refer Fig. 1). The atom numbering follows respective PDB structures.

For NSP12d1 (Fig. 2B) the average RMSD (0.3 nm), RMSF (0.16 nm) and Rg (3 nm) were similar to that of the apo-protein. The ∆RMSF was lower for the phytochemical and repurposed drug interacting residues. For NSP12d2 (NSP12-NSP7 interface) the average RMSD (0.3 nm) were similar with the apo-protein for N12d1P2, N12d1P3, N12d1P4 and N12d1R6 compared to the slightly higher value obtained for N12d1P1 (0.36 nm). Similar trend in average RMSF and Rg were obtained N12d1 as evident from the ∆RMSF and Rg plots illustrated in Fig. 2C. In the case of N12d3 binding (Fig. 2D), the average RMSD (0.27 nm), RMSF (0.14 nm) of N12d3P1 was slightly lower than the apo-protein, the other phytochemicals and the repurposed drug. The compactness of all the complexes were similar with average Rg (0.3 nm). The compactness of all the simulated structures were similar with average Rg of ~1.44 nm.

In the case of NSP16 (Fig. 2B) the lowest average RMSD was obtained for N16P2 (0.28 nm) compared to the apo-protein (0.33 nm). The average RMSF fluctuations of the phytochemicals were like that of the apo-protein (~0.15 nm). It is to be noted here that the repurposed drug (N16R4) did not form stable complex with NSP16 after 20 ns of the MD simulation and thus was not analysed further. Lower fluctuations at the inhibitor interacting residues were obtained for the phytochemicals w.r.t. the apo-protein as shown by the ∆RMSF plot.

In the case of NSP15 (Fig. 2E), N15P1 showed comparable average RMSD, RMSF and Rg with the apo-protein. N15P2 and N15R2 did not form stable complex with NSP15 with remarkably high average RMSD (>0.7 nm) and thus were not further analysed. N15P1 also showed decrease in fluctuations at the inhibitor binding site as evident from the ∆RMSF plot.

In the case of NSP9 (Fig. 2F), the phytochemicals showed comparable RMSD with an average of ~0.3 nm, average RMSF (0.17 nm) and average Rg (1.4 nm). The ∆RMSF plot showed lower fluctuations at the inhibitor interacting residue. Interestingly a lower dip in ∆RMSF was obtained for Phe76 for N9P2 and N9R4.

In the case of NSP3 (Fig. 2G) the average RMSD of N3P1 was comparable to the apo-protein. The other phytochemicals and the repurposed drugs showed slightly higher RMSD. The ∆RMSF plot revealed that lowest fluctuations were obtained for residues 126–132. N3P4 and N3R5 showed higher fluctuations at residues 38–50. The average Rg was 2.4 nm. It can be noted that the compactness of N3R4 started increasing after 60 ns.

Trajectory analyses for relative free energy of binding and hydrogen bond occupancy of the NSP complexes

The relative binding free energy were calculated for each complex of NSP using MM/PBSA. The standard errors were calculated by 500 steps of bootstrap analysis. It was observed that the relative binding free energies obtained for the complexes were in agreement with the RMSD, RMSF and Rg of their respective trajectories (Fig. 2A). The computed lowest relative free energy of binding among all the complexes with respect to each NSP has been presented in Table 2. Interestingly repurposed drugs were computed to possess high relative free energy of binding with NSP3, NSP15 and NSP16. These complexes were thus not analysed further. Among all the complexes the phytochemical with the lowest relative binding energy (< ~ −120 kJ/Mol) was obtained for N12d3P1, N3P1, N3P4 and N16P1. For the repurposed drug lowest relative binding free energy was obtained for N10R3 (−120 kJ/Mol) and N12d3R2 (~ −223 kJ/Mol). Interestingly among the NSP12 complexes the D3 site (NSP12-NSP8 interface) showed better binding with phytochemicals and the repurposed drugs than the D1 (active site) and D2 (NSP12-NSP7 interface) sites. This reveals that NSP12-NSP8 interface is a good target site for drug targeting. NSP10-NSP16 interface and NSP3 also showed good binding energies with the phytochemicals. Whereas NSP15 showed the highest relative binding energies (worst performer in binding to phytochemicals and repurposed drugs) compared to the rest of the NSP analysed through this work.

Additionally, in order to understand the contribution of each residue of the NSPs towards the total relative binding free energy, decomposition of the energy for each residue was computed (Supplementary Table S1). The top 3 residues with the lowest relative free energy of binding have been presented in Table 2. Most of the residues in all the NSP complexes with lowest binding free energy were observed to be non-polar and aromatic.

Further, hydrogen bond (H-bond) occupancy was calculated over 100 ns trajectory for the complexes (Supplementary Table S2). Interestingly the residues with low relative free energy of binding were also with highest H-bond occupancy. In the case of NSP10 highest H-bond occupancy for N10P1 was obtained for Ile4308-CD atom (31.43%). In N10R3, the atom-His4336-CA was with 24% occupancy. In the case of NSP12d1, occupancy of 21.4% was obtained for Leu437-N in N12d1P1 and Phe249 in N12d1R2 (23%). In the case of N12d2, Ala550-N and Trp800-CD1 were with 30% and 53% occupancies for N12d2P4 and N12d2R6 respectively. In the case of N12d3P1 and N12d3R2, Tyr273-CB and Phe326-CB were with 31% and 57% occupancies, respectively. From relative free energy and H-bond analysis of NSP12, it was realized that phytochemicals and repurposed drugs do not bind with the catalytic site residues S759-D761 for N12d1 and the NSP7 interacting interface (NSP12d2). However, it was realized that NSP12-NSP8 interface is a better target site for these phytochemicals and repurposed drugs as these complexes had comparatively the lowest relative free energy of binding. The residues involved in binding to NSP8 were also observed to interact with the phytochemicals and the repurposed drugs. For NSP16P1, Tyr6950-CZ was with 33.2% occupancy. The residues which are involved in high H-bond occupancy (Supplementary Table S2) are near the NSP16-NSP10 interface. In the case of N9P2 and N9R4, Ala108-CA (23%) and Leu104-O (39%) were with the highest H-bond occupancies, respectively. These residues are present at the NSP9 dimerization interface (PDB ID: 6W4B). In the case of NSP15, His235-CB atom was with 31% H-bond occupancy. Thus, the active site residues were involved in H-bond with the phytochemicals. Finally, in the case of NSP3, N3P1 and N3P4, Leu126-CB (31%) and Ile131-CG2 (40%) were with the highest H-bond occupancies. These residues interact with Adenosine-5-Diphosphoribose in the crystal structure of NSP3 (6W02).

Principal component and dynamic cross correlation matrix analysis of the NSP-phytochemical and repurposed drug complexes

To extract the structural variations in detail for the complexes with the lowest relative free energy of binding (Table 2), principal component analysis (PCA) on the Cα atom was performed and compared to their respective apo-proteins. The results of PCA corroborates with the RMSD, RMSF, Rg and relative binding free energy obtained for the phytochemicals and the repurposed drugs. The first two eigen vectors captured around ~50% of the motions (Fig. 3 ). In the case of NSP10 much larger conformational space was explored by N10P1 and N10R3 (Fig. 3A). Their corresponding residue fluctuation plots showed overall lower RMSF for N10P1 and N10R3 compared to the apo-protein. In the case of NSP9 comparable conformational space were explored by the apo-protein and the complexes (Fig. 3B). Residue fluctuation were lower for N9P2 compared to the apo-protein. However, increase in fluctuations were observed for N9R4 in comparison to the apo-protein. This indicates that binding of the phytochemical was more stabilizing than the repurposed drug. In the case of NSP15, N15P1 explored more conformational space (Fig. 3C). Residue fluctuations were reduced in comparison to the apo-protein, showing creation of a stable environment upon ligand binding. In the case of N12d1 (Fig. 3D), more conformational space was explored by N12d1R2 in PC1. Comparable conformational space was explored by the apo-protein and the N12d1P1. Residue fluctuations were similar for N12d1P1 and N12d1R2 compared to the apo-protein (residues 500–800). Residue fluctuations were also similar in the case of N12d2 (residues 300–400) and N12d3 (200–500) at the ligand binding site in comparison to the apo-protein. In the case of NSP16 (Fig. 3E), larger conformational space was explored by N16P1 compared to the apo-protein. Overall decrease in residue fluctuations were observed for N16P1. In the case of NSP3 similar conformational space was explored by N3P1 and N3P4 compared to the apo protein (Fig. 3F). Slight decrease in residue fluctuations at the C-terminal was observed for N3P1. Residue fluctuations increased for N3P4 compared to N3P1 and apo-protein. Therefore, PCA showed overall reduction in residue fluctuations upon ligand binding for all the NSP complexes.

Fig. 3

Comparative PCA of the 8 NSP binding sites with the phytochemicals and the repurposed drugs along with their apo-proteins (OP). The phytochemicals with the most negative free energy of binding have been illustrated here. The projection of PC1 on PC2. The continuous color spectrum from blue to whit to red represents simulation time. The initial timescale represented by blue, intermediate by white and final by red. The RMSF of residue contribution to PCA. (black: PC1, blue: PC2). RES COBTRIB: Residue contribution to PCA (RMSF). For abbreviations please refer main text and Table 1. A) NSP10, B) NSP9, C) NSP15, D) NSP12d1, NSP12d2 and NSP12d3, E) NSP16, F) NSP3.

Further, dynamic cross-correlation matrix was computed by using the coordinates of Cα atoms from the trajectories. Dynamic cross correlation plots of the NSP complexes reveal that anti-correlation (indicative of ligand binding) slightly increased and correlation (indicative of stability upon ligand binding) increased at the binding sites on binding of phytochemicals and repurposed drugs in comparison to the apo-proteins. In case of NSP10 complex (Fig. 4 A), increase in correlation was more for N10P1 than N10R3.In case of NSP 16 (Fig. 4B), N16P1 showed increased anti-correlation and correlation in comparison to the apo-protein. In the case of NSP15 (Fig. 4C), overall increase in correlated motion was observed. In case of NSP12- N12d1, N12d2 and N12d3 there was overall decrease in anti-correlated motions. In the case of NSP9 (Fig. 4D) slight increase in anti-correlation was observed on ligand binding. In the case of NSP3 (Fig. 4F), both anti-correlation and correlation increased on ligand binding compared to the apo-protein. Conclusively it can be inferred that ligand binding created a more stable environment in all the NSPs with respect to the apo-proteins.

Fig. 4

C) DCCM plots for 8 NSP binding sites with the phytochemicals and the repurposed drugs along with their apo-proteins (OP). The positive value represents the positively correlated motions (cyan), while negative values represent the anti-correlated motions (pink). For abbreviations please refer main text and Table 1. A) NSP10, B) NSP16, C) NSP15, D) NSP12d1, NSP12d2 and NSP12d3, E) NSP9, F) NSP3.

To perform pathway enrichment analysis for the 22 phytochemicals, SwissTarget Database was used to retrieve the top 15 human protein targets for each phytochemical. The SARS-CoV-2 human protein targets were retrieved from literature (Chakrabarty et al., 2020; Zhou et al., 2020). Enrichr and REVIGO analysis (Supek et al., 2011; Kuleshov et al., 2016) of Gene Ontology enrichment of biological processes for phytochemicals and SARS-CoV-2 with respect to p-value calculated with Fisher's exact test has been presented in Fig. 5 A-B respectively. Among the top 10 GO biological processes cytokine mediated pathway and phosphorylation were with the most significant p-value for the SARS-CoV-2 and phytochemicals interacting with human protein targets, respectively. Among the top 100 enriched GO terms- the biological processes common to both SARS-CoV-2 and phytochemicals were cytokine-mediated signaling pathway (GO:0019221), cellular response to cytokine stimulus (GO:0071345), positive regulation of macromolecule metabolic process (GO:0010604), interleukin-12-mediated signaling pathway (GO:0035722), cellular response to interleukin-12 (GO:0071349), MAPK cascade (GO:0000165), response to organic cyclic compound (GO:0014070), negative regulation of phosphoprotein phosphatase activity (GO:0032515) and JAK-STAT cascade (GO:0007259).

Fig. 5

Enrichr analysis of biological processes enriched by the human proteins targets of A) phytochemicals and the B) SARS-CoV-2. The processes are colored based on the dispensability score obtained from REVIGO- higher scores are represented by darker shades of red (Supek et al., 2011).C) Top 10 KEGG pathway and VirusMint enrichment for phytochemicals. D) Top 10 KEGG and VirusMint enrichment for SARS-CoV-2. The bars in panel C and D represents the p-values computed using the Fisher's exact test. The longer and lighter colored bars reflects that the term is more significant.

Interestingly, Enrichr analysis of KEGG pathway revealed that Influenza A and Epstein-Barr virus infection pathways were enriched for SARS-CoV-2 and the neuroactive ligand-receptor interactions were enriched for the phytochemicals. Again through VirusMint database enrichment analysis it was observed that the human protein targets of both SARS-CoV-2 and the phytochemicals interacted with viral proteins from Human Immunodeficiency Virus I with the most significant p- value calculated with Fisher exact test (Fig. 5C-D).

Discussions

COVID-19 pandemic causative agent SARS-CoV-2 is an RNA virus which codes for 16 non structural proteins (NSPs) from its ORF1a/b (da Silva et al., 2020). Among the 16 NSPs, 6 NSP have been well characterized (NSP3, NSP9, NSP10, NSP12, NSP15 and NSP16) to play pivotal role in viral replication (Astuti and Ysrafil, 2020). Therefore, these NSPs qualify as potential drug targets of SARS-CoV-2. It is believed that combination of drugs that can simultaneously target these NSPs will have higher possibility in combating COVID-19 infection. In our previous work we identified 18 different steroidal lactones from Withania somnifera and one each phytochemical from Solanum nigrum, Nigella sativa and Azadirachta indica, by virtual screening from 55 Indian medicinal plants, can be used to target these NSPs (Parida et al., 2020a). In the present work binding potential of (22 phytochemicals and 10 repurposed drugs docked to 6 NSPs from SARS-CoV-2) were evaluated by employing 100 ns molecular dynamics simulation analysis. The results revealed the phytochemicals with the lowest binding energy against the 6 NSP targets were majorly from Withania somnifera. These phytochemicals thus have the potential for being developed into a concoction of natural drugs which can be used as prophylactic therapy for containing COVID-19 infections. Interesting observations were made while exploring the dynamic binding of the phytochemicals and repurposed drugs by RMSD, RMSF, Rg, PCA, DCCM, relative free energy of binding and hydrogen bond occupancy analyses. RMSD and RMSF were analysed as they reflect dynamic stability of the trajectory, necessary for obtaining good binding affinities (Dubey et al., 2013). Rg analysis reflects protein compactness which generally increase if binding with ligand forms a stable environment. High hydrogen bond occupancy is another parametric indicator of strong binding affinity (Khan et al., 2020). PCA and DCCM reflect overall conformational variance and atomic motions correlation upon ligand binding respectively (Hünenberger et al., 1995; Islam et al., 2019).

It was reported that in NSP3 (ATP ribose phosphatase of SARS-CoV-2), targeting the Adenosine-5-Diphosphoribose substrate binding site may result in enzyme inhibition (Wu et al., 2020). This work showed by MD simulation analysis of NSP3 complexes, 24,25-dihydrowithanolide D and 27-Deoxy-14-hydroxywithaferin A formed more stable complex among the 4 phytochemicals and one repurposed drug (Baricitinib) analysed. Phe132, Ile131 in its active site pocket which forms hydrogen bond (<~3 Å distance) with the Adenosine-5-Diphosphoribose substrate, were with the lowest relative free energy of binding. Lowest RMSF fluctuations were also obtained for residues 126–132. Leu126 and Ile131 were also with the highest H-bond occupancies. Leu126 is also involved in non-bonded contacts with Adenosine-5-Diphosphoribose C4′ and C5′ atoms at a distance of ~3 Å. Both PCA and DCCM analysis revealed that ligand binding created a more stable environment.

It was also reported that inhibition of dimerization of NSP9 (important in viral replication) can lead to therapeutic benefits by affecting its RNA binding and SARS-CoV viral proliferation (Sutton et al., 2004; Frieman et al., 2012). It was reported that 97% sequence identity is shared by NSP9 from SARS-CoV and SARS-CoV-2 (Littler et al., 2020). Interestingly, analysis of NSP9 complexes revealed that the phytochemical 12-Deoxywithastramonolide and the repurposed drug Baricitinib were with the lowest relative free energy of binding. RMSD, RMSF, Rg were also stable for the complexes. PCA showed that the phytochemical complex has lower residue fluctuations compared to the repurposed drug. DCCM also showed slight increase in anti-correlation was observed on ligand binding. Interestingly a lower dip in ∆RMSF was obtained for Phe76. As Phe76 of one chain in NSP9 is an important residue involved in inter-protein contact with Pro7 of another chain (PDB ID: 6w4b). Thus, its interaction with the inhibitors may affect protein dimerization. Further, interaction with the interface residues, Ala108 and Leu104 with the high H-bond occupancies and Leu113, Phe76 and Leu104 high affinity binding reflected by their binding free energy can affect dimerization with these compounds are primed with NSP9. Leu 104 is also part of the conserved protein-binding motif Gly100-Gly105 (Littler et al., 2020).

It was also reported that targeting the NSP10-NSP16 binding interface may result in attenuation of SARS-CoV-2 infection. NSP 10 is known to stimulate methyl transferase activity of NSP 16 in SARS-CoV (Bouvet et al., 2014). Therefore, similar role is expected in SARS-CoV-2 due to its shared sequence homology. It can be postulated that inhibiting the complex formation of NSP10-NSP16 will aid in containing the SARS-CoV-2 infection. 27-Hydroxywithanolide B from Withania somnifera and Solvanol from Solanum nigrum were with the lowest relative free binding energies with NSP10 and NSP16 respectively. The repurposed drug Baricitinib also showed low relative free binding energy with NSP10. The RMSD, RMSF, Rg, PCA and DCCM analysis also corroborated with the relative free binding energy results reflecting on the stability of the compounds over 100 ns trajectories. The residues, with low relative binding free energy and involved with H-bond with the compounds, lie at the NSP10-NSP16 binding interface and are also involved in zinc coordination (His4336). Zinc coordination was reported to stimulate NSP16 activity in SARS-CoV (Bouvet et al., 2014). Similarly, residues in NSP16 at the NSP10-NSP16 interface were observed to form H-bonds with the compounds and have low relative free energy of binding. Thus, phytochemicals interacting with these residues may interfere with the NSP10-NSP16 complex formation.

Further, targeting three sites of NSP12- RNA-dependent RNA polymerase (RdRp) involved in viral replication, revealed interesting results. The three sites targeted were the catalytic pocket, the NSP7-NSP12 and NSP8-NSP12 interfaces. The NSP12-NSP8 interface was realized to be the best binding site by lowest free energy obtained for Withaferin A and Lopinavir. Docking of Nsp8 to NSP12 was reported to be important in synthesis of the primer required for Nsp12-RdRp RNA synthesis (Wu et al., 2020). All the other trajectory analyses could be positively correlated to the computed relative free energy of binding for the NSP12 complexes.

MD simulation results of NSP15 (uridylate specific endoribonuclease) (Kim et al., 2020), showed that somniferine was with the lowest free energy of binding. The other trajectory analyses agreed with the obtained relative free energy of binding. The active site residues (His235) were involved in H-bond with the phytochemical and positively contributed to lowering the relative free energy of binding. However, NSP15 showed the highest relative binding energies compared to the rest of the NSP analysed through this work. Thus, it can be assumed to be a poor SARS-CoV-2 target in binding to these phytochemicals and repurposed drugs compared to the other NSPs.

Interestingly, through trajectory analysis for H-bond occupancies, lower occupancy of hydrogen bonds was observed in all the trajectories (~ <40%). This reveals that van der Waal's interactions were more pronounced in maintaining residue-residue contacts in the NSP complexes rather than electrostatic interactions. This is also reflected by the much lower values of van der Waal's (vdW) energy than the electrostatic energy, contributing to the overall relative free binding energy. Therefore, it is seen that hydrophobic interactions are dominant in all complexes. High negative value of vdW energy represents greater hydrophobic interaction between the NSP complexes. This can be attributed to the presence of high number of non-polar residues at the binding sites of the 6 NSPs. Further, due to the presence of ergostane framework of steroidal lactone ring structures in the phytochemicals from Withania somnifera, that formed stable complexes with the NSPs, can also be related to lower hydrogen bond occupancy but higher prevalence of hydrophobic interactions. Interestingly, the repurposed drugs with good binding affinity also resembles the phytochemicals in having steroidal rings in their structures and a greater number of hydrogen bond acceptors than donors.

Additionally, enrichment of human biological pathway analysis was performed for the phytochemicals in comparison to the SARS-CoV-2 human protein targets. Interestingly the results revealed that pathways related to immunity, RNA virus mediated disease pathways and apoptosis were observed for both SARS-CoV-2 and the phytochemicals. This strongly implies that the action of phytochemicals in these pathways will be therapeutic in case of controlling COVID-19 infections. For example, biological process enrichment analysis revealed that COVID-19 targets the cytokine mediated pathway. This corroborates with previous report of COVID-19 association with “cytokine storms” including interleukin-6 (IL-6), interleukin 12 (IL-12), tumor necrosis factor alpha (TNFα) and chemokines- causing damage to the host (Hirano and Murakami, 2020). The biological process of phosphorylation associated with the phytochemicals corroborates with research where phytochemicals were reported to be protein modulators by potentiating signal transduction cascades (Frigo et al., 2002). Again, the MAPK cascade was a common target for both SARS-CoV-2 and phytochemical associated human protein targets. The MAPK pathway is expressed in response to external stress such as viral infection and results in activation of p38 MAPK. This pathway also takes part in cell death and was reported to be activated by SARS-CoV (Mizutani, 2010). Phytochemicals have been previously reported to modulate MAPK signaling pathway mediated apoptosis (Kaur et al., 2006). Further inhibition of anus-kinase/Signal transducer and activator of transcription (JAK/STAT) cascade by drugs like Baricitinib and Ruxolitinib has been reported to control SARS-CoV-2 infection (Magro, 2020). Thus, phytochemicals acting on this pathway can result in containing the SARS-CoV-2 infection. Furthermore, KEGG pathways enriched were Influenza A andneuroactive ligand-receptor interactions for SARS-CoV-2 targets and the phytochemical targets, respectively. Association of SARS-CoV-2 protein targets with Influenza virus disease pathway was expected as both lead to severe respiratory diseases (van den Brand et al., 2014) reflecting the analysis accuracy. Interestingly, neuroactive ligand-receptor interactions are responsible to transmit signals for cell growth and survival. This pathway was reported to interact with Ebola virus transcription (Yu et al., 2018). Therefore, phytochemicals interacting with this pathway can be strongly assumed to control SARS-CoV-2 infections. Further, enrichment of VirusMint database resulted in human proteins that are targeted by Human Immunodeficiency Virus I for both the identified SARS-CoV-2 and phytochemicals associated human protein targets. "Cloaked similarity" was reported for SARS-CoV and HIV-I infections (Kilger and Levanon, 2003). Also there have been sporadic reports of phytochemicals in the control of HIV-I infections (Kashiwada et al., 2005). Further, anti‐HIV drugs have been proposed to be effective against SARS‐CoV‐2 (Martinez, 2020). This reflects that these phytochemicals can act as antivirals in not only controlling SARS-CoV-2 but with potential in acting against RNA viruses. Another interesting observation was negative regulation of phosphoprotein phosphatase activity being enriched for both SARS-CoV-2 and phytochemicals. Viruses use phosphoprotein phosphatase activity especially Phosphoprotein Phosphatase 2A to specifically weaken key survival pathways of their hosts (Guergnon et al., al.,2011). Thus, involvement of this biological process with SARS-CoV-2 and phytochemicals indicates that Phosphatase 2A can be a potential drug target for controlling SARS-CoV-2 infections. A novel observation about phytochemical-SARS-CoV-2–human systems was that the human proteins targeted were enriched in the regulation of metabolic processes. Metabolic processes like glycolysis have been reported to be targeted by viruses (DurmuşTekir et al., 2012). Therefore, metabolic pathways should be studied further to understand SARS-CoV-2 infection process in order to comprehend the effects of the phytochemicals in controlling the SARS-CoV-2 infection.

Conclusions

The ongoing SARS-CoV-2 pandemic makes us painfully realize that treating the infection is still a far-fetched goal. Therefore, any contribution in this direction may lead to development of a treatment regime for controlling the infection at an early stage. In this work stability of phytochemicals and repurposed drug complexes with 6 NSPs from SARS-CoV-2 were analysed by 100 ns MD simulations. Overall, the simulation run confirms the stability of interaction of phytochemicals from Withania somnifera with the 8 NSPs binding sites (three sites belonging to NSP12. Further, the important residues interacting with the phytochemicals showed that phytochemicals and repurposed drugs with steroidal moieties in their chemical structures formed stable interactions with the NSPs. Again, Human protein pathway analysis revealed that pathways related to immunity, RNA virus mediated disease pathways and apoptosis were enriched for both SARS-CoV-2 and the phytochemicals. Thus it can be said that the multi-potency of these phytochemicals can be used to regulate COVID-19 infection by targeting these pathways along with the NSPs. Conclusively this work warrants further study from the perspective of experimental validation to test the efficacy of these phytochemicals and the repurposed drug for the prevention and treatment of COVID-19.

CRediT authorship contribution statement

Pratap Kumar Parida: Conceptualization, Methodology, Software, Formal analysis. Dipak Paul: Data curation, Writing - original draft. Debamitra Chakravorty: Conceptualization, Visualization, Investigation, Supervision, Validation, Writing - review & editing.

Declaration of Competing Interest

No potential conflict of interest was reported by the authors.