Identification of potential inhibitors against SARS-CoV-2 by targeting proteins responsible for envelope formation and virion assembly using docking based virtual screening, and pharmacokinetics approaches.

WHO has declared the outbreak of COVID-19 as a public health emergency of international concern. The ever-growing new cases have called for an urgent emergency for specific anti-COVID-19 drugs. Three structural proteins (Membrane, Envelope and Nucleocapsid protein) play an essential role in the assembly and formation of the infectious virion particles. Thus, the present study was designed to identify potential drug candidates from the unique collection of 548 anti-viral compounds (natural and synthetic anti-viral), which target SARS-CoV-2 structural proteins. High-end molecular docking analysis was performed to characterize the binding affinity of the selected drugs-the ligand, with the SARS-CoV-2 structural proteins, while high-level Simulation studies analyzed the stability of drug-protein interactions. The present study identified rutin, a bioflavonoid and the antibiotic, doxycycline, as the most potent inhibitor of SARS-CoV-2 envelope protein. Caffeic acid and ferulic acid were found to inhibit SARS-CoV-2 membrane protein while the anti-viral agent's simeprevir and grazoprevir showed a high binding affinity for nucleocapsid protein. All these compounds not only showed excellent pharmacokinetic properties, absorption, metabolism, minimal toxicity and bioavailability but were also remain stabilized at the active site of proteins during the MD simulation. Thus, the identified lead compounds may act as potential molecules for the development of effective drugs against SARS-CoV-2 by inhibiting the envelope formation, virion assembly and viral pathogenesis.

Introduction

On 31st December 2019, China revealed to the world health organization (WHO) and the rest of the world, the occurrence of symptoms of unexplained pneumonia in a cluster of cases from Wuhan city (Rodríguez-Morales et al., 2020; Zhou et al., 2020a). The causative agent was later identified as a novel strain of coronavirus, named as 2019-nCOV and the disease as COVID-19 (Zhou et al., 2020a). On 30th January 2020, WHO declared the outbreak of COVID-19 as a public health emergency of international concern and also called a pandemic. The 2019-nCoV shared a 79.5% sequence identity to SARS-CoV. Recently, the Coronaviridae Study Group (CSG) of the International Committee on Taxonomy of Viruses (ICTV) renamed it SARS-CoV-2 (Coronaviridae Study Group of the International Committee on Taxonomy of Viruses, 2020). The ever-growing infections and the mortality rate across the globe have called an urgent emergency for specific anti-COVID-19 therapeutics and extensive screening of presently available drugs for the treatment and prevention of SARS-CoV-2. Coronaviruses (CoVs) are enveloped positive-stranded RNA viruses and Coronaviridae can be subdivided into four groups- alpha-, beta-, gamma- and delta- CoV (Perlman and Netland, 2009; Fehr and Perlman, 2015). Members of this virus family infect the mammalian respiratory organ from the upper respiratory tract (URTIs) to the lower respiratory tract (LRTIs) and gastrointestinal tract by incompletely understood mechanisms (Fehr and Perlman, 2015; Cong and Ren, 2014). SARS-COV-2 is the seventh-known SARS virus that will infect people after 229E, NL63, OC43, HKU1, MERS-CoV and the original SARS-CoV (Zhu et al., 2020). SARS-CoV-2 is a member of the subgenus Sarbecovirus (beta-CoV lineage B) with genomic sizes approx 30,000 bases in length and distinctive among other beta coronaviruses in its corporation of a polybasic cleavage site (Wong et al., 2019; Coutard et al., 2020; Hoffmann et al., 2020; Zhou et al., 2020b).

CoVs are quite structurally complex, along with a variety of minor components like non-structural and host cell-derived proteins. The purified viruses particles consist of four structural proteins, namely- Membrane (M) protein, Spike (S) protein, Nucleocapsid (N) protein and Envelope (E) protein (Dent et al., 2015). Recent studies state that SARS-CoV strongly depends on the M, E and N protein for the envelope formation, virion assembly and pathogenesis, i.e., to form a complete infectious virion (DeDiego et al., 2007). Apart from their primary roles, they are also involved in other aspects of the replication cycle.

The M protein is the crucial player in CoV assembly. When expressed alone, it accumulates in the Golgi complex in the homo-multimeric complex and is insufficient for virion formation (Klumperman et al., 1994; Neuman et al., 2011; Lim and Liu, 2001). However, in combination with the E protein, virus-like particles (VLPs) are assembled, indicating the combination of M and E proteins as the minimal requirements for the envelope formation or the production and release of VLPs (Bos et al., 1996; Vennema et al., 1996).

The N protein is the only structural protein that directly links up with the replicase-transcriptase complexes (RTCs) (Hurst et al., 2010; Verheije et al., 2010). This protein binds to the CoV RNA genome, which is essential for incorporating the viruses' genetic material into CoV particles (Ma et al., 2010; Kuo et al., 2016). Its function involves entering the host cell, binding to the viral RNA genome and forming the ribonucleoprotein core (McBride et al., 2014; Zeng et al., 2020). Though, for SARS-CoV, the interaction of N with M was found to be independent of viral RNA (He et al., 2004; Luo et al., 2006). The localization of N to the ER-Golgi region has proposed a function in assembly and budding (Tooze et al., 1984; Klumperman et al., 1994). Studies suggest that binding of N protein to M protein stabilizes the nucleocapsid as well as the internal core of the virions, promoting the completion of viral assembly (Escors et al., 2001; Narayanan et al., 2000). Apart from their primary function, N protein is also involved in virus replication, transcription, translation and host pathological response to viral infection (McBride et al., 2014; Hsin et al., 2018).

The E protein contains a single hydrophobic domain and is a minor component of the virion envelope (Venkatagopalan et al., 2015). Though expressed abundantly, only a small amount is incorporated into the viral envelope; the rest of the E protein is localized to ER, Golgi and ERGIC (Nieto-Torres et al., 2011). Studies have shown that the recombinant CoVs lacking E protein exhibit significantly reduced viral titer or incompetent virus propagation (Ortego et al., 2007). Based on literature three roles have been proposed for E protein, i.e., assembly process based on its interaction with M protein, resulting in driving the production of VLPs; release of infectious viruses that require the hydrophobic domain (HD) of the E protein and lastly, form ion channels, which is mainly associated with pathogenesis (Opstelten et al., 1995; Ye and Hogue, 2007; Ruch and Machamer, 2011; Nieto-Torres et al., 2014; Gupta et al., 2020).

Due to the high infection rate of the COVID-19 patient, it calls for an urgent need for effective therapies or existing treatment strategies (Bruning et al., 2018; Gretebeck and Subbarao, 2015). Due to the extensive membrane perturbation caused by the CoVs, it is not surprising that the membrane-bound proteins have been considered here as a potential target for the development of drugs that could act as anti-viral (Dyall et al., 2017). Since an extensive amount of research and time is a prerequisite for discovering a compound as potent inhibitors or drugs, in-silico evaluation can bridge that gap extensively.

The present study aims to identify potential molecules against SAR-CoV-2 proteins, responsible for envelope formation, virion assembly and pathogenesis. Both natural and synthetic anti-viral compounds were selected for virtual screening against the three structural proteins, i.e., envelope (E), membrane (M) and nucleocapsid (N) protein. Molecular docking and MD simulation results suggested that the natural compounds rutin, caffeic acid, ferulic acid, synthetic anti-virals doxycycline, grazoprevir and simeprevir may be explored as promising drug candidates in the therapy of COVID-19.

Materials and methods

Sequence retrieval

SARS-CoV-2 envelope (E) and membrane (M) protein sequences were taken from Genbank ® of the National center for biotechnology information (NCBI) (Sayers et al., 2019). The E protein and M protein sequences were assembled in FASTA format from the NCBI database with GenBank accession number MT308700.1 and MT093631.2.

Homology modeling

Homology n class="Gene">modeling for both E and M protein was accomplished by I-TASSER online platform for protein structure and function predictions (Yang and Zhang, 2015). 3D models of proteins were built based on multi-threading alignments by LOMETS (Wu and Zhang, 2007) in I-TASSER itself. I-TASSER only uses the template of the highest significance with the best normalized Z-score (>1) that indicates a proper alignment and vice versa (Wu and Zhang, 2007). In the prediction of the 3D structure by threading, the protein with PDB ID: 2MM4 with Z-score 7.01 and PDB ID: 4f91B with Z-score of 1.15 was operated as a template for E and M protein respectively (Table S1). The crystal structure (3D structure) of SARS-CoV-2 nucleocapsid (N) protein was downloaded from the protein database (PDB ID: 6M3M) and saved in PDB format (Kang et al., 2020).

Energy minimization and model validation

Energy minimization of E, M and N protein structure was carried out by YASARA Energy Minimization Server (Krieger et al., 2009) to obtain an energy-minimized and highly stable protein structure validated by PROCHECK (Laskowski et al., 1993). Structural quality and reliability of E and M protein structures were validated through ERRAT, Verify3D and ProSa (Colovos and Yeates, 1993; Eisenberg et al., 1997; Wiederstein and Sippl, 2007).

Prediction of the active or binding site

Binding site residues were anticipated through literature study and different pocket-binding site-recognition web servers such as the CASTp server, and the HotSpot Wizard 3.0 server (Pal et al., 2020; Tian et al., 2018; Lahiri et al., 2019). CASTp 3.0 provides dependable, inclusive and global topological identifications and dimensions of protein designating the identification of residues in the binding site pocket and its volume, cavities and channels. The binding pocket size with the greater surface area was considered the active site and the amino acid residues in it were also generated and shown. HotSpot Wizard 3.0, on the other hand, is a semi-automated process for determining the pocket binding site or hotspots improving the protein stability, catalytic activity, substrate specificity and enantioselectivity. HotSpot Wizard 3.0 server comprises sequence, structural and evolutionary information obtained from 3 databases and 20 computational tools. The functional hotspots depict the functional residues in the binding pocket hotspot (Tian et al., 2018; Lahiri et al., 2019).

Ligand selection and ligand file preparation

Both natural anti-viral compounds and synthetic anti-viral drugs were selected as ligands against E, M and N protein. More than 200 natural compounds, including alkaloids, flavonoids, quinone, tannins, terpenes, steroids, thiophenes, polyacetylenes, lactones, butenolide and lectins were selected as ligands (El Sayed, 2000; Bekhit and Bekhit, 2014). In contrast, an anti-viral compound library with a unique collection of 348 anti-viral compounds used for investigating novel anti-viral drugs was selected from Selleckchem.com (https://www.selleckchem.com/). All the natural compounds were obtained from the PubChem database in smi format (SMILES format) and stored in a single file, whereas the anti-viral compounds from Selleckchem.com were downloaded and saved in sdf (MDL MOL format) file. The SMILES and MDL MOL format was converted to pdbqt (AutoDock PDBQT format) file by OPENBABEL, to molecular docking.

Molecular docking and interaction analysis

In molecular docking, validation of docking protocol is necessary to ensure that ligands bind within the binding pocket in the correct conformation, which is done by validating the size and center of the grid box's coordinates across the binding pocket (Lim et al., 2011). PyRx, a virtual screening software for computational drug discovery, was used to screen the ligand files against E, M and N protein (Dallakyan and Olson, 2015). PyRx uses AutoDock 4 and AutoDock Vina as docking software implying the Lamarckian Genetic Algorithm and Empirical Free Energy Scoring Function. PyRx was carried out using our choice of inhibitors on the predicted energy minimized E, M and N protein structures. Using PyRx platform, the macromolecular structure of all the three proteins (E, M and N) and inhibitors were prepared and then docking was performed into the binding site residues inside a grid box with X, Y and Z axis and dimensions adjusted to 20.74 Å × 20.96 Å × 11.75 Å, 14.93 Å × 62.33 Å × 24. 06 Å and 18.19 Å × 15.05 Å × 13.62 Å for E, M and N protein respectively. The docking protocol was then run at exhaustiveness of 8 and set to only output the lowest energy pose. The interactions between our targeted proteins and the ligands were studied using Ligplot (Wallace et al., 1995), and the figures were processed and prepared using Pymol Molecular Visualization Software (Lill and Danielson, 2011).

Pharmacokinetics studies

The pharmacological profiling of the selected ligands was carried out by analyzing Lipinski's violation of 5, which was evaluated through Molsoft L.L.C.: Drug-Likeness and molecular property prediction for drug-likeness (http://www.molsoft.com/mprop/). Apart from evaluating Lipinski violation, Molsoft L.L.C labeled all the selected inhibitors with an individual drug score.

Molinspiration (https://molinspiration.com/cgi-bin/properties) was used to check the bioactivity of the selected inhibitors, while the ADMET studies were carried by the admetSAR prediction tool (http://lmmd.ecust.edu.cn/admetsar2/).

Molecular dynamics (MD) simulation

All-atoms MD simulations were performed on the atomic coordinates of the best-docked complex of SARS-CoV-2 target proteins, E, M and N, using Gromacs v5.1.4 with force field CHARMM27 and water model TIP3P (Abraham et al., 2015) and the ligands parameters were defined from Zoete et al. (Zoete et al., 2011). The simulation box was defined with buffer distance (10 Å) from the centrally placed protein-ligand complex. The prepared system was solvated with water molecules and neutralized with the addition of 0.15 M counter ions (Na+ and Cl−) (Joung and Cheatham 3rd, 2008). The energy minimization process involves 50,000 steps for each steepest descent, followed by conjugant gradients. PBC condition was defined for x, y and z directions (Darden et al., 1993) simulations were performed at physiological temperature, 300 K. SHAKE algorithm was applied to constrain all bonding involved hydrogen and long-range electrostatic forces treated with PME (Particle mesh Ewald). The system was equilibrated in two steps, NVT and NPT at 300 K for a period of 500 ps. During the simulation, Berendsen thermostat (Berendsen et al., 1987) and Parrinello-Rahman pressure (Parrinello and Rahman, 1980) were used to maintain pressure and temperature. LINC algorithm was used to constrain the bonds and angles (Hess et al., 1997). The van der Waals interactions are taken care of by LJ potential with a cutoff of 0.10 nm. Using the NPT ensemble, production runs were performed for the period of 100, with time integration. The energy, velocity and trajectory were updated at the time interval of 10 ps. All production runs were done on CUDA enabled Tesla GPU machine (DELL T640 with V100 GPU) and OS Centos 7 (Singh et al., 2019; Prakash et al., 2018) and the Gromacs utilities were used for the analyses of obtained MD trajectories.

Results

I-TASSER evaluates the full-length model using paran class="Gene">meters such as C-score, which is the confidence score, TM-score (Template Modeling Score) and RMSD (root mean square deviation) of the estimated models (Wu and Zhang, 2007). The C-score value is usually in the range of −5 to 2 with a higher value signifies a model with a higher confidence score and vice-versa. The TM-score measures the structural similarity between the query and template structures. TM-score for a model with value < 0.17 implies random similarity, while values > 0.5 signify a model of correct topology and structure (Yang and Zhang, 2015). I-TASSER platform generated five different protein models for both E and M protein. We selected the best models, which were Model 1 in the case of both E and M protein with a confidence score of −0.81 and −3.35 and TM-score of 0.61 ± 0.14 and 0.34 ± 0.12, respectively (Table S1). The C-score for both the models was adequate, signifying the proper homology modeling for both the proteins and acceptable TM-score and RMSD values (Table S1). The crystal structure of the SARS-CoV-2 nucleocapsid (N) protein N-terminal RNA binding domain (PDB ID: 6M3M) was downloaded from the protein database and saved in PDB format (Kang et al., 2020).

YASARA was used for energy minimization for E, M and N protein (Krieger et al., 2009). The energy minimized structures for all three proteins with low energy and high stability was validated using PROCHECK (Laskowski et al., 1993). The energy minimized structure of both E and M protein validated through PROCHECK for Ramachandran plot analysis, which showed 99.9% and 99% residues in favored, additionally allowed and generously allowed regions, respectively (Figs. S2A & S3A). No residues were assigned in the disallowed regions for E protein, whereas only 1% of the residues were in the disallowed regions for M protein. The energy minimized N protein (PDB ID: 6M3M) showed that all residues were in favored, additionally allowed and generously allowed regions (Fig. S4A). Thus, we obtained good quality models for E, M and N protein. Errat, Verify3D and ProSa (Colovos and Yeates, 1993; Eisenberg et al., 1997; Wiederstein and Sippl, 2007) further evaluated the energy minimized structure with the calculation of Errat value of 96.92% and 95.93% for E and M protein structure, respectively (Figs. S2B & S3B). The Z-score and energy plots by ProSa indicated an overall good quality model for both E and M protein.

Prediction of the binding site

The active site of protein was generated using a literature study and cross verified with CASTp and Hotspot Wizard 3.0 web servers (Pal et al., 2020; Tian et al., 2018; Lahiri et al., 2019; Kang et al., 2020). The functionality of binding pockets and its residues evaluated from literature & web server were compared and the residues, which superimposed each other, were considered the binding site residues for E, M and N protein in our docking study (Fig. 1A–C & Table 1 ). Binding pockets of E, M and N protein are shown in Fig. S1.

Fig. 1

(A) Structure of envelope (E) protein predicted by I-Tasser with active site residue shown in green. (B) Structure of the membrane (M) protein predicted by I-Tasser with active site residue shown in red. (C) Structure of the Nucleocapsid (N) protein (pdb id: 6M3M) with active site residue shown in orange. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

Summary of predicted binding sites of the protein models E and M from literatures and those predicted for the model N via CASTp 3.0 and HotSpot Wizard 3.0.

Protein	Residues at the binding site
Envelope (E)	Ala32, Thr35, Tyr42, Tyr57
Membrane (M)	Lys50, Leu51, Leu54, Leu93, Ala98
Nucleocapsid (N)	Asn48, Thr49, Ala50, Ser51, Phe53, Arg88, Ala90, Thr91, Arg92, Arg93, Asp103, Leu104, Ser105, Arg107, Tyr109, Tyr111, Pro117, Arg149

The docking study carried out with the natural as well as synthetic anti-viral agents as ligands against the E, M and N protein of SARS-CoV-2 provided the structural insight of molecular interactions. We find that most of the natural ligands interact with all three target proteins more or less efficiently, but the ligands with the lowest dock energy against a particular protein were assigned as a potential target and their interaction studies in detail. Natural bioflavonoid rutin and anti-bacterial as well anti-viral agent doxycycline (Rothan et al., 2014; Wu et al., 2015) showed the best result with E protein, having a docking score of −9.3 kcal/mol and −9.1 kcal/mol, respectively (Table 2 ). The two hydrocinnamic acids, i.e., caffeic acid and ferulic acid is belonging to phenylpropanoids, obtained as best docked against the M protein of SARS-CoV-2 with docking score of −8.4 kcal/mol and −8.3 kcal/mol, respectively (Table 2). Simeprevir and grazoprevir showed an excellent docking result against N protein with docking energy of −8.7 kcal/mol and −8.5 kcal/mol, respectively, in PyRx (Table 2). The best docking pattern of compounds are shown in Fig. 2, Fig. 3, Fig. 4 , respectively.

Table 2

Binding energy values and interactions of the selected inhibitors with the key residues of E, M and N protein of SARS CoV-2 evaluated by PyRx docking.

Ligands	Targets
	Envelope (E)		Membrane (M)		Nucleocapsid (N)
	Binding energy (kcal/mol)	Key residues interraction	Binding energy (kcal/mol)	Key residues interraction	Binding energy (kcal/mol)	Key residues interraction
Rutin	−9.3	Thr35	−7.7	Lys50	−7.9	Thr91, Arg92, Tyr109
Doxycycline	−9.1	Thr35	+47.60	No H-bond interaction	+23.10	No H-bond interaction
Caffeic acid	−6.1	Thr35	−8.4	Lys50	−7.4	Ala90, Ser105
Ferulic acid	−6.0	Thr35	−8.3	Lys50	−7.2	Thr91, Tyr111
Simeprevir	+14.20	No H-bond interaction	+63.60	No H-bond interaction	−8.7	Thr91, Arg93
Grazoprevir	+22.70	No H-bond interaction	+66.40	No H-bond interaction	−8.5	Tyr109, Tyr111

Significance of bold: Indicate Lowest binding energy for respective ligand with target protein.

Fig. 2

The docking results of (A) Rutin and (B) Doxycycline within binding pocket of the E protein of SARS-CoV-2. Hydrogen bonded interactions are shown as yellow dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

The docking results of (A) Caffeic acid and (B) Ferulic acid within binding pocket of the M protein of SARS-CoV-2. Hydrogen bonded interactions are shown as yellow dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

The docking results of (A) Simeprevir and (B) Grazoprevir within binding pocket of the N protein of SARS-CoV-2. Hydrogen bonded interactions are shown as yellow dotted lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Ligplot analysis reveals H-bond formation between selected ligands with their target protein (Fig. 5 A-F). PyRx docking was carried out with grid box formation to include all the active site residues for all three proteins, i.e., E, M and N, respectively. Interaction studies mainly focused on the efficient binding of our selected ligands with binding site residues within the grid box, thus evaluating a proper and efficient docking method. Results show that rutin involved in H-bond with the active site residues of E protein, such as Leu31, Thr35, Ser55, Arg69 and Pro71 (Fig. 5A). In contrast, doxycycline engaged in H-bond interaction with residues Leu31, Thr35, Val52 and Ser55 of E protein (Fig. 5B), indicating their strong tendency for efficient binding towards E protein and can inhibit envelope formation, virion assembly and pathogenesis. Thus, Ligplot analysis of rutin and doxycycline with E protein shows a frequent H-bond interaction with the active site residue, Thr35.

Fig. 5

Diagrammatic sketch illustrating the interactions between (A) E protein and Rutin. (B) E proten and Doxycycline. (C) M protein and Caffeic acid. (D) M protein and Ferulic acid. (E) N protein and Simeprevir and (F) N protein and Grazoprevir by LigPlot. Ligand is shown in purple and: green dashed lines indicate hydrogen bonds with distance in angstrom (Å), spoked red arcs indicate hydrophobic contacts, atoms are shown in black for carbon, blue for nitrogen, red represents oxygen and yellow represents sulfur. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Caffeic acid and ferulic acid formed H-bond with the active site residue Lys50 of M protein and having frequent interaction with Tyr47 (Fig. 5 C–D). Results of compounds with N protein show the best docking score with simeprevir and grazoprevir. Ligplot analysis shows that simeprevir formed H-bond with the active site residues Thr91 and Arg93 (Fig. 5E), whereas, grazoprevir shows interaction with two tyrosine residues Tyr109 and Tyr111 through H-bonding (Fig. 5F).

Adding to our excellent docking result and interaction analysis which assessed rutin to be the best inhibitor against E protein, it was also observed that rutin interact and form H-bond with active site residues of both the M and N protein with a docking score of −7.7 kcal/mol and −7.9 kcal/mol respectively in PyRx (Table 2). While caffeic acid and ferulic acid were found to form H-bond with active site residues of both E and N protein. The docking score of caffeic acid against E and N protein was found to be −6.1 kcal/mol and −7.4 kcal/mol respectively (Table 2), whereas, the docking scores of ferulic acid was −6.0 kcal/mol and −7.2 kcal/mol against E and N protein respectively (Table 2). The synthetic anti-viral agent's grazoprevir and simeprevir, on the other hand, displayed a weak affinity towards the other two proteins (E and M) apart from their interaction with N protein (Fig. S5 C–F). Doxycycline, on the other hand, displayed a weak affinity towards the other two proteins (M and N) apart from their interaction with E protein. This observation was validated by Ligplot interaction analysis (Fig. S5 A–B).

The pharmacological studies first assessed the selected ligands against E, M and N protein for a good oral administration based on the Lipinski rule of five (Lipinski, 2004), which was evaluated through Molsoft L.L.C.: Drug-Likeness and molecular property prediction. Lipinski's “rule of five” is a heuristic approach for predicting drug-likeness stating that molecules had Molecular weight (M.W. < 500 Da), high lipophilicity expressed as LogP (LogP < 5), hydrogen bond donors (HBD < 5) and hydrogen bond acceptors (HBA < 10) have good absorption or permeation across the cell membrane. Doxycycline, caffeic acid and ferulic acid did not show any Lipinski violation, while rutin was found to violate Lipinski rule as observed from Table 3 in terms of its high molecular weight, more than 5 HBD and more than 10 HBA. Apart from possessing a high molecular weight, both simeprevir and grazoprevir did not violate any other Lipinski rule. Drug molecules of higher molecular weight are not readily transported and have poor absorption compared to molecules with low molecular weight (Srimai et al., 2013). The absorption percentage (AB%) was calculated using the formula (Zhao et al., 2002; Husain et al., 2016).

Table 3

Evaluation of Lipinski's rule of 5 along with drug likeness score through Molsoft L.L.C.: Drug-Likeness and molecular property prediction of the selected inhibitors after docking.

Ligands	Mass (<500)	Hydrogen bond donor (<5)	Hydrogen bond acceptor (<10)	LOGP (<5)	TPSA (≤140)	Absorption percentage (AB%) (>50%)	Drug likeness Score (>0)
Rutin	610.15	10	16	−1.55	213.63	35	0.91
Doxycycline	444.15	6	9	−0.44	138.73	61	1.44
Caffeic acid	180.04	3	4	1.27	61.72	87	−0.35
Ferulic acid	194.06	2	4	1.61	52.80	90	−0.61
Simeprevir	749.29	2	9	4.97	131.74	63	1.17
Grazoprevir	766.34	3	10	4.64	165.83	51	0.88

The AB% (>50%) of a drug is a good indication of its excellent bioavailability, distribution and circulation by oral route a (Ghose et al., 1999; Hebert, 2013; Mokhnache et al., 2019). Only rutin was found to possess an AB% of 35 (<50%), which was due to its high polar surface area (PSA) of 213.63 Ǻ2 while PSA ≤ 140 Ǻ2 was proposed optimal (Veber et al., 2002). All the remaining inhibitors apart from rutin possessed AB% greater than 50% (Table 3), an excellent evaluation of these selected ligands to have proper bioavailability for an oral route (Ghose et al., 1999). Only caffeic acid and ferulic acid showed a negative drug-likeness score as assessed by Molsoft L.L.C. software (Table 3). Ligands with negative drug score are generally not considered as a potential drug candidate (Chandrasekaran and Thilak Kumar, 2016). However, in vitro anti-viral effect of caffeic acid on type 1 human immunodeficiency virus, hepatitis B & C virus, Herpes Simplex Virus (HSV), VSV-Ebola pseudotyped and vaccinia viruses have already been demonstrated (Lai et al., 1992; Wang et al., 2009; Langland et al., 2018; Shirasago et al., 2019). At the same time, ferulic acid was inhibitory to murine macrophage cell line in response to the respiratory syncytial virus (RSV), a potential influenza H1N1 Neuraminidase Inhibitor, possessed excellent antioxidant properties for its remedial use against various diseases (Sakai et al., 1999; Hariono et al., 2016; Srinivasan et al., 2007).

The bioactivity of the selected ligands checked through molinspiration and calculated the activity against GPCR ligand, ion channel modulator, a kinase inhibitor, nuclear receptor ligand, protease inhibitor and enzyme inhibitor (Mokhnache et al., 2019). The bioactivity values were interpreted as: active (bioactivity score > 0), moderately active (bioactivity score: −5.0–0.0) and inactive (bioactivity score < −5.0) (Ungell, 1997). Rutin and doxycycline were evaluated as active enzyme inhibitors with values 0.12 and 0.59, respectively. Caffeic acid and ferulic acid were moderately active against enzyme inhibition as well as nuclear receptor ligand. Simeprevir and grazoprevir were evaluated as active protease inhibitors by molinspiration, which supports the fact that these two compounds are an HCV NS3/4A protease inhibitor (Fried et al., 2013; Talwani et al., 2013; Hayashi et al., 2014; Summa et al., 2012; Yeh et al., 2018). The bioactivity prediction by molinspiration is tabulated in Table 4 .

Table 4

Bioactivity prediction of the selected inhibitors against SAR CoV-2 by molinspiration.

Ligands	GPCR ligand	Ion channel modulator	Kinase Inhibitor	Nuclear receptor ligand	Protease Inhibitor	Enzyme Inhibitor
Rutin	−0.05	−0.52	−0.14	−0.23	−0.07	0.12
Doxycycline	−0.47	−0.43	−0.68	−0.38	−0.07	0.59
Caffeic acid	−0.48	−0.23	−0.81	−0.10	−0.79	−0.09
Ferulic acid	−0.47	−0.30	−0.72	−0.14	−0.81	−0.12
Simeprevir	−0.53	−1.96	−1.36	−1.59	0.10	−0.87
Grazoprevir	−0.50	−2.03	−1.66	−1.84	0.19	−1.00

The pharmacodynamic study for the selected ligands was carried out through admetSAR to understand the drug's action inside a host's body. The ADMET study in this work focused on the parameters that can define absorption, distribution, metabolism, excretion, toxicity, solubility (LogS), human intestinal absorption (HIA), CaCO-2 permeability, P-glycoprotein substrate inhibition, cytochrome substrate/inhibitor, AMES toxicity and acute rat toxicity (LD50). Due to the emerging and effectiveness of remdesivir, dexamethasone, chloroquine and its derivative hydroxychloroquine in treating COVID-19 patients in clinical trials (Wang et al., 2020), the ADMET properties of the selected ligands were thus compared with it. All the selected ligands showed optimal solubility with values higher than −4 (>−4) (Mokhnache et al., 2019) and immensely satisfying results with some of the values even better than the control drugs, as shown in Table 5 .

Table 5

Pharmacodynamics profile of the selected inhibitors with the Hydroxy-chloroquine as control through admetSAR.

Ligands	Log S (>−4)	Blood Brain Barrier (BBB)	Human intestinal absorption (HIA)	Caco2 permeability	CYP substrate/Inhibitor	Ames toxicity	Carcinogenicity	LD50 (rat acute toxicity) (mol/kg)
Rutin	−2.77	0.85	0.80	0.91	Non-substrate/Non-inhibitor	Non-toxic	Non-carcinogen	2.49
Doxycycline	−3.07	0.98	0.85	0.87	Substrate/Non-inhibitor	Non-toxic	Non-carcinogen	2.31
Caffeic acid	−1.69	0.63	0.93	0.56	Non-substrate/ Non-inhibitor	Non-toxic	Non-carcinogen	1.40
Ferulic acid	−2.47	0.56	0.96	0.71	Non-substrate/Non-inhibitor	Non-toxic	Non-carcinogen	1.43
Simeprevir	−3.69	0.72	0.99	0.84	Substrate/Inhibitor	Non-toxic	Non-carcinogen	2.64
Grazoprevir	−3.75	0.84	0.92	0.51	Substrate/Inhibitor	Non-toxic	Non-carcinogen	2.62
Hydroxy-chloroquine	−3.56	0.56	0.98	0.51	Non-substrate/Non-inhibitor	Ames-toxic	Non-carcinogen	2.63
Chloroquine	−4.34	0.74	0.99	0.58	Substrate/Non-inhibitor	Non-toxic	Non-carcinogen	2.95
Remdesivir	−3.47	0.74	0.88	0.65	Substrate/Non-inhibitor	Non-toxic	Non-carcinogen	2.71
Dexamethasone	−3.70	0.97	0.99	0.83	Substrate/Non-inhibitor	Non-toxic	Non-carcinogen	2.14

Pharmacodyene">namics profile of the selected inhibitors with the Hydroxy-chloroquine as control through admetSAR.

The selected inhibitors against E, M and N protein of SARS CoV-2 have excellent human intestinal absorption (HIA), good Caco2 permeability values and were expected to be highly absorbed in the intestine upon oral administration (Table 5) (Gautret et al., 2020). Simeprevir showed the highest HIA at 0.99, which is indicative of its better absorption rate. All the selected inhibitors also possessed better values for the blood-brain barrier (BBB) and Caco2 permeability (Table 5). In terms of metabolism, rutin, doxycycline, caffeic acid and ferulic acid were found to be non-inhibitor of cytochrome P450 (CYP 450), indicating its proper metabolism by CYP450 (Table 5). A non-inhibitor of CYP450 means that the drug molecule will not hamper in the biotransformation of drug metabolism by CYP 450 enzyme. All the selected inhibitors were non-AMES toxic and non-carcinogens, while hydroxychloroquine was evaluated as AMES toxic in admetSAR (Table 5).

A compound with lower LD50 doses is n class="Gene">more lethal than the compound having higher LD50. admetSAR analysis revealed all the selected inhibitors with optimal LD50 doses indicating they are nonlethal (Table 5).

Molecular dynamic (MD) simulation studies

Based on the docking and pharmacological analysis, the best and effective ligands against proteins E, M and N were considered for further molecular dynamic (MD) simulation studies. MD simulation is a useful computational tool that provides a better understanding of protein-ligand interactions, atomic description of the spatial orientation of ligands at the active site of a protein, conformational dynamics of protein and active site residues and the molecular stability (Kumar et al., 2019; Luthra et al., 2009; Le Rouzic et al., 2002; Yan et al., 2020). All atoms MD simulations were performed on the protein-ligand complexes of rutin and doxycycline with protein E, caffeic and ferulic acids with M protein and simeprevir and grazoprevir with N protein, for the period of 100 ns at 300 K. Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuations (RMSF) were analyzed for understanding the deviation of Cα atoms of the protein from its backbone and also the fluctuations associated with the amino acid residues of the protein during the simulation (Mishra et al., 2018; Malathi et al., 2019). The structural order parameter, RMSD result of E protein complex with rutin, shows that the structure of the protein-ligand complex stabilized around 40 ns and remains stable till the simulation finished at 100 ns (Fig. 6A). Comparatively, the RMSD trajectory of doxycycline achieves dynamics equilibrium at ~25 ns, and the structural stability of E protein- doxycycline complex seen stable for the remaining period of simulation (Fig. 6A). The docked structure of caffeic and ferulic acids with M protein shows quite stable conformational dynamics during the simulation (Fig. 6A). The complex structure of caffeic with M protein quickly attains equilibrium at RMSD ~0.5 nm during the initial 0–10 ns, which is continued till 100 ns, whereas ferulic acid shows a marginal increase in RMSD during 0–20 ns which dropped at ~20 ns (Fig. 6A).

Fig. 6

The 100 ns MD results of four protein–ligand complexes (E-Rutin, E-Doxycycline, M-Caffeic acid, and M-Rerulic acid). (A) RMSD values of backbone atoms. (B) RMSF of carbon alpha of complex structure. (C) Area per residue over trajectory. (D) SASA of the ligands. (E) Rg of backbone atoms. (F) Potential energies.

Moreover, the stable equilibrium can be seen for a period of 20–100 ns. With the marginal perturbation, the structure of simeprevir and grazoprevir with N protein achieved equilibrium ~20 ns, and the conformational dynamics of complex structures find consistent up to 100 ns (Fig. 7A). All the fluctuations for the protein-ligand complexes altogether indicated excellent structural stability of the complexes (Li et al., 2019; Jeyaram et al., 2019).

Fig. 7

The 100 ns MD results of two protein–ligand complexes (N-Simeprevir, and N- Grazoprevir). (A) RMSD values of backbone atoms. (B) RMSF of carbon alpha of complex structure. (C) SASA of the ligands. (D) Rg of backbone atoms. (E) Potential energies.

The RMSF plot of all Cα-atoms of E, M and N protein docked complexes with ligands are shown in Figs. 6B and 7B. We find restricted fluctuation in the active site residues of E and M proteins as compared to adjacent and N-terminal residues. RMSF complies that the average fluctuation of residues belonging to stable secondary conformations remains below 0.25 nm, which is also observed consistent with binding pocket residues and Thr35 for E and Lys50 for M, which was found in H-bond interactions in docking results (Fig. 6B). Whereas RMSF fluctuations for N protein residues against simeprevir and grazoprevir also showed fewer fluctuations (<1.0 nm) in the binding site residues Thr91 and Arg93 for simeprevir and residues Tyr109 and Tyr111 for grazoprevir (Fig. 7B). The binding site residues for all the three proteins in protein-ligand complexes from the RMSF study are favorable in molecular interactions; thus, a stable conformational dynamic of protein-ligand interaction is observed during simulation, which was in excellent agreement with docking analyses.

Furthermore, we also examine the solvent-accessible surface area (SASA) of protein-ligand complexes. SASA ensures the contribution of hydrophobic interactions of the nonpolar amino to the conformation stability of proteins in the solvent environment (Zhang et al., 2020; Bhowmik et al., 2020; Rocco et al., 2008; Kumar et al., 2018). Results show quite a stable SASA plot of E protein complex with rutin and doxycycline, which is stabilized throughout simulation with 42 nm2 (Fig. 6D). Similarly, the stable structure of M protein with caffeic acid and ferulic acid is maintained with SASA value 135 nm2 (Fig. 6D). Whereas, the SASA plots of N protein docked complex with simeprevir and grazoprevir show that an average SASA value of 167 and 167.48 is maintained during the simulation (Fig. 7C). Thus, SASA results indicate no significant changes in the proteins solvent accessible, which again provided the stable conformational dynamics of ligand, bound structures.

To determine the conformational stability of protein-ligands, we also computed the time evolution plot of R g (Figs. 6E and 7D). Results show that the R g trajectory of rutin optimizes at ~40 ns and remains stable during 40–100 ns, signifying that ligand is well accommodated at the active site of the E protein (Fig. 6E). The R g trajectory of doxycycline reaches equilibrium ~20 ns and a stable equilibrium is seen during 20–100 ns, which indicates that ligand is well fitted in the binding pocket. We do not find any significant change in the R g trajectory of Caffeic acid (Fig. 6E). R g achieves equilibrium quickly in a few nanoseconds, and an unperturbed structure is sustained until simulation finished at 100 ns with the average change in R g value 1.83 ± 0.01(Fig. 6E). We find an initial rise in R g of Ferulic acids; however, complex structure smoothly settles ~20 ns, and a compact structure is propagated for the remaining period of simulation (Fig. 6E). The R g plots of simeprevir and grazoprevir with N protein reveals that the structural compactness of complex structures remains stable, with the average change in R g value 15. 3 ± 0.02 and 15. 25 ± 0.02 throughout the simulation period, thus, favoring the spatial interaction at the active site of N protein (Fig. 7D). Figs. 6F & 7E showed; stable potential energy for all complexes, which confirmed the stability of all six complexes.

The structural fold of the protein is largely stabilized by the various inter and intra H-bond interactions. These H-bond interactions play a crucial role to accommodate the ligand at the active site of a protein. Thus, we also calculated the evolution plot of H-bond interactions with the respective ligands, shown in Fig. 8 . The result shows the maximum propensity of five H-bonds between E protein and rutin; however, four H-bonds remains consistent during the simulation (Fig. 8A). Whereas, four H-bonds observed between doxycycline and E protein, remains stable, which can be seen during the simulation, time 20–100 ns (Fig. 8A). Caffeic and ferulic acid formed a maximum of three H-bond with M protein; however, only two remains stable during the progression of simulation (Fig. 8 A & B). Simeprevir with N protein shows the possibility of four H-bond interactions at the active site, out of which three remains consistent throughout the simulation period (Fig. 8C). However, grazoprevir shows five H-bonding with active site residues, but three remain stably bond until the simulation finished at 100 ns (Fig. 8C). These results found consistent, as we observed in molecular docking. Thus, MD simulation of the protein-ligand complexes for all the three proteins with selected ligands leads to the establishment of structurally compact and stable conformation dynamics protein-ligand molecular interactions.

Fig. 8

Hydrogen bond analysis. Hydrogen bond formation between the target protein and ligand over 100 ns MD simulation. (A) E-Rutin, E-Doxycycline, and M-Caffeic acid. (B) M-Ferulic acid. (C) N-Simeprevir, and N-Grazoprevir.

Hydrogen boene">nd analysis. n class="Chemical">Hydrogen bond formation between the target protein and ligand over 100 ns MD simulation. (A) E-Rutin, E-Doxycycline, and M-Caffeic acid. (B) M-Ferulic acid. (C) N-Simeprevir, and N-Grazoprevir.

The H-bond analysis through the gmx hydrogen bond showed stable and efficient H-bond interaction of rutin and doxycycline against E protein without leaving the pocket (Fig. 8A). Caffeic acid and ferulic acid showed quite stable H-bond against M protein (Fig. 8A). Ferulic acid does leave the pocket during the 35–40 ns run and then again maintained a stable H-bond during the 100 ns simulation run with M protein while caffeic acid leaves the pocket at the last 85-90 ns but again formed stable and reliable H-bond afterward as shown in Fig. 8 A and B respectively. Simeprevir and grazoprevir, on the other hand, showed a strong H-bonding throughout the 100 ns simulation run against protein N without leaving the pocket indicating their strong tendency towards N protein and probably inhibit viral assembly and further pathogenesis (Fig. 8C). Thus, MD simulation of the protein-ligand complexes for all the three proteins against their selected ligands leads to the establishment of our complexes to be basically and structurally stable and reliable.

Discussions

In this study we evaluated the envelope (E), membrane (M) and nucleocapsid (N) protein as a potential drug target and molecular docking was carried out to find a potential inhibitor against these three proteins, which are essential for coronavirus assembly, envelop formation, pathogenesis and replication. We found that rutin, a bioflavonoid and doxycycline, a known anti-bacterial as well as anti-viral compound, possesses a higher binding affinity with E-protein. Similarly, the phenylpropanoids, caffeic acid and ferulic acid showed a higher affinity with M protein, while the anti-viral agent's simeprevir and grazoprevir demonstrated higher affinity against the N protein. Moreover, the natural compounds, i.e., rutin, caffeic acid and ferulic acid, were found to interact with the other target proteins, revealing a high tendency of these compounds to interact with the potential targets of SARS-CoV-2. Doxycycline interacted efficiently only with E protein, while simeprevir and grazoprevir had more affinity towards the N protein in terms of the lowest binding energy (Table 2).

Pharmacological studies of the selected inhibitors for the Lipinski rule of 5 indicated rutin, simeprevir and grazoprevir with a high molecular weight (Srimai et al., 2013). Although simeprevir and grazoprevir did not violate any other Lipinski parameters, as shown in Table 3. However, earlier studies described rutin as a potential inhibitor against rotavirus and enterovirus with IC50 values of 98 μM and 109.63 μM, respectively (Bae et al., 2000). Also, rutin has powerful antioxidants and different pharmacological properties (Lin et al., 2012; Sharma et al., 2013). Doxycycline, caffeic acid and ferulic acid did not violate the Lipinski rule and possessed excellent pharmacological properties for a better oral drug (Table 3). Generally, compounds with a negative drug score are not considered as a potential drug (Yeh et al., 2018), but the anti-viral activity of caffeic acid and ferulic acid against different viruses as discussed above, cannot be neglected (Lai et al., 1992; Wang et al., 2009; Langland et al., 2018; Shirasago et al., 2019), and both can be considered as potential inhibitors against M protein. Simeprevir and grazoprevir with high molecular weights have been studied and well demonstrated for their anti-viral activity against hepatitis C genotype 1 (Sakai et al., 1999; Hariono et al., 2016; Srinivasan et al., 2007). Safety and tolerability of simeprevir and grazoprevir were also well documented. Simeprevir with a single oral administration combined with ribavirin demonstrated potent antiviral activity against hepatitis C genotype 1b. It was well tolerated, with no new or further adverse effects other than the mild, reversible and asymptomatic bilirubin increase without any hepatic damage (Fried et al., 2013; Talwani et al., 2013; Hayashi et al., 2014). Grazoprevir was well tolerable with a dose of once-daily or in combination with Elbasvir against Hepatitis C replication and possessed good pharmacokinetic profiles (Yeh et al., 2018). The bioactivity of the selected inhibitors against SARS CoV-2 structural E, M and N protein was carried out, which was also supportive of our selected inhibitors. The ADMET studies for the selected inhibitors were checked and evaluated through the admetSAR. All the compounds were compared with the control/ reference drugs taken in the study.

Lastly, results from MD simulation supported our choice of inhibitors evaluated against their particular protein target. The protein-ligand complexes for E, M and N were stable throughout the simulation period as interpreted by RMSD, RMSF, SASA, Rg, PE and H-bond graphs (Fig. 6, Fig. 7, Fig. 8).

Thus molecular, docking, pharmacokinetic profiling and MD simulation outcome supports our choice of inhibitors as potential inhibitors against E, M and N protein of SARS CoV-2 and may inhibit envelope formation, virion assembly and viral pathogenesis.

Conclusion

The study evaluated the bioflavonoid rutin and the anti-bacterial and anti-viral agent doxycycline as a potential inhibitor for E protein; Caffeic acid and ferulic acid of the class phenylpropanoids as a potential inhibitor for M protein while the anti-viral agent's simeprevir and grazoprevir to have high affinity against the N protein. Besides that the natural anti-viral agent's rutin, a bioflavonoid and the two hydrocinnamic acids, i.e., caffeic acid and ferulic acid which are a class of phenylpropanoids were found to interact with all the three structural proteins of SARS-CoV-2 as well, considered under this study revealed a strong tendency and efficiency of the natural anti-viral compounds against SARS-CoV-2 proteins. The anti-bacterial and anti-viral agent doxycycline docked efficiently only against the E protein, whereas the protease inhibitors simeprevir and grazoprevir were particularly efficient against N protein only. The protein-ligand complexes for E, M and N proteins were stable throughout MD simulation period. The pharmacological profiling evaluated our selected inhibitors as a potent drug, indicating our stronghold on the choice of inhibitors and their potential interaction with the SARS-CoV-2 structural proteins that might lead to the interruption in the interaction between these surface structural proteins that are responsible for viral assembly and pathogenesis. However, more in vivo work may be required to authenticate the effectiveness and safety of these inhibitors against SARS-CoV-2.

Author contribution statements

DK, DB, RJ, RN, NK and AP, collected data. DB, RN, AP and DK prepared original draft of manuscript and edited the manuscript. DK designed the study, edited the final version of manuscript and conceptualized the study. All authors read the manuscript and approved the final manuscript.

Declaration of Competing Interest

The authors declare that no conflict of interest exists.