In silico screening of hundred phytocompounds of ten medicinal plants as potential inhibitors of nucleocapsid phosphoprotein of COVID-19: an approach to prevent virus assembly.

Currently, there is no specific treatment to cure COVID-19. Many medicinal plants have antiviral, antioxidant, antibacterial, antifungal, anticancer, wound healing etc. Therefore, the aim of the current study was to screen for potent inhibitors of N-terminal domain (NTD) of nucleocapsid phosphoprotein of SARS-CoV-2. The structure of NTD of RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2 was retrieved from the Protein Data Bank (PDB 6VYO) and the structures of 100 different phytocompounds were retrieved from Pubchem. The receptor protein and ligands were prepared using Schrodinger's Protein Preparation Wizard. Molecular docking was done by using the Schrodinger's maestro 12.0 software. Drug likeness and toxicity of active phytocompounds was predicted by using Swiss adme, admetSAR and protox II online servers. Molecular dynamic simulation of the best three protein- ligand complexes (alizarin, aloe-emodin and anthrarufin) was performed to study the interaction stability. We have identified three potential active sites (named as A, B, C) on receptor protein for efficient binding of the phytocompounds. We found that, among 100 phytocompounds, emodin, aloe-emodin, anthrarufin, alizarine, and dantron of Rheum emodi showed good binding affinity at all the three active sites of RNA binding domain of nucleocapsid phosphoprotein of COVID-19.The binding energies of emodin, aloe-emodin, anthrarufin, alizarine, and dantron were -8.299, -8.508, -8.456, -8.441, and -8.322 Kcal mol-1 respectively (site A), -7.714, -6.433, -6.354, -6.598, and -6.99 Kcal mol-1 respectively (site B), and -8.299, 8.508, 8.538, 8.841, and 8.322 Kcal mol-1 respectively (site C). All the active phytocompounds follows the drug likeness properties, non-carcinogenic, and non-toxic. Theses phytocompounds (alone or in combination) could be developed into effective therapy against COVID-19. From MD simulation data, we found that all three complexes of 6VYO with alizarin, aloe-emodin and anthrarufin were stable up to 50 ns. These phytocompounds can be tested further for in vitro or in vivo and used as a potential drug to cure SARS-CoV-2 infection. Communicated by Ramaswamy H. Sarma.

Introduction

During the end of year 2019, a venomous corona disease was emerged from the city of Wuhan, China (Baloch et al., 2020). In a short span of time, corona has spread all across the globe. This newly emerged viral disease was named as coronavirus disease (COVID-19 and also known as SARS-CoV-2) and World Health Organization (WHO) declared COVID-19 as global health emergency. The COVID-19 is transmitted through person to person and by so far has affected 213 countries across the World. Like other coronaviruses, SARS-CoV-2 has large genome (positive-sense, single-stranded RNA) size of around 30 kb, subdivided into 14 open reading frames (ORFs) (Gordon et al., 2020). The 5′-ORF 1a/ORF 1ab encodes for a polyprotein, which is further processed into 16 non-structural proteins (nsp1-16) by proteolytic cleavage (Gordon et al., 2020). The 3′- end of the RNA genome encodes for 13 ORFs for four structural proteins, such as the matrix (M), small envelope (E), spike (S) and nucleocapsid phosphoprotein (N), and facilitate entry of SARS-CoV-2entry to human cells through human angiotensin converting enzyme 2 (ACE2) receptor (Gordon et al., 2020). Both structural and non-structural proteins ensure virus entry into the cells and its replication inside the host cell and therefore serve as potential as drug targets. The current antiviral drugs developed to treat coronavirus (CoV) infections primarily target S protein, the 3 C-like (3CL), and papain-like (PLP) proteases (Ramajayam et al., 2011; Wrapp et al., 2020). The mutations in the S-protein could lead to escape of viral target by the drugs (Wrapp et al., 2020) whereas antiviral proteases could act on off targets and leads to potential side effects. Therefore, it is important to identify novel antiviral targets that could inhibit assembly of the virus particles. Viral nucleocapsid protein (N) is a potential antiviral drug target, serving multiple critical functions during the viral life cycle, including packaging the viral RNA genome (Baric et al., 1988; Cong et al., 2020; McBride et al., 2014; Nelson et al., 2000). The N-protein binds viral RNA genome and results in the formation of ribonucleoprotein (RNP) complex (Masters & Sturman, 1990; Nelson et al., 2000). The coronavirus N protein has distinct N-terminal RNA-binding domain (NTD), a C-terminal dimerization domain (CTD), and an intrinsically disordered central Ser/Arg (SR)-rich linker (Chen et al., 2013; Lo et al., 2013). It is been proposed that NTD is responsible for RNA binding protein and very crucial for transcription and replication of viral RNA (Grossoehme et al., 2009; Saikatendu et al., 2007). Many critical amino acid residues have been identified for RNA binding and virus infectivity in the N-terminal domain (NTD) of coronavirus N proteins (Tan et al., 2006). Therefore, NTD of nucleocapsid (N) protein RNA binding domains could be one of the essential drug targets to prevent assembly of virus particles (Kang et al., 2020).

Currently, there are really no approved drugs for the specific treatment of COVID-19. The COVID-19 is spreading at a much higher rate and the speed of drug discovery is quite slow and lengthy process. In order to speed up the process of drug development against COVID-19, there is an urgent need to take parallel and multidirectional approach to counter the spread of the virus. The approach to develop a vaccine from the scratch is a time taking process and not very good option for stopping an ongoing pandemic. The alternative approach is repurposing existing drugs such as remedesivir and chloroquine against COVID-19 could be a big step forward to curtail the growing threat (Sanders et al., 2020), but not yet approved by US-FDA for the treatment of COVID-19. This approach to new drug discovery against COVID-19 could be further boosted by in silico (hypothesis-based drug designing) approach and repurposing more existing FDA approved drugs (Ganjhu et al., 2015; Gupta, 2020). Moreover, phytocompounds of many medicinal plants have been traditionally used to cure many viral and respiratory diseases (Ben-Shabat et al., 2020; Ganjhu et al., 2015). Schwarz et al. (2011), specifically reported that emodin (major phytocompound of Rheum emodi) act as inhibitor of 3a ion channel of corona virus SARS-CoV and HCoV-OC43 as well as virus release from HCoV-OC4. Ho et al. (2007) identified emodin as an effective blocker of interaction of the SARS-CoV S protein with the ACE2 and the infection by S protein-pseudo-typed retrovirus. It was shown that the major phytocompounds of Thymus serpyllum, Berberis aristate, Moringa oleifera, Zanthoxylum armatum, Piper nigrum, Myristica fragrans, Thalictrum foliolosum, Cymbopogon citrates and Allium sativum are responsible for the antimicrobial, antiviral, antioxidant, anti-inflammatory, anti-HIV anti-carcinogenic properties (Aja et al., 2014; Fenwick & Hanley, 1985; Ganjhu et al., 2015; Jaric et al., 2015; Meng et al., 1993; Nai-Lan et al., 1993; Potdar et al., 2012; Sharma et al., 2020; Singh & Shikha, 2017; Takooree et al., 2019; Tsai et al., 1985; Oladeji et al., 2019; Weber et al., 1992).

