Novel adenosine derivatives against SARS-CoV-2 RNA-dependent RNA polymerase: an in silico perspective.

BACKGROUND: SARS-CoV-2 is a newly emerged human coronavirus that severely affected human health and the economy. The viral RNA-dependent RNA polymerase (RdRp) is a crucial protein target to stop virus replication. The adenosine derivative, remdesivir, was authorized for emergency use 10 months ago by the United States FDA against COVID-19 despite its doubtful efficacy against SARS-CoV-2.
METHODS: A dozen modifications based on remdesivir are tested against SARS-CoV-2 RdRp using combined molecular docking and dynamics simulation in this work.
RESULTS: The results reveal a better binding affinity of 11 modifications compared to remdesivir. Compounds 8, 9, 10, and 11 show the best binding affinities against SARS-CoV-2 RdRp conformations gathered during 100 ns of the Molecular Dynamics Simulation (MDS) run (- 8.13 ± 0.45 kcal/mol, - 8.09 ± 0.67 kcal/mol, - 8.09 ± 0.64 kcal/mol, and - 8.07 ± 0.73 kcal/mol, respectively).
CONCLUSIONS: The present study suggests these four compounds as potential SARS-CoV-2 RdRp inhibitors, which need to be validated experimentally.

Introduction

SARS-CoV-2 is a coronavirus identified as the cause of severe acute respiratory syndrome that spread worldwide, leading to an outbreak of a global pandemic. The virus emerged at the end of 2019 in Wuhan, China, infecting more than 170 million people with mortality exceeding 3.5 million till now [1-3]. The virus also harmed the worldwide economy [4]. With all these massive losses in lives and money, there is a desperate need for a new drug that can eradicate this life-threatening severe virus. Since the beginning of the pandemic, many researchers have focused on drug repositioning, where approved drugs or active molecules that were used in other medical indications can be used to treat COVID-19 patients [5]. Although many medications were mistakenly believed to help overcome the disease as hydroxychloroquine [6], others, as remdesivir, showed an advantage in decreasing the severity of the disease [7, 8]. Remdesivir is an adenosine derivative that was previously used against SARS-CoV, a human coronavirus similar to COVID-19. Remdesivir acts as an adenosine nucleoside triphosphate analog that inhibits the RNA-dependent RNA polymerase enzyme (RdRp), which is a crucial enzyme for virus replication and proliferation [9, 10]. Remdesivir is a prodrug, where it is converted by the action of esterases and phosphoramidases to GS-441524, which is then phosphorylated to the active triphosphate form (Fig. 1). The potency of remdesivir is found not enough to treat all COVID-19 patients, and consequently, drug optimization is required.

Fig. 1

A diagram showing the activation of Remdesivir drug

In this study, our strategy is carrying out reasonable modifications on the general features of the remdesivir scaffold and modeling the effect of such alterations on the COVID-19 RdRp binding. The modifications are categorized into two groups: sugar modifications and nucleobase modifications. Sugar modifications include studying the impact of changing the stereogenic configuration at positions 2 and 3 of the ribosyl ring and converting it to arabinose and xylose, respectively (structures 2 and 3, Fig. 2) [11]. Substitution of the hydroxyl group at the ribosyl 5-position with an amino or substituted amino groups can also mimic the triphosphate’s binding (structures 4–8, Fig. 2). On the other hand, the nucleobase modifications of remdesivir are based mainly on substituting the pyrrolotriazine ring with the isosteric ring pyrazolopyrimidine. Other changes involve substitutions at the 4-amino group with an attached catechol ring (structures 9–10, Fig. 2) [12, 13]. The effect of connecting a hydrazine group or a benzylthio group at 6′ position is also explored (structures 11–12, Fig. 2) [14]. The rationale beyond this approach is to identify optimized candidates with improved selectivity.

Fig. 2

A schematic diagram showing both the sugar modifications and nucleobase modifications to the remdesivir drug

Materials and methods

Structural retrieval

The ligands’ structures are sketched and prepared to be in their active (triphosphate) form with the aid of Gaussian 09 software [15]. After that, geometry optimization is performed with the help of the classical Universal Force Field (UFF) [15, 16]. This is followed by the semi-empirical force field Parameterization Method 6 (PM6) [17]. Remdesivir and cinnamaldehyde are used as positive and negative controls, respectively.

The structure of the protein SARS-CoV-2 RNA-dependent RNA polymerase (RdRp) (PDB ID: 7BTF) is downloaded from the Protein Data Bank [18]. PyMOL software is utilized to prepare the protein structure for docking [19]. Water molecules and ligands are removed while missing hydrogen atoms are added to the structures. AutoDock Tools software is used to prepare the ligands and the protein structures for docking (adding charges and missing hydrogen atoms) [20].

Molecular dynamics simulation

SARS-CoV-2 RdRp is subjected to a round of molecular dynamics simulation run for a period of 100 ns at Normal Volume and Temperature (NVT) ensemble [21]. This step is started after a round of equilibration of the protein solution (NaCl of 154 mM) for 1.1 ns in the NPT ensemble. After volume adjustment, NAMD software is used to run the MDS production run, while the graphical software VMD is used to prepare the system and analyze the data [22, 23]. CHARMM 36 force field is used during the simulation, while the TIP3P water model is used to represent the water molecules [24, 25]. Temperature is raised gradually during the equilibration period to the physiological temperature of 37 °C. After the MDS, cluster analysis for the trajectories is performed utilizing the Chimera software tools default parameters [26]. Eight different clusters are found in the 100 ns MDS trajectories; each cluster’s representative structure is used to test the ligands’ binding affinities.

Molecular docking

AutoDock Vina software is utilized in the docking experiments using eight parallel processors [27]. The two consecutive active site aspartates (D760 and D761) are treated as being flexible during all the docking experiments. The grid box used in the docking experiment is 20 × 20 × 20 Å3 in size and centered at the active site aspartates (D760 and D761). Protein–Ligand Interaction Profiler (PLIP) software is utilized to analyze the docking complexes [28].

CHARMM-GUI webserver [29, 30] was used to prepare the necessary files for each cluster for MDS using NAMD. The same conditions used in the primary MDS calculations are set for the best four complexes of the SARS-CoV-2 RdRp and the compounds 8, 9, 10, and 11. After that, a production step of 20 ns was performed in the constant number of atoms, constant volume, and constant temperature (NVT) ensemble. To calculate the binding free energy for each complex, Molecular Mechanics Generalized Born Surface Area (MM-GBSA) implemented in MMPBSA.py in Ambertools 20 was utilized [31, 32]. All frames (200 frames) were used to calculate the free binding energy, saltcon variable was set to 0.154 M, and the generalized born (igb) method was set to five. In addition, decomposition of the free energy was obtained to determine the contribution of each amino acid toward the binding for each complex.

Results

Binding affinities of the adenosine derivatives against SARS-CoV-2 RdRp

Figure 3 shows the average binding energies (in kcal/mol) for the 12 novels adenosine derivatives in addition to the positive and negative controls, remdesivir (red column), and cinnamaldehyde (orange column). The docking experiments are performed on eight different conformations for SARS-CoV-2 RdRp gathered during the 100 ns MDS. Each conformation resembles a cluster from the trajectories. Molecular dynamics is performed to equilibrate the protein system and ensure the protein’s different possible conformations.

Fig. 3

The average binding affinities (in kcal/mol) for the 12 adenosine derivatives and the positive and negative controls (red and orange columns, respectively). Error bars resemble the standard deviation (SD) of the mean. Eight different conformations of the protein are used in the docking experiments after MDS cluster analysis. The best compounds are depicted in green, while blue columns show the compounds that give better average binding affinities than Remdesivir

The PLIP software is used to check the established interactions to analyze the compounds’ binding modes against SARS-CoV-2 RdRp. Table 1 shows the interactions established for each ligand after docking to the active site of SARS-CoV-2 RdRp. Based on the binding energy values, the best four compounds are 8, 9, 10, and 11 (Fig. 2), where the blue sticks represent the interacting residues, while the yellow sticks depict the ligands. H-bonds are shown in blue lines, while the hydrophobic contacts in dashed-gray lines. The π–π stacking interactions are represented in a dashed-green line with gray balls at the interacting ring center. Additionally, the salt bridge interactions are described in a dashed-yellow line with yellow balls at the positions of interacting atoms.

Table 1

The interactions that were established after docking the 12 adenosine derivatives, positive and negative controls into the SARS-CoV-2 RdRp

Compound	AutoDock score (kcal/mol)	H-bonding		Hydrophobic interaction
		Number	Amino acids involved	Number	Amino acids involved
Remdesivir	− 7.5	6	I548, S549, A550, Q815, D845 and R858a	2	F441 and F441b
1	− 7.5	10	W617, Y619(2), K621, N691, D760, D761(2), W800, and E811	0	N/A
2	− 7.2	10	W617, Y619, D761(4), A762, K798, K798a and W800	2	D618 and Y619
3	− 7.2	11	R555, D618, T619, K621, D623, R624, T680(2), T687, N691, and D760	0	N/A
4	− 6.4	7	D618, D761(3), W800, C813, and S814	2	E811(2)
5	− 7.1	11	S759, D761(2), W800, E811(2), C813, S814(3), and R836	0	N/A
6	− 7.6	11	K551, R569(2), R624(2), A685, T687, A688, S759(2), and D760	1	T687
7	− 7.6	8	R555, W617, Y619, K621, C622, D623(2), and D761	0	N/A
8	− 8.4	8	K551, C622, D623, N695, D761, E811, C813, and S814	3	Y619c, D623, and D760
9	− 8.2	13	R555, D618, K621, C622, D623, D761(2), W800, E811, C813, and S814(3)	2	Y619 and D760
10	− 8.3	13	R555, Y619, R624, D760(2), D761(2), W800, E811(2), C813, and S814(2)	1	Y619
11	− 8.1	12	K551, K551a, R555, R555a, D623(2), R624(3)a, N691, D760, and D761	3	D623, N691, and F694
12	− 7.7	15	R555(2), Y619(3), K621, D623(2), D760, D761(3), E811, C813, and S814	1	Y619
Cinnamaldehyde	− 5.4	2	R858(2)	4	F441(2), A550, and Y546b

Bold residues are the two consecutive active site aspartates (D760 and D761). The best four compounds (8, 9, 10, and 11) are bold and underlined

aResidues that interact with salt bridges

bResidues are that interact with π-stacking

cResidues represent the residues interacting through π–π stacking

After docking MDS

MDS for the best four complexes (8, 9, 10, and 11) were performed for 20 ns using NAMD and VMD software, and then the MM-GBSA was performed using Amber tools. Figure 5 shows the Root Mean Square Deviation (RMSD), Radius of Gyration (RoG), Surface Accessible Surface Area (SASA), and the number of H-bonds in the four complexes. The active site aspartates (D760 and D761) are shown underlined in Table 2 as well.

Fig. 5

The RMSD (in Å), RoG (in Å), SASA (in Å2), and the number of formed H-bonds versus the simulation time (in ns) for the best four complexes (8, 9, 10, and 11 with SARS-CoV-2 RdRp)

Table 2

The MM-GBSA calculations for the best four compounds (8, 9, 10, and 11) calculated after 20 ns MDS

Complex	RdRp-compound 8 complex		RdRp-compound 9 complex		RdRp-compound 10 complex		RdRp-compound 11 complex
	Residue	Binding energy (kcal/mol)	Residue	Binding energy (kcal/mol)	Residue	Binding energy (kcal/mol)	Residue	Binding energy (kcal/mol)
	Y619	− 2.58	K798	− 1.35	W800	− 3.31	D623	− 3.58
	C622	− 1.84	C813	− 1.33	D623	− 2.40	Y619	− 3.34
	N695	− 1.18	W800	− 1.27	F812	− 2.27	D760	− 2.09
	D761	− 1.10	R555	− 1.07	A762	− 1.85	R624	− 1.54
	S814	− 0.93	S814	− 1.06	H810	− 1.81	N691	− 1.35
	K621	− 0.90	H810	− 1.01	G616	− 1.56	R555	− 1.23
	N691	− 0.78	W617	− 1.00	Y619	− 1.50	C622	− 0.82
	D760	− 0.70	D623	− 0.92	K621	− 1.04	P620	− 0.60
	R624	− 0.46	C799	− 0.90	S814	− 0.99	K551	− 0.54
	F694	− 0.41	F812	− 0.89	R624	− 0.98	T687	− 0.49
	E811	− 0.34	R553	− 0.47	W617	− 0.98	S759	− 0.43
			N552	− 0.38	V763	− 0.67	I757	− 0.40
			G616	− 0.31	C799	− 0.60	N695	− 0.37
			R836	− 0.30	R555	− 0.58	T680	− 0.32
					M615	− 0.52
					E811	− 0.50
					C813	− 0.49
					C622	− 0.37
Total binding energy	− 10.65 kcal/mol		− 11.08 kcal/mol		− 19.03 kcal/mol		− 17.75 kcal/mol

Bold residues are that have binding energy greater than or equal to 1.0 kcal/mol

At least one of these two active site residues has a high contribution in binding compounds 8 (D761) and 11 (D760) with binding energy values of − 1.10 and − 2.09 kcal/mol, respectively. The residue Y619 has a high contribution in the binding of the RdRp to compounds 8 (− 2.58 kcal/mol), 10 (− 1.50 kcal/mol), and 11 (− 3.34 kcal/mol). In addition, the residues G616, W617, W800, H810, and F812 are main contributors for RdRp binding with compounds 9 (− 0.31, − 1.00, − 1.27, − 1.01, and − 0.89 kcal/mol), and 10 (− 1.56, − 0.98, − 3.31, − 1.81, − 2.27 kcal/mol). On the other hand, R555 and D623 contribute highly in binding compounds 9 (− 1.07 and − 0.92 kcal/mol), 10 (− 0.58 and − 2.40 kcal/mol), and 11 (− 1.23 and − 3.58 kcal/mol) to RdRp. In addition, C622 and R624 contribute in binding compounds 8 (− 1.84 and − 0.46 kcal/mol), 10 (− 0.37 and − 0.98 kcal/mol), and 11 (− 0.82 and − 1.54 kcal/mol) to RdRp, while N691 and N695 both contribute in binding compounds 8 (− 0.78 and − 1.18 kcal/mol), and 11 (− 1.35 and − 0.37 kcal/mol) to RdRp of SARS-CoV-2. The residue S814 contributed with − 0.93, − 1.06, and − 0.99 kcal/mol to bind the RdRp to compounds 8, 9, and 10, while the residue K621 contributed with − 0.90 and − 1.04 kcal/mol in binding the protein to compounds 8 and 10.

Figure 6A shows the ligands (8, 9, 10, and 11) RMSD in Å calculated for the four ligand-RdRp complexes during the 20 ns MDS period. The snapshot of each complex before and after the MDS is shown in Fig. 6B, with the ligands represented in colored sticks

Discussion

The anti-Ebola Virus (EBOV) drug, remdesivir, is one of the new drugs that give exciting results in silico, in vitro, and in vivo, and is recently approved by the Food and Drugs Administration (FDA) of the United States to treat COVID-19 emergency patients [33-37]. Remdesivir can be used as a seed for a more efficient drug against the SARS-CoV-2 pandemic as it targets one of the crucial viral proteins, the RNA-dependent RNA polymerase [38, 39].

The current study shows the results of the novel 12 adenosine analogs against the SARS-CoV-2 RdRp structure (PDB ID: 7BTF) using combined molecular docking and MDS. As reflected from Fig. 3, the best four compounds based on their average binding affinities to SARS-CoV-2 RdRps are compounds 8, 9, 10, and 11 (shown in green), ranging from − 8.07 ± 0.73 kcal/mol (11) down to − 8.13 ± 0.45 kcal/mol (8). In addition, compounds 1, 2, 3, 5, 6, 7, and 12 (blue columns) show comparable results regarding the average binding affinities, with remdesivir (red column), the positive control. Their average binding energies range from − 7.26 ± 0.66 kcal/mol (5) down to − 7.71 ± 0.65 kcal/mol (7), while the average binding energy of remdesivir is − 7.47 ± 0.28 kcal/mol. On the other hand, compound 4 (yellow column) shows higher average binding energy of − 6.3 ± 0.45 kcal/mol (lower affinity) compared to the positive control, while it still lower (better) than the negative control, cinnamaldehyde (− 5.39 ± 0.32 kcal/mol).

The primary type of interaction that is established after docking the 12 adenosine derivatives is the H-bonding. At least 7 bonds (4 and 9) and up to 15 H-bonds (12) are reported for the adenosine derivatives. In addition, few hydrophobic contacts are also reported in Table 1 for some compounds [maximum of three contacts (11)]. The active site residues D760 and D761 (bold residues in Table 1) are common and form H-bonding with all of the adenosine derivatives. This is common for the nucleotide inhibitors in different RNA viruses as the ligand will be added to the newly formed RNA [40-47]. The residues that are most common interacting with the adenosine derivatives ranked by the number of reported interactions are D761 (21), Y619 (15), D623 (11), S814 (11), E811 (10), D760 (9), R555 (8), and W800 (6). Compound 11 forms five salt bridges through K551, R555, and R624 (3), while compound 8 forms π–π stacking with Y619. On the other hand, for the positive control remdesivir, a salt bridge is formed with residue R858, while both positive and negative control show π-stacking interactions with F441 and Y546, respectively.

As reflected in Table 1 and Fig. 4, the active site aspartates D760 and D761 form H-bonds and few hydrophobic contacts with the ligands. The four adenosine derivatives 8, 9, 10, and 11 show the best binding affinity to the SARS-CoV-2 RdRp; hence they could present possible inhibitors of the polymerase.

Fig. 4

The interaction patterns for the best four adenosine derivatives (8, 9, 10, and 11) after docking into the active site of SARS-CoV-2 RdRp. The residues from SARS-CoV-2 RdRp that form contacts to the ligand are depicted in blue sticks, while the ligands are in orange. H-bonds are shown in blue lines, while hydrophobic interactions are in dashed-gray lines. Salt bridges are in dashed-yellow lines with yellow balls, while π–π stacking is in dashed-green lines with gray balls

As reflected in Fig. 5, the systems are equilibrated during the simulation period. In Table 2, the residual contribution for the binding of SARS-CoV-2 RdRp and the best four compounds (8, 9, 10, and 11) are listed with the bold residues for the highest contributed residues in the binding with binding energy greater than 1 kcal/mol. The compounds binding affinity against RdRp ranked are 10 > 11 > 9 > 8 with calculated binding affinity of − 19.03, − 17.75, − 11.08, and − 10.65 kcal/mol, respectively.

Calculation of the four ligands’ RMSD was done using TK console in VMD for each trajectory. Figure 6 shows the change in RMSD values for each compound along the trajectories. Compound 8 shows the smallest changes for most of the simulation duration with an increase in the last 4 ns, as shown in Fig. 6A (blue line). Compound 9 (orange line) shows a small RMSD for approximately 6 ns before increasing to an average value of 3 Å. Compound 10 (gray line) shows an increasing RMSD for 5 ns before stabilizing at an average value of 2.6 Å. Compound 11 (yellow line) offers the most stable fluctuation for approximately 10 ns before rising to an average value of 2.3 Å. Figure 6B shows the differences between the first and last frames of the simulations for compounds 8, 9, 10, and 11.

Fig. 6

A The ligands-RMSD (in Å) versus the simulation time (in ns) for the best four complexes (8, 9, 10, and 11 with SARS-CoV-2 RdRp). B The superposition of the first and last frames for the complexes where the ligands are represented in colored sticks as shown in the figure

The presented data show the favorable binding affinity of 11 adenosine derivatives against the SARS-CoV-2 RdRp compared to remdesivir. Four derivatives (8, 9, 10, and 11) show the best binding affinity to the RdRp and can be potential inhibitors against SARS-CoV-2 RdRp after in vitro activity testing and toxicity parameters checking experimentally.

Conclusion

Remdesivir shows an excellent result in terms of its efficacy in reducing the hospitalization time for COVID-19 patients. Twelve novel derivatives of adenosine triphosphates are tested against SARS-CoV-2 RdRp using the in silico techniques. The results reveal the binding potential of 11 derivatives against the viral polymerase. These adenosine modifications can be potential inhibitors of SARS-CoV-2 RdRp. Further, in vitro and other experimental validations are suggested for the best four compounds (8, 9, 10, and 11) as future work.