
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12817-4
10.1016/j.heliyon.2024.e36786
e36786
Research Article
Computational investigation of novel synthetic analogs of C-1′β substituted remdesivir against RNA-dependent RNA-polymerase of SARS-CoV-2
Cardoza Savio a
Singh Anirudh b
Sur Souvik c
Singh Mintu a
Dubey Kshatresh D. d
Samanta Sintu Kumar b
Mandal Ajay e
Tandon Vibha vtandon@mail.jnu.ac.in
vtandon@iicb.res.in
af⁎
a Special Centre for Molecular Medicine, Jawaharlal Nehru University, New Delhi, 110067, India
b Department of Applied Sciences, Indian Institute of Information Technology Allahabad, Allahabad, 211012, Uttar Pradesh, India
c Research and Development Center, Teerthanker Mahaveer University, Moradabad, Uttar Pradesh, 240001, India
d Department of Chemistry, Shiv Nadar University, Gautam Buddha Nagar, Uttar Pradesh, 201314, India
e Symbol Discovery Ltd, ASPIRE-TBI, School of Life Sciences, University of Hyderabad, Gachibowli, Hyderabad, 500046, India
f CSIR- Indian Institute of Chemical Biology (IICB), 4, Raja S C Mullick Road, Jadavpur, Kolkata, 700032, India
⁎ Corresponding author. Special Centre for Molecular Medicine, Jawaharlal Nehru University, New Delhi, 110067, India. vtandon@mail.jnu.ac.invtandon@iicb.res.in
23 8 2024
15 9 2024
23 8 2024
10 17 e3678618 1 2024
22 8 2024
22 8 2024
© 2024 Published by Elsevier Ltd.
2024

https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Remdesivir, a C-nucleotide prodrug binds to the viral RNA-dependent-RNA polymerase (RdRp) and inhibits the viral replication by terminating RNA transcription prematurely. It is reported in literature that interaction between the C-1’β–CN moiety of Remdesivir (RDV) and the Ser861 residue in RdRp enzyme, causes a delayed chain termination during the RNA replication process and is one of the important aspect of its mechanism of action. In the pursuance of increasing the biological activity of RDV and enhancing the SAR studies, against RNA viruses, we have designed its fourteen C1’β substituted analogs, 10 –23 bearing 4/5-membered heterocyclic rings. The docking and 100 ns molecular dynamics (MD) simulations of 10-23 to the RdRp protein (PDB ID: 7L1F) revealed important interactions between 2’,3’-diol, oxo group of phosphoramidate, nitrogen residues of heterocyclic rings of synthetic molecules with Arg555, Arg553, Ser759, Cys622, Asn691, Asp623 amino acid residues of protein. The docking score of 2-ethylbutyl ((S)-(((2R,3S,4R,5R)-5-(4-aminopyrrolo[2,1-f][1,2,4]triazin-7-yl)-3,4-dihydroxy-5-(1H-1,2,3-triazol-4-yl)tetrahydrofuran-2-yl)methoxy)(phenoxy)phosphoryl)-L-alaninate, 11 was found to be the higher than RDV among 14 new compounds i.e. -5.20 kcal/mol. Out of 3 compounds, 10, 12 and 13 submitted for MD simulations and Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) analysis, trifluoro-oxadiazole derivative, 13 showed higher binding energy as compared to Remdesivir. The predicted ADMET properties of 14 compounds showed their potential for being drug candidates. The present study suggests that substitution at the C1’β position by 4/5-membered rings plays an important role in the interactions between nucleoside/tide and target protein.

Keywords

Remdesivir analogs
Molecular docking
MD simulations
==== Body
pmc1 Introduction

New viruses threaten human survival as had been encountered in the recent COVID-19 pandemic [1,2]. Exploring small-molecule drugs presents a viable and prompt approach to counter such threats. Nucleoside/tide derivatives have been promising in the search of anti-viral drug candidates [3,4]. WHO has recommended the usage of Remdesivir, a C-nucleotide, for mild to moderate COVID-19 cases for patients at risk of hospitalization [5,6]. Inspired by Remdesivir's efficacy, there has been an extensive effort in the development of small-molecule inhibitors of SARS-CoV-2 since the outbreak in 2019 [[7], [8], [9]]. Computational tools have played a pivotal role in search of small-molecules targeting relevant proteins important in the SARS-CoV-2 machinery [10]. Molecular docking complemented by MD simulations reveal repurposed molecules targeting SARS-CoV-2 proteins.

The RNA-dependent RNA polymerase (RdRp) is a highly studied target of various nucleotide drugs due to its conserved binding site across different viruses [11,12]. It is responsible for the replication of RNA strand using an RNA template strand, which must be avoided to decelerate viral replication. Following the cryo-EM elucidation of SARS-CoV-2 RdRp structure [13], Cramer's group reported that a steric clash between the C-1′β -CN substituent and S861 amino acid residue of nsp12 causes a translocation barrier to the RNA strand formation [14]. Recently, Taylor and Johnson et al. observed that the addition of four covalently bonded Remdesivir monophosphate molecules to the reconstructed cryo-EM complex of RdRp, RNA primer, and template strands, exhibits a delayed chain termination [15]. A comparative computational study between Remdesivir, Ribavirin, and Favilavir portrayed a more favourable binding interaction of Remdesivir and RdRp which mimicked interactions between ATP molecule and RdRp [16]. Additional reports by various theoretical chemists point to the hypothesis that Remdesivir binds strongly to the RdRp protein [17,18]. Due to its essential role in viral replication, RdRp is an attractive target for anti-SARS-CoV-2 drug discovery and various molecules are known to effective target this protein (Fig. 1) [[19], [20], [21], [22], [23]].Fig. 1 Small-molecules targeting RdRp enzyme.

Fig. 1

Remdesivir is a C-nucleotide containing a C-C bond between the nucleobase and the ribose ring, thus making the C-glycosyl bond stable towards enzymatic degradation [24,25]. C-nucleosides have been extensively researched in the last decade with the majority of modifications on the nucleobase [[26], [27], [28], [29], [30], [31], [32]]. In an elaborative work, Gilead Sciences led the drug discovery program to identify potential 1,1′-disubstituted C-nucleosides/tides showing moderate inhibitory activity against the HCV virus [33]. These compounds contained a nitrile group at the 1′β position and a methyl group at the 2′ position [34]. Mackman et al. conducted a study wherein they screened ∼1000 molecules out of which, novel C-nucleosides bearing 1′-CN group and their phosphoramidate derivatives showed potent antiviral activity against Ebola virus, Hepatitis C virus and Respiratory synctical virus (Fig. 2) [35].Fig. 2 1,1′-disubstituted anti-viral C-nucleotides synthesized by Mackman et al.

Fig. 2

A number of researchers investigated RdRp inhibitors as one of the target to inhibit virus growth, but till date not much success was achieved. Our research focuses on proposing novel derivatives featuring C-1,1′ disubstitution on ribonucleoside/tide and understanding its impact on ligand-receptor interactions. In the present manuscript, we focussed on select 4/5-membered heterocyclic rings, which have not been studied previously on RdRp. These heterocycles have a rigid framework and can form significant H-bonds or aromatic pi-pi interactions in the binding pocket of RdRp. Previously our lab synthesized C-1,1′ substituted ribonucleoside/tide analogs containing tetrazole, triazole, and 1,2,4-oxadiazole, a few showing inhibition of SARS-CoV-2 infection at micromolar concentrations in Calu-3 cell lines [36]. In the present work, computational investigations were conducted on the binding of proposed C-1,1′ disubstituted ribonucleotides to the RdRp enzyme to elucidate the effects of moieties at the C-1,1’ β position of the RDV framework (Fig. 3). The current research aims to uncover the potential benefits of these modifications in different nucleoside analogs for combating RNA viruses. The results obtained with modified remdesivir analogs clearly suggests that designed compounds showed enhanced binding affinity, inhibitory constant (Ki), reduced toxicity and comparable anti SARS-CoV-2 activity with remdesivir.Fig. 3 Proposed derivatives for this work.

Fig. 3

2 Result and discussion

We docked the proposed compounds 10–23 on the RdRp protein (PDB id: 7L1F) at the ligand binding site, containing the critical catalytic residues of nsp12 i.e.: Asp618, Phe753, Asn767, Ser759, Asp760, and Asp761 which is conserved in most of the RNA viruses and compared their docking scores with that of Remdesivir (RDV) (Table 1). The 14 compounds showed a docking score in the range of −3.18 to −5.32 kcal/mol, of which compound 11 showed docking score of −5.20 kcal/mol close to that of RDV. All the compounds were observed to interact with the residues at the catalytic sites of nsp12 i.e. between Asn611-Met626 and between Phe753 - Asn767 [37]. The docking scores and binding scores were also calculated using Autodock and the docking scores of proposed derivatives were found to be better than RDV. The binding constant, Ki is an important parameter which determines the affinity of the molecule for the binding cavity and that for compounds 10, 11, 12, 14 and 17 were found to be better than RDV. Thus, the proposed derivatives have the potential of increased binding affinity compared to that of RDV.Table 1 Docking scores of synthesized C1′β analogs of Remdesivir. aDocking score of compounds 10–23, RDV obtained from Glide.

Table 1

Our results suggest that RDV forms hydrogen bonds via its 2′,3′-diol group and –NH2 of the nucleobase with the side chain's hydroxy group of Ser759 and Tyr619 of RdRp enzyme respectively (Fig. 4A). In addition to it, we observed the H-bond interaction between the oxygen (on the P=O group) and the side chain of NH2 group of Arg555 which is part of the crucial catalytic diad (Fig. 4). Whereas, compound 10 interacts via the tetrazole -NH group to the carboxylate group of Asp623 along with the carboxylate of Asp760, hydroxy group of Ser759, and carboxamide group of Asn691 (Fig. 4B). Triazole derivative, 11, was seen to have a characteristic π-π interaction and hydrogen bond with the side chain of NH2 group of diad Arg555 and Arg553 respectively and the triazole ring is exposed to the solvent-accessible part of the protein (Fig. 4C). The oxadiazole 12 interacts with the side chain of NH2 group of Arg553, Arg555, and SH group of Cys622 via its phosphoramidate group (Fig. 4D). The amino group on the nucleobase forms hydrogen bonds with the hydroxy group of Ser682 residue. 13 forms fewer hydrogen bonds with the hydroxy group of Ser759 and carboxamide group of Asn691 residues (Fig. 4E). The tetrazole derivative (10), methyl-substituted 1,2,4-oxadiazole (12), and trifluoromethyl-substituted 1,2,4-oxadiazoles (13) were chosen as representative compounds for MD simulation studies using GROMACS 2019.4. Present study aims to provide a detailed understanding of how these specific compounds interact with the RdRp protein at the molecular level.Fig. 4 2D representation of docked poses and binding interactions of A) RDV, B) 10, C) 11, D) 12, E) 13.

Fig. 4

Various structural order parameters, including Root mean square deviation (RMSD), Root mean square fluctuation (RMSF), Radius of gyration (Rg), and Solvent accessible surface area (SASA), were assessed for RdRp and its complexes (Fig. 5). RMSD measures the root mean square deviation, crucial for understanding macromolecular conformational changes during simulation. RdRp protein's backbone-RMSD reached equilibrium with an average deviation of approximately 0.4 nm. RdRp complexed with RDV stabilized within 25 ns and remained stable until 100 ns (Fig. 5A). Complexes 10, 12, and 13 with RdRp exhibited average RMSD values ranging from ∼0.31 to ∼0.39 nm, reaching equilibrium within 20 ns. These RMSD values indicate that the proposed models occupy favourable spatial positions within the receptor's binding pocket.Fig. 5 A) Root Mean Square Distance (RMSD) of RdRp and complexes of RDV, 10,12 AND 13 with RdRp, B) Radius of gyration (Rg) of RdRp and complexes of RDV, 10,12 and 13 with RdRp C) Solvent accessible surface area (SASA) of RdRp and complexes of RDV, 10,12 and 13 with RdRp.

Fig. 5

The Rg parameter determines the compactness of the complex, aiding in understanding any significant changes in the structural integrity. RdRp and its complex with RDV maintained consistent Rg values (RdRp: 3.16 ± 0.02 nm; RdRp-RDV: 3.15 ± 0.01 nm) (Fig. 5B). Complexes with ligands 10, 12, and 13 exhibited Rg values ranging from 3.14 ± 0.01 nm to 3.17 ± 0.01 nm. No significant differences were observed in comparison to the apo protein and RDV-RdRp complex.

SASA provides critical information regarding the surface area of a protein that is accessible to solvent molecules in molecular dynamics (MD) simulations. It offers insights into the protein's interactions with its surrounding environment, such as water molecules. Changes in SASA values over the course of a simulation can indicate alterations in the protein's conformation, folding, or interactions with ligands, which are crucial for understanding its stability and behaviour [38]. All complexes and unbound RdRp protein showed nominal deviations in SASA, suggesting stable dynamics (Fig. 5C). Average SASA values observed were: RdRp (445.93 ± 6.24 nm), RdRp-RDV (454.50 ± 8.50 nm), Complex 10 (455.26 ± 6.50 nm), Complex 12 (446.57 ± 5.61 nm), Complex 13 (451.46 ± 4.93 nm).

The molecular dynamics simulations of protein-ligand complexes and unbound protein maintained stable conformations throughout the simulation. Minimal structural changes or denaturation occurred under the simulation conditions.

The RMSF analysis provides insights into the average fluctuations of amino acid residues within the protein structure. High RMSF values suggest certain residues undergo significant fluctuations (i.e. deviations from their average positions), often observed in terminal residues and loops, indicating lower structural stability. Conversely, lower RMSF values indicate rigidity and stability, with residues experiencing minimal deviations from their average positions. In our study of the RdRp protein, RMSF quantification revealed an average value of ∼0.15 nm, indicative of a stable secondary structure (Fig. 6). However, complexes bearing ligands exhibited RMSF values ranging from 0.16 to 0.34 nm, suggesting varying degrees of flexibility within the complexes.Fig. 6 RMSF of RdRp and its complexes with ligands.

Fig. 6

The time evolution plots of no. of hydrogen bonds between protein and RDV analogs were computed, showing the no. of H-bonds between RDV analogs proposed in present work and RdRp within each time frame. RDV formed an average of ∼2 hydrogen bonds from 0 to 40 ns, increased to ∼3 bonds from 40 to 70 ns, and then decreased back to ∼2 bonds from 70 to 100 ns (Fig. 7). Specifically, the binary complex of compound 10-RdRp, carrying a tetrazole-moiety showed an average of ∼6 hydrogen bonds from 0 to 75 ns and ∼2 bonds from 75 to 100 ns. Similarly, the oxadiazole bearing compound 12-RdRp complex, maintained consistent hydrogen bonding throughout the simulation, while compound 13-RdRp complex carrying -CF3 substituted 1,2,4-oxadiazole exhibited fluctuating hydrogen bond formation. Overall, the analysis suggests that tetrazole-incorporated compounds tend to form more hydrogen bonds within equilibration states.Fig. 7 Time-evolution plot of hydrogen bonds formed between ligands and the protein.

Fig. 7

MD simulations generate vast amounts of data regarding the positions and motions of atoms over time. PCA (Principal component analysis) reduces high-dimensional data into a smaller set of principal components (PCs), which represent the most significant modes of motion in the system. This simplification makes it easier to analyze and interpret the data. The Free energy surface (FES) plots allows for the simultaneous exploration of both the structural variability (captured by PCA) and the associated free energy landscape.

In Fig. 8A-E, the x and y axes typically represent the principal components (PC1 and PC2, principal component 1 and 2), which capture the most significant modes of motion. The color or contour lines represent the free energy landscape, with lower energy regions corresponding to more stable or favored conformations, and higher energy regions indicating less stable or less populated states. By overlaying the FEL onto the PCA plot, researchers can visualize how the energy landscape correlates with the structural dynamics captured by PCA. The energy basins represent clusters of similar conformations with low potential energy. Narrow energy barriers between basins indicate easy transitions, while higher barriers suggest less frequent transitions requiring more energy. The RdRp explores a broad conformational space, effortlessly transitioning between states within energy basins along PC2 with a small energy barrier of <3 kcal/mol (Fig. 8A). In contrast, the RDV complex with RdRp navigates a narrower space, hinting at a strongly bound ligand-protein interaction due to a single minimum and a 5 kcal/mol transition barrier (Fig. 8B). Similarly, the 10-RdRp complex displays stable transitions within a confined conformational space, with a transition energy of <2 kcal/mol (Fig. 8C). Conversely, the 12-RdRp complex showcases a broader space than RDV-RdRp, with multiple energy basins and <4 kcal/mol transition barriers (Fig. 8D). The 13-RdRp complex exhibits an even wider space, with <3 kcal/mol barriers (Fig. 8E). These findings imply stronger binding energetics of tetrazole and oxadiazole substituted compounds with the protein similar to Remdesivir. Further analysis, including binding affinity measurements and biological assays, are crucial to confirm these predictions and assess drug-like properties for potential development.Fig. 8 The 2D plot of the Free Energy Landscape of proposed complexes and protein A) RdRp model B) RDV-RdRp complex, C) 10-RdRp complex, D) 12-RdRp complex and E) 13-RdRp complex.

Fig. 8

The MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) analysis provides a detailed understanding of the energetics of protein-ligand binding at the atomic and residue level. It helps to identify key residues that are crucial for ligand binding and can guide further research into drug design, optimization, and the study of protein-ligand interactions [38,39]. Remdesivir (RDV) and 13 have comparable binding energies during the entire duration, whereas 10 and 12 have higher binding energy (Fig. 9). The average binding energy in the last 50 ns for RDV, 10, 12, and 13 were found to be −60.90 kcal/mol, 8.40 kcal/mol, 0.090 kcal/mol, −75.89 kcal/mol.Fig. 9 Plots of binding energy (kcal/mol) and simulation time (ps) for prodrugs.

Fig. 9

The 1′β substitution was seen to have a favourable impact on the binding energy for the prodrugs. MMPBSA calculations do not take into effect the entropy factor and positive ΔG may be observed due to it.

QikProp 3.1-Schrodinger, a widely utilized tool, facilitates the prediction of absorption, distribution, metabolism, and excretion (ADME) properties, making it a good choice in pharmaceutical research for high throughput screening. Predictions of ADME properties were conducted for the proposed molecules (Table 2), aiding in the assessment of their druggability [40].Table 2 Pharmacokinetics properties of 14 analogs (10–23) of Remdesivir.

Table 2Compounds	SASA	FOSA	FISA	PISA	WPSA	volume	logBB	logKhsa	% HOA	
10	911.96	339.41	305.09	266.42	1.04	1760.80	−3.74	−0.83	9.72	
11	936.54	348.12	285.46	300.70	2.28	1782.72	−3.64	−0.67	16.78	
12	956.18	441.87	244.39	269.05	0.88	1834.01	−3.18	−0.63	39.45	
13	972.78	339.10	237.79	269.36	126.54	1870.18	−2.81	−0.56	44.84	
14	912.05	337.63	249.50	324.13	0.79	1758.11	−3.13	−0.76	36.39	
15	979.61	347.18	314.86	315.32	2.25	1846.39	−4.20	−0.94	21.80	
16	1024.96	535.74	222.62	265.70	0.90	1994.76	−3.053	−0.33	49.41	
17	843.86	311.44	227.14	305.27	0	1719.26	−2.58	−0.49	43.63	
18	805.58	329.19	161.78	314.58	0.04	1692.68	−1.82	−0.55	55.91	
19	893.94	427.92	199.38	266.63	0	1758.73	−2.48	−0.41	51.82	
20	861.73	391.81	196.94	272.45	0.53	1736.87	−2.34	−0.56	49.12	
21	926.82	411.79	199.34	259.79	55.90	1804.76	−2.43	−0.25	54.88	
22	904.66	440.06	234.35	230.25	0	1774.24	−2.89	−0.25	32.46	
23	896.11	393.07	248.55	248.74	5.74	1817.51	−2.89	−0.72	36.45	
RDV	874.89	346.16	252.23	275.63	0.87	1692.62	−3.13	−0.68	36.75	
[Range: SASA (300.0–1000.0), FOSA (0.0–750.0), FISA (7.0–330.0), PISA (0.0–450.0), WPSA (0.0–175.0), volume (500.0–2000.0), logBB (−3.0 to 1.2), log Khsa (−1.5 to 1.5), %HOA (>80%-high; <25%-poor)].

The solvent-accessible surface area (SASA), in this case is indicative of the extent of contact between solvent and molecule, larger SASA often correlates with better absorption and distribution of the drug in the body, as it indicates increased surface area for interaction with biological membranes, which typically ranges from 300.0 to 1000.0 Å2. All compounds show values greater than 800 Å2 which is indicative of good profile. Fully Occupied Solvent Accessible Surface Area (FOSA) represents the maximum surface area that solvent molecules can occupy around the compound and assess the solvation behaviour of the compounds, falls within 311–535 Å2. Except compound 16, all other RDV analogs FOSA values suggest that these compounds will be able to cross cell membrane. Additionally, hydrophilic component (FISA) which indicates drugs affinity for water molecules suggesting its solubility and pi component (PISA) that is indicative of the pi-electron surface area of the molecule which helps in aromatic interactions were predicted to be within acceptable ranges (7.0–33.0 Å2 and 0.0–450.0 Å2, respectively). The weakly polar component of SASA (WPSA) and volume enclosed by the solvent-accessible molecular surface of compounds were found to be in acceptable range as well. These findings collectively suggest promising drug profiles for the compounds in solvation behaviour. However, brain-blood partition coefficients for most compounds slightly exceeded the acceptable range of −3.0 to 1.2, indicating potential central nervous system entry. Prediction results indicated acceptable human intestinal absorption capacity for all compounds except 10 and 11. Overall, the predicted pharmacokinetic data supports the viability of these compounds for further development as antiviral agents.

The toxicity profiles were also estimated for the proposed derivatives using an open-source software ADMETSAR (Table 3) [41]. ADMETSAR serves as an invaluable resource for researchers in pharmaceutical science, offering a comprehensive platform for the prediction of crucial drug properties such as absorption, distribution, metabolism, excretion, and toxicity (ADMET). Drawing from an extensive compilation of over 200,000 meticulously curated data points sourced from scientific literature, ADMETSAR provides comprehensive insights into the ADMET profiles of nearly 96,000 unique compounds and new compounds based on the structural similarites. We reported here the probability of proposed derivatives being AMES toxic, carcinogenic and acute oral toxic. Log HERG was predicted using Qikprop. The models used have been given a probability output and based on it proposed derivatives have been predicted to be non AMES-toxic and non carcinogenic. Proposed derivatives including Remdesivir have been predicted to be slightly toxic and classified under class III for acute oral toxicity. "Log HERG K″ denotes the logarithm of the half-maximal inhibitory concentration (IC50) values concerning a compound's interaction with the human Ether-à-go-go-Related Gene (hERG) potassium ion channel. This channel is vital for regulating the heart's electrical activity, and compounds that obstruct hERG channels may induce cardiac arrhythmias and other cardiovascular issues. Evaluation of HERG K+ channel blockers' potential toxicity through forecasted IC50 values indicated no increase in toxicity due to newly added fragments in the RDV structure.Table 3 Predicted toxicity of 14 analogs (10–23) of Remdesivir.

Table 3Compounds	Probability of AMES mutagenicity	Probabilty of being Non-Carcinogen	Probability of Acute Oral Toxicity	Log HERGa	
10	0.60	0.83	0.54	−6.47	
11	0.60	0.83	0.54	−6.86	
12	0.61	0.81	0.54	−6.73	
13	0.58	0.77	0.55	−6.79	
14	0.60	0.81	0.53	−6.76	
15	0.61	0.81	0.53	−7.19	
16	0.62	0.81	0.53	−6.88	
17	0.60	0.81	0.54	−5.86	
18	0.63	0.83	0.53	−5.42	
19	0.61	0.87	0.55	−6.23	
20	0.62	0.82	0.54	−5.89	
21	0.61	0.81	0.54	−6.44	
22	0.60	0.77	0.54	−6.13	
23	0.58	0.64	0.55	−5.89	
RDV	0.62	0.83	0.54	−6.32	
a logHerg (<−5; concern).

3 Conclusion

In summary, docking score = −5.20 kcal/mol of compound 11 further supports its potential as an inhibitor of RdRp, which is significant in the context of antiviral drug development, particularly for targeting RNA viruses like SARS-CoV-2. The Ki values of 10–14, 17,18, 20 and 21 were observed better than remdesivir. Further experimental and computational studies may be conducted to validate these findings and evaluate the potency of compounds as antiviral agents against RNA virus. MD simulations studies of tetrazole 10, -CH3 substituted oxadiazole 12 and -CF3 substituted oxadiazole 13 proved that binary complex between these compounds and RdRp got stabilized within first 20 ns with an RMSD range of 0.31–0.38 nm. Radius of gyration for the complexes of above 3 compounds with RdRp did not show major deviations revealing a stable binding interaction. The solvent accessible surface area did not deviate much from the apo enzyme and RDV-RdRp complex. Therefore, the structural order parameters of the complex were found to be consistent with that of RDV bound RdRp complex. The Free energy surface (FES) analysis too were found to show a stable energetics for the binding of compounds. The MM-PBSA calculations of complexes reveal comparable binding energies of RDV and 13. The complex between compounds 10, 12 and RdRp showed relatively higher binding energy than RDV during the 100 ns simulations respectively. In conclusion, the positive outcomes from computational studies and analysis of drug likeliness properties of proposed compounds provide strong motivation to continue exploring the potential of 1’ substitutions on RDV as a strategy for designing effective inhibitors of RdRp.

4 Experimental section

4.1 Preparation of RdRp and compounds

The structure of the RdRp protein was retrieved from Protein Data Bank (PDB ID: 7L1F). The compounds including Remdesivir, were built using Chemdraw Professional 16.0 by PerkinElmer and Chemdraw 3D. The docking studies were done in Schrodinger, LLC using the Maestro interface [42]. The 7L1F was pre-processed followed by energy optimization and minimization using Schrodinger's protein preparation wizard. The receptor grid was prepared with the size of 72 Å × 72 Å x 88 Å around the 4 covalently bonded Remdesivir monophosphates. 14 proposed compound structures were minimized using the LigPrep suite in Schrodinger.

4.2 Molecular docking of 14 proposed analogs of RDV to RdRp

Molecular docking was done using the Glide suite with extra precision in Schrodinger package (v 31207), retaining the chiralities from the input 3D structures with possible ionizations at pH = 7.0 ± 2.0. Besides, molecular docking was also performed using AutoDock which involved several key steps to predict the preferred orientation of a ligand bound to a receptor. First, the receptor's 3D structure was cleaned by removing non-standard residues and water molecules, with polar hydrogens added and Kollman charges assigned. Ligands were prepared by adding Gasteiger charges and converting the structure to PDBQT format. The docking grid box was defined using AutoDock Tools (ADT) to encompass the receptor's binding site. The Autogrid module generated grid maps, while the Autodock module performed the docking simulation using specified parameters, typically involving the Lamarckian Genetic Algorithm.

4.3 molecular dynamics of complexes of RdRp and compounds

Gromacs 2019.4 was used to build systems and conduct MD simulations for selected complexes using CHARMM36 forcefield and transferable intermolecular potential 3P (TIP3P) [43] as water models. The parameters of selected compounds were generated from the Swissparam website by Swiss Institute of Bioinformatics [44]. CHARMM36 forcefield was used for MD simulation of RdRp protein without bound ligand. Water molecules were added within the box shape of a dodecahedron with dimensions ∼1 nm. 10 Na+ ions were added to the system to neutralize the charges. The periodic boundary conditions were defined in x, y, and z directions, and the electrostatic interactions were evaluated using Particle-mesh Ewald methods [45]. The system was minimized using a steep descent algorithm for 50000 steps. Equilibration was first performed for 100 ps in the NVT ensemble followed by 100 ps in the NPT ensemble. The simulations were run for 100 ns on a Graphics Processing Unit (GPU) server.

4.4 Molecular dynamics analysis of complexes of RdRp and compounds

The structural order parameters measured using GROMACS were root mean square distance (RMSD) of protein and ligand, root-mean-square fluctuations (RMSF) for protein in complexes, solvent accessible surface area (SASA), Radius of gyration (Rg), and hydrogen bond interaction (H-Bonds). The hydrogen bonds were defined with the distance cut-off = 3.5 Å angle cut-off = 30°. Two components PC1 and PC2 were taken into consideration for principal component analysis [46,47] and FES analysis was done using the gmx-sham algorithm [48].

The principal component analysis (PCA) is the study of the collective motion of all atoms in the simulation [49]. It is known that the function of proteins is dependent upon the collective motion of all their atoms. This collective motion is used as a parameter to observe the stability of the protein. The large number of trajectories defined during MD simulations can be reduced to be visualized in an orthogonal space. This space is defined by eigenvectors which are produced by diagonalization of the covariance matrix, Cij.Cij = (xi − xj) (i,j = 1,2,3, …,3N)

where i and j are the pairs of the 3 N cartesian coordinates and xi and xj are instantaneous values of the i-th and j-th Cα atom, respectively. This reduction method allows one to understand an overall change in protein conformation that occurs during the simulation time. The 3-dimensional cartesian coordinates generated at each trajectory are mapped to the linear combination of orthogonal vectors called principal components (PC1, PC2).

The Free energy surface provides details on the minimum conformations achieved by the entire system. This method helps in identifying conformational transitions due to the protein-ligand binding interactions. The FES is plotted using the Boltzmann inversion formulaG = -RTlnP

Where P is the probability distribution of PC1 and PC2.

4.5 Binding free energy estimation of complexes of RdRp and compounds

Molecular Mechanics Generalized Born surface area (g_MMPBSA) in GROMACS was used to estimate the binding energy of complexes derived from MD trajectories [50]. It defines the molecular stability of the compound in the active site of the receptor. ΔGbind, the binding energy of complex is:ΔGbind. = ˂GCom > - ˂GR˃ - ˂GC˃

Where, GCom, GR, and GC are free energy of RdRp-compound complex, receptor, and compound respectively and terms within <, > denote them as ensemble average.

4.6 ADMET properties of 14 proposed compounds (remdesivir analogs)

The properties and descriptors of analyzed structures were predicted using QikProp 3.1 from Schrodinger package (v 31207) [40]. QikProp was run in normal mode. The three-dimensional structures of compounds were prepared in LigPrep 2.2 using settings recommended in the QikProp's user manual [51]. Surface area components and other ADMET values were predicted by the QikProp module in Schrodinger for the following descriptors: total solvent accessible surface area (SASA: 300–1000), hydrophobic component of the SASA (FOSA: 0.0–750.0), hydrophilic component of SASA (FISA: 7.0–300.0), π component of SASA (PISA: 0.0–450.0) and weakly polar component of SASA (WPSA: 0.0–175.0), predicted IC50 value for blockage of HERG K+ channels, logHerg (<−5; concern), predicted brain/blood partition coefficient, logBB (−3.0 to 1.2), prediction of binding to human serum albumin, log Khsa (−1.5 to 1.5) and predicted human oral absorption on 0 %–100 % scale, %HOA (>80%-high; <25%-poor). ADMETSAR was used to predict the toxicity profile of the proposed derivatives of remdesivir [41].

CRediT authorship contribution statement

Savio Cardoza: Writing – original draft, Validation, Methodology, Investigation. Anirudh Singh: Visualization, Investigation. Souvik Sur: Writing – review & editing. Mintu Singh: Writing – review & editing. Kshatresh D. Dubey: Supervision. Sintu Kumar Samanta: Writing – review & editing, Resources. Ajay Mandal: Conceptualization. Vibha Tandon: Writing – review & editing, Resources, Project administration, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Vibha Tandon reports financial support was provided by 10.13039/501100001843 Science and Engineering Research Board . If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The work was supported by the DST-SERB-CRG sponsored project no. SPG/2021/004179, and DBT-BIRAC-AMR grant no. BT/PR31944/MED/29/1408/2019 sanctioned to Prof. Vibha Tandon, SCMM, JNU. SC and MS acknowledge the fellowship support from UGC, CSIR-SRF and ICMR project no 2021–13818/SCR/ADHOC-BMS.
==== Refs
References

1 Ryu W.S. New Emerging Viruses Molecular Virology of Human Pathogenic Viruses 2017 289 302
2 Shang Z. Chan S.Y. Liu W.J. Li P. Huang W. Recent insights into emerging coronavirus: SARS-CoV-2 ACS Infect. Dis. 7 6 2021 1369 1388 33296169
3 Yates M.K. Seley-Radtke K.L. The evolution of antiviral nucleoside analogues: a review for chemists and non-chemists. Part II: complex modifications to the nucleoside scaffold Antivir. Res. 162 2019 5 21 30529089
4 K. L.Seley-Radtke Yates M.K. The evolution of nucleoside analogue antivirals: a review for chemists and non-chemists. Part 1: early structural modifications to the nucleoside scaffold Antivir. Res. 154 2018 66 86 29649496
5 WHO WHO Recommends Highly Successful COVID-19 Therapy and Calls for Wide Geographical Distribution and Transparency from Originator.
6 Consortium W.S.T. Repurposed antiviral drugs for Covid-19—interim WHO solidarity trial results New Engl. J. Med. 384 6 2021 497 511 33264556
7 C D. Savi Hughes D.L. Kvaermo L. Quest for a COVID-19 cure by repurposing small-molecule drugs: mechanism of action, clinical development, synthesis at scale, and outlook for supply Org. Process Res. Dev. 24 6 2020 940 976 37556267
8 Owen D.R. Allerton C.M.N. Anderson A.S. Aschenbrenner L. Avery M. Berritt S. Boras B. Cardin R.D. Carlo A. Coffman K.J. Dantonio A. Di L. Eng H. Ferre R. Gajiwala K.S. Gibson S.A. Greasley S.E. Hurst B.L. Kadar E.P. Kalgutkar A.S. Lee J.C. Lee J. Liu W. Mason S.W. Noell S. Novak J.J. Obach R.S. Ogilvie K. Patel N.C. Pettersson M. Rai D.K. Reese M.R. Sammons M.F. Sathish J.G. Singh R.S.P. Steppan C.M. Stewart A.E. Tuttle J.B. Updyke L. Verhoest P.R. Wei L. Yang Q. Zhu Y. An oral SARS-CoV-2 Mpro inhibitor clinical candidate for the treatment of COVID-19 Science 374 6575 2021 1586 1593 34726479
9 Urvashi Senthil Kumar J.B. Das P. Tandon V. Development of azaindole-based frameworks as potential antiviral agents and their future perspectives J. Med. Chem. 65 9 2022 6454 6495 35477274
10 Gao K. Wang R. Chen J. Cheng L. Frishcosy J. Huzumi Y. Qiu Y. Schluckbier T. Wei X. Wei G. Chem Rev. 122 2022 11287 11368 35594413
11 Hillen H.S. Kokic G. Farnung L. Dienemann C. Tegunov D. Cramer P. Structure of replicating SARS-CoV-2 polymerase Nature 584 7819 2020 154 156 32438371
12 Venkataraman S. Prasad B. Selvarajan R. RNA dependent RNA polymerases: insights from structure, function and evolution Viruses 10 2 2018 76 99 29439438
13 Yin W. Mao C. Luan X. Shen D.D. Shen Q. Su H. Wang X. Zhou F. Zhao W. Gao M. Chang S. Xie Y.C. Tian G. Jiang H.W. Tao S. Shen J. Jiang Y. Jiang H. Xu Y. Zhang S. Zhang Y. Xu H.E. Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by remdesivir Science 368 6498 2020 1499 1504 32358203
14 Kokic G. Hillen H.S. Tegunov D. Dienemann C. Seitz F. Schmitzova J. Farnung L. Siewert A. Höbartner C. Cramer P. Mechanism of SARS-CoV-2 polymerase stalling by remdesivir Nat. Commun. 12 1 2021 279 33436624
15 Bravo J.P.K. Dangerfield T.L. Taylor D.W. Johnson K.A. Remdesivir is a delayed translocation inhibitor of SARS-CoV-2 replication Mol. Cell. 81 7 2021 1548 1552 33631104
16 Byléhn F. Menéndez C.A. Perez-Lemus G.R. Alvarado W. de Pablo J.J. Modeling the binding mechanism of remdesivir, Favilavir, and ribavirin to SARS-CoV-2 RNA-dependent RNA polymerase ACS Cent. Sci. 7 1 2021 164 174 33527086
17 Nguyen H.L. Thai N.Q. Truong D.T. Li M.S. Remdesivir strongly binds to both RNA-dependent RNA polymerase and main protease of SARS-CoV-2: evidence from molecular simulations J. Phys. Chem. B 124 50 2020 11337 11348 33264025
18 Eweas A.F. Alhossary A.A. Abdel-Moneim A.S. Molecular docking reveals ivermectin and remdesivir as potential repurposed drugs against SARS-CoV-2 Front. Microbiol. 11 2021 3602
19 Chien M. Anderson T.K. Jockush S. Tao C. Li X. Kumar S. Russo J.J. Kirchdoerfer R.N. Ju J. J. Proteome Res. 19 11 2020 4690 4697 32692185
20 Mondal S.K. Mukhoty S. Kundu H. Ghosh S. Sen M.K. Das S. Brogi S. In silico analysis of RNA-dependent RNA polymerase of the SARS-CoV-2 and therapeutic potential of existing antiviral drugs Comput. Biol. Med. 135 2021 104591 104604 34216889
21 Elfiky A.A. Ribavirin Remdesivir Sofosbuvir Galidesivir Tenofovir against SARS-CoV-2 RNA dependent RNA polymerase (RdRp): a molecular docking study Life Sci. 253 2020 117592 117598 32222463
22 Kumar R. Gideon D.A. Mariadasse R. Nirusimhan V. Castin J. Jeyakanthan J. Dhayabaran V.V. In Silico Evaluation of Isatin-Based Derivatives with RNA-dependent RNA Polymerase of the Novel Coronavirus SARS-CoV-2 2021 10.31219/osf.io/rkz9n
23 Wakchaure P.D. Ghosh S. Ganguly B. Revealing the inhibition mechanism of RNA-dependent RNA polymerase (RdRp) of SARS-CoV-2 by remdesivir and nucleotide analogues: a molecular dynamics simulation study J. Phys. Chem. B 124 47 2020 10641 10652 33190493
24 Štambaský J. Hocek M. Kočovský P. C-nucleosides: synthetic strategies and biological applications Chem. Rev. 109 12 2009 6729 6764 19761208
25 Bzowska A. Kulikowska E. Shugar D. Purine nucleoside phosphorylases: properties, functions, and clinical aspects Pharmacol. Ther. 88 3 2000 349 425 11337031
26 Hosoya M. Shigeta S. Nakamura K. De Clercq E. Inhibitory effect of selected antiviral compounds on measles (SSPE) virus replication in vitro Antivir. Res. 12 2 1989 87 97 2480744
27 Morrey J.D. Smee D.F. Sidwell R.W. Tseng C. Identification of active antiviral compounds against a New York isolate of West Nile virus Antivir. Res. 55 1 2002 107 116 12076755
28 De Clercq E. Ebola virus (EBOV) infection: therapeutic strategies Biochem. Pharmacol. 93 1 2015 1 10 25481298
29 Kitaoka S. Konno T. De Clercq E. Comparative efficacy of broad-spectrum antiviral agents as inhibitors of rotavirus replication in vitro Antivir. Res. 6 1 1986 57 65 3008646
30 Warren T.K. Wells J. Panchal R.G. Stuthman K.S. Garza N.L. Van Tongeren S.A. Dong L. Retterer C.J. Eaton B.P. Pegoraro G. Protection against filovirus diseases by a novel broad-spectrum nucleoside analogue BCX4430 Nature 508 7496 2014 402 405 24590073
31 Huggins J.W. Robins R.K. Canonico P.G. Synergistic antiviral effects of ribavirin and the C-nucleoside analogs tiazofurin and selenazofurin against togaviruses, bunyaviruses, and arenaviruses Antimicrob. Agents Chemother. 26 4 1984 476 480 6151377
32 Walker J.A. Liu W. Wise D.S. Drach J.C. Townsend L.B. Synthesis and antiviral evaluation of certain novel pyrazinoic acid C-nucleosides J. Med. Chem. 41 8 1998 1236 1241 9548814
33 Eastman R.T. Roth J.S. Brimacombe K.R. Simeonov A. Shen M. Patnaik S. Hall M.D. Remdesivir: a review of its discovery and development leading to emergency use authorization for treatment of covid-19 ACS Cent. Sci. 6 5 2020 672 683 32483554
34 Mish M.R. Cho A. Kirschberg T. Xu J. Sebastian Zonte C. Fenaux M. Park Y. Babusis D. Feng J.Y. Ray A.S. Kim C.U. Preparation and biological evaluation of 1′-cyano-2′-C-methyl pyrimidine nucleosides as HCV NS5B polymerase inhibitors Bioorg. Med. Chem.Lett. 24 14 2014 3092 3095 24907145
35 Siegel D. Hui H.C. Doerffler E. Clarke M.O. Discovery and synthesis of a PhosphoramidateProdrug of a pyrrolo[2,1-f][triazin-4-amino] adenine C-nucleoside (GS-5734) for the treatment of Ebola and emerging viruses J. Med. Chem. 60 5 2017 1648 1661 28124907
36 Cardoza S. Shrivash M.K. Riva L. Chatterjee A.K. Mandal A. Tandon V. Multistep synthesis of analogues of remdesivir: incorporating heterocycles at the C-1′ position J. Org. Chem. 88 13 2023 9105 9122 37276453
37 Padhi A.K. Shukla R. Saudagar P. Tripathi T. High-throughput rational design of the remdesivir binding site in the RdRp of SARS-CoV-2: implications for potential resistance iScience 24 1 2021 101992
38 Kumari R. Kumar R. Consortium O.S.D.D. Lynn A. g_mmpbsa - a GROMACS tool for high-throughput MM-PBSA calculations J. Chem. Inf. Model. 54 7 2014 1951 1962 24850022
39 Genheden S. Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities Expert Opin. Drug Discov. 10 5 2015 449 461 25835573
40 Schrödinger Release 2023-3: QikProp 2023 Schrödinger, LLC New York, NY
41 Cheng F. Li W. Zhou Y. Shen J. Wu Z. Liu G. Lee P.W. Tang Y. admetSAR: a comprehensive source and free tool for evaluating chemical ADMET properties J. Chem. Inf. Model. 52 11 2012 3099 3105 23092397
42 Halgren T.A. Murphy R.B. Friesner R.A. Beard H.S. Frye L.L. Pollard W.T. Banks J.L. Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening J. Med. Chem. 47 7 2004 1750 1759 15027866
43 Jorgensen W.L. Chandrasekhar J. Madura J.D. Impey R.W. Klein M.L. Comparison of simple potential functions for simulating liquid water J. Chem. Phys. 79 2 1983 926 935
44 Zoete V. Cuendet M.A. Grosdidier A. Michielin O. SwissParam: a fast force field generation tool for small organic molecules J. Comput. Chem. 32 11 2011 2359 2368 21541964
45 Essmann U. Perera L. Berkowitz M.L. Darden T. Lee H. Pedersen L.G.A. A smooth particle mesh Ewald method J. Chem. Phys. 103 19 1995 8577 8593
46 Laberge M. Yonetani T. Molecular dynamics simulations of hemoglobin A in different states and bound to dpg: effector-linked perturbation of tertiary conformations and HbA concerted dynamics Biophys. J. 94 7 2008 2737 2751 18096633
47 Doharey P.K. Singh V. Gedda M.R. Sahoo A.K. Varadwaj P.K. Sharma B. In silico study indicates antimalarials as direct inhibitors of SARS-CoV-2-RNA dependent RNA polymerase J. Biomol. Struct. Dyn. 40 12 2022 5588 5605 33475021
48 Prakash A. Kumar V. Banerjee A. Lynn A.M. Prasad R. Structural heterogeneity in RNA recognition motif 2 (RRM2) of TAR DNA-binding protein 43 (TDP-43): clue to amyotrophic lateral sclerosis J. Biomol. Struct. Dyn. 39 1 2021 357 367 31914861
49 Papaleo E. Mereghetti P. Fantucci P. Grandori R. De Gioia L. Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: the myoglobin case J. Mol. Graph. Model. 27 8 2009 889 899 19264523
50 Wang E. Sun H. Wang J. Wang Z. Liu H. Zhang J.Z.H. Hou T. End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design Chem. Rev. 119 16 2019 9478 9508 31244000
51 Schrödinger Release 2023-3 2023 LigPrep, Schrödinger, LLC New York, NY
