Computational and network pharmacology studies of Phyllanthus emblica to tackle SARS-CoV-2.

Background: Since December 2019, SARS-CoV-2 had been a significant threat globally, which has accounted for about two million deaths. Several types of research are undergoing and have reported the significant role of repurposing existing drugs and natural lead in the treatment of COVID-19. The plant Phyllanthus emblica (Synonym-Emblica officinalis) (Euphorbiaceae) is a rich source of vitamin C, and its use as an antiviral agent has been well established. Purpose: The present study was undertaken to investigate the potency of the several components of Phyllanthus emblica against three protein targets of 2019-nCoV viz. NSP15 endoribonuclease, main protease, and receptor binding domain of prefusion spike protein using molecular docking and dynamics studies.
Methods: The docking simulation studies were carried out using Schrödinger maestro 2018-1 MM share version, while dynamics studies were conducted to understand the binding mechanism and the complexes' stability studies.
Results: Out of sixty-six tested compounds, Chlorogenic acid, Quercitrin, and Myricetin were most effective in showing the highest binding energy against selected protein targets of SARS-CoV-2. The network pharmacology analysis study confirmed these compounds' role in modulating the immune response, inflammatory cascade, and cytokine storm through different signaling pathways.
Conclusion: Current pharmacoinformatic approach shows possible role of Phyllanthus emblica in the treatment and management of COVID-19.

Introduction

In December 2019, a unique and unprecedented worldwide pandemic arose from the Wuhan city of China, known as the ‘Severe Acute Respiratory Syndrome Coronavirus-2’ (SARS-CoV-2) and had become a significant peril to the world to high mortality and morbidity rates. The coronaviruses are single-stranded RNA (genome size-26-32-kilobases), which plays a dangerous role in the initial RNA synthesis of the infectious cycle and acts as a substrate packaging into the progeny virus (Patil et al., 2020). In coronaviruses, two-thirds of the genome encodes a replicase polyprotein processed by viral protease, which cleaves into 16 non-structural proteins (NSPs) and is involved in transcription and replication (Chikhale et al., 2020a). The non-structural protein 15 of SARS-CoV-2 (Nsp15) is a nidoviral RNA uridylate-specific endoribonuclease (NendoU), and the exact functional relevance of Nsp15 remains unknown. Recent studies claimed that the NendoU activity of NSP15 is mainly attributed to protein interference with innate immune response and is considered an essential segment for the natural progression of coronaviruses (Patil et al., 2020). Another NSP that helps enter 2019-nCoV into the host cells is spike (S) protein composed of two subunits S1 and S2. This S protein plays a vital role in receptor recognition and the host cell membrane's fusion process. During viral infection in the host cells, the target cell proteases activate the S protein by cleaving it into two subunits, which are essential to activate the membrane fusion domain after viral entry into target cells (Hoffmann et al., 2020).

The current focus is the development of novel therapeutics including natural antivirals and specific vaccines. In-silico tools have employed a diverse set of computational approaches to understand the relative performance of the predictions for SARS-CoV-2 by repurposing existing drugs and the alleged role of natural lead in treating COVID-19 (Chikhale et al., 2020a, 2020b, 2021; Patil et al., 2020; Khanal et al., 2020, 2021a, 2021b). However, it may take months or years to develop such effective treatments, and hence the quest for prompt treatment is the utmost need of the hour. During the 2003 severe acute respiratory syndrome (SARS) pandemic, traditional medicinal plants' efficacy has been established (Li et al., 2005). Consequently, based on historical and traditional evidence, researchers have started investigations on medicinal plants for SARS-CoV-2 (Chikhale et al., 2020a, 2020b, 2021).

The medicinal plants, with their therapeutic prominence are a gift to humankind to attain a healthy life. Phyllanthus emblica L. (syn. Emblica officinalis Gaertn.) (Euphorbiaceae), commonly-known as Indian Gooseberry or Amla, is a remarkable tropical south-east Asian shrub, cultivated throughout India. Moreover, it is reported for its anti-inflammatory, antipyretic, antioxidant, anticancer, anti-hyperlipidemic, adaptogenic, anti-diabetic, nootropic, antimicrobial and immunomodulatory, anti-bacterial potential (Variya et al., 2016;). Antiviral efficacy of P. emblica in the treatment of HIV (El-Mekkawy et al., 1995), Coxsackie VB3 (Liu et al., 2009), hepatitis B virus (Xiang et al., 2010), and herpes simplex virus (Xiang et al., 2011) has evidenced with scientific reports. Out of various parts of E. officinalis, the fruits have a significant place in Rasayana to treat several infectious and non-infectious diseases (Variya et al., 2016).