In silico screening of such phytocompounds against specific targets of COVID-19 could provide a continuous feed for in vitro testing and clinical trials. Recently, Kang et al. (2020), solved the crystal structure of SARS-CoV-2 nucleocapsid protein RNA binding domain and proposed as unique target for developing therapeutics against COVID-19. Therefore, the current study was undertaken to test in silico inhibitory potential of 100 major phytocompounds of ten different medicinal plants including major phytocompounds of R. emodi against NTD of nucleocapsid protein RNA binding protein. The phytocompounds identified through in silico studies could become a feeder line for in vitro testing and clinical trials and therefore speed up the process of drug discovery against COVID-19 Figure 1.

Figure 1.

Proposed Binding mechanism of phytocompounds with NTD of nucleocapsid protein RNA binding protein of COVID-19 (Prepared using BioRender).

Methods

Structure reliability of NTD of nucleocapsid phosphoprotein of SARS coronavirus 2

The structure of N-terminal domain (NTD, involved in RNA binding) of nucleocapsid phosphoproteinof SARS coronavirus 2 (N protein PDB ID: 6VYO; MMDB ID: 186265) was retrieved (http://www.rcsb.org/structure/6VYO; https://www.ncbi.nlm.nih.gov (Figure 3A) it based on the previous study by Yadav et al. (2020). The reliability of RNA binding domain of nucleocapsid phosphoprotein was done using he PROCHECK server, which displays Ramachandran plot with distinctive properties of protein such as most favored region, additional allowed regions, generously allowed regions, Disallowed regions, The End residues, Proline, Glycine and total number of residues. Multiple sequence alignment of SARS-N-terminal domain (NTD) of nucleocapsid phosphoprotein (6VYO) was done to find out the conserved regions.

Figure 3.

(A) Secondary structure of NTD (amino acid residues 50-173) of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO); (B) Topology of secondary structure; and C) The Ramachandran plot describes the favored and disfavored regions of residues. The maximum (92.5%) residues are located in an allowed region.

Preparation of NTD of nucleocapsid phosphoprotein of SARS coronavirus 2 for docking

The three dimensional structure of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) (global quality) was prepared using Schrodinger’s Protein Preparation Wizard, by rectifying bond orders, addition of hydrogen bonds, making zero-order bonds to metal atoms, make disulphide bonds, making selenomethionine to methionine conversion, fill in lost side chains prime, erase waters from het groups, and generate het states using EpikPh 7 to +-2 (Georgiev et al., 2005). Protein hydrogen bonds were optimized to remodel by the overlying hydrogen’s and minimized using OPLS3e force field (Lanka et al., 2019).

Receptor grid generation and active sites prediction of NTD of nucleocapsid phosphoprotein of SARS coronavirus 2 for docking

Identification of the potential ligand binding sites or active site residues involved in interaction is a foremost step in docking to discover novel chemical substances (Kalyaanamoorthy & Chen, 2011). The active site residues provide a pocket containing a variety of hydrogen donors, acceptors, hydrophilic regions and hydrophobic regions (Lionta et al., 2014). The binding cavity in NTD of nucleocapsid protein was identified using computational tools. Grid was generated using Glide module in Maestro 12.0 Schrodinger suite. The atoms of protein were settledinside using the default parameters of the radii of Vander Waal’s scaling factor of 1 Å with partial charge cut-off of 0.25 Å using OPLS3e force field. Grid of the protein was generated by Yadav et al. (2020).

Preparation of phytocompounds (ligands)

