Exploring RdRp-remdesivir interactions to screen RdRp inhibitors for the management of novel coronavirus 2019-nCoV.

A novel coronavirus recently identified in Wuhan, China (2019-nCoV) has resulted in an increasing number of patients globally, and has become a highly lethal pathogenic member of the coronavirus family affecting humans. 2019-nCoV has established itself as one of the most threatening pandemics that human beings have faced, and therefore analysis and evaluation of all possible responses against infection is required. One such strategy includes utilizing the knowledge gained from the SARS and MERS outbreaks regarding existing antivirals. Indicating a potential for success, one of the drugs, remdesivir, under repurposing studies, has shown positive results in initial clinical studies. Therefore, in the current work, the authors have attempted to utilize the remdesivir-RdRp complex - RdRp (RNA-dependent RNA polymerase) being the putative target for remdesivir - to screen a library of the already reported RdRp inhibitor database. Further clustering on the basis of structural features and scoring refinement was performed to filter out false positive hits. Finally, molecular dynamics simulation was carried out to validate the identification of hits as RdRp inhibitors against novel coronavirus 2019-nCoV. The results yielded two putative hits which can inhibit RdRp with better potency than remdesivir, subject to further biological evaluation.

Introduction

SARS-CoV-2 is an enveloped, positive-sense, single-stranded RNA β-coronavirus similar to the Severe Acute Respiratory Syndrome (SARS) and Middle East Respiratory Syndrome (MERS) viruses [1]. Potential antiviral targets encoded by the viral genome include non-structural proteins (e.g. 3-chymotrypsin-like protease, papain-like protease, RNA-dependent RNA polymerase (RdRp) and its helicase), structural proteins (e.g. the capsid spike glycoprotein) and accessory proteins [2]. Kaletra (lopinavir/ritonavir) is thought to inhibit the 3-chymotrypsin-like protease of the SARS and MERS coronaviruses and was associated with improved clinical outcomes in a trial against SARS [3]. Ascletis, a biotechnology company, also reported that a patient with COVID-19 improved rapidly when administered with an HIV protease inhibitor combination [4].

From the start of the COVID-19 outbreak, medical practitioners have followed China’s guidelines set up in January and treated hospitalized patients with α-interferon combined with the repurposed drug Kaletra, an approved cocktail of the HIV protease inhibitors ritonavir and lopinavir [5]. The World Health Organization has noted that this combination could provide some clinical benefits. Kaletra is also being tested in other combinations, for instance, with the guanosine analogue and RNA synthesis inhibitor ribavirin, with reverse transcriptase inhibitors (emtricitabine/tenofoviralafenamide fumarate) or with membrane fusion inhibitor umifenovir. Umifenovir is also in trials as a single agent [6].

Many research labs have been working on repurposing studies using existing drugs, approved for other viruses, as treatments in the coronavirus outbreak. Pharmaceutical companies such as Ascletis Pharma are also testing two HIV protease inhibitors (ritonavir and ASC09) to treat COVID-19 [7], while Gilead Sciences is investigating remdesivir (GS-5734), a broad-spectrum antiviral originally developed to treat Ebola virus and then dropped, which has shown significant results against coronavirus infection [8]. Academic research groups are also focusing on utilizing reported inhibitors for virtual screening analysis against different viral targets such as 3-C like protease [9,10], papain-like protease (PLpro) [11], etc. However, recently researchers have favoured targeting a virus-specific protein such as RdRp, noting that coronaviruses do not contain or use a reverse transcriptase [12]. Supporting this hypothesis, remdesivir, a nucleotide analogue antiviral that blocks the RNA polymerase of the Ebola virus and so prevents replication, has also shown positive signs of being effective against SARS-CoV-2 [13]. A recent study attempted to study possible mechanism of anti-RdRp drugs (e.g. favipiravir, sofosbuvir, ribavirin and galidesivir) that may inhibit the SARS-CoV-2 RdRp. They reported that sofosbuvir, ribavirin, galidesivir, remdesivir, favipiravir, cefuroxime, tenofovir, and hydroxychloroquine can bind to the RdRp active site tightly, and are thought to be good candidates for clinical trials [14].

In vitro studies have shown remdesivir to be an active agent against a clinical isolate of SARS-CoV-2 [15]. Experimental data in animal models with the related MERS virus also showed that the drug was better than a combination of lopinavir/ritonavir and interferon beta in improving lung function [16]. Several patients with confirmed COVID-19 have been reported to improve after being treated for 1 day with remdesivir, although this could not be directly attributed to the drug’s effect [17]. Since then, remdesivir has been shown to reduce the severity of disease, virus replication and damage to the lungs in a non-human primate model of MERS.