The present research investigation highlights phytoconstituents from the Indian gooseberry using an in-silico approach targeting three different SARS-CoV-2 proteins viz. NSP15 endoribonuclease, main protease, and receptor-binding domain (RBD) of the prefusion spike protein.

Material and methods

Molecular docking studies

Schrödinger Glide SP module was used for Molecular docking studies of a total of sixty-six phytoconstituents of P. emblica with the experimentally solved crystal structures of NSP15 endoribonuclease (PDB: 6W01), SARS-CoV-2 spike RBD (PDB: 6M0J), and SARS-CoV-2 main protease (PDB: 6WNP) having the resolution of 1.9, 2.45 and 1.45Å, respectively (Supplementary File) (Chikhale et al., 2020a).

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

The AMBER18 software package was used for MD simulations, and ligands were parameterized with ANTECHAMBER. The prepared protein-ligand complexes were subjected to 100 ns MD simulations on Nvidia V100-SXM2-16GB 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).

Network construction and analysis

The down-regulated and up-regulated protein-based targets of phytoconstituents were retrieved from DIGEP-Pred at the pharmacological activity (Pa)>0.5. The complete proteins were then queried in the STRING (https://string-db.org/) to predict the protein-protein interaction and GO analysis (Gene Ontology Consortium, 2004) for cellular and molecular functions and the biological spectrum. Similarly, the probably regulated pathways were identified concerning KEGG pathways. Also, the network interaction of bioactive with proteins and regulated pathways was constructed using Cytoscape (Shannon, 2003). The network was treated as directed and evaluated by setting map node size from “low values to small sizes” and map node color from “low values to bright colors” based on edge count for both.

Results and discussion

Molecular docking

To assess the possible potential and understand the possible mechanism of sixty-six phytoconstituents molecular docking simulations were carried out on three proteins i.e. NSP15 endoribonuclease (PDB: 6W01), main protease (PDB: 6M0J), and receptor-binding domain (RBD) of prefusion spike protein (PDB:6WNP). Based on the docking score (Supplementary Table S1), only the best five tested ligands and the reference drug with respect to different proteins have been discussed in the present investigation.

NSP15 endoribonuclease bears a catalytic C-terminal domain, a hexamer, responsible for the specific cutting of double-stranded RNA substrates through endoribonuclease. Every monomeric unit comprises ~345 amino acids folded into three domains: N-terminal, middle, and nidoviral RNA-specific endoribonuclease (NendoU) C-terminal catalytic domain. The two anti-parallel β-sheets of NendoU C-terminal catalytic domain contain six key amino acids viz. His 235, His 250, Lys 290, Thr 341, Tyr 343, and Ser 294. Among them, His 235, His 250, Lys 290 is a catalytic triad, His 235 functions as a general acid, His 250 works as a basis, while Ser 294 and Tyr 343 are found to regulate U specificity (Kim et al., 2020, p. 2). The docking study was carried out on three proteins to check the potency and possible 62 herbal ligands' possible mechanism.

The docking results of crystal structure of NSP15 endoribonuclease from SARS CoV-2 revealed that chlorogenic acid was adequately positioned into the binding pocket assembled by polar Gln 245, Thr 341, Ser 294, Asn 278 charged Hip 235, Hip 250, Lys 290, Lys 345 hydrophobic Leu 346, Val 292, Tyr 343 and Gly 248 amino acids with docking score -8.397 kcal/mol. The hydroxyl group of the terminal phenyl ring of chlorogenic acid exhibited bi-furcated H-bonding with Leu 346 and Asn 278, carboxyl group of cyclohexane ring also exhibited bi-furcated H-bonding with (protonated His) Hip 235, Hip 250, Lys 290, Thr 341 and carboxyl group of ester linkage showed H-bonding with Tyr 343. The carboxyl group of cyclohexane ring formed a salt bridge with Lys 290 (Supplementary Fig. S1a). The hydroxyl group of the terminal phenyl ring of phyllanemblinin D exhibited bi-furcated H-bonding with Leu 346, Asn 278, oxygen of oxane ring showed H-bonding with Tyr 343, and the carboxyl group depicted H-bonding interaction with Gln 245, Gly 248, Hip 235, Thr 341 amino acids with docking score-8.341 kcal/mol. The isochromene ring of phyllanemblinin D exhibited π-cation interaction, π-π stacking with Hip 235, and carboxyl group formed a salt bridge with Hip 250 (Supplementary Fig. S1b). Chebulinic acid showed H-bonding with Val 292, Gln 245, and bi-furcated H-bonding with Leu 346, Asn 278 amino acids with docking score-8.334 kcal/mol. The carboxyl group depicted H-bonding interaction with Thr 341, Hip 250, Hip 235, and Gly 248 and formed a salt bridge with Lys 290 (Supplementary Fig. S1c). The carboxyl group of punigluconin exhibited H-bonding interaction with Thr 341 and Hip 250, carboxyl of ester linkage depicted H-bonding interaction with Hip 235, Gly 248, Lys 290, and hydroxyl group of phenyl ring showed H-bonding interaction with Gln 245, Glu 340 amino acids with docking score-8.014 kcal/mol. The carboxyl group formed a salt bridge, and the phenyl ring exhibited π-cation interaction with Hip 235 (Supplementary Fig. S1d). The carboxyl group of isochromene ring of phyllanemblinin F exhibited bi-furcated H-bonding interaction with Hip 250, Thr 341and the oxygen of isochromene ring showed H-bonding interaction with Hip 235 with docking score -7.709 kcal/mol (Supplementary Fig. S1e). The hydroxyl group of tetrahydrofuran ring of Remdesivir showed H-bonding with Val 292, and oxygen of phosphoramide exhibited bi-furcated H-bonding with Hip 250, Thr 341 amino acids with docking score-5.94 kcal/mol (Supplementary Fig. S1f). Glu 267 amino acid revealed H-bonding with a hydroxyl group of a terminal amino chain. It formed a salt bridge with a quaternary amine group of hydroxychloroquine with a docking score -4.45 kcal/mol (Supplementary Fig. S1g). The amide group of Lopinavir exhibited H-bonding with Val 292, Hip 235, and hydroxyl showed H-bonding with Lys 290. The phenyl ring demonstrated π-cation interaction with Hip 235 amino acid with docking score -4.382 kcal/mol (Supplementary Fig. S1h).

The second docking study was done on the crystal structure of the SARS-CoV-2 spike receptor-binding domain bound with ACE2. SARS-CoV-2 uses a homotrimeric glycoprotein spike to get entry into the host cells ACE-2 receptor. Such a binding relationship is enhanced by attaching the subunit S1 to the host cell receptor and converting the subunit S2 to a highly stable postfusion conformation. Besides, the S1 subunit receptor-binding domain (RBD) consists of 5 twisted β sheets β1, β2, β3, β4, and β7 that are anti-parallel to each other. There was an extended inclusion of β4 to β7 containing some α loops called receptor binding motif (RBM). This RBM includes most of the residues necessary to link n-COVID-19 to ACE-2. Recent research reported that of all the residues Arg 319 to Phe 541, only 17 residues Lys 417, Gly 446, Tyr 449, Tyr 453, Leu 455, Phe 456, Ala 475, Phe 486, Asn 487, Tyr 489, Gln 493, Gly 496, Gln 498, Thr 500, Asn 501, Gly 502, Tyr 505 are essential for the attachment of ACE-2 to the 20 residues of the ACE-2 N-terminal peptidase domain. Among these 17 residues Gln 493, Asn 501, Tyr 449, Tyr 489, and Tyr 505 are strongly linked by H-bonding and Lys 417 by salt bridge interaction (Lan et al., 2020).

Myricetin was found to interact with catalytic site constructed by polar Thr 376, Ser 375 charged Lys 378, Asp 405, Arg 408 hydrophobic Val 433, Tyr 380, Val 407, Tyr 508, Val 503, Ala 411, Ile 410, and Gly 404 amino acids with docking score -6.782 kcal/mol. The hydroxyl group of chromenone ring exhibited bi-furcated H-bonding with Ile 410, Tyr 380, and the hydroxyl group of phenyl ring showed H-bonding with Ser 375, Gly 404. The chromenone ring revealed two π-cation interactions with Lys 378 (Supplementary Fig. S2a). Quercetin also exhibited bi-furcated H-bonding with Tyr 380, Ile 410, and the hydroxyl group of phenyl ring showed H-bonding with Gly 404, Tyr 508 amino acids with docking score -6.449. The chromenone ring exhibited two π-cation interactions with Lys 378 (Supplementary Fig. S2b). Phyllaemblicin F indicated bi-furcated H-bonding with Thr 415, Glu 406. The other hydroxyl group of phyllaemblicin F showed H-bonding with Asp 420, Asp 405. In contrast, a carboxyl group showed H-bonding interaction with Tyr 421, Asn 460, Arg 403 amino acids with docking score -6.127 kcal/mol (Supplementary Fig. S2c). Naringenin has chromenone ring substituted by hydroxyl group showed H-bonding interaction with Tyr 380, Ile 410 amino acids with docking score -5.708 kcal/mol. The chromenone ring exhibited π-cation interactions with Lys 378 (Supplementary Fig. S2d). The chromenone ring of quercetin-3-L-rhamnoside substituted by hydroxyl group exhibited H-bonding interaction with Glu 406, Arg 403. The hydroxyl group of oxane ring showed bi-furcated H-bonding interaction with Ser 494, Tyr 453 with docking score-5.607 kcal/mol. The chromenone ring of quercetin-3-L-rhamnoside demonstrated π-cation interaction with Arg 403 and π-π stacking with Tyr 505 (Supplementary Fig. S2e). The hydroxyl group of tetrahydrofuran ring of Remdesivir showed H-bonding with Glu 406, oxygen of phosphoramide exhibited bi-furcated H-bonding with Arg 403, Tyr 505, and carboxyl group showed H-bonding with Tyr 453 amino acids with docking score-4.685 kcal/mol (Supplementary Fig. S2f). The carboxyl group and nitrogen of tetrahydro-pyrimidinone ring of Lopinavir exhibited H-bonding with Gln 493, Ser 494, and hydroxyl group depicted H-bonding with Tyr 449 amino acids with docking score-4.62 kcal/mol (Supplementary Fig. S2g). The hydroxyl group of the terminal amino chain of quinoline exhibited H-bonding with Asn 450. The quinoline ring's quaternary nitrogen showed H-bonding with Glu 484 amino acids with a docking score-4.6 kcal/mol (Supplementary Fig. S2h).

The third docking study was performed on the X-ray structure of SARS-CoV-2 main protease bound to Boceprevir. The SARS-CoV-2 main protease (3CLpro) comprises approximately 306 amino acids, responsible for coronavirus replication and polypeptide processing into functional proteins. Each 3CLpro consists of three domains, with domains I (residues 8-101) and II (residues 102-184) having an anti-parallel β-barrel structure. Domain III (201-303 residues) has five α-helices, arranged in a mostly anti-parallel globular cluster and further connected to Domain II through an extended loop region (residues 185-200). 3CLpro has a Cys–His catalytic dyad, and the substrate-binding position is in a cleft between Domain I and Domain II (Anson and Mesecar, 2020). Quercitrin was found to interact with catalytic site assembled by polar Thr 26, Ser 144, Thr 25, Thr 190, Asn 142, Hie 172, Hie 164, Hie 41, His 163, Gln 192, Gln 189 charged Asp 187, Arg 188, Glu 166 hydrophobic Leu 27, Leu 141, Cys 145, Phe 140, Met 165, Met 49 and Gly 143 amino acids with docking score -9.043 kcal/mol. The hydroxyl group of chromenone ring exhibited H-bonding with Arg 188. The hydroxyl group of the oxane ring showed bi-furcated H-bonding with Thr 26, Gly 143, Asn 42, and the hydroxyl group of phenyl ring showed H-bonding with Glu 166. The ether linkage between chromenone ring and the oxane ring exhibited H-bonding with Cys 145 (Supplementary Fig. S3a). The pyran and chromenon ring compounds were also reported for their activity against HIV (De Clercq, 2000; Park et al., 2008). A US patent (US005843990A) suggests using pyran-chromenone compounds to prevent viral growth or replication (Baker, 1998). The hydroxyl group of oxane ring of Phyllaemblicin C exhibited H-bonding with Glu 166, Asn 142, Hie 41, hydroxymethyl group of oxane ring showed H-bonding with Hie 41, the phenyl hydroxyl group of Phyllaemblicin C showed H-bonding with Glu 166, and oxygen of oxane ring showed H-bonding interaction with Cys 145 amino acids with docking score-8. 881 kcal/mol (Supplementary Fig. S3b). The hydroxyl group of Rutin's oxane ring exhibited H-bonding with Glu 166, Gln 189, Ser 46, and phenyl hydroxyl group showed H-bonding with Thr 26 Asn 142 amino acids with docking score-8. 831 kcal/mol (Supplementary Fig. S3c). The chromenone ring of Avicularin exhibited H-bonding with Arg 188, the hydroxyl group of phenyl ring showed H-bonding with Glu 166, the hydroxyl group of tetrahydrofuran ring showed H-bonding with Thr 26, Asn 142, and an ether linkage between chromenone ring and tetrahydrofuran ring exhibited H-bonding with Cys 145 amino acids with docking score-8.69 kcal/mol (Supplementary Fig. S3d). Chromenone ring of Quercetin-3-L-rhamnoside substituted by hydroxyl group demonstrated H-bonding interaction with Arg 188. The hydroxyl group of phenyl ring showed H-bonding interaction with Glu 166. The tetrahydrofuran ring substituted with the hydroxyl group showed H-bonding interaction with Thr 26, Asn 142 amino acids with docking score-8.644 kcal/mol and ether linkage between chromenone ring and tetrahydrofuran ring exhibited H-bonding with Cys 145 (Supplementary Fig. S3e). The two secondary amine and hydroxyl groups of co-crystalized ligand boceprevir exhibited H-bonding with Glu 166. His 164 and terminal hydroxyl group showed bi-furcated H-bonding with Gly143, Cys 145 amino acids with docking score-8.405 kcal/mol (Supplementary Fig. S3f). Remdesivir has a pyrrolotriazine ring substituted with an amine group that showed H-bonding with Thr 26, and phosphoramide's oxygen exhibited H-bonding with Glu 166 amino acids with docking score-7.766 kcal/mol. The phenyl ring of Remdesivir depicted π-π interactions with Hie 41 (Supplementary Fig. S3g). The amide group of Lopinavir showed H-bonding with Glu 166, Cys 145 amino acids with docking score -6.793 kcal/mol (Supplementary Fig. S3h). The hydroxyl group of a terminal amino chain of hydroxychloroquine revealed H-bonding with Hie 41. The quaternary nitrogen of the quinoline ring showed H-bonding with Arg 188 amino acids with a docking score -5.411 kcal/mol. The terminal quaternary nitrogen depicted π-cation interactions with Hie 41 (Supplementary Fig. 3h).

Network analysis

Target mining identified the apigenin to regulate the highest number of proteins i.e. 32, in which 14 (NFE2L2, MMP7, NOS2, MMP2, CAT, CTNNB1, KLK3, MDM2, EGLN1, MMP3, CASP8, FLT1, PPARA, and TOP2A) are down-regulated, and 18 are upregulated (KRT1, PLAT, NFE2L2, PLAU, CYP1A1, PRDX2, HMOX1, PRDX4, TIMP1, AR, TP73, SIRT1, EPAS1, TP53I3, VDR, GADD45B, GSS, and PGR). Similarly, Survival of motor neuron 2 (SMN2) was predicted to be majorly upregulated by 30 phytoconstituents. The gene ontology of interacted targets (Supplementary Fig. S4) was identified for 61 cellular components, 83 molecular functions, and 951 biological processes. The top 25 hits of gene ontology are presented in Fig. 1 . Further, KEGG identified 87 pathways to be primarily regulated in which pathways in cancer was majorly modulated at the false discovery rate of 2.68E-10 by regulating 19 proteins (GADD45B, HMOX1, MMP2, IFNG, RARA, MDM2, CCND2, EPAS1, KLK3, NOS2, ESR2, CTNNB1, RAC1, CASP8, EGLN1, AR, NFE2L2, PRKCA, and RXRA) against 515 background proteins (Supplementary Table S2). Similarly, the interaction of phytoconstituents, regulated proteins, and associated pathways are represented in Fig. 2 .

Fig. 1

GO analysis (a) Cellular components, (b) Molecular function and (c) Biological process.

Fig. 2

Interaction of phytoconstituents with their targets and regulated pathways.

The gene ontology of regulated proteins identified the regulation of the multiple pathways that are concerned with infectious and non-infectious diseases. Some pathways involved with viral infections, like human papillomavirus, herpesvirus, and Epstein-Barr virus infection are also modulated. The subjects with lower immunity are more prone to the COVID-19 infection. In the present study, the majority of the regulated pathways like HIF-1, p53, IL-17, PI3K-Akt, FoxO, Wnt, NF-kappa B, TNF, MAPK, Rap1, and AMPK signaling pathways are directly or indirectly concerned with immunity manipulation and regulating inflammation (Johnson and Chen, 2012; Li et al., 2019; Muñoz-Fontela et al., 2016; Palazon et al., 2014; Šedý et al., 2014; Zenobia and Hajishengallis, 2015). Further, the secondary phytoconstituents from P. emblica may also be involved in managing the cytokine storm as the multiple pathways related to chemokines are also regulated (Supplementary Table S1).

During the bioactive-target-pathway network analysis, the apigenin was identified as the lead hit by regulating proteins in which 14 were downregulated (NFE2L2, MMP7, NOS2, MMP2, CAT, CTNNB1, KLK3, MDM2, EGLN1, MMP3, CASP8, FLT1, PPARA, TOP2A) and 18 were upregulated (KRT1, PLAT, NFE2L2, PLAU, CYP1A1, PRDX2, HMOX1, PRDX4, TIMP1, AR, TP73, SIRT1, EPAS1, TP53I3, VDR, GADD45B, GSS, PGR). However, the binding of apigenin with selected SARS-CoV-2 proteins was less than three leads of docking results i.e., chlorogenic acid, Quercitrin, and Myricetin. Although, chlorogenic acid regulated 12 genes in which 4 were down-regulated (CHEK1, FLT1, MDM2, CASP8), and 8 were up-regulated (RARA, NPPB, SMN2, PLAU, RAC1, PLAT, HMOX1, CD14). Quercitrin regulated 19 targets in which 5 were down-regulated (MMP7, CHEK1, NOS2, NFE2L2, CTNNB1), and 14 were up-regulated (TIMP1, PLAU, PLAT, TNFRSF1A, PRDX4, HMOX1, SMN2, TP73, KRT1, NPPB, CBR1, EPAS1, AR, NFE2L2). In contrast, myricetin regulated 29 targets in which 11 were downregulated (NFE2L2, MMP7, CAT, NOS2, CTNNB1, EGLN1, KLK3, MMP2, MMP3, CASP8, MDM2) and 18 were upregulated (HMOX1, NFE2L2, KRT1, PLAT, PLAU, PRDX2, SIRT1, TP73, TIMP1, PRDX4, CYP1A1, TP53I3, EPAS1, AR, TP63, GADD45B,VDR, and SMN2). These results somehow reflect that although apigenin is directly involved in the prime modulation of the multiple proteins in the network and reported for the anti-viral activity for various viral strains (Lv et al., 2014), it may not be directly involved in the inhibition of the coronavirus infection. However, it could play an essential role in immune regulation and possible anti-inflammatory activity in COVID-19 infection as reported for both properties (Hosseinzade et al., 2019). Though chlorogenic acid, quercitrin, and myricetin reflected the comparatively less modulatory effect of multiple proteins in the network, they may not be directly involved in regulating multiple pathways as apigenin does. It suggests the apigenin from E. officinalis in COVID-19 infection could majorly contribute to the immune regulation and possess anti-inflammatory action in the infective tissue, whereas, chlorogenic acid, quercitrin, and myricetin could directly act over the coronavirus to amplify its anti-viral activity.

Molecular dynamics

The molecular dynamics studies (MDS) of target protein-ligand complexes were performed to understand the binding mechanism and the complexes' stability over the period, explaining and modeling its in-vitro or in-vivo efficacy. The top-scoring ligands in complex proteins from the molecular docking experiments were selected for the MDS. Remdesivir, a drug approved by the US-FDA for emergency use in the COVID-19 pandemic, was also studied as the control for these simulation studies. The NSP15 endoribonuclease bound to Chlorogenic acid and Remdesivir were simulated for 100 ns in the explicit solvent model at physiological salt concentration. The NSP15 endoribonuclease protein RMSD converged with lower fluctuations between 1.5 and 2 Å for the first 40 ns, which later fluctuated with a steep fall and rose to 3 Å and remained between 2 and 3 Å for the remaining simulations period (Fig. 3 .1a). The RMSF for individual residues show a higher fluctuation for amino acids (aa) between 35 and 50 of 3.5 Å, but the rest of aa had lower changes (Fig. 3.1b). The bound ligand, chlorogenic acid, had a low RMSD of 0.5 Å for the first 50 ns, which gradually rose for 20 ns to reach and stabilize at 3 Å (Fig. 3.1c). It suggests the protein-ligand complex's stability. The gradual rise in Ligand RMSD means a smooth positional change from its original binding site to a more stabilized position (Fig. 3.1d and 1e). This shifting of chlorogenic acid was stabilized by forming hydrogen bonds between the hydroxy groups and the residue Leu347 (Fig. 3.1e).

Fig. 3

1(a) RMSD of NSP15 endoribonuclease; 1(b) RMSF of NSP15 endoribonuclease; 1(c) RMSD of Chlorogenic acid; 1(d) Initial binding pose of Chlorogenic acid with NSP15 endoribonuclease and; 1(e) Binding pose of Chlorogenic acid with NSP15 endoribonuclease towards the end of the simulation. 2(a) RMSD of NSP15 endoribonuclease; 2(b) RMSF of NSP15 endoribonuclease; 2(c) RMSD of Remdesivir; 2(d) Initial binding pose of Remdesivir with NSP15 endoribonuclease and; 2(e) Binding pose of Remdesivir with NSP15 endoribonuclease towards the end of the simulation.

In the MDS of the NSP15 endoribonuclease-Remdesivir complex, the protein RMSD was more stable with fluctuations between 1 and 3 Å (Fig. 3.2a) compared to the NSP15 endoribonuclease-chlorogenic acid complex. The RMSF for amino acid residues was very similar for the NSP15 endoribonuclease-Remdesivir and NSP15 endoribonuclease-chlorogenic acid complexes (Fig. 3.2b). The Ligand RMSD for Remdesivir converged between 4 and 6 Å from 5 ns to the end of the simulation with fewer fluctuations between 25 and 50 ns (Fig. 3.2c), which could be attributed to the conformational changes in the Remdesivir ligand (Fig. 3.2d and 3.e).

The SARS-CoV-2 main protease-Quercitrin and SARS-CoV-2 main protease-Remdesivir complexes were simulated for 100 ns, and their binding modes were analyzed. The protease in complex with Quercitrin converged slowly to 3 Å during the initial 20 ns and remained stable throughout the simulation with few fluctuations around 30 and 60 ns (Fig. 4 .1a). The RMSF lower than 2.5 Å reflects the high stability of the amino acid residues (Fig. 4.1b). The ligand RMSD for Quercitrin converged and stabilized between 2 and 3.5 Å with occasional fluctuations (Fig. 4.1c). The visual inspection of the trajectory shows the breaking of the hydrogen bond with Arg188 around 20 ns of the simulation, which is also reflected from the ligand RMSD. The formation of a hydrogen bond with Glu166 stabilizes the complex (Fig. 4.1d and e). The protein RMSD of SARS-CoV-2 main protease-Remdesivir shows a low RMSD, between 2 and 3 Å (Fig. 4.2a), the receptor RMSF was also low between 0.5 and 2 Å (Fig. 4.2b), suggesting high stability of the complex. The ligand RMSD of Remdesivir bound to the main protease was between 2 and 3 Å with fluctuations of 1 Å. The visualization of trajectories and the simulation frames shows a hydrogen bond between the ‘cyano’ group of Remdesivir. The Cys145 remains intact with a short bond length keeping the complex intact (Fig. 4.2c, 4.2d, and 4.2e).

Fig. 4

1(a) RMSD of SARS-CoV-2 main protease; 1(b) RMSF of SARS-CoV-2 main protease; 1(c) RMSD of Quercitrin; 1(d) Initial binding pose of Quercitrin with SARS-CoV-2 main protease and; 1(e) Binding pose of Quercitrin with SARS-CoV-2 main protease towards the end of the simulation. 2(a) RMSD of SARS-CoV-2 main protease; 2(b) RMSF of SARS-CoV-2 main protease; 2(c) RMSD of Remdesivir; 2(d) Initial binding pose of Remdesivir with SARS-CoV-2 main protease and; 2(e) Binding pose of Remdesivir with SARS-CoV-2 main protease towards the end of the simulation.

The SARS-CoV-2 receptor-binding domain (RBD) was docked with the compounds under investigation, and remdesivir was included as a standard drug for comparing the binding on the natural products. Remdesivir is a prodrug that acts by interfering with the RNA-dependent RNA polymerase function. However, recently the alternate mechanism of action reported the concentration-dependent effect of Remdesivir. These reports provide for further investigation in the possible alternate binding sites and mechanism of action for Remdesivir. The top-scoring compound, Myricetin bound to RBD, and Remdesivir docked on the RBD were simulated for 100 ns, and the binding modes were analyzed. The RBD domains RMSD was between 1.5 and 2.5 Å, and the RMSD for its aa was below 2.5 Å, which is highly acceptable for a domain region of a protein (Fig. 5 .1a and b). The RBD bound ligand Myricetin had RMSD below 2 Å for the initial 25 ns of the simulation, which later fluctuated between 2 and 5 Å for the rest of the simulation. The visualization of the trajectory and different frames from the trajectory shows the intact hydrogen bond between Ile78 and Myricetin throughout the period. This complex gets more stabilized towards the end of the simulation by forming another hydrogen bond with Gly72 (Fig. 5.1c, d and e).

Fig. 5

1(a) RMSD of SARS-CoV-2 RBD; 1(b) RMSF of SARS-CoV-2 RBD; 1(c) RMSD of Myricetin; 1(d) Initial binding pose of Myricetin with SARS-CoV-2 RBD and; 1(e) Binding pose of Myricetin with SARS-CoV-2 RBD towards the end of the simulation. 2(a) RMSD of SARS-CoV-2 RBD; 2(b) RMSF of SARS-CoV-2 RBD; 2(c) RMSD of Remdesivir; 2(d) Initial binding pose of Remdesivir with SARS-CoV-2 RBD and; 2(e) Binding pose of Remdesivir with SARS-CoV-2 RBD towards the end of the simulation.

The SARS-CoV-2 RBD bound Remdesivir showed initial stability in its RMSD but has several very high fluctuations of about 5-10 Å. The aa residues showed higher RMSF of 0.5 to 3.5 Å over the simulation period suggesting lowered stability or no effect of Remdesivir binding on the RBD (Fig. 5.2a and b). The ligand had very high RMSD due to structural conformation changes and the absence of initial stage interactions. However, towards the end of the simulation, it forms hydrogen bond interaction with Glu74 and Arg71, and possibly these interactions keep the ligand attached to the receptor (5.2c, 5.2d, and 5.2e).

The Molecular Mechanics/Generalized Born Surface Area calculations were performed to calculate the binding free energy of the complexes. The complete trajectories for 100 ns were used for the study, and the results are tabulated in Table 1 . The individual analyses suggest higher binding energy of the Quercitrin-SARS-CoV-2 main protease complex (ΔGbind = -36.27) over the Remdesivir-SARS-CoV-2 main protease complex (ΔGbind = -27.59). It indicates a more substantial and stable drug-receptor complex and better efficacy of the Quercitrin to bind with SARS-CoV-2 main protease. The Myricetin-SARS-CoV-2 RBD complex has higher binding energy (ΔGbind = -17.41) than the Remdesivir-SARS-CoV-2 RBD (ΔGbind = -0.91). This observation shows the low affinity of the Remdesivir towards the SARS-CoV-2 RBD, which would suggest the alternative mechanism of action for remdesivir, as discussed earlier. Similarly, chlorogenic acid has week binding energy towards the target NSP15 endoribonuclease (ΔGbind = -0.42).

Table 1

MMGBSA results for binding energies of selected complexes.

Compounds-Protein	Glide Score Kcal/mol	MM-GBSA*
		ΔEVDW	ΔEELE	ΔGGB	ΔGSurf	ΔGgas	ΔGSol	ΔGbind
Chlorogenic acid (PDB: 6W01)	-8.397	-17.91 (3.72)	-113.92 (26.11)	128.98 (20.69)	-2.90 (0.47)	-126.50 (24.28)	126.08 (20.55)	-0.42 (4.98)
Remdisvir (PDB: 6W01)	-5.94	-36.63 (4.56)	-21.07 (6.89)	39.26 (5.96)	-4.99 (0.61)	-47.52 (8.60)	34.26 (5.77)	-13.31 (5.05)
Quercitrin-(PDB: 6WNP)	-9.043	-35.24 (4.20)	-33.14 (16.06)	46.81 (10.81)	-5.20 (0.23)	-77.88 (15.64)	41.60 (10.57)	-36.27 (5.73)
Remdesivir-(PDB: 6WNP)	-7.766	-46.83 (3.82)	-44.35 (7.82)	58.47 (5.18)	-5.99 (0.30)	-80.07 (7.29)	52.47 (5.08)	-27.59 (4.32)
Myricetin-(PDB: 6M0J)	-6.782	-24.28 (2.46)	-21.62 (7.12)	32.66 (6.10)	-3.89 (0.26)	-46.18 (7.09)	28.77 (6.01)	-17.41 (2.62)
Remdesivir-(PDB: 6M0J)	-4.685	-22.49 (4.33)	-23.38 (11.05)	40.82 (9.18)	-3.59 (0.53)	-38.14 (12.29)	37.23 (4.80)	-0.91 (4.80)

Conclusion

In the present research investigation, a total of sixty-six phytoconstituents of Phyllanthus emblica were subjected for their plausible antiviral properties using in-silico approaches against SARS-COV-2- proteins i.e. NSP15 endoribonuclease, main protease, and receptor binding domain of prefusion spike protein. The docking studies and MDS confirmed the promising potential of Chlorogenic acid, Quercitrin, and Myricetin to inhibit the n-CoV-2 key viral proteins. The network pharmacology analysis demonstrated the involvement of selected phytoconstituents in modulating multiple signaling pathways that could play a significant role in immunomodulation, regulating inflammation, and managing the cytokine storm. The current pharmaco-informatics approach would give the researchers new insights for further experimental work on E. officinalis phytoconstituents to investigate the lead molecules demonstrated in the present study to inhibit SARS-COV-2 vital proteins.

CRediT authorship contribution statement

Rupesh V. Chikhale: Investigation, Software, Methodology, Validation. Saurabh K. Sinha: Resources, Software, Methodology, Validation. Pukar Khanal: Investigation, Software, Methodology, Validation. Nilambari S. Gurav: Data curation, Formal analysis, Visualization, Writing – review & editing. Muniappan Ayyanar: Data curation, Formal analysis, Visualization, Writing – review & editing. Satyendra K. Prasad: Validation, Visualization, Data curation, Formal analysis. Manish M. Wanjari: Conceptualization, Project administration, Supervision, Writing – review & editing. Rajesh B. Patil: Writing – original draft, Writing – review & editing. Shailendra S. Gurav: Conceptualization, Project administration, Supervision, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests.