Identification and characterization of novel RdRp and Nsp15 inhibitors for SARS-COV2 using computational approach.

The World Health Organization has declared COVID-19 as a global health emergency. COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and highlights an urgent need for therapeutics. Here, we have employed a series of computer-aided drug repurposing campaign to discover inhibitors of RNA dependent RNA polymerase (RdRp) and Nsp15/EndoU. Subsequently, MD simulation has been performed to observe dynamic behavior of identified leads at the active site of RdRp and Nsp15. We successfully identified novel lead molecule such as Alectinib for RdRp while Naldemedine and Ergotamine for NSP15. These lead molecules were accommodated in the active site of the enzyme and stabilized by the networks of the hydrogen bond, pi type and hydrophobic interaction with key residues of either target. Interestingly, identified compounds show molecular mimicry in terms of molecular interactions with key residues of RdRp and Nsp15 essential for catalysis and substrate interaction. Previously, Alectinib, Naldemedine and Ergotamine were used as drug in different diseases might be repurposed against selected protein targets of COVID19. Finally, we propose that the identified inhibitors represent a novel lead molecule to design a more effective inhibitor to stop the progress of pathogen.Communicated by Ramaswamy H. Sarma.

Introduction

Coronaviruses (CoVs) is an RNA virus belonging to subfamily Coronavirinae within the family Coronaviridae, which infects humans, birds and mammals (Gorbalenya et al., 2004; Snijder et al., 2003). It was first identified in the Wuhan city of Hubei province of China during the epidemic outbreak in 2019 and is now designated as COVID-19. Electron microscopic image shows a spike-like projection on the outer surface of the virus, giving it crown-like shape, hence the name (Richman et al., 2016). In China, over the last two decades, β genera coronavirus of an animal was a crossover to humans through the intermediate host, which resulted in severe respiratory syndrome causing the death of 916 patients with the fatality rate of 11% (Chan-Yeung & Xu, 2003). In 2012, bat originated MERS coronavirus was also crossover to humans via camels as an intermediate host affecting 2494 people and caused 858 death with the mortality rate of 34% in Saudi Arabia (https://www.who.int/emergencies/mers-cov/en/). Wuhan outbreak showed an increased number of pneumonia patients for unknown reasons. However, high-throughput screening of samples obtained from the lower respiratory tract showed novel coronavirus and was designated as SARS-CoV-2 and 2019-nCoV (Hasan et al., 2020; Zhu et al., 2020). This novel coronavirus showed more than 70% similarity with severe acute respiratory syndrome coronavirus (SARS-CoV), while more than 95% homology was observed with bat coronavirus (Singhal, 2020). According to the world health organization, the fatality rate of COVID-19 was 2%. Till now seven different strains of HCoVs have been found, which includes SARS, OC43, COVID-19, HKU1 and MERS strains of HCoVs belongs to β genera coronavirus while NL63 and 229E strains of HCoVs belongs to α genera coronavirus (Elfiky et al., 2017; Hui et al., 2020). To date, no drugs, monoclonal antibodies or vaccines used or approved in treatment against COVID-19 infections (Colson et al., 2020; Liu et al., 2020). Therefore, the pressing need to develop or design novel therapeutic strategies to combat against COVID-19 infections.

HCoVs proteins classified into two groups, namely structural proteins and non-structural proteins. Spike is a structural protein present in a homo-trimeric state and is involved in recognition, attachment and entry of virus in the host cell (Belouzard et al., 2012; Ibrahim et al., 2019). The replicase-transcriptase complex consists of 16 non-structural proteins (Nsp1-16) that are auto-proteolytically cleaved by a single Orf encoded polyprotein. RdRp/Nsp12 is a unique class of nucleic acid polymerases encoded by RNA viruses. RdRp is a catalytic core forming a human right-hand design consisting of thumb, fingers and palm domains (Ago et al., 1999; Bressanelli et al., 1999; Hansen et al., 1997; Lesburg et al., 1999). The fingers and palm domains of RdRP are the most conserved parts and consist of 7 conserved polymerase catalytic motifs, which are responsible for the catalysis of only RNA (Elfiky, 2020; Elfiky & Azzam, 2020; Jia & Gong, 2019). The functional part of RdRP is fused with the catalytic part resulting in the diversity of the whole structure and its regulatory mechanism.

On the other hand, Nsp15 is also a non-structural protein of coronavirus and is recognized as an endoribonuclease (EndoU; Bhardwaj et al., 2004; Ivanov et al., 2004). Coronavirus does not require EndoU/Nsp15 activity for the synthesis of RNA (Deng & Susan, 2018; Deng et al., 2019). It has been proved that viruses with catalytic mutant versions of Nsp15 trigger elevated interferon response and PKR-dependent apoptosis (Deng et al., 2017; Kang et al., 2007; Kindler et al., 2017). Nonetheless, Nsp15 activity plays an essential role in evading host sensing of viral dsRNA (Deng et al., 2017; 2019; Kindler et al., 2017). However, inactivating Nsp15 catalytic activity showed rapid clearance of infectious virus from spleen and liver of infected mice, and macrophage showed early induction of type-1 interferon in the virus-infected mice. Protective immune response was also generated in the infected mice, where Nsp15 catalytic activity was inactivated. Thus, RdRps and Nsp15 are attractive targets for the development of antiviral therapy (Machitani et al., 2020; Deng & Susan, 2018; Wu et al., 2020). The present study aims to identify novel inhibitor molecules against RdRp and Nsp15 targets. In this respect, we predicted 3D structures of RdRps and Nsp15 of COVID-19, followed by stereochemical validation to check the quality of the predicted model. Subsequently, we screened chemical libraries of 3277 compounds (1355 FDA approved drugs and 1922 world not FDA) obtained from ZINC database to identify novel inhibitory molecules against RdRp and Nsp15 using a virtual screening approach. Subsequently, MD simulation has been performed on top two complexes of RdRp and Nsp15. We successfully identified novel lead molecule such as Alectinib for RdRp while Naldemedine and Ergotamine for NSP15. These lead molecules will provide a detailed framework to understand the nature and selectivity of ligand and may help in the development of safe and effective inhibitors against RdRp and Nsp15.

Methods

Homology modeling of RNA-dependent RNA polymerases (RdRPs) and Nsp15

The complete sequences of RdRP (Acc. No: YP_009725307) and Nsp15 (Acc. No: YP_009725310) were retrieved from the NCBI coronavirus2 data hub. The template search was performed using NCBI’s BLASTp program against a protein data bank (Altschul et al., 1997). The blast search gave hits with significant identity against query protein sequences. The SARS-coronavirus structure of NSP12 complexed with NSP7 and NSP8 has a 96% identity with RdRp (PDB ID: 6NUR; Kirchdoerfer & Ward, 2019). Whereas, the crystal structure of NSP15 of SARS-COV shows 88% identity with Nsp15 (PDB ID: 2OZK; Joseph et al., 2007). These crystal structures were used as a template to predict homology models of RdRp and Nsp15. In total, 20 homology models were generated for RdRp and Nsp15 using Modeller 9.23 as per previous studies (Barage et al., 2018; Elmezayen et al., 2020; Šali & Blundell, 1993). Among 20 models, the top model of each structure was selected based on the DOPE score (Shen & Sali, 2006). The structural validation and characterization of the predicted models were performed using PROCHECK (Laskowski et al., 1993) and PROSA-Web Z-Score (Wiederstein & Sippl, 2007).

Preparation of ligand dataset

The FDA approved, and the world not FDA datasets of compounds were retrieved from the ZINC database. Only purchasable compounds from FDA approved (1355 compounds), and the world not FDA (1922 compounds) were extracted from the ZINC database (Sterling & Irwin, 2015). The world not FDA drug dataset consist drug approved by major national regulatory agencies other than FDA. Both datasets were merged and submitted to the FROG2 server for 3D conformation generation (Miteva et al., 2010). In total, 3277 compounds 3D conformations were minimized using the MMFF94 force field and the steepest descent optimization algorithm implemented in the PyRx tool (Dallakyan & Olson, 2015). The minimized structures of all the ligands were converted into pdbqt format compatible for autodock.

Target preparation and virtual screening

The predicted models of RdRP and Nsp15 were prepared for docking using AutoDock Tool 1.5.6 (Morris et al., 2009). A well-established docking protocol was used in the present study as per previous studies (Barage et al., 2017; Meshram et al., 2020). The protonation states and Kollman united atom charges were assigned to receptor atoms using the Autodock tool 1.5.6. Then, both the structures were loaded in PyRx and converted into pdbqt format using make macromolecule option of PyRx (Dallakyan & Olson, 2015). The residues are known to play an essential role in catalytic activity, and substrate interaction of RdRp and Nsp15 were treated as flexible residues. The RdRp residues Val557, Asp623, Thr680 and Asn691, are the active site residues and are found conserved in the CoV family as per the previous literature (Kirchdoerfer & Ward, 2019; Jia & Gong, 2019). Whereas, the residues His 234, His 242, His 249 and Lys 289 contributes in the catalytic activity of Nsp15 (Joseph et al., 2007). These active site residues of RdRp and Nsp15 were treated as flexible, and the grid box was set to 50 × 52 × 62 and 74 × 66 × 50, respectively. The grid box was centered on selected flexible residues and leaving adequate space for translation (0.25 Å) and rotation (5 Å) of the ligand with 0.375 Å grid spacing. The 10 docking runs were performed for each ligand. The maximum number generation and individual population size were set to 27,000 and 150, respectively. All other parameters were kept to their default values. The generated docked conformation was ranked by predicted binding free energy. The top 10 ligand molecules were shortlisted based on their binding free energy. The binding mode and molecular interaction of the top 10 ligands were analyzed. Subsequently, the conformational clustering was performed on the top 2 ligand conformation with root mean square deviation (RMSD) tolerance 1.5 Å. The top two ligand molecules for RdRp and Nsp15 were selected form the highest populated cluster with the lowest binding energy. The obtained docked complexes were analyzed for its molecular interactions and subjected for MD simulation.

Molecular dynamics simulations

Molecular dynamics (MD) simulation was performed to study the dynamic behavior of protein ligand complex by applying CHARMm27 (Van Der Spoel et al., 2005) force field implemented in GROMACS 5.0.7 package. SwissParam was used to generate the topology files for ligands (Zoete et al., 2011). The system was solvated in the cubic box of TIP3P water model. Na+ counter ions were added to neutralize the negatively charged system. The system was energy minimized by steepest descent algorithm to eliminate the possible hard contacts from the initial structure until tolerance of 2500 kJ/mol. The energy-minimized system was subjected to 500 ps equilibration under NVT ensembles by using a V-rescale thermostat at 300 k temperature. In the next step, 1000 ps of NPT equilibration was done at 1 bar using Parrinello-Rahman barostat. Finally, the equilibrated systems were subjected to 30 ns production run at 300 K temperature and 1 bar pressure. LINCS algorithm was applied to constrain the bonds (Kumar et al., 2017) while; Particle Mesh Ewald was used to calculate the long-range electrostatic interactions with a cutoff distance of 1.2 nm (Steindl & Langer, 2004). The MD simulations were done in periodic boundary conditions. GROMACS tools, DS and visual molecular dynamics software were used to analyze the trajectories.

Binding free energy calculations

Molecular mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) was used to compute the binding free energy as described previously (Bavi et al., 2016, 2017; Kumar et al., 2019). For computing binding free energy, 30 snapshots of protein-ligand complex were chosen evenly from 0 to 30 ns of molecular dynamics simulation trajectories as per the previous studies (Bavi et al., 2019; Hou et al., 2011). MM/PBSA algorithm calculates different energy parameters by using the same snapshots (Hou et al., 2011; Vorontsov & Miyashita, 2011). The binding interaction between protein and ligand is calculated in three terms such as electrostatic contribution (ΔEele), van der Waals contribution (ΔEvdw) and the solvation contribution (ΔEsol).

