key: cord-0848601-e4j2k7gq
authors: Elfiky, Abdo A.; Mahran, Hanan A.; Ibrahim, Ibrahim M.; Ibrahim, Mohamed N.; Elshemey, Wael M.
title: Molecular dynamics simulations and MM-GBSA reveal novel guanosine derivatives against SARS-CoV-2 RNA dependent RNA polymerase
date: 2022-01-20
journal: Rsc Advances
DOI: 10.1039/d1ra07447d
sha: 3ec373de42066ba130c1f5db1ed3e623c5f4dc32
doc_id: 848601
cord_uid: e4j2k7gq

According to the World Health Organization (WHO), SARS-CoV-2 is responsible for more than 5 M deaths and is reported in 223 countries infecting 250+ M people. Despite the current vaccination momentum, thousands of people die every day by COVID-19. Suggesting possible blockers of the viral RNA-dependent RNA polymerase is highly needed for potential effective therapeutics against SARS-CoV-2. This study utilizes combined molecular dynamics simulation and molecular docking to test novel guanosine derivatives against SARS-CoV-2 RdRp. Results reveal the binding potency of nineteen guanosine derivatives against SARS-CoV-2 solved structures. The bulky moieties (hydroxyl or fluorated phenyl moieties) added to the 2′ position of the ribose ring positively impacted the binding affinity to RdRp. The current in silico study represents a one-step-ahead for suggesting new possible blockers of SARS-CoV-2 RdRp that are yet to be verified in the wet lab. It offers new potential binders or blockers of RdRp that bind to the protein active site tighter than remdesivir. The latter was approved by the food and drug administration (FDA) for emergency use against COVID-19 last year.

The viral pandemic Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has killed over 5 million people in 223 countries since it was rst reported in Wuhan, China, in December 2019. According to the WHO website, the infection and death numbers are continuing to grow (https:// www.who.int/emergencies/diseases/novel-coronavirus-2019, accessed 12 November 2021). 1, 2 The SARS-CoV-2 virus is a new strain of betacoronaviruses that causes the deadly COVID-19 disease. Betacoronaviruses are a family of coronaviruses (crown-like virus) that can cause respiratory and gastrointestinal illnesses in humans and animals, unlike gammacoronaviruses and deltacoronaviruses, which infect birds and possibly mammals but have never been linked to human illness. 3 

Syndrome coronavirus (MERS-CoV) (2012) are two other betacoronaviruses that were reported to cause human epidemic diseases in the past two decades. 1, [4] [5] [6] The SARS-CoV-2 genome contains two large genes, ORF1a and ORF1b, which code for 16 non-structural proteins (NSP1-NSP16), including NSP12, which codes for RNA-dependent RNA polymerase (RdRp). 4 RdRp is a crucial enzyme for RNA viruses, including the human coronaviruses. It promotes viral gene transcription and replication in collaboration with other viral (nsp7 and nsp8) and host-cell factors. 7, 8 It has no human homolog, so this paves the way for antiviral drug production and lowers the chances of a protein being harmed in human cells. 9, 10 RdRp was targeted in different RNA viruses such as Hepatitis C Virus (HCV), Ebola Virus (EBOV), Zika Virus (ZIKV), and human coronaviruses (HCoVs). [11] [12] [13] [14] [15] [16] [17] Many directly acting antiviral drugs are officially approved against different viruses, such as sofosbuvir and ribavirin against RdRp of HCV and remdesivir against EBOV. 11, 18 These drugs are modied nucleotides that compete with physiological nucleotides for the RdRp active site. Ribavirin is a wide-acting antiviral drug used to combat a wide range of viruses in combination with immunomodulators or other antiviral medications, while remdesivir is a nucleotide analog used against EBOV, Marburg virus, MERS-CoV, and SARS-CoV. [19] [20] [21] [22] Ribavirin reported a half-maximal Effective Concentration (EC50) against SARS-CoV-2 of 109.5 mM, while remdesivir reported an EC50 of 1.76 mM. 23, 24 This study focuses on nucleotide inhibitors because of the strong evidence of inhibition against viral RdRps. Aer Molecular Dynamics Simulation (MDS) for up to 50 ns for each structure, seventy-six guanosine triphosphate derivatives are tested against ten different solved structures of SARS-CoV-2 RdRp. Aer cluster analysis of the RdRp structures' trajectories, molecular docking is utilized to test the binding affinity and the guanosine derivatives' compulsory mode against SARS-CoV-2 RdRp. Finally, the top hits are subjected to additional 50 ns MDS followed by binding free energy calculations.

Ten different SARS-CoV-2 RNA-dependent RNA polymerase structures were downloaded from the Protein Data Bank (PDB) database. 25 The PDB ID codes for the structures were 6M71,  6XEZ, 6XQB, 6YYT, 7BTF, 7BV1, 7BV2, [26] [27] [28] [29] [30] [31] The structures were solved experimentally using the cryo-EM technique, yet they maintained their high resolution. The PyMOL soware was utilized to prepare the structures for docking and molecular dynamics simulation studies. 32 Water and ligands were removed, while any missing hydrogen atoms were added to the structures. The structures that contain other cofactors were prepared by removing any cofactors or nucleic acids before the simulation studies.

The structures of the ligands used in this study were drawn using SCIGRESS (version FJ 2.9.1) soware. 33 The ligands were modied guanosine derivatives based on the guanosine triphosphates (GTP) modied at position 2 0 while the hydroxyl group was replaced by a bulky group, as shown in Fig. 1 . The structures were geometry optimized rst using the classical mechanical force eld (MM3) then using the semi-empirical parameterization method 6 (PM6) before the docking study. 34, 35 Molecular dynamics simulation The ten SARS-CoV-2 RdRp structures were examined using 50 ns molecular dynamics simulation runs with the help of the Nanoscale Molecular Dynamics (NAMD) 2.14 soware (CHARMM 36 force eld). 36, 37 Visualizing Molecular Dynamics (VMD) 1.9.3 soware was utilized to build the input les and analyze the trajectories along with some in-house codes. 38, 39 Before the production run, the structures were solvated using the TIP3P water model, while Na and Cl ions were added to keep the system's charge neutral. 40 Additionally, the system was minimized for 10 000 steps using the steepest descent algorithm, followed by gradually increasing the system temperature until the physiological temperature of 310 K at the equilibration step. While the number of atoms, the pressure, and temperature were kept constant (NPT) during the equilibration period (1.1 ns), the volume was adjusted and kept constant (NVT) during the production runs (50 ns).

Aer MDS, cluster analyses of the trajectories were performed using TTClust soware. 41 Two clusters were dened for each protein, and a representative conformation was randomly selected from each cluster for the docking study. A total of twenty different conformations of the RdRps were used to test the guanosine derivatives.

Aer MDS trajectory clustering, the docking experiment was performed using the twenty structural conformations of SARS-CoV-2 RdRps with the help of AutoDock Vina soware. 42 Seventy-six guanosine derivatives were tested against the different conformations of the SARS-CoV-2 RdRps. Guanosine triphosphate (GTP), sofosbuvir, ribavirin, and remdesivir were used as positive controls, while cinnamaldehyde and thymoquinone were selected as negative controls. 8,15 AutoDock Tools 1.5.6 was utilized to prepare ligands and the different protein conformations (merging Gasteiger charges and generating dockable PDBQT format for molecular docking). 43 The receptor grid boxes for the twenty conformations of SARS-CoV-2 RdRp were adjusted to cover the binding site represented by D760 and D761, which were treated as exible during the docking simulations. The box size was set to 30 Â 30 Â 30Å 3 centered between the D760 and D761 residues that have different coordinates in each protein conformation. The calculations were done on eight cores, and the default exhaustiveness value of 8 was set in all docking calculations.

Aer docking, SARS-CoV-2 RdRp structures were mined for the established interactions between ligands and proteins using the Protein-Ligand Interaction Proler (PLIP) web server (Technical University of Dresden) and then tabulated. 44, 45 MDS of the best two complexes and free energy calculations

The best two complexes in binding RdRps (C13 and C14) were subjected to 100 ns MDS run following the same protocol using NAMD 2.13 soware. This was followed by Molecular Mechanics Generalized Born Surface Area (MM-GBSA) calculation using MMPBSA.py script implemented in amber tools 17. 46 All frames (1000) were used in the calculations. The salt concentration was set to 0.154 M, while all other settings were kept as default.

2 0 analog as a potent nucleotide derivative

The nucleotide modied at its 2 0 position of the ribose ring was found successful as an inhibitor against different viral polymerases, including Hepatitis C virus (HCV), Zika virus, and human coronaviruses (such as MERS-CoV, SARS-CoV, and SARS-CoV-2). 8, 12, 47, 48 Sofosbuvir, for example, was approved by the Food and Drugs Administration (FDA) against HCV NS5B RdRp eight years ago. 15, 49, 50 In addition, it showed excellent in silico binding affinity against other viruses and pathogenic fungi. 17,51-53 In Fig. 1 , the guanosine derivatives giving good binding energies are listed with the modications added at the ribose 2 0 position instead of the hydroxyl group in the parent compound Guanosine Triphosphate (GTP). The added bulky groups (phenyl or OH, SH, F substituted phenyl moieties) will interfere with RNA recognition aer the addition of the guanosine derivative to the primer RNA at the viral RdRp active site.

The ten RdRp solved structures are equilibrated for 50 ns MDS as shown from the average Root Mean Square Deviation (RMSD) curve of Fig. 2A . The blue curve is the average over the ten RMSD curves for individual protein structure versus time, while the error bars represent the standard deviation values. The average RMSD reects that the protein systems are equilibrated during the rst 30 ns, and the average RMSD curve reaches a saturation value of 2.75Å. The protein systems are also stable during the dynamics simulation as deduced from the average per-residue Root Mean Square Fluctuations (RMSF) curve in Fig. 2B . Again, the blue curve is the average RMSF, while the error bars are the standard deviation values. The majority of the protein residues are stable during the dynamic's simulation except for the highly uctuating N and C terminals (about 50 residues at the beginning and 50 at the end of the RMSF curve) which reach up to <5Å. The active site residues D760 and D761 exhibit minor uctuations (average RMSF up to 1.13Å), while the regions F334-V338 (red), V359-S364 (yellow), Q822-D825 (blue), and D846-M855 (magenta) show moderate average RMSF values (<2 A). These four regions are apart from the active site D760 & D761 (black sticks) and represent movable turns and loops as shown in the structure of SARS-CoV-2 RdRp (PDB ID: 6XEZ) in Fig. 2B .

The 50 ns trajectories of the ten RdRps are subjected to cluster analysis using TTClust, and two clusters are dened for each structure. The binding affinities of the GTP and its 76 derivatives against SARS-CoV-2 RdRps conformations are calculated using the AutoDock Vina soware. Fig. 3 shows the average binding affinities for the best 19 derivatives in blue alongside the GTP in green, the positive controls (sofosbuvir, remdesivir, and ribavirin) in yellow, and the negative controls (cinnamaldehyde and thymoquinone) in red. The colored bars represent the average values (in kcal mol À1 ), while the error bars represent the standard deviations. One can notice that the guanosine derivatives have average binding affinities comparable to that of the approved drug remdesivir (À7.70 AE 0.58 kcal mol À1 ) and to those of the other positive controls sofosbuvir (À7.55 AE 0.43 kcal mol À1 ) and ribavirin (À7.35 AE 0.48 kcal mol À1 ). This indicates that these guanosine derivatives could be successful therapeutic agents against SARS-CoV-2 RdRp as its binding is better than or equal to the approved drug remdesivir. We discussed in a previous study the appropriate use of anti-HCV, anti-HIV, and anti-EBOV drugs against SARS-CoV-2 RdRp, but our study was based on in silico built models of the RdRp. 9,10 Fortunately, the current research reproduced almost the same binding affinity for remdesivir and sofosbuvir using ten different solved structures of the RdRp. This indicates the effectiveness of our previous model of RdRp to test anti-COVID-19 therapeutics. Additionally, we compared the remdesivir and sofosbuvir against RdRp of SARS-CoV-2 in an earlier study, and it gave the same binding affinity to the current research. 8 Therefore, even the interacting residues are conserved in the case of SARS-CoV-2 RdRp against remdesivir and sofosbuvir.

As shown in Table 1 , the nineteen guanosine derivatives successfully establish strong contacts with the active site of the RdRp of SARS-CoV-2. Hydrogen bonds or salt bridges' formation is the primary type of interaction found between the protein and the guanosine derivatives, with few hydrophobic (7), N691 (5), and D618 (4). Simultaneously, the most reported residue that forms hydrophobic contact with the guanosine derivatives is the active site residue D760 (8 connections). Additionally, other residues from RdRps contribute to developing four (Y619 and D761) or three (D623 and N691) hydrophobic contacts with the guanosine derivatives.

The best two guanosine derivatives are the C13 (3,5-dihydroxy phenyl substitution) and C14 (3-hydroxyphenyl substitution) with average binding affinities of À7.94 AE 0.58 and À7.96 AE 0.43 kcal mol À1 , respectively. Fig. 4 shows a graphical representation of the interactions between the best two compounds and selected RdRps. The RdRp residues are labeled by their oneletter codes and represented in blue sticks, while the guanosine derivatives are in yellow sticks. H-Bonds and salt bridges are depicted by blue lines and dashed-yellow lines, respectively, Table 1 The interactions that were established upon docking the guanosine triphosphate and its derivatives into the SARS-CoV-2 RdRp active site. Docking scores are calculated using AutoDock Vina software, while PLIP is used to get the interactions established. The interactions of positive (sofosbuvir, remdesivir, and ribavirin) and negative (cinnamaldehyde and thymoquinone) control compounds are listed as well in the (2),

À7.9 10 Y619, K621, R624(2), S682, N691, S759(2), D760, D761 while dashed-gray lines represent the hydrophobic contacts.

Notably, the added bulky moiety at the 2 0 position of the ribose contributes to the binding against the protein's active site pocket residues.

MDS for the best two complexes (C13 and C14 against RdRp) and the GTP-RdRp complex are performed for 100 ns using NAMD 2.13 soware, followed by MM-GBSA calculations for the same complexes using Amber tools. The residual contribution for the binding of RdRp against the GTP-RdRp complex and the best two compounds (C13 and C14) are listed in Table 2 Thus, the total binding energy for the RdRp-C14 complex is lower (À46.95 AE 6.3 kcal mol À1 ) than RdRp-C13 and the RdRppositive control (GTP) (À5.14 AE 6.6 and À13.82 AE 6.2 kcal mol À1 , respectively); hence it is suggested to be the best SARS-CoV-2 RdRp potential inhibitor in the current study. Therefore, we recommend testing this compound experimentally against SARS-CoV-2 RdRp. Fig. 5 supports the docking data as reected from the attened RMSD curves (Fig. 5A ) of the GTP and the two complexes (C13 and C14) versus SARS-CoV-2 RdRp. The three complexes were equilibrated during the 100 ns as reected from the attened RMSD (A), the equilibrated Radius of Gyration (RoG) (B), and Surface Accessible Surface Area (SASA) (C) curves. The systems RMSD, RoG, and SASA mean values are equilibrated at around 2.5Å, 29.7Å, and 40 000Å 3 , respectively. The GTP-RdRp complex possesses elevated RoG and SASA values compared to the C13 and C14 complexes. The total number of the formed Hbonds is almost stable around 1300 during the simulation period for C13 and C14 complex but 1500 for the GTP complex, as shown in Fig. 5D . The per-residue RMSF is shown in Fig. 5E for the two complexes, where C13-RdRp is represented by the orange line, while C14-RdRp is represented by the gray one. Apart from the N and C-terminals, the RMSF show high uctuations (RMSF < 3) in the D469-G480 region in the case of the C13-RdRp complex compared to the C14-RdRp. Additionally, the C14-RdRp complex shows lower RMSF values in 460 and 560 residues and is more stable than the C13-RdRp complex. These regions are apart from the active site aspartates (D760 & D761) that are shown in black sticks in the structure at the top of the RMSF curves. So, C14-RdRp is a more stable complex than C13-RdRp complex, which supports our calculated binding energies. Table 2 The MM-GBSA calculations for the best two complexes (C13 and C14 complexes with RdRp) after 100 ns MDS. Red coloured residues represent the residue have negative contribution on the binding (positive binding energies). The average binding free energies and its individual terms are shown at the bottom for each complex with its standard deviations Fig. 5B , C, and D reect that the GTP complex has higher RoG, SASA, and H-bonds. This, in fact, is due to the protein size used in this complex which is 925 residues (PDB ID: 6XEZ), while for C13-RdRp and C14-RdRp complexes, the protein size is 833-residues (PDB ID: 7BV1). In general, all the nineteen guanosine derivatives presented in Fig. 3 could be successful candidates against SARS-CoV-2 RdRp, while compound C14 is the best in binding SARS-CoV-2 RdRp based on the calculated MM-GBSA and the RMSF curves. Therefore, C14 should be tested experimentally to prove its binding affinity and binding mode against the SARS-CoV-2 RdRp active site.

The pandemic SARS-CoV-2 is still causing everyday life casualties since its emergence in December 2019. The availability of solved structures for SARS-CoV-2 RdRp allows us to test a bunch of new compounds, based on guanosine triphosphate, as possible RdRp blockers. Based on the current computational study, nineteen of the guanosine derivatives show promising binding affinities against SARS-CoV-2 RdRp solved structures (À7.83 AE 0.08 kcal mol À1 ) comparable to that of remdesivir (À7.70 AE 0.58 kcal mol À1 ), sofosbuvir (À7.55 AE 0.43 kcal mol À1 ), and ribavirin (À7.34 AE 0.49). The added bulky groups at the 2 0 position of the ribose ring positively impact the binding affinity. C14 (3-hydroxyphenyl substitution at 2 0 position) shows the tightest binding affinity to the RdRp using MM-GBSA calculations and hence is suggested as a potent binder to the viral protein. However, this is yet to be conrmed experimentally. 

All the authors declare that there is no conict of interest for this work.

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health -The latest 2019 novel coronavirus outbreak in Wuhan, China

Pneumonia of unknown aetiology in Wuhan, China: potential for international spread via commercial air travel

A Review of Human Coronaviruses' Receptors: The Host-Cell Targets for the Crown Bearing Viruses

Quantitative structure-activity relationship and molecular docking revealed a potency of anti-hepatitis C virus drugs against human corona viruses

Middle East respiratory syndrome coronavirus: another zoonotic betacoronavirus causing SARS-like disease

East respiratory syndrome coronavirus (MERS-CoV)

RNA-dependent RNA polymerase: Structure, mechanism, and drug discovery for COVID-19

The anti-HCV, Sofosbuvir, versus the anti-EBOV Remdesivir against SARS-CoV-2 RNA dependent RNA polymerase in silico

Ribavirin, Remdesivir, Sofosbuvir, Galidesivir, and Tenofovir against SARSCoV-2 RNA dependent RNA polymerase (RdRp): A molecular docking study

Molecular dynamics and docking reveal the potency of novel GTP derivatives against RNA dependent RNA polymerase of genotype 4a HCV

Novel Guanosine Derivatives as Anti-HCV NS5b Polymerase: A QSAR and Molecular Docking Study

Applications of Computer-Aided Approaches in the Development of Hepatitis C Antiviral Agents

Molecular Modeling and Docking revealed superiority of IDX-184 as HCV polymerase Inhibitor

IDX-184 is a superior HCV direct-acting antiviral drug: a QSAR study

Molecular dynamics simulation revealed binding of nucleotide inhibitors to ZIKV polymerase over 444 nanoseconds

Zika virus: novel guanosine derivatives revealed strong binding and possible inhibition of the polymerase

Anti-HCV drugs in the pipeline

The clinical pharmacology of ribavirin

Peginterferon alfa-2b plus ribavirin compared with interferon alfa-2b plus ribavirin for initial treatment of chronic hepatitis C: a randomised trial

Nucleotide polymerase inhibitor sofosbuvir plus ribavirin for hepatitis C

Therapeutic efficacy of the small molecule GS-5734 against Ebola virus in rhesus monkeys

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Broadspectrum antiviral activity of the IMP dehydrogenase inhibitor VX-497: a comparison with ribavirin and demonstration of antiviral additivity with alpha interferon

Announcing the worldwide Protein Data Bank

The cost and benets of segregation (by age)

Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by remdesivir

Structure of the RNA-dependent RNA polymerase from COVID-19 virus

Structure of replicating SARS-CoV-2 polymerase

Ligand docking and binding site analysis with PyMOL and Autodock/Vina

Structural properties of metal-free apometallothioneins

Molecular mechanics. The MM3 force eld for hydrocarbons. 3. The van der Waals' potentials and crystal data for aliphatic and aromatic hydrocarbons

Application of the PM6 semiempirical method to modeling proteins enhances docking accuracy of AutoDock

Scalable molecular dynamics with NAMD

CHARMM36 all-atom additive protein force eld: validation based on comparison to NMR data

VMD: visual molecular dynamics

Host-cell recognition through GRP78 is enhanced in the new UK variant of SARS-CoV-2, in silico

Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K

TTClust: A Versatile Molecular Simulation Trajectory Clustering Program with Graphical Summaries

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

AutoDock4 and AutoDockTools4: Automated docking with selective receptor exibility

PLIP: fully automated protein-ligand interaction proler

PLIP 2021: expanding the scope of the protein-ligand interaction proler to DNA and RNA

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

Hepatitis C RNA-dependent RNA polymerase inhibitors: a review of structure-activity and resistance relationships; different scaffolds and mutations

Molecular docking revealed the binding of nucleotide/side inhibitors to Zika viral polymerase solved structures

Genotype and subtype proling of PSI-7977 as a nucleotide inhibitor of hepatitis C virus

Mechanism of activation of PSI-7851 and its diastereoisomer PSI-7977

The clinically approved antiviral drug sofosbuvir inhibits Zika virus replication

The FDA-approved drug sofosbuvir inhibits Zika virus infection

The antiviral Sofosbuvir against mucormycosis: an in silico perspective

Shaheen supercomputer of King Abdullah University of Science and Technology (KAUST) is utilized to perform the aer dockings MDS calculations (under Project Number k1482). Additionally, the Bibliotheca Alexandrina HPC is used to perform the core MDS calculations. This work is supported by the Cairo University COVID-19 fund received by AE.