The RdRp complex, a protein complex responsible for viral RNA transcription and replication, represents a primary target for the antiviral drug development [18]. Importantly, in April 2020, the high-resolution crystal structure of RdRp was released in a paper in Science, describing the structure of the polymerase protein [19]. Therefore, in the current study, we tried to utilize the information gained via the interaction of remdesivir and RdRp complex, obtained via molecular docking analysis of the same, to further screen a library of other previously reported RdRp inhibitors, to identify possible inhibitors more potent than remdesivir which could be then utilized for the management of COVID-19. The protocol used in this in silico analysis involved analysis of the binding pocket in the 3D structure of RdRp of novel corona virus, followed by molecular docking analysis of remdesivir into the identified binding pocket, which led to the identification of key residues involved in the binding of inhibitors. Further considering these interactions, a library of already reported RdRp inhibitors was subjected to structure-based drug design analysis via molecular docking. The hits, selected on the basis of maintaining conserved interactions, were then subjected to clustering analysis to filter out structurally similar hits. Followed by scoring refinement using AutoDock tools, the hits maintaining significantly high binding affinity were then subjected to molecular dynamics simulation analysis, leading to identification of two hits showing higher putative potency than remdesivir against RdRp, which, subject to biological evaluation, can be utilized for the management of COVID-19.

Material and methods

Database preparation and optimization

For the current study, a database was created by downloading the list of RNA-dependent RNA polymerase (RdRp; EC 2.7.7.48) inhibitors reported to date, from BRENDA, a comprehensive enzyme information system database [20]. Overall, the database consisted of more than 350 molecules reported as RdRp inhibitors. The structures were prepared using MOE (Molecular Operating Environment) as follows: (a) explicit hydrogen atoms were added; (b) partial charges were added to the structures; (c) energy minimization was then carried out.

Protein preparation

The crystal structure of inhibitor-free RdRp of SARS-CoV-2 was retrieved from the RCSB Protein Data Bank (PDB ID: 7BV2) [21]. It consisted of multiple chains; chain A, containing the binding pocket, was chosen for this study. The protein target was prepared using AutoDock 4.25 [22]. Briefly, water molecules outside the binding pocket and sulphate ions were removed, and hydrogen atoms were added using the ADT module implemented in AutoDock. Atom type was modified into ADT type and charges were adjusted using the Gasteiger charges module for proteins implemented in AutoDock.

Docking-based virtual screening (MOE docking procedure)

The MOE docking protocol was applied to the receptor model, which was optimized by selecting AMBER99 as force field and fixing hydrogens and charges [23]. Then, the dockable space was set by selecting the specific binding cavity and putting dummy atoms inside it. Keeping this selection as docking site, the previously prepared database of RdRp inhibitors was used as Ligand. The Triangle Matcher protocol was used as placement feature for the compounds in the binding pocket. Finally, the London dG parameter was utilized as scoring function to assess the docking results. Remdesivir was also docked following the same protocol to obtain a cut-off score to filter the obtained hits.

Clustering analysis

Clustering is an invaluable cheminformatics technique for subdividing a typically large compound collection into small groups of similar compounds. Molecular fingerprints based on structural fragments are utilized to cluster structurally similar molecules. In the current work, MACCS keys fingerprints [24] were calculated to assess structure similarity, and on the basis of these fingerprints clustering of the hits obtained from the docking analysis was performed.

Scoring refinement (AutoDock docking refinement)

To refine the hit search, an optimized docking procedure was used. Briefly, docking was performed with AutoDock version 4.2, using the empirical free energy function and the Lamarckian protocol [22]. The atomic charges for the protein were assigned using the Gasteiger–Marsili method. Mass-centred grid maps were generated with 80 grid points for every direction and with 0.375 Å spacing by the AutoGrid program. Ten independent docking runs were carried out for each ligand. The docking results were analysed for the binding mode and conserved interactions such as hydrogen bond, hydrophobic and п-п interactions between the hits and the active site of the protein. The common interactions in all the complexes were scrutinized. Finally, Remdesivir was also docked following the same protocol and the obtained pose was utilized to screen the obtained hits.

Molecular dynamics simulation analysis