Results and discussion

The emergence of highly pathogenic coronavirus, causing a severe acute respiratory infection is rapidly spreading throughout the world (Dhama et al., 2020; Sharma et al., 2020). Unfortunately, to date, no drug or vaccine is available to inhibit the pathogenic effect of coronavirus. Thus, the development of safe and effective drug molecules is important to control the disease and its spreading (Li & De Clercq, 2020). Therefore, non-structural viral proteins act as an attractive target for drug discovery (Zumla et al., 2016). The normal drug discovery is time-consuming (Stromgaard et al., 2016). Thus, drug repositioning, also known as repurposing strategies, is effective in identifying an antiviral agent to combat this disease (Li et al., 2016). The RdRp and Nsp15 are well studied non-structural proteins of viruses and acts as an attractive target for antiviral drug discovery. In the present study, we used these targets for the identification of potential inhibitors using a drug repurposing strategy.

Structural modeling of RdRp and Nsp15

The homology models of RdRp and Nsp15 were predicted using Modeller 9.23, as discussed in the method section. The stereochemical quality of the RdRp and Nsp15 model was performed using various tools. The Ramachandran plot analysis of the RdRP model shows that 92.1% of residues were in the favored region, and 7.6% residues were in the allowed region (Laskowski et al., 1993). Whereas 82.4% of residues of the Nsp15 model falls in the favored region and 15.5% residues are in the allowed region (Supporting Information Fig. S1). Also, the overall G-factor and main-chain parameters indicate a better quality model, as depicted in Supporting Information Fig. S2. Furthermore, the Z-score of the model is calculated using the ProSA-Web server. The Z-score is indicative of the model quality and measures the deviation of the total energy of structure (Wiederstein & Sippl, 2007). The Z-score of RdRp and Nsp15 was –12.79 and –7.6, respectively, as showed in Supporting Information Fig. S3. These results revealed that predicted models of RdRp and Nsp15 are in native conformation with well validated geometry (Laskowski et al., 1993; Wiederstein & Sippl, 2007). The structural superposition of the predicted model was also performed with a corresponding template structure using PDBefold. The structural alignment is a powerful tool for the identification of model residues occupying an equivalent geometric position similar to template structure (Fiser, 2010; Krissinel & Henrick, 2004). The overall structure of RdRp and Nsp15 are very similar to the corresponding templates with the backbone RMSD of 0.01 and 0.39, respectively (Krissinel & Henrick, 2004). Overall, structural analysis indicates that the predicted models of RdRp and Nsp15 are good in terms of stereochemical quality and well-validated geometry. The RdRp and Nsp15 models were further used for the docking study.

Figure 1.

Predicted 3D structures rendered as tube (A) RdRp structure composed of Finger (Green), Palm (Blue) and thumb (Orange) domain with interface domain (Grey) adjacent to polymerase domain. (B) Nsp15 structure with important structural element labelled and highlighted in different color.

Figure 2.

(A) Left panel: RdRp surface representation showing front view of long narrow active site tunnel (yellow) with buried inhibitor (Orange stick) molecule. Right panel: close-up view of bound conformation of top two ligand (upper-Chlorohexidine and lower-Alectinib) in active site tunnel. (B) Left Panel: RdRp tube representation showing binding mode of top two ligand (stick) at active site. Right panel: close-up view showing comparison of binding mode of top two ligands Chlorohexidine (red) and Alectinib (yellow).

Figure 3.

Molecular interactions between ligand and RdRp residues, (A) Left Panel: RdRp rendered as flat ribbon with key interacting residues in ball and stick (Green) and Chlorohexidine (Orange sticks) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of RdRp residues with Chlorohexidine (Orange). (B) Left Panel: RdRp rendered as flat ribbon with key interacting residues in ball and stick (Green) and Alectinib (magenta sticks) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of RdRp residues with Alectinib (magenta).

Virtual screening of compound against RdRp and Nsp15

The ligand dataset of 3277 compounds was subjected to virtual screening against RdRp and Nsp15 protein targets of COVID-19. The virtual screening was performed using the autodock tool implemented in PyRx (Dallakyan & Olson, 2015).

Interaction analysis and selection of top leads as RdRp inhibitors

The screening of the ligand dataset was performed against a structural model of RdRp of COVID-19. The RdRp catalyzes RNA synthesis from the RNA template crucial for viral genome replication and transcription process. Furthermore, the 3D architecture of RdRp was found conserved and similar to the polymerase family of viruses. It is composed of three sub-domains includes finger (residues S397-A581 and K621-G679), palm (residues T582-P620 and T680-Q815) and thumb (residues H816-E920; Figure 1(A); Bruenn2003).

The sequence and structural analysis of RdRp from RNA viruses revealed that the RdRp active site is surrounded by seven catalytic motifs A-G distributed within palm and finger sub-domain (Jia & Gong, 2019; Wu et al., 2010). The obtained docked conformations were ranked based on binding free energy and top 10 lowest binding energy ligand conformations were extracted. The binding free energy of the top 10 ligand conformation ranges from –10 to –9 Kcal/mol. This indicates that these compounds have a good binding affinity with RdRp. The resulting top 10 ligand conformations were evaluated for its binding mode and molecular interaction with active site residues of RdRP (Table 1). The significant differences have been observed in molecular interactions, binding mode and side-chain orientation between shortlisted top 10 compounds with RdRp (Table 1). Among them, the two compounds having the lowest binding energy observed in the active site of RdRp as compared with other lead molecules as given in Figure 2.

Table 1.

Molecular docking results of top10 ligand molecules of RdRp with binding free energy and their molecular interactions.

Compound structure and ZINC ID	Binding Free Energy	H-Bonding interactions (Conventional H Bond, Carbon-H Bond)	Hydrophobic interactions (Pi-Sigma, Pi-Alkyl, Alkyl)	Other interactions (Pi-Cation, Pi-Anion, Pi-Lone Pair, Salt Bridge, Halogen)
ZINC000014768621 (Chlorohexidine)	–10.11	ASP618,SER759,ASP760	ALA688,LYS621,CYS622	ASP684
ZINC000066166864 (Alectinib)	–9.6	TYR619,THR680,ARG624,THR540	TYR456,CYS622,ALA558	MET542,TYR456,ARG624
ZINC1481815 (Exjade)	–9.2	ARG553,ARG555,ARG624,LYS676,THR680	VAL667,MET542,ALA558,TYR456	ASP623,THR556
ZINC000011677837 (Eliquis)	–8.91	ARG553,THR556,ASP452,ALA554,ASP760	CYS622,ALA558	ARG624,ASP623,MET542
ZINC000043206370 (Niraparib)	–8.85	LYS676,THR680	MET542,MET626,ALA558	ARG624,ASP623,ASP760
ZINC000011679756 (Eltrombopag)	–8.74	ARG631,SER682,GLY679	MET542,MET626,CYS622,VAL667,LYS676,ALA558,THR687	ASP623
ZINC000001530886	–8.66	LYS621,ASP452,THR680	CYS622,TYR455	ARG553,ARG624
ZINC000012503187 (Conivaptan)	–8.55	ARG624,THR556,ASP623	CYS622,ALA690,MET626,ARG555,THR687	ARG553,ASP452
ZINC000003914810 (Degalol)	–8.48	ARG553,ARG555,ALA554,	LYS621,ARG624,TYR455	–
ZINC000019796087 (Nicardipine)	–8.47	LYS545,ARG555,THR680,ASP623	LEU758	ASP760

The first compound, namely, Chlorohexidine with binding free energy –10.11 Kcal/mol efficiently accommodated in the active site of RdRp and do molecular mimicry of nucleotides in terms of substructure interactions with RdRp residues (Figure 2(A)). Similarly, the second compound Alectinib with binding free energy –9.6 Kcal/mol was observed in the catalytic groove of RdRp (Table 1 and Figure 2). The resulting conformation of Chlorohexidine and Alcteinib was evaluated for molecular interactions with active site residues of RdRp.

The Chlorohexidine consists of two 4-chlorophenyl-guanido groups separated by hexane. The 4-chlorophenyl moiety interacts with Asp684 and Ala688 through hydrogen bonding and Pi-sigma interactions, respectively. While the guanido group forms direct hydrogen-bonding interactions with Ser759, Asp760 and Asp 618. Moreover, the Cys622 and Lys621 make pi-alkyl interactions with hexane and 4-chlorophenyl group, respectively, as shown in Figure 3(A) and Table 1.

Additionally, the RdRp residues Tyr619, Pro620, Asp623, Thr680, Ser682, Thr687 and Asn691 shows van der Waals interaction with Chlorohexidine (Figure 3(A)). The Chlorohexidine occupies active site pocket of RdRp and forms hydrogen bonding, Pi-alkyl and van der Waals interaction with RdRp residues. Further, the bound conformation of Alectinib was analyzed. The Alectinib consist of a substructure of four planer ring structure of benzo[b]carbazole-3-carbonitrile shows a slight deviation when compared with Chlorohexidine. Whereas, 4-pyrrolidine-1-yl-piperidine moiety occupies relevant space similar to Chlorohexidine (Figure 2). The variation in binding mode results in altered RdRp residue interaction with Alectinib. The 4-pyrrolidine-1-yl-piperidine moiety forms hydrogen bonding and Pi-alkyl interactions through Tyr619 and Cys622, respectively (Figure 3(B) and Table 1), whereas benzo[b]carbazole-3-carbonitrile core exists as planer conformation interacts with residues Thr680, Arg624, Thr540 through hydrogen bonding. Additionally, RdRp residues Tyr456, Ala558, Met542, Tyr456 and Arg624 show pi interaction with a substructure of Alectinib (Figure 3(B)). Finally, the RdRp complexed with Chlorohexidine and Alectinib were subjected to MD simulation.

Interaction analysis and selection of top leads as Nsp15 inhibitors

Nsp15 is another important non-structural protein specifically cleave RNA at 3′ of uridylates, forming a 2′, 3′-cyclic phosphate product (Fung & Liu, 2019; Bhardwaj et al., 2004). Recently, it has been reported that Nsp15 suppress the activation of host dsRNA sensor, including PKR kinase (Deng et al., 2017). The active site of Nsp15 is composed of three histidine at positions 234, 242 and 249 located on the activation loop and Lys 289 opposite to histidine on supporting loop (Figure 1(B)). The residue His 234, His242, His 249 and Lys 289 were treated as flexible during docking of Nsp15 (Joseph et al., 2007; Ricagno et al., 2006). The prepared ligand datasets were screened against the active site of Nsp15. The selection of the top 10 ligands from docked conformations was based on the binding free energy. The binding free energy of the top 10 ligand ranges from –10.95 to –9.4 Kcal/Mol, as depicted in Table 2. The binding modes of the top 10 compounds were evaluated. The top eight compounds occupied a similar binding site and last two compounds bound different sites than the catalytic region of Nsp15. Even though similar binding mode of top eight leads, but significant variation was observed in the orientation of ligand substructure concerning active site pocket. Subsequently, the top two compounds were evaluated for its binding mode, and molecular interaction. Naldemedine, the top lead molecule, accommodated in the active site pocket of Nsp15 (Figure 4). The octahydromethanobenzofuro isoquinoline moiety was observed near to the activation loop pocket while benzene-oxadiazol in the pocket formed by the supporting loop. The molecular interaction of Naldemedine sub-structures with the key residues of Nsp15 was located on the activation and supporting loop.

Table 2.

Molecular docking results of top10 ligand molecules of Nsp15 with binding free energy and their molecular interactions.

Compound structure and ZINC ID	Binding energy	H-Bonding interactions (Conventional H Bond, Carbon-H Bond)	Hydrophobic interactions (Pi-Sigma, Pi-Alkyl, Alkyl)	Other interactions (Pi-Cation, Pi-Anion, Pi-Lone Pair, Salt Bridge, Halogen)
ZINC000100378061 (Naldemedine)	–10.95	SER261,ARG257,LYS289	PHE263,LEU254	ASP282,ASP239
ZINC000052955754 (Ergotamine)	–10.12	CYS290,ARG257	LEU254,PHE263,LYS289	HIS249
ZINC000003780893 (Simvastatin)	–9.93	LEU254,ARG257	CYS292,VAL294,ILE280,CYS290,LYS289	–
ZINC000001530886	–9.93	ARG257,LYS289,CYS290,VAL291,GLY246	CYS292,LEU254,	HIS249
ZINC000003830500 (Cephalexin)	–9.9	PRO262	LYS289,PHE263,HIS249,ILE280,LEU254,CYS292,CYS290	ASP282
ZINC000000538621 (Olmesartan)	–9.7	TYR237,ARG257, LYS289	HIS249,PHE263,LEU254,ILE280	CYS292
ZINC000004215736 (Dihydrocodeine)	–9.67	LYS289,GLY238,CYS290	CYS292,LEU254,HIS249	–
ZINC000003813003 (Klonopin)	–9.62	ARG257,TYR237	CYS290,CYS292,LEU254,VAL294	ASP239
ZINC000255977092	–9.41	ASP239,GLY238	CYS290	–
ZINC000257362202	–9.4	VAL291,LYS289	CYS290,MET330,ILE235,	–

Figure 4.

(A) Left panel: Nsp15 surface representation showing front of view active site pocket (yellow) with bound inhibitor (Red stick) molecule. Right panel: close-up view of bound conformation of top two ligands (upper-Naldemedine and lower-Ergotamine) in active site pocket. (B) Left Panel: Nsp15 tube representation showing binding mode of top two ligand (stick) at active site. Right panel: close-up view showing comparison of binding mode of top two ligands Naldemedine (red) and Ergotamine (green).

Naldemedine formed hydrogen bonding interactions with Ser261, Arg257 and Lys289 of Nsp15 (Joseph et al., 2007; Ricagno et al., 2006). The residues Asp282, Phe263, Asp239, Leu254, Arg 257 and Lys289 show pi interaction with Naldemedine, as shown in Figure 5(A) and Table 1. The active site surrounding residues such as Gly238, Gy247 Ser287, Ser 288, Cys290, Ser293 shows van der Waals interaction with the ligand molecules (Figure 5(A) and Table 1; Joseph et al., 2007; Ricagno et al., 2006). The second lead molecule, Ergotamine bound conformation, showed slight variation in the side-chain orientation with respect to Naldemedine binding in the active site of Nsp15(Figure 4). The hexahydroindolo quinoline structure accommodates in the pocket of a supporting loop while hexahydrooxolopyrrolo-pyrazine moiety was observed in the space between activation and supporting loop (Figure 4(B)). The resulting bound conformation of Ergotamine makes hydrogen bonding interaction with residues Cys290 and Arg257. The residues Lys289, Phe263, His249, Leu254 and Lys289 show pi type interactions with Ergotamine. Furthermore, Val291, Cys292, Ile280, Pro262, shows van der Waals interactions with the sub-structure of Ergotamine (Figure 5(B) and Table 2). Subsequently, MD simulation has been performed on NSP15 complexed with Naldemedine and Ergotamine.

Figure 5.

Molecular interactions between ligand and Nsp15 residues, (A) Left Panel: Nsp15 rendered as flat ribbon with key interacting residues in ball and stick (Green) and Naldemedine (dark purple) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of Nsp15 residues with Naldemedine (magenta). (B) Left Panel: Nsp15 rendered as flat ribbon with key interacting residues in ball and stick (Green) and Ergotamine (tyrian purple) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of Nsp15 residues with Ergotamine (pink).

MD Simulation of selected docked complexes

MD simulation is essential to check stability of predicted interaction between identified ligand molecule and target proteins, selected docked complexes were subjected to atomistic simulation. We selected top 2 complexes for each target protein RdRp and Nsp15 and subjected to 30 ns MD simulation. We determined structural deviation of protein and ligand molecule during simulation time by analysis of each complex with set of standard parameters such as RMSD and root mean square fluctuations (RMSF). Moreover, binding free energy of ligands with protein molecule has been calculated throughout simulation using MM/PBSA method to validate binding affinity of ligand with protein target.

Conformational flexibility and stability analysis of RdRp during MD

To examine the changes conformational stability of protein-ligand complexes, we calculated RMSD of protein backbone and ligand with respect to initial structure for all simulations. The average RMSD value is relatively low (2.0 Å) throughout simulation time for protein backbone. The average RMSD value of ligand Alectinib and Chlorohexidine is 1.0 Å and 2.0 Å, respectively, as depicted in Figure 6(A).

Figure 6.

(A) RMSD plots of RdRp backbone and ligands during 30 ns MD simulations. (B) RMSF plot of RdRp complexed with Chlorohexidine and Alectinib for entire simulation time. (C) Estimated binding free energy using MM/PBSA method of RdRp complexes throughout simulation time.

It indicates that the bound conformation of Alectinib is more stable and not deviated from docked region whereas initially Chlorohexidine shows higher value due to conformation rearrangement in active site followed by stable for entire simulation time (Figure 6(A)). On the other hand, RMSD of protein backbone found stable throughout simulation time with small drift and plateaus (Figure 6(A)). The RMSF calculation is helpful to understand protein residue fluctuation during simulation time. There are no significant differences observed in residue fluctuation in both complexes. The higher residues fluctuation observed only at N-terminal and C-terminal regions and loop region (Figure 6(B)). Also, the residues at position 137, 432 and 590 shows major fluctuation due to located on extended loops when compared RdRp-Chlorohexidine with RdRp-Alectinib complex.

Analysis of molecular interaction of RdRp with ligands

The representative conformation from highest populated cluster was extracted for analysis of binding mode and molecular interaction. The binding mode and interaction of Chlorohexidine and alcteinib with RdRp compared with predicted in docking pose. The Chlorohexidine consists of two 4-chlorophenyl-guanido groups separated by hexane. The noticeable difference has been observed in binding mode and molecular interaction of Chlorohexidine in catalytic groove of RdRp. The hexane chain observed folded in docked pose while it is found as straight chain in simulated structure (Supporting Information Fig. S4A). This conformational rearrangement result into highest plateau observed in ligand RMSD plot at 1200 PS to 2000 PS followed by stable for entire simulation time (Figure 6(A)). The Chlorohexidine conformational variation leads to alteration in molecular interaction depicted in Figure 7(A). The RdRp residues Arg553, Tyr455 and Lys 500 shows hydrogen bonding and Pi-Alkyl interactions while Asn497, Lys545, Val557 and Ala558 make van der Waal interactions with Chlorohexidine (Figure 7(A)). The significant variation in interaction is due to altered conformation Chlorohexidine in catalytic groove of RdRp.

Figure 7.

Interactions between ligands and RdRp residues after MD simulation, (A) Left Panel: RdRp rendered as flat ribbon with key interacting residues in ball and stick (Green) and Chlorohexidine (Orange sticks) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of RdRp residues with Chlorohexidine (Orange). (B) Left Panel: RdRp rendered as flat ribbon with key interacting residues in ball and stick (Green) and Alectinib (magenta sticks) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of RdRp residues with Alectinib (magenta).

Similarly, the representative structure of RdRp-Alectinib extracted from simulation for interaction analysis. The binding mode of Alectinib in simulated structure found similar with predicted docked pose (Supporting Information Fig. S4). It is also evident from RMSD of Alectinib found stable for entire simulation time (Figure 6(A)). Furthermore, the RdRp residues Asn695 and Thr680 forms direct hydrogen bonding and Arg624 make Pi-Cation interactions which are preserved as observed in docked pose (Figure7(B)). The Cys622 and Ala558 forms alkyl interaction with Alectinib (Figure 7(B)).

Further, the binding free energy of both ligands has been calculating using MM/PBSA method. The snapshot of RdRp and ligand complexes has chosen evenly form MD trajectories for binding free energy evaluation. The favorable ΔG value has been observed for both complexes range from –20 to –140 kJ/mol as depicted in Figure 6(C). However, noticeable difference has been observed in binding free energy values of Chlorohexidine and Alectinib. The average binding free energy for RdRp-Chlorohexidine and RdRp-alcteinib is –25 kJ/mol and –101 kJ/mol, respectively (Figure 6(C)). The observed difference in binding free energy might be due to altered binding mode and conformational variation in Chlorohexidine during simulation time.

Taken together, the Alectinib shows excellent interaction in terms of hydrogen bonding, Pi stacking and van der Waals with RdRp residues. The docked pose and simulated structure show similar binding mode of Alectinib along with conserved molecular interactions. Alectinib binding mode and interactions found conserved and similar to nucleotide-binding in HCV. Structurally catalytic chamber of HCV and COVID-19 is highly conserved as reported in previous studies (Appleby et al., 2015; Elfiky & Azzam, 2020; Jia & Gong, 2019). Conversely, Chlorohexidine docked pose observed in the catalytic groove of RdRp but after simulation significant difference in Chlorohexidine orientation results into altered residue interaction. Thus, based on stability and molecular interaction and binding mode, Alectinib acts as lead molecule to design novel RdRp inhibitor to block the viral replication.

Conformational flexibility and stability analysis of NSP15 during MD

The RMSD value of protein backbone and ligand molecules has been calculated for entire simulation. The average RMSD value of protein backbone is 2.0 Å indicate stability of protein structure (Figure 8(A)). The ligand RMSD values are 1.0 Å and 1.5 Å for Naldemedine and Ergotamine, respectively (Figure 8(A)). The RMSF plot shows highest fluctuation of terminal regions and loop region of NSP15 (Figure 8(B)). The similar residue fluctuation observed in NSP15-Naldemedine and NSP15-Ergotamine complexes except position 150 to 153 (Figure 8(B)). This region shows more flexibility due to loop region.

Figure 8.

(A) RMSD plots of Nsp15 backbone and ligands during 30 ns MD simulations. (B) RMSF plot of Nsp15 complexed with Naldemedine and Ergotamine for entire simulation time. (C) Estimated binding free energy using MM/PBSA method of Nsp15 complexes throughout simulation time.

Analysis of molecular interaction of NSP15 with ligands

The conformational clustering was performed to select representative structure of NSP15 with ligand from simulation run. The NSP15-Naldemedine structure evaluated for binding mode of Naldemedine in simulated structure with reference to its predicted docked pose. The Naldemedine observed in active site pocket of NSP15 formed by activation and supporting loop as depicted in Supporting Information Fig. S4B. However, Naldemedine is slightly displaced and flipped from predicted docked pose results into altered residue interactions in simulated structure as compared with docked pose (Supporting Information Fig. S4B). The NSP15 residue Arg257 and His234 shows hydrogen bonding, pi-stacking and Pi-cation interaction while Cys290 and Lys289 forms Pi-sulfur and Pi-alkyl interactions, respectively (Figure 9(A)). In addition, Gly238, His249 and Tyr237 form van der Waals interactions with Naldemedine (Figure 9(A)). The differences in NSP15 residues interaction with Naldemedine in simulated structure and predicted docked pose is due to flipped orientation of Naldemedine (Supporting Information Fig. S4B and Figure 4(B)).