Hundred major phytocompounds of ten medicinal plants (Rheum emodi, Thymus serpyllum, Cymbopogon citrates, Moringa oleifera, Thalictrum foliolosum, Berberis aristata, Piper nigrum,Allium sativum, Myristica fragrans and Zanthoxylum armatum) were utilized for virtual screening and molecular docking study against NTD of nucleocapsid (N-protein) of COVID-19 target proteins. All the hundred phytocompounds were retrieved from the PubChem Database (https://pubchem.ncbi.nlm.nih.gov/). Ligands were prepared by using the ligand preparation module of molecular modeling package with suitable parameters like optimization, ring conformation, 2 D to 3 D conversion, and determination of protomers, tautomers and ionization states at pH 7.0; along with partial atomic charges using OPLS3e force field with root mean square deviation (RMSD) value of 0.11 Å.

Molecular docking

Molecular docking study of selected 100 phytocompounds of 10 medicinal plants for High throughput virtual screening was performed against the NTD of nucleocapsid N-protein (PDB ID: 6YVO). HTVS with flexible docking was performed on Schrodinger’s maestro 12.0 between Ligands and active sites of target protein (6YVO). All the hundred ligands were filtered out by Standard Precision docking (SP) using ligand-docking process. Among hundred Ligands, the best configuration with highest docking score of virtual screening were sorted by Extra Precision (XP) with Flexible docking.

Binding free energies calculation using prime MM/GBSA approach

The binding free energies of ligand and receptor complex were analyzed by prime MM/GBSA (Molecular Mechanics Generalized Born Surface Area) module of Schrodinger suite with the OPLS3e force field. The prime MM GBSA−1 approach is based on the docking complex and is used to calculate ΔGbind of each ligand using the below mentioned equationand as described earlier (Ramatenki et al., 2015; Tripathi et al., 2013).

Where, ΔEMM is the difference in the minimized energiesbetween the 6VYO-ligand complex and the sum of the energies of the unliganded6VYO and inhibitor. ΔGSolv isthe difference in the GBSA solvation energy of 6VYO inhibitor complex and sum of the solvation energies for the free 6VYO receptor and inhibitor complex and ΔGSA is the difference in surface area energies for the complex and sum of the surface area energies for the unliganded 6VYO receptor and inhibitor. MM GBSA−1 approach is used as a rescoring function to prioritize the lead inhibitors. The binding free energies obtained from MM GBSA1 calculations were considered along with docking score to optimize the identified drug molecules.

Admet and toxicity prediction of active phytocompounds

Absorption, distribution, metabolism, excretion, and toxicity (ADMET) screening was done to determine the absorption, toxicity, and drug-likeness properties of ligands.The 3 D structures of phytocompounds (aloe-emodin, anthrarufin, alizarine, dantron and emodin) were saved in .smiles format and drug were uploaded on SWISSADME (Molecular Modeling Group of the SIB (Swiss Institute of Bioinformatics), Lausanne, Switzerland), admetSAR (Laboratory of Molecular Modeling and Design, Shanghai, China), and PROTOX webservers (Charite University of Medicine, Institute for Physiology, Structural Bioinformatics Group, Berlin, Germany) for ADMET screening. SwissADME is a web tool used for the prediction of ADME and pharmacokinetic properties of a molecule. The predicted result consists of lipophilicity, water solubility, physicochemical properties, pharmacokinetics, drug-likeness, medicinal chemistry, and Brain or Intestinal Estimated permeation method (blood–brain barrier and PGP ± prediction). admetSAR provides ADMET profiles for query molecules and can predict about fifty ADMET properties. Toxicity classes are as follows: (i) Category I contains compounds with LD50 values ≤ 50 mg kg−1, (ii) Category II contains compounds with LD50 values > 50 mg kg−1 but 500 mg kg−1 but 5000 mg kg−1. PROTOX is a Rodent oral toxicity server predicting LD50 value and toxicity class of query molecule. The toxicity classes are as follows: (Class 1: fatal if swallowed (LD50 ≤ 5), Class 2: fatal if swallowed (55000) (Banerjee et al., 2018).

MD Simulation

Molecular dynamic simulation was performed to analyze the stability of ligand-protein interactions with respect to the physical transition of the structural aspect of macromolecules to the functional relevance of the complex. In brief, MD simulation demonstrates strength, pattern, dynamic conformational changes and intermolecular properties of the interactions. Among all the complexes, best three complexes with binding site A were selected for MD simulation up 50 ns by using AMBER 18 software. The root means square deviation (RMSD), root mean square fluctuation (RMSF), intermolecular hydrogen bonds interactions, binding free energy and radius of gyration (RoG) were evaluated for the elucidation of conformational, structural and compactness of the protein-ligand complexes (Nandekar et al., 2013). Molecular receptor topology was created by leap module of software.

Subsequently, following established protocol such as neutralization and submerging complex into the water molecule rectangular box (TIP3P) where the protein atoms were 10 Å away from the nearest edge of the box. The minimization of the solvent system was achieved by freezing the protein and removing the bad contact with restraint on the heavy atoms, first through 2500 steepest descent method and then conjugation gradient for 2000 further steps. The system was gradually heated from 0 to 300 K temperature at 1 atm pressure (NPT conditions). It was achieved for 50 ps followed by 50ps of density equilibrium and 1 ns of constant equilibrium to exchange the potential and kinetic energies. Afterward, the temperature was kept constant by using the Berendsen algorithm. The MDS was run for 50 ns to evaluate the stability of the docked ligand-receptor complexes preceded by the system equilibration for 500ps on the canonical NPT ensemble. The covalently bound hydrogen atoms were constrained by the SHAKE algorithm with 2 fs time and temperature control by Langevin dynamics. Finally, the recording was made for every 5ps using the particle-mesh Ewald summation method for treating electrostatic interaction. The RMSD, RMSF and RoG values were recorded for 50 ns.

Results

Protein structure reliability by ramachandran plot

The scheme for the current study has been summarized in Figure 2. The secondary structure of NTD (amino acid Ala 50 to Ala 173), the RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2 (6VYO) is shown in Figure 3A and secondary structure topology is shown in Figure 3B. The PROCHECR server presents Ramachandran plot displaying allowed and disallowed regions with regard to back bone dihedrals of protein residues. More than 90% score of the most favored regions ensures that model has good quality. RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2 (6VYO) protein showed 92.5% of amino acid residues in most allowed regions (Table 1 and Figure 3C). Multiple sequence alignment of SARS-N-terminal domain (NTD) of nucleocapsid phosphoprotein (6VYO) was done by blast to find out the conserved region. The domain hits showed similarity with Coronavirus nucleocapsid protein with E-value of 3.14e-70 and with N-terminal domain of nucleocapsid (N) protein of coronavirus; The coronavirus nucleocapsid (N) with E-value of 1.23e-54 (Figure 4A and B). The phylogenetic tree of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 is shown in and Multiple sequence alignment of nucleocapsid phosphoprotein from SARS coronavirus 2 in .

Figure 2.

Schematic of the experiments followed in the current study.

Table 1.

Ramachandran Plot statistics showing favored regions in RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) (Yadav et al., 2020).

Properties	RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2(6YVO)
	Residues	%
The Most favored regions	358	92.5
Additional allowed regions	29	7.5
Generously allowed regions	0	0
Disallowed regions	0	0
The End residues (excluding Glycine and Proline)	8	–
Proline	64	–
Glycine	40	–
Total number of residues	499	100

Figure 4.

(A) Conserved domains of (6 VYO_1) chains A, B, C, D nucleoprotein of Severe acute respiratory syndrome coronavirus 2. (B) List of domain hits with E-values.

Preparation of NTD, the RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2 (6VYO)

The three dimensional homo-tetrameric structure of RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2 (PDB ID: 6VYO) was retrieved from RCSB PDB (http://www.rscb.org/pdb) in .sdf format. The protein structure was prepared by correcting bond orders, addition of hydrogen bonds, creating zero-order bonds to metal atoms, creating disulphide bonds, creating selenomethionine to methionine conversion, fill in the missing side chains prime, deleting waters from het groups, and generate het states using EpikPh 7 to +-2 by using Schrodinger’s Protein Preparation Wizard (Figure 5A-D).

Figure 5.

Active site prediction. 3-D structure representation of tetrameric RNA binding domain (N-terminal domain) of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) showing three predicted binding sites shown as Active site A, Active site B and Active site C (A). Panel B, C, and D represent tetrameric RNA binding domain showing each active site as indicated.

Receptor grid generation and active site prediction

To predict the active site of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO), protein grid was generated by using glid module in maestro 12.0 schrodinger suit. In RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) protein, we found three different active sites, namely A, B and C with D score 1.039, 0.99 and 1.045 respectively. Active site of target protein is based on previous study Yadav et al. (2020). The dimensions of the grid box and receptor setup were x = 80 Å, y = 80 Å, z = 80 Å and x = 32 Å, y = 32 Å, z = 32 Å respectively with a grid space of 1 Å for all the three binding sites. Active site A was predicted between chain A and B. The chain A contributes 26 amino acid and chain B contributes 21 amino acids. The Active site B was predicted between chain A (25 amino acids) and chain B (26 amino acids). Active site C was predicted between chain C and D with 15and 21 amino acids respectively (Table 2 and Figure 5 A-D). This is important to mention that Thr 76, Ser 78, Ser 105, Thr 166 which are potential phosphorylation sites are involved in formation of active sites A, B and C (Table 2; Yadav et al., 2020).

Table 2.

Amino acids involved in active sites of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) (Yadav et al., 2020).

Active Sites	D Score	Size	Site Score	Chains and amino acid residues numbers
Site-A	1.039	216	1.009	CHAIN A: ALA 50, SER 51, PHE 53, THR 54, ALA 55, LEU 56, THR 57, HIP 59, ARG 92, SER 105, ARG 107, TYR 109, ARG 149, ALA 156, ILE 157, VAL 158, GLN 160, LEU 161, PRO 162, GLN 163,GLY 164, THR 165, LEU 167, PHE 171, TYR 172,ALA 173CHAIN B: TRP 52, PHE 53, PRO 73, ILE 74, ASN 75, THR 76, ASN 77, SER 78, ASP 82, HIE 145, ILE 146, THR148, ARG 149, ASN 150, ASN 153, ASN 154, ALA 155, VAL 158, LEU 159, GLN 160, ILE 157
Site-B	0.99	196	0.98	CHAIN A: TPR 52, PHE 53, PRO 73, ILE 74, ASN 75, THR 76, ASN 77, TYR112, GLY 114, THR 115, GLY 116, ALA 119, ASP 144, HIE 145, ILE 146, GLY 147, THR 148, ARG 149, ASN 150, ASN 153, ASN 154, ILE 157, VAL 158, LEU 159, GLN 160CHAIN D: THR 54, ALA 55, THR 57, HIE 59, ALA 90, THR 91, ARG 92, ILE 94, GLY 96, GLY 97, ASP 98, LYS 102, LEU 104, SER 105, PRO 106, ARG 107, TYR 109, VAL 158, LEU 159, GLN 160, LEU 161, GLN 163, THR 165, THR 166, LEU 167, ALA 173
Site-C	1.045	183	1.039	CHAIN C: ALA 50, SER 51, PHE 53, THR 54, ALA 55, LEU 56, THR 57, HIP 59, ARG 92, ARG 107, TYR 109, ARG 149, ALA 156, ILE 157, VAL 158CHAIN D: TRP 53, PRO 73, ILE 74, ASN 75, THR 76, ASN 77, SER 78, ASP 82, HIE 145, ILE 146, GLY 146, THR 148, ARG 149, ASN 150, ASN 153, ASN 154, ALA155, ILE157, VAL158, LEU159, GLN160

Each amino acid residue is numbered followed by three letter code. Amino acid residues in boldface letters are the phosphorylation sites on Serine and Threonine were predicted by Davidson et al. (2020). Nonstandard amino acids are CL, GOL, HOH, MES, ZN.

Ligand preparation

Hundred major phytocompounds of ten selected medicinal plants (Rheum emodi, Thymus serpyllum, Cymbopogon citrates, Moringa oleifera, Thalictrum foliolosum, Berberis aristata, Piper nigrum, Allium sativum, Myristica fragrans and Zanthoxylum armatum) were used for the docking studies. The 3-dimensional structures of all the phytocompounds were obtained from pubchem (www.pubchem.com) in .sdf format. The .sdf file of phytocompounds was converted into PDB format by using the tool Open Babel. Table 3 summarizes the names of the medicinal plants, pharmacological properties and major phytocompounds. Phytocompounds are numbered as 1-100. Figure 6 showed five phytocompounds of R. emodi, which showed very efficient binding with NTD of nucleocapsid phosphoprotein of SARS coronavirus 2.

Table 3.

List of hundred major phytocompounds (numbered 1-100) of 10 medicinal plants with pharmacological properties.

S.No	Medicinal Plants	Pharmacological properties	Major phytocompoundsSelected for the current study
1.	Rheum emodi	Treatment of Severe Acute Respiratory Syndrome (SARS) (Ho et al., 2007; Schwarz et al., 2011). Antioxidant, antibacterial antifungal and Bioenhancer of antibiotics (Dev et al., 2017; Rolta et al., 2018a; 2020).	(1) Alizarin, (2) Aloe-emodin, (3) Anthraquinone, (4) Anthrarufin, (5) Anthrone, (6) Chryophanol, (7) Clicoemodin, (8) Dantron, (9) Emodin, (10) Gallic acid, (11) Juglone, (12) Physicon, (13) Piceatannol, (14) Quinizarin, (15) Rubiadin
2.	Thymus serpyllum	Antimicrobial, antiviral, antioxidant, anti-inflammatory and anti-carcinogenic (Jaric et al., 2015)	(16) 3-Octanone, (17) 8-Prenylnaringenin, (18) Apigenin, (19) Apigenin-7-O-glucoside, (20) Bornyl acetate, (21) Camphene, (22) Camphor, (23) Carvacrol, (24) Caryophyllene oxide, (25) Catechin, (26) Eriodictyol, (27) Eriodictyol-7-glucuronide, (28) Geranyl acetate, (29) Linalool, (30) Myrcene, (31) Rutin, (32) Spathulenol, (33) Taxifolin, (34) Thymol, (35) trans-Nerolidol, (36) α-Pinene, (37) β-Bisabolene, (38) β-Caryophyllene, (39) β-Pinene, (40) δ-Cadinene.
3.	Cymbopogon citrates	Anti-obesity, anti-bacterial, anti-fungal, anti-nociceptive, anti-oxidants anti-diarrheal, and anti-inflammatory (Singh & Shikha, 2017).	(41) Citral alpha, (42) citronellal, (43) Geraniol, (44) lemonene, (45) nerol, (46) terpinolene
4.	Moringa oleifera	Antimicrobial (Pal et al., 1995) antiviral potential (Aja et al., 2014).	(47) Glucomoringin, (48) Niazinicin A, (49) Niazinin, (50) Pterygospermin
5.	Thalictrum foliolosum	Antiplasmodial, antiviral, antibacterial and anti-gastric ulcer activity (Sharma et al., 2020)	(51) 8-Oxyberberine, (52) Berberine, (53) Jatrorrhizine, (54) Noroxyhydrastinine, (55) Palmatine, (56) Thalicarpine, (57) Thalidasine, (58) Thalirugidine, (59) Thalirugine, (60) Thalisopine, (61) Thalrugosaminine, (62) Thalrugosidine
6.	Berberis aristata	Anti-diabetic, anti-microbial, anti-cancer, antilipidemic, anti-HIV, anti-pyretic (Potdar et al., 2012)	(63) Berbamine, (64) Columbamine, (65) Isocorydine, (66) Lambertine, (67) Lupeol, (68) Oleanolic acid, (69) Oxyberberine, (70) Oxycanthine, (71) Stigmasterol
7.	Piper nigrum	Antimicrobial, antioxidant, anticancer, analgesic, anticonvulsant, neuroprotective, hypoglycemic, hypolipidemic, and anti-inflammatory activities (Takooree et al., 2019)	(72) Sabinene, (73) Piperolein B, (74) Piperolein A, (75) Piperine, (76) Piperazine, (77) Guaiol, (78) b-caryophyllene
8.	Allium sativum	Anti-influenza A and B (Fenwick & Hanley, 1985), antiviral (Meng et al., 1993; Nai-Lan et al., 1993; Tsai et al., 1985)	(79) Vinyllithium, (80) Alliin, (81) Allicin, (82) (Z)-Ajoene, (83) (E)-Ajoene
9.	Myristica fragrans	Antioxidant, antimicrobial, aphrodisiac, anticancer, hepatoprotective, anti-inflammatory, antidepressant, and cardioprotective activity (Ganjhu et al., 2015)	(84) Trimyristin, (85) safrole, (86) Sabinene, (87) myristicin, (88) methoxyeugenol, (89) meso-dihydroguaiaretic acid, (90) isoelemicin, (91) elemicin, (92) camphene, (93) 9-epoxylignan, (94) 4-terpinol
10.	Zanthoxylum armatum	Antifungal, antibacterial, antiviral, anti-inflammatory, anti-oxidant, cardiovascular and immunomodulatory (Singh & Shikha, 2017)	(95) β-Sitosterol-β-D-Glucoside, (96) Xanthyletin, (97) Tambulin, (98) Sesamin, (99) Bergapten, (100) Armatamide

Figure 6.

Structures of potential phytocompounds of R. emodi: A) Anthrarufin, (B) Aloe- emodin, (C) Alizarine, (D) Dantron, (E) Emodin.

Molecular docking analysis of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) target protein and phytocompounds

Molecular docking of all the 100 phytocompounds with three binding sites (A, B and C) of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) protein was done by using Schrodinger's Software. Active phytocompounds were analyzed by scoring function and binding free energies MM GBSA−1. Out of the hundred phytocompounds, only five phytocompounds showed the best interactions with all three-binding site of RNA binding domain of nucleocapsid phosphoprotein from SARS coronavirus 2 (PDB ID: 6VYO) protein (Tables 4–6 and Figures 7–12). Molecular structures of five best phytocompounds are shown in Figure 6. The remaining 95 phytocompounds did not show efficient binding (data not shown) and hence were not pursued further.

Molecular docking of phytocompounds with active site A

Active site A was predicted between chain A and Chain B with D score 1.039. Among all the selected phytocompounds, only aloe-emodin, anthrarufin, alizarine, dantron and emodin showed the best interaction with docking score −8.508, −8.456, −8.441, −8.322 and −8.299 Kcal mol−1respectively. Prime MM GBSA−1 (Kcal mol−1), hydrogen bonding and interactive amino acids are shown in Table 4, Figures 7 and 8. Apart from hydrophobic interactions, H-bonds were also detected. Ser 78 of chain B, which has been predicted as phosphorylation site showed interaction with all the five phytocompounds.

TABLE 4.

Interacting amino acids and docking score of phytocompounds with binding site A (1.039) of 6VY0.

Compound CID no.	Name of Ligand (Plant source)	ΔG Bind (kcal mol-1)	Docking Score	H-bond	Interactive Residues
10207	Aloe- emodin (Rheum emodi)	−25.45	−8.508	Chain A: ARG 107Chain B: HIE 145, ASN 75, ASN 77	Chain A: LEU 159, VAL 158, THR 54, ALA 55, THR 57, HIP 59, ARG 92, ARG 107Chain B:SER 78, ILE 157, ASN 154, TRP 52, ILE 146, HIE 145, ASN 75, ASN 77
8328	Anthrarufin (Rheum emodi)	−38.099	−8.456	Chain B: ASN 75	Chain A: THR 54, ALA 55, THR 57, HIP 59, ARG 92, ARG 107, VAL 152Chain B: TRP 52, ASN 154, ILE 157, ASN 77, SER 78, ILE 146, HIE 145, ASP 82, ASN 75
6293	Alizarine (Rheum emodi)	−31.2891	−8.441	Chain A: ASN 75Chain B: ASN 77, ASN 75, HIE 145	Chain A: ARG 107, HIP 59, ARG 92, ALA 55, THR 54, LEU 159, VAL 158, ASN 75Chain B: TRP 52, TYR 112, ILE 157, ASN 154, SER 78, ILE 146, HIE 145, ASN 77, ASN 75
2950	Dantron (Rheum emodi)	−39. 8129	−8.322	Chain A: THR 57, ALA 55Chain B: ASN 75, ASN 77	Chain A: HIP 59, LEU 56, THR 54, ARG 107, LEU 159, VAL 158, THR 57, ALA 55Chain B: ASP 82, HIE 145, ILE 146, SER 78, TRP 52, ILE 157, ASN 75, ASN 77
3220	Emodin (Rheum emodi)	−43.87	−8.299	Chain A: ALA 55Chain B: ASN 77, ASN 75	Chain A: VAL 158, LEU 159, THR 54, LEU 56, HIP 59, ARG 107, ALA 55Chain B: ASP 82, ILE 157, TRP 52, SER 78, HIE 145, ILE 146, ASP 82, ASN 77, ASN 75

Each amino acid residue is numbered followed by three letter code. Amino acid residues in boldface letters are the phosphorylation sites on Serine and Threonine were predicted by Davidson et al. (2020).

Figure 7.

The 2D molecular interactions of best 5 phytocompounds with active site A of 6VYO: (A) Alizarin, (B) Aloe- emodin, (C) Anthrarufin (D), Dantron, (E) Emodin.

Figure 8.

Surface view of 6VYO (active site A) complex with best 5 compounds of selected medicinal plants: (A) Alizarin, (B) Aloe- emodin, (C) Anthrarufin (D), Dantron, (E) Emodin. Phytocompounds are colored in green.

Molecular docking of phytocompounds with active site B

Active site B was predicted between chain A and chain With D score 0.99. Among all the selected phytocompounds, only emodin, dantron, alizarine, aloe-emodin and anthrarufin showed the best interaction with docking score −7.714, −6.99, −6.598, −6.433 and −6.354 Kcal mol−1respectively. Prime MM GBSA−1 (Kcal mol−1), hydrogen bonding and interactive amino acid is shown in Table 5, Figures 9 and 10. Both hydrophobic and H-bond interactions were observed. Ser 78 of Chain A and Ser 105 of chain B, the potential phosphorylation residues were also involved in hydrophobic interaction. Exceptionally, SER 78 of chain A also showed H-bond formation with Aloe-emodin.

Table 5.

Interacting amino acids and docking score of phytocompounds with binding site B (0.99).

Compound CID no.	Nameof Ligand (Plant source)	ΔG Bind (kcal mol-1)	Docking Score	H-bond	Interactive Residues
3220	Emodin (Rheum emodi)	−40.60	−7.714	Chain A: HIE 145, ASN 77, ASN 75, ASN 154Chain D: ARG 92	Chain A: ILE 146, SER 78, TRP 52, ILE 157, HIE 145, ASN 77, ASN 75, ASN 154Chain D: HIE 59, ALA 173, THR 57, LEU 104, SER 105, ARG 107, ARG 92
2950	Dantron (Rheum emodi)	−45.48	−6.99	Chain A: ASN 154, ILE 146, ASN 77	Chain A: TRP 52, ASN 75, ULE 157, HIE 145, ASN 75, SER 78, ASN 154, ILE 146, ASN 77Chain D: ARG 107, SER 105, LEU 104, ARG 92
6293	Alizarine (Rheum emodi)	−33.59	−6.598	Chain A: ASN 77, ASN 75	Chain A: ILE 146, HIE 145, ARG 92, SER 78, TRP 52, ILE 157, ASN 154, ASN 77, ASN 75Chain D:SER 105, ARG 107, HIE 59, THR 57
10207	Aloe- emodin (Rheum emodi)	−36.92	−6.433	Chain A: ASN 154, ILE 157, SER 78, ILE 146, HIE 145, TRP 52, ASN 75, ASN 77	Chain A: ASN 154, ILE 157, SER 78, ILE 146, HIE 145, TRP 52, ASN 75, ASN 77Chain D: ARG 107, SER 105, LEU 104, ARG 92, HIE 59, THR 57
8328	Anthrarufin (Rheum emodi)	−34.9567	−6.354	Chain A: ASN 77	Chain A:SER 78, ASN 75, TRP 52, ILE 146, ILE 157, ASN 77Chain D: ARG 107, SER 105, ARG 92, THR 57, HIE 59

Each amino acid residue is numbered followed by three letter code. Amino acid residues in boldface letters are the phosphorylation sites on Serine and Threonine predicted by Davidson et al. (2020).

Figure 9.

The 2 D molecular interactions of best 5 phytocompounds with active site B of 6VYO: (A) Alizarin, (B) Aloe- emodin, (C) Anthrarufin (D) Dantron (E) Emodin.

Figure 10.

Surface of 6VYO (active site B) with best 5 compounds: (A) Alizarin, (B) Aloe- emodin, (C) Anthrarufin (D) Dantron (E) Emodi. Phytocompounds are colored in green.

Molecular docking of phytocompounds with active site C

Active site C was predicted between chain A and chain d with D score 1.045. Among all the selected phytocompounds, anthrarufin, aloe-emodin, alizarine, dantron and emodin showed the best interaction with docking score −8.538, −8.508, −8.841, −8.322 and −8.299 Kcal mol−1respectively. Prime MM GBSA−1 (Kcal mol−1), hydrogen bonding and interactive amino acid is shown in Table 6, Figures 11 and 12. Hydrophobic interactions and H-bonds were detected all the phytocompounds. Ser 78 of chain D, the potential phosphorylation residue was also involved in hydrophobic interaction.

Table 6.

Interacting amino acids and docking score of phytocompounds with binding site-C (1.045).

Compound CID no.	Name of Ligand (Plant source)	ΔG Bind (kcal mol-1)	Docking Score	No. of H-bond	Interactive Residues
8328	Anthrarufin (Rheum emodi)	−39.48	−8.538	Chain D: ASN 75, ASN 77	Chain C: TTR 54, ALA 55, THR 57, HIP 59, ARG 107, ARG 92Chain D: TRP 52, ASN 154, ILE 157, HIS 145, ILE 146, SER 78, ASP 82, ASN 75, ASN 77
10207	Aloe- emodin (Rheum emodi)	−32.61	−8.508	Chain C: THR 57, LEU 159Chain D: ASN 75	Chain C: HIP 59, LEU 56, ALA 55, THR 54, VAL 158, LEU 159, THR 57, LEU 159Chain D: ASN 77, SER 78, ASN 154, ILE 157, VAL 158, LEU 159, ILE 146, HIS 145, ARG 107, ASN 75, ASN 75
6293	Alizarine (Rheum emodi)	−33.84	−8.441	Chain C: ALA 55, LEU 56, THR 57Chain D: ASN 75	Chain C: THR 54, LEU 159, VAL 158, ALA 55, LEU 56, THR 57Chain D: HIS 145, ILE 146, ASN 77, SER 78, ASN 154, ILE 157, TRP 52, TYR 112, ASN 75
2950	Dantron (Rheum emodi)	−35.459	−8.322	Chain C: THR 57, ALA 55	Chain C: HIP 59, THR 57, LEU 56, ALA 55, THR 54, ARG 107, VAL 158, LEU 159Chain D: HIS 145, ILE 146, TRP 52, SER 78, ASN 77, ASN 154, ILE 157, ASP 82
3220	Emodin (Rheum emodi)	−32.78	−8.299	Chain C: ARG 107Chain D: ILE 146, ASN 175	Chain C: TYR 109, ARG 107, ARG 92, HIP 59, THR 57, ALA 55, THR 57Chain D: ASN 153, ASN 154, ILE 146, HIS 145, TRP 52, ASP 82, SER 79, SER 77, ASN 75

Each amino acid residue is numbered followed by three letter code. Amino acid residues in boldface letters are the phosphorylation sites on Serine and Threonine predicted by Davidson et al. (2020).

Figure 11.

The 2 D molecular interactions of best 5 phytocompounds with active site C of 6VYO: (A) Alizarin, (B) Aloe- emodin, (C) Anthrarufin (D) Dantron (E) Emodin. Phytocompounds are colored in green.

Figure 12.

Surface of 6VYO (active site C) and 5 best compounds complex: (A) Alizarin, (B) Aloe-emodin, (C) Anthrarufin (D) Dantron (E) Emodin. Phytocompounds are colored in green.

Admet prediction and toxicity analysis of active phytocompounds

The ADME properties predicted by SwissADME of selected phytocompounds are summarized in Table 7. Consensus Log Po/w value of ˂5 indicates good aqueous solubility, which means that an adequate amount of drug can reach the target and can be maintained inside the body through oral administration. All the five phytocompounds have consensus Log Po/w value of <5 (Table 7). TPSA indicates a compound ability to permeate into cells. A TPSA value of <140 Å2 is required for good permeation of compound into the cell membrane and value <90 Å2 is required to permeate through blood–brain barrier. All the selected phytocompounds have TPSA value <90 Å2, indicating good permeability of selected phytocompounds through blood–brain barrier. Lipinski’s rule of five helps to determine drug‑likeness of the phytocompounds; an orally active drug should not violate the rule. Both ligands showed 0 violations for Lipinski’s rule of five.

Table 7.

ADME properties of selected phytocompounds predicted by SwissADME.

Compounds	SwissADME
	Water Solubility	GI Absorption	Drug likeness Rules
			Lipinski’s rule of 5					Ghose	Veber	Egan	Muegge
			Logp	nON	nOHNH	TPSA (Å2)	MW
Emodin	Soluble	High	3.01	5	3	94.83	270.24	Yes	Yes	Yes	Yes
Anthrarufin	Moderately soluble	High	3.13	4	2	74.60	240.21	Yes	Yes	Yes	Yes
Alizarine	Soluble	High	2.90	4	2	74.60	240.21	Yes	Yes	Yes	Yes
Aloe- emodin	Soluble	High	2.42	5	3	94.83	270.24	Yes	Yes	Yes	Yes
Dantron	Soluble	High	3.13	4	2	74.60	240.21	Yes	Yes	Yes	Yes

Toxicity predicted by PROTOX‑II and admetSAR is summarized in Table 8. Data in Table shows that all the selected phytocompounds are non-carcinogenic and non-cytotoxic in nature and are safe to administer.   However, LD50 value of emodin and artemisinin calculated from Protox II was higher than that of all other phytocompounds and chloroquine, indicating that these natural phytocompounds are safer than that of chemically synthesized chloroquine.

Table 8.

Toxicity prediction of active phytocompounds predicted by admetSAR and PROTOX-II software.

Phytocompounds	admet SAR		Protox II
	Carcinogens	Rate Acute toxicity (LD50) kg mol-1	LD50, (mg kg-1)	Cytotoxicity
Emodin	non-carcinogen	2.5826 (III)	5000 (Class 5)	Inactive
Anthrarufin	Non-carcinogen	2.8601 (III)	5000 (class 5)	Inactive
Alizarine	Non-carcinogen	2.5292 (III)	7000 (class 5)	Inactive
Aloe- emodin	Non-carcinogen	2.9280 (II)	5000 (class 5)	Inactive
Dantron	Non-carcinogen	2.8601 (III)	7000 (class 6)	Inactive

MD Simulations of protein-ligand complexes

Among all the 100 phytocompounds, Alizarine, Aloe- emodin and Anthrarufin showed the good binding affinity with binding site NTD of nucleocapsid protein RNA binding protein of covid-19 (6VYO) and these phytocompounds follows the all parameters (RMSD, RMSF, Intermolecular hydrogen bonding, RoG and Binding free energy through MDS) for molecular dynamic simulation study.

RMSD analysis:

RMSD analysis showed structural stability of protein-ligand complex with in particular time period. Protein-ligand complex behavior upto 50 ns by MDS. The RMSD values of the 6VYO with Alizarine is gradually stabilized at 3 Å and 30 ns and values are same for 50 ns. Similarly, in case of aloe-emodin RMSD value is stabilized at 4 Å and 35 ns and values are same for 50 ns. In case of Anthrarufin RMSD value is stabilized at 3.5 − 4 Å and 45 ns and values are same for 50 ns (Figure 13A).