Finally, to validate the results obtained after following the in silico approach, molecular dynamics simulations was performed on the selected best binding pose of the hits docked in the catalytic domain of the target protein [25]. The simulations were performed to analyse the stability of the protein–hit complex and to study the most stable interactions which are retained after the simulation time period, by observing its 3D-interaction diagram [26]. This analysis was carried out using MOE software [27] with AMBER99 force field. Partial charges were calculated and energy minimizations were performed. The protocol of molecular dynamics simulations involved solvation of the protein–ligand complex using SPC water in a spherical box. Molecular dynamics simulation was carried for the time duration of 20 ns. The NPT statistical ensemble was conducted at 310 K with constant pressure. The Nose–Hoover–Anderson equations were used to solve the equations of motion. The simulation was carried out using the NPT ensemble and a time step of 0.002 fs and the coordinate data were stored in the database. The temperature was fixed at 310 K using the Nose–Hoover method as the thermostat and pressure of 1 bar using Berendsen barostat. The root-mean-square-deviation (RMSD) value was calculated to determine the stability of the complex after the time period of 20 ns [28].

Results and discussion

Following the in silico protocol, first of all binding cavity analysis of the 3D structure the RdRp of SARS-CoV-2 co-crystallized with remdesivir (PDB ID: 7BV2; 2.5 Å) available at RCSB was performed. Following the identification of the binding cavity, the first step was re-docking analysis of remdesivir – a nucleotide analogue reported as successful for the management SARS-CoV-2 – to validate and verify the key amino acid residues involved in its binding in the catalytic domain of RdRp. After re-docking analysis, remdesivir retained its co-crystallized pose with the RMSD < 1. This analysis also disclosed that remdesivir maintains key H-bond interactions with binding pocket amino acid residues such as Arg553, Arg555, Thr556 and Asn691, which are essential for the inhibitory potential of small molecule heterocycles (Figure 1). Remdesivir also maintained key H-bond interaction with U-20 nucleotide of RNA, highlighting the crucial parameter required for RdRp inhibitors, i.e. to maintain stabilizing interactions with both RNA and the catalytic domain of RdRp. This information was crucial, and was utilized as a screening parameter in further analysis. In the next step, a database of RdRp inhibitors previously reported in the literature provided by the BRENDA library was prepared. This database of RdRp inhibitors was then subjected to docking-based virtual screening utilizing the catalytic domain of RdRp. Special attention was given to the binding mode and binding affinity of the hits obtained after screening. In this primary screening, the docking score obtained for remdesivir was considered as the cut-off value. Docking results (Table 1) revealed that 42 RdRp inhibitors, out of the complete database, possessed a higher docking score than remdesivir. These hits occupied well within the binding pocket of RdRp, maintaining key interactions essential for the inhibitory potential, as suggested by the binding of remdesivir. Although the selected hits showed comparably high docking scores, there was variation in the binding mode due to the shape and size of the hits. Therefore, the hits obtained from this preliminary screening were then subjected to clustering analysis on the basis of their structural attributes to remove the structurally similar hits.

Figure 1.

3D interaction diagram of remdesivir in the binding pocket of RdRp

Table 1.

Results after MOE-based molecular docking, MACCS key fingerprint analysis-based clustering and scoring refinement analysis via AutoDock

Compound ID	Docking score (London dG)	Cluster ID	Docking score (kcal/mol)	Predicted activity	H-bond interactions
IN-1	−14.8021	1	−3.29	3.9 mM	Arg553, Lys621
IN-2	−15.3838	2	−4.25	764 µM	Lys551, Lys621,Ser795
IN-3	−14.7440	3	−6.64	13.6 µM	-
IN-4	−15.0525	4	−7.19	5.36 µM	-
IN-5	−16.0244	5	−7.59	2.74 µM	Arg553, Thr556
IN-6	−14.5860	6	−6.53	16.32 µM	Arg555, Ala550, U10 (RNA)
IN-7	−15.0054	7	−4.21	802 µM	Lys621
IN-8	−14.9487	8	−4.34	733 µM	Lys591
IN-9	−15.1409	9	*	*	*
IN-10	−14.1418	10	−6.38	21.1 µM	U10 (RNA)
IN-11	−16.8280	9	−4.76	375 µM	Arg553, Thr556
IN-12	−18.4770	12	−6.19	28.8 µM	Lys551, Lys621, Ser795
IN-13	−17.2873	12	*	*	*
IN-14	−18.1005	14	−5.28	108 µM	Lys551, Lys621
IN-15	−15.5923	15	−7.4	3.77 µM	-
IN-16	−14.0843	16	−4.56	324 µM	Arg553, Thr556
IN-17	−14.7859	17	−8.61	491 nM	Lys551,U20 (RNA)
IN-18	−17.3240	18	−6.1	33.9 µM	Arg553, Thr556, Ser682
IN-19	−18.0230	19	−7.4	3.76 µM	Arg553, Thr556
IN-20	−14.1306	20	−4.97	226 µM	Arg555,U10 (RNA)
IN-21	−14.3338	21	−4.85	197 µM	Arg555
IN-22	−18.8620	22	−2.74	9.79 mM	Ser549, Lys551, Arg553, Arg555
IN-23	−14.0651	23	−5.76	59.8 µM	U20 (RNA)
IN-24	−14.1158	24	*	*	*
IN-25	−16.1310	25	−4.14	923 µM	Lys551, Ser795
IN-26	−19.4613	26	−3.84	1.34 mM	-
IN-27	−14.3212	27	*	*	*
IN-28	−14.2804	28	*	*	*
IN-29	−14.3415	29	−4.73	341 µM	-
IN-30	−14.8812	24	−4.43	697 µM	Ser549
IN-31	−14.7497	24	*	*	*
IN-32	−14.7420	24	*	*	*
IN-33	−15.0047	33	−4.57	525 µM	Arg555
IN-34	−14.1123	34	−5.21	111 µM	-
IN-35	−14.0488	24	*	*	*
IN-36	−14.0786	24	*	*	*
IN-37	−14.0306	24	*	*	*
IN-38	−14.2075	24	*	*	*
IN-39	−15.1382	39	−5.11	123 µM	-
IN-40	−15.8917	27	−6.24	26.5 µM	Arg553,Ser549
IN-41	−14.2085	28	*	*	*
IN-42	−16.6071	28	−5.31	129 µM	U20 (RNA)
Remdesivir	−14.0953	-	−5.97	42 µM	Arg553, Arg555,U20 (RNA)

*Docking Score (kcal/mol), Predicted activity and H-bond interactions, determined during scoring refinement analysis, were only calculated for hits obtained after cluster analysis

For clustering analysis, MACCS key fingerprints for each hit were calculated and then all 42 hits were clustered. The top hit on the basis of the docking score from each cluster was selected, which resulted in 30 structurally diverse hits representing each cluster (Table 1). To validate these results, the 30 structurally diverse hits were selected for AutoDock-based scoring refinement. Through molecular docking studies, it was found that 11 molecules exhibited excellent binding energy scores (AutoDock score), higher than the standard remdesivir, and showed higher predicted inhibitory potential (Table 1). Interestingly enough, three hits, IN-3, IN-4 and IN-15 did show higher binding affinity than remdesivir but did not maintain the key H-bond interactions with crucial amino acid residues of the binding pocket, establishing them as false positives. Another key point observed during the analysis was that, out of all the top hits, six hits, IN-5, IN-10, IN-12, IN-18, IN-19 and IN-40, did show more binding affinity than remdesivir and did maintain key H-bond interaction with binding pocket amino acid residues (Lys551, Arg553, Arg555 and Thr556), but they did not maintain the interaction with RNA. Among all, two ligands, IN-6 and IN-17, maintained the key conserved H-bond interaction with Lys551, Arg553, Arg555 and Thr556, similar to that of remdesivir, which is expected to be an essential requirement for RdRp inhibitory activity. Moreover, these two hits also maintained the crucial H-bond interaction with U20 of the RNA, confirming the potential of these hits to maintain stabilizing complexes. Interestingly, IN-17 exhibited the highest binding energy among these top hits (−8.61 kcal/mol), a lot higher than remdesivir (−5.97 kcal/mol). Thus these top two can be considered as potent putative RdRp inhibitors (Figure 2). The 3D interaction diagram of the top hits identified as novel hits for putative SARS-CoV-2 inhibitory potential is shown in Figure 2.

Figure 2.

3D interaction diagrams of top hits after AutoDock scoring refinement. (a) IN-17; (b) IN-6

In the last step, molecular dynamics simulations were run to analyse and validate the interactions, stability and binding of the retrieved hits with their proteins. They were run to analyse and validate the interactions, stability and binding mode of the hits within the catalytic domain of RdRp. Thus, previously docked complexes of two putative inhibitors with RdRp were considered for the molecular dynamics simulations. These complexes were exposed to molecular dynamics simulations for a time period of 20 ns. Careful investigation of interaction diagrams revealed that all three molecules retained the conserved interactions in the binding pocket of the protein which are expected to be essential for inhibitory activity, thus justifying the claim of putative inhibitors. Furthermore, complexes were analysed for their stability by calculating the RMSD of the hits in the binding pocket to validate the obtained results (Figure 3). In case of the best hit, IN-17, the values were found in the range of 1.3–1.7 Å for the ligand. Initial variations in RMSD can be justified by the fact that slight adjustment occurs at the beginning of a simulation study. However, RMSD and interactions within the active site of both the protein molecules conclude that the complex was fairly stable after initial adjustment. The best retrieved molecule (IN-17) possessed conserved interactions with the catalytic domain amino acid including Arg553, Arg555, Thr556, Lys551, Lys621 and Ser795of RdRp and,U10/U20 of RNA in complex with RdRp, suggesting potential inhibitory potential of the identified hit.

Figure 3.

3D interaction diagram of best hit (IN-17) in the catalytic domain of RdRp, along with the RMSD plot

Finally, the ADMET evaluation [29,30] of IN-17 was also performed using a freely available web interface SwissADME (http://www.swissadme.ch/) [31], as given in Table 2. Results suggest that the molecule possesses the required number of H-bond acceptors and H-bond donors, and therefore fulfils all the drug-like criteria as per Lipinski rule, except for the molecular weight, which is 517.57, slightly higher than 500. Considering the lipophilicity of the molecule, SwissADME provides multiple variants of log P utilizing different methodologies; however, the Consensus log Po/w, the average of the log P values calculated by different methods, was found to be 5.51, which is still slightly high for drug-like molecules. Similarly, the solubility parameter also suggested that the molecule is poorly water soluble. However, both these issues can be managed via formulation-based optimizations. One approach, if the biological validation confirms the potency of the molecule, could be developing a prodrug of the hit molecule which would not alter the structural integrity of the lead but can definitely improve the physicochemical properties. Also considering the acidic functional group in the side chain, the molecule is suitable for the development hydrolysable prodrugs which could manage the solubility and permeability criteria of the molecule.

Table 2.

Various predicted ADME properties of IN-17

S.No.		Properties	Values
1.	Physicochemical	Molecular weight	517.57 g/mol
2.		Num. H-bond acceptors	6
3.		Num. H-bond donors	2
4.		Molar Refractivity	150.74
5.		TPSA	105.32 Å2
6.	Lipophilicity	log P (iLOGP)	4.03
7.		log P (XLOGP3)	6.49
8.		log P (WLOGP)	7.36
9.		log P (MLOGP)	3.57
10.		log P (SILICOS-IT)	6.10
11.		Consensus Log Po/w	5.51
12.	Water Solubility	log S (ESOL)	−7.19 (Poorly soluble)
13.		log S (Ali)	−8.50 (Poorly soluble)
14.		log S (SILICOS-IT)	−10.71 (Poorly soluble)
15.	Pharmacokinetics	GI absorption	Low
16.		BBB permeant	No
17.		P-gp substrate	No
18.		CYP1A2 inhibitor	No
19.		CYP2C19 inhibitor	Yes
20.		CYP2C9 inhibitor	No
21.		CYP2D6 inhibitor	Yes
22.		CYP3A4 inhibitor	No
23.		log Kp (skin permeation)	−4.85 cm/s
24.	Druglikeness	Lipinski	1 violation: MW > 500
25.		Ghose	3 violations: MW > 480, WLOGP > 5.6, MR > 130
26.		Veber	Yes
27.		Egan	1 violation: WLOGP > 5.88
28.		Muegge	1 violation: XLOGP3 > 5
29.		Bioavailability Score	0.56
30.	Medicinal Chemistry	PAINS	0 alert
31.		Brenk	0 alert
32.		Leadlikeness	2 violations: MW > 350, XLOGP3 > 3.5
33.		Synthetic accessibility	4.06

Further, the pharmacokinetic predictions regarding P-gp substrate and blood–brain barrier permeant was found to be negative. Also the molecule was found to fully comply with Veber rules of drug-likeness. Finally, the hit was found not to be a PAIN molecule, establishing the importance of further exploration.

Conclusion

SARS-CoV-2 has been wreaking ongoing global havoc. Out of all potential targets, researchers have favoured targeting a virus-specific protein such as the RdRp. Therefore, in the present study, we have performed an in silico analysis to identify previously reported RdRp inhibitors as potential agents to inhibit RdRp of the SARS-CoV-2. Initial analysis of the binding pocket of RdRp and interaction pattern of remdesivir with this pocket laid grounds for the detailed analysis. This was followed by a structure-based virtual screening protocol to screen a library of already reported RdRp inhibitors to determine their potential in the management of SARS-CoV-2. Overall, the analysis disclosed two putative hits which could possibly inhibit RdRp at around 1 µM concentration. However, this is simply an in silico analysis, and even though virtual screening makes it possible to discover molecules relatively quickly, these compounds still need to be experimentally tested.

Click here for additional data file.