Figure 9.

Interactions between ligand and Nsp15 residues after MD simulation, (A) Left Panel: Nsp15 rendered as flat ribbon with key interacting residues in ball and stick (Green) and Naldemedine (dark purple) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of Nsp15 residues with Naldemedine (magenta). (B) Left Panel: Nsp15 rendered as flat ribbon with key interacting residues in ball and stick (Green) and Ergotamine (tyrian purple) of the predicted poses. Right panel: Schematic representation of 2D interaction plot of Nsp15 residues with Ergotamine (pink).

On the other hand, the molecular interactions in NSP15-Ergotamine complex the binding position of Ergotamine is similar to docked pose with slight variation. The Ergotamine hexahydroindolo quinoline moiety accommodated in pocket of a supporting loop similar to docked pose. Whereas, hexahydrooxolopyrrolo-pyrazine substructure is observed near activation loop (Supporting Information Fig. S4B) as it present between activation and supporting loop of NSP15 in docked pose (Figure 4(B)). The interaction of residues Arg257, Lys289, Cys290 and Phe263 with hexahydroindolo-quinoline moiety of Ergotamine similar to docked pose (Figure 9(B)). Also, the hexahydrooxolopyrrolo-pyrazine substructure shows Pi-T shaped interaction with residue His234 located on activation loop of NSP15 (Figure 9(B)). This interaction not observed in docked pose due to change in hexahydrooxolopyrrolo-pyrazine substructure orientation of Ergotamine in simulated structure (Figure 9(B)).

Furthermore, we calculated binding free energy of NSP15-Naldemedine and NSP15- Ergotamine using MM/PBSA method. The ΔG value ranges from –300 to –400 kJ/mol as depicted in Figure 8(C). The average binding free energy for Naldemedine and Ergotamine is –330 and –326 kJ/mol, respectively. Thus, no significant energy difference has been observed in Naldemedine and Ergotamine due to similar binding mode and molecular interaction with NSP15.

Overall, the identified top lead molecule Naldemedine and Ergotamine was observed in the catalytic groove of Nsp15. The Naldemedine and Ergotamine shows excellent interaction with active site residues along with the surrounding residues of Nsp15. In addition, only slight variation in binding free energy has been observed between both ligand molecules. Thus, substructure properties and interaction specificities of Naldemedine and Ergotamine useful to design novel NSP15 inhibitor which will interfere with RNA cleavage and trigger an elevated interferon response and PKR dependent apoptosis. Interestingly, most of identified interacting residues contribute to the interaction with another monomer of Nsp15. Thus, identified ligands may play a dual role as an inhibitor and also interfere with its interaction with other monomers.

In summary, the two non-structural proteins, the RdRp and Nsp15 play a crucial role in viral replication and RNA processing used in our study. We have performed homology modeling, virtual screening technique to selected potential lead molecule inhibitors. Subsequently, selected top complexes of RdRp and Nsp15 with identified lead molecules subjected to MD simulations. The promising compound such as Alectinib for RdRp while Naldemedine and Ergotamine for NSP15 were successfully identified which will act as a lead inhibitor for RdRp and Nsp15. The identified lead molecules lodged in the catalytic groove of enzymes and form hydrogen bonding as well as pi-pi, pi-sigma, pi-alkyl and vdW interactions with non-polar residues of the targets. As discussed earlier, most of the identified interactions are conserved and similar to the previous substrate and inhibitor interaction with these targets. Our study gives an important structural explanation of the promiscuity and specificity of ligand substructure necessary to inhibit the activity of RdRp and Nsp15. The identified lead inhibitors open up possibilities for designing potent and selective inhibitor against RdRP and Nsp15.

Conclusion

Considering the global health emergency and lack of effective treatment available, drug repurposing is the only option of choice. Under the current adverse situation, we employed a series of in silico drug repurposing screening campaigns against important non-structural proteins RdRp and Nsp15. The Alectinib, Naldemedine and Ergotamine, an FDA approved drug, shows high affinity and crucial molecular interaction with key residues involved in catalysis and interaction of substrate in RdRp and Nsp15. Our MD simulation analysis showed excellent stability of these drugs with active site of target proteins. Thus, considering the severity of coronavirus, the current study is in-line with the concept of finding the new inhibitors against the replication-associated proteins of the coronavirus to expedite the process of drug discovery. Finally, we propose that the identified inhibitors represent a novel lead molecule to design a more effective inhibitor to stop the progress of pathogen.

SHB and A. Karthic are thankful to Amity Institute of Biotechnology, Amity University, Mumbai for infrastructure support. We would like to thank Amity University administrative team for their support to carry out this research in a safe and healthy environment during CoVID-19 pandemic. Rohit Bavi is thankful to the University Grants Commission for Dr D. S. Kothari postdoctoral fellowship.

Click here for additional data file.