Figure 13.

A) RMSD graph of three protein-ligand complexes for 50 ns. (B) RMSF graph of three ligand-protein complexes for 50 ns.

RMSF analysis:

The RMSF determines the deviation of the particle from the original position of macromolecules of three dimensional structure of protein. It describes the conformational flexibility of the protein structure. We found the low fluctuation in loop of Alizarine and 6VYO complex; aloe- emodin and 6vyo complex and high fluctuation anthrarufin and 6vyo complex (Figure 13B).

Intermolecular hydrogen bonding analysis

It measures the binding affinity of ligand. Therefore more no. of hydrogen bonds between the protein and Ligands showed strong binding affinity. In this study we found maximum 27 hydrogen bond at the end of 50 ns simulation 23 hydrogen bonds in Alizarine and 22 in aloe- emodin and 24 in Anthrarufin (Figure 14A).

Figure 14.

(A) Number of Intermolecular hydrogen bonds between the ligands and amino acid residues of the target protein for 50 ns. (B) Radius of gyration (RoG) results of the ligand-protein complex for 50 ns.

RoG analysis:

It describes the folding competence of the protein- ligand complexes. Through RoG analysis we confirm the conformational changes in the three-dimensional structure of protein after interaction of ligand. The high RoG value showed the loose packing and folding behavior of the protein ligand complex. All the complexes showed the high value at approximately 27.2 Å throughout 50 ns (Figure 14B).

Discussion

Molecular docking studies explored the interaction of 100 major phytocompounds of ten selected medicinal plants (Rheum emodi, Thymus serpyllum, Cymbopogon citrates, Moringa oleifera, Thalictrum foliolosum, Berberis aristata, Piper nigrum, Allium sativum, Myristica fragrans and Zanthoxylum armatum) for interaction with RNA binding domain (NTD) of nucleocapsid phosphoprotein of SARS coronavirus 2.

The molecular docking study showed that out to 100 different phytocompounds,five phytocompounds such as emodin, anthrarufin, alizarine, aloe- emodin and dantron of Rheum emodishowed efficientbinding affinity and pharmacokinetic properties with NTD of RNA binding domain of nucleocapsid phosphoprotein of SARS coronavirus 2. It is very interesting to note that all the all the five phytocompounds (emodin, aloe-emodin, anthrarufin, alizarine, and dantron) of R. emodi bind at three different active sites. Thebinding energies of emodin, aloe-emodin, anthrarufin, alizarine, and dantron was −8.299, −8.508, −8.456, −8.441, and −8.322 Kcal mol−1 respectively (site A), −7.714, −6.433, −6.354, −6.598, and −6.99 Kcal mol−1 respectively (site B), and −8.299, 8.508, 8.538, 8.841, and 8.322 Kcal mol−1 respectively (site C). Similar to our study, Kumar et al (Kumar et al., 2020) reported the binding affinities of Nelfinavir, Rhein, Withanolide D, Withaferin A, Enoxacin, and Aloe-emodin as −8.4, −8.1, −7.8, −7.7, −7.4, and −7.4 Kcal mol−1 respectively against main protease (6LU7) of COVID 19. Latha and Pandit (Latha & Pandit, 2020) reported the binding affinity of phytocompounds derived from Silybum marianum (Silybin), Withania somnifera (Withaferin A), Tinospora cordifolia (Cordioside) and Aloe barbadensis (Catechin and Quercetin) with SARS-CoV-2 target. The binding energies of phytocompounds was more than the widely used hydroxychloroquine and other repurposed drugs used for treatment of COVID–19 infection. Tatar and Turhan (2020) reported the binding affinity of nafamostat, rapamycin, saracatinib, Imatinib and Camostat with binding energy −10.24, −9.28, −9.66, −9.23 and 9.07 respectively with N-Terminal Domain of SARS-CoV-2 Nucleocapsid. Ubani et al. (2020)reported best binding affinity of scopodulcic acid and 249 Dammarenolic acid with the Spike glycoprotein (6VSB) and the Mpro 250 (6FLU7) respectively. The in-silico studies were supported by the fact that Schwarz et al. (2011) reported emodin as inhibitor of 3a ion channel of corona virus SARS-CoV and HCoV-OC43 as well as virus release from HCoV-OC43. In another study, Ho et al. (2007) identified emodin as an effective to block the interaction of the SARS-CoV S protein with the ACE2 and the infection by S protein-pseudo-typed retrovirus.

Similar to our study, Ahmed and Shohael (Ahmed & Shohael, 2019), reported that the chrysophanol, aloe-emodin, rhein and emodin fulfill all the criteria of Lipinski’s rule of five and ADME/T. Nucleocapsid N protein of SARS coronavirus 2 has been predicted to contain many phosphorylation sites at Serine, Threonine and Tyrosine residues (Davidson et al., 2020). In our study, we found that some of the potential Serine and Threonine residues contribute for the formation of active site A, B and C. Serine at position 78, 79 and 105 also showed interaction with all the five phytocompounds of R. emodi. Involvement of potential phosphorylation residues in interaction with phytocompounds could further ensure that NTD of Nucleocapsid N protein of SARS coronavirus 2 render inactive and prevents the assembly of the virus and hence stop the infection.

Conclusion

Hundred phytocompounds of ten medicinal plants were screened on the bases of binding energy against NTD of Nucleocapsid protein of COVID-19 target. Among all the 100 phytocompounds, emodin, Anthrarufin, Alizarine, Aloe- emodin, and Dantron showed the good binding affinity at three different active sites of N-Terminal Domain of SARS-CoV-2 Nucleocapsid Protein and potentially prevent the assembly of virus particles and stop the infection, ADMET prediction revealed that emodin, Anthrarufin, Alizarine, Aloe- emodin, and Dantron are potential drugs for COVID-19 treatment. Molecular dynamic simulation of the three protein- ligand complexes (alizarin, aloe-emodin and anthrarufin) showed stable interaction for 50 ns.

Click here for additional data file.