Novel guanosine derivatives against MERS CoV polymerase: An in silico perspective

Abstract The Middle East Respiratory Syndrome Coronavirus (MERS CoV), also termed camel flu, is a new viral infection that first reported in the year 2012 in the Middle East region and further spread during the last seven years. MERS CoV is characterized by its high mortality rate among different human coronaviruses. MERS CoV polymerase shares more than 20% sequence identity with the Hepatitis C Virus (HCV) Non-structural 5b (NS5b) RNA dependent RNA polymerase (RdRp). Despite the low sequence identity, the active site is conserved between the two proteins, with two consecutive aspartates that are crucial in the nucleotide transfer reaction. In this study, seven nucleotide inhibitors have been tested against MERS CoV RdRp using molecular modeling and docking simulations, from which four are novel compounds. Molecular Dynamics Simulation for 260 nanoseconds is performed on the MERS CoV RdRp model to test the effect of protein dynamics on the binding affinities to the tested nucleotide inhibitors. Results support the hypothesis of using the anti-polymerases (Anti-HCV drugs) against MERS CoV RdRp as a potent candidates. Besides four novel compounds are suggested as a seed for high performance inhibitors against MERS CoV RdRp. Communicated by Ramaswamy H. Sarma


Introduction
The Middle East Respiratory Syndrome coronavirus (MERS CoV) is a new viral infection that was reported in the Kingdom of Saudi Arabia in the year 2012 for the first time (Zaki et al., 2012). MERS CoV belongs to a group of viruses that are termed human coronaviruses (Raj et al., 2014;van den Brand et al., 2015). MERS CoV has a flat spread rate (total number of infections is 2494 until today) but with a high mortality rate (34.3%) (Sharif-Yakan & Kanj, 2014).
Severe Acute Respiratory Syndrome (SARS) was the latest human coronavirus reported before MERS CoV characterized by severe acute pneumonia (Stadler et al., 2003). SARS CoV was reported for the first time in 2002 in China (Y. Guan et al., 2003). The main differences between SARS and MERS CoVs are the mortality and spread rates (Coleman & Frieman, 2014). SARS is characterized by a low mortality rate (10%), although it has a higher spread rate compared to MERS CoV (Fouchier et al., 2004;Qinfen et al., 2004). HKU1, OC43, NL63, 229E, SARS, MERS CoV, and SARS-CoV-2 are the seven human coronaviruses strains that were recorded in the last 70 years (Zumla et al., 2015). Except for SARS CoV, MERS CoV, and SARS-CoV-2 all human coronaviruses are of low epidemic importance outbreaks (Fouchier et al., 2004;Zumla et al., 2015). Human coronaviruses are zoonotic viruses. MERS CoV infects human through the dromedary camel while SARS CoV is hosted in palm civet cat before it is transmitted to humans. The newly evolved SARS-CoV-2 is suggested to be hosted by an unknown animal before jubing to humans as well. Infection can occur from animal to animal and from animal to human in the case of close contact with infected animals (Azhar et al., 2014;Han et al., 2016). Additionally human to human transmission was reported for human coronaviruses (Elfiky, 2020a(Elfiky, , 2020bIbrahim et al., 2020;Yang, 2020). SARS, MERS, and SARS-CoV-2 caused diseases are characterized by a lower respiratory ailment like bronchitis, bronchiolitis, and pneumonia, which leads to death with different ratios (Bogoch et al., 2020;Graham et al., 2013;van den Brand et al., 2015;Wu et al., 2020).
During the past two decades, HCV has been profoundly studied, and several anti-HCV drugs have been either approved or under clinical trials (Gonzalez-Grande et al., 2016;Mayhoub, 2012;Yang et al., 2011). Computational methods such as molecular docking and molecular dynamics simulations represent a powerful tool to mimic the properties of biological molecules (Ganesan & Barakat, 2017;Leach, 2001). In this study, three of the anti-polymerase drugs used as inhibitors for the HCV NS5B RdRp are tested against MERS CoV RdRp using a computational approach. The tested compounds include sofosbuvir (approved against HCV in 2013), ribavirin (wide acting antiviral), and IDX-184 (tested in clinical trials). Besides, four novel guanosine derivatives are also tested and compared to the three drugs and the parent guanosine nucleotide. All compounds have been tested against a MERS CoV RdRp model built in silico and equilibrated by Molecular Dynamics Simulation (MDS).

Sequence retrieval and alignment
MERS RdRp structure has not been solved experimentally yet. Therefore, we utilized a molecular modeling approach to construct the all atoms 3D structure of MERS CoV RdRp. The protein database of the National Center for Biotechnology Information (NCBI) was used to retrieve the sequences for the polymerases of all human coronaviruses (MERS, SARS, HKU1, OC43, NL63, and 229E) (NCBI, 2020). Multiple sequence alignment was done using CLUSTAL Omega web server (Sievers et al., 2011) to reveal the sequence conservations among the downloaded sequences for human coronaviruses and HCV polymerase sequences (PDB ID: 2XI3). ESPript 3.0 software is utilized to prepare the multiple sequence alignment (Robert & Gouet, 2014). Structural alignment of the MERS CoV RdRp model and HCV polymerase structure (PDB ID: 2XI3) was done by the aid of Chimera software (Pettersen et al., 2004) (Root Mean Square (RMS) difference of 2.7 Å).

Structure prediction and docking study
I-TASSER web server was used in this study to build the allatoms 3D structure of MERS CoV polymerase from the sequence (ID AHY61336.1) (Yang et al., 2015). Different protein modelling webservers were used to build the 3D structure of MERS HCoV RdRp, while the model built by I-I-TASSER was the best model based on structural validation servers . The structure was validated using the Ramachandran plot, ERRAT, PROVE, and verify-3D software from Structural Analysis and Verification Server (SAVES) (Hooft et al., 1996;Laskowski et al., 1996;SAVES, 2020). Guanosine triphosphate (GTP), Uridine triphosphate (UTP), IDX-184 (GTP derivative), sofosbuvir (UTP derivative), ribavirin (wide acting antiviral drug), and four suggested guanosine derivatives (Elfiky & Elshemey, 2018) were sketched using SCIGRESS 3.0 tools (Summers et al., 2012). The structures were optimized classically using the MM3 force field (Lii & Allinger, 1989) then were further optimized using semi-empirical parameterization methods 6 (PM6) (Stewart, 1991). Finally, the quantum mechanical density functional theory (DFT) was used to optimize the ligands' structure (Becke, 1993). The quantum mechanical functional B3LYP was also used to calculate the infrared transition spectra of the optimized ligands to ensure reality (Saleh et al., 2014).
AutoDock Vina (Morris et al., 2009;Trott & Olson, 2009) was employed in this study to assess the binding affinities and possible binding modes of the interactions between the ligands and MERS CoV RdRp. Four nucleotide inhibitors (based on anti-HCV drugs (guanosine inhibitors)) are utilized in this study. Sofosbuvir, IDX-184, and ribavirin were also tested against MERS CoV polymerase. AutoDock Tools (ADT) software is used to prepare both the small molecules and the protein 3D-structures for the docking experiment. The grid box was set to be 30 Â 30 Â 30 Å and its center is selected to be between the residues, D255 and D256. Flexible ligand in a flexible active site docking approach is used in this study. Moreover, the Vina scoring function is applied to score the resulting complexes. The docking study is conducted using different conformations of the protein corresponding to the protein at different dynamical states (every 10 ns) during the Molecular Dynamics Simulation (MDS) run (Leach, 2001).

Molecular dynamics simulation study
To ensure the binding of the ligands within the MERS CoV RdRp, we used molecular dynamics simulation for 260 nanoseconds to ensure the equilibration of the protein system since any changes in the structure can affect the small molecule binding. NAMD software (Phillips et al., 2005), installed in the Cyprus Institute of Science supercomputing facility, is utilized applying CHARMM force field (Small and MacKerell Jr, 2015). Before the simulation, MERS CoV RdRp was solvated using a TIP3P water model at pH 7. Two Magnesium ions were fixed to the active site to resemble the active site conformation. The coordinates of these two ions were taken from the solved polymerase structure (PDB ID 2XI3). The total charge of the protein system was neutralized by adding ten chlorine ions (Noorbatcha et al., 2010). So, the simulation mimics the protein in the aqueous environment of the host cells.
Before the equilibration, the water and ions were minimized for 10000 steps followed by 100 ps MDS. After that, another 10000-step minimization of the whole system (water, ions, and the protein) was performed. Equilibration of the system was performed for 5 ns at normal pressure (1 atm) and temperature (310 K) (NPT ensemble). Periodic Boundary Conditions (PBC) were used with a simulation box of size 85.3 Â 85.3 Â 85.3 Å and the box center is (64.63, 64.00, 65.78) Å as calculated from the equilibration simulation period at NPT ensemble.
A production run for 260 ns in constant volume and temperature (NVT ensemble) was conducted on the Cy-Tera supercomputing facility of the Cyprus Institute of science (Project number pro15b114s1). NAMD and VMD software (Humphrey et al., 1996;Phillips et al., 2005) of the University of Illinois, NIH Center for Macromolecular Modeling and Bioinformatics, Theoretical and Computational Biophysics group, were used to prepare the structure, run the simulation and analyze the trajectories. Every ten ns of the MDS, the protein coordinates were extracted from the trajectory file to be used in the docking experiment. The first 44 ns part of the production run is excluded (to ensure that the information is taken from a full equilibration system). A total of 23 different conformations of MERS CoV RdRp were employed in the docking study.

Multiple sequence and structural alignment
Polymerase active site is conserved in many organisms (Doublie & Ellenberger, 1998;Ferrer-Orta et al., 2015;Gonzalez-Grande et al., 2016;J acome et al., 2015;Mayhoub, 2012). Figure 1A shows the multiple sequence alignment of the MERS CoV RdRp along with other human coronavirus' polymerases. The hepatitis C Virus polymerase sequence and secondary structure (PDB ID 2XI3) were also included for comparison. From the alignment, we can deduce the polymerase sequences conservation, particularly around the active site aspartates (between b9 and b10). Structural conservation is also apparent from the superposition of MERS CoV (comparatively modeled structure ) and HCV (PDB ID: 2XI3) polymerases, as shown in Figure 1B. The b-hairpin structural fold (b9 and b10 of Figure 1A) is conserved structurally. The two consecutive aspartate residues D255 and D256 (red-colored) of the MERS CoV polymerase protrude from the beta-turn with the same orientation of that for HCV polymerase (green colored D318 and D319) which suggests the possibility of inhibiting MERS CoV with anti HCV drugs.

Molecular dynamics simulation
The MERS CoV comparative model that was generated using I-TASSER web server is valid based on the Ramachandran plot (Laskowski et al., 1996) (90.6% in the most favored region, while only 3.2% in the disallowed area) and ERRAT software (Hooft et al., 1996) (overall error factor 75.6%). This values are very good in terms of the low sequemce identity between the target protein and the best homologous solved structure. Additionally, we performed long MDS run to ensure its equilibration and the reliability of the binding affinity data. Dramatic changes in the dynamics of protein subdomains, even apart from the active site, can alter small molecule binding. It was reported that the finger domain can move close to the palm domain during nucleotide addition to the RNA primer and hence we have to simulate the dynamics for a longer time (Elfiky & Ismail, 2019). Figure 2A shows the Root Mean Square Deviation (RMSD in Å) (blue line), Radius of Gyration (Rodrigues & Bonvin, 2014) (RoG in Å) (orange line), Surface Accessible Surface Area (SASA in Å 2 ) (gray line) and the number of H-bonds (yellow line) versus time in nanoseconds. The 260 nanoseconds were enough for the system to be equilibrated, as shown by the saturation of the RMSD curve, which reached after ̴ 30 ns. Besides, RoG, SASA, and number of H-bonds are stable during the simulation period with distributed values, as shown in panel B of the figure (the same coloring scheme of Figure 2A). The RMSD values of the equilibrated system were around 7.75 Å while it was a little bit smaller at the beginning of the simulation (before 100 ns). On the other hand, RoG and SASA values were around 22 Å and 20000 Å 2 at the beginning of the simulation but these values are reduced during the simulation to ̴ 21 Å and 18000 Å 2 , respectively as shown in Figure  2B. This indicates a slight reduction of the surface accessible area and protein radius of gyration at the end of the simulation compared to the first 100 ns of the simulation. Figure 2C shows the per residue-Root Mean Square Fluctuations RMSF. Four regions of highly fluctuating residues (plus the N and C-terminal regions) are reported here (RMSF value as much as 8 Å). The highly fluctuating regions (red-colored cartoons in the right panel and red surface in the left panel structure, Figure 2C) are S73-G91, C140-Y148, V162-G173, and N206-K216, in addition to the N-terminal region (F1-E18) and the C-terminal region (A292-F307). All the highly fluctuating regions are loops connecting the secondary structural motifs (S73-G91 connecting a2-a3, C140-Y148 connecting a4-a5, V162-G173 connecting a5-a6, while N206-K216 connecting a6-b3). On the other hand, the active aspartates, D255 and D256, (magenta sticks and magenta surface in the right and left panel structures in Figure 2C) show a low level of fluctuations (RMSF > 2 Å). Almost all the highly fluctuating regions are accessible surface loops (see the red colored surface representation in Figure 2C) while the embedded loops (D113-N123, C260-S263 and the active site b-turn (L253-D256)) show a low level of fluctuations (RMSF > 2.5 Å).

Molecular docking
MERS CoV polymerase conformations after 44 ns were used in the docking experiment to ensure the system's equilibration. A total of 23 different conformations were used with time intervals of 10 ns (see Figure 2D). As shown in the figure, all the ligands were able to bind to the polymerase active site with good binding affinities. The worst value was recorded for UTP at 94 ns (-5.9 kcal/mol) while, the best value was reported for compound 3 at 174 ns (-8.8 kcal/mol). Figure 3A shows the 2D structures for the four novel guanosine derivatives used in this study. Figure 3B shows the average binding energies calculated by AutoDock Vina for the different small molecules to the different conformations of MERS CoV polymerase. GTP and UTP (blue) are used here as a positive control. Sofosbuvir, IDX-184, and ribavirin (red) are compared to the four suggested compounds (green). As expected, based on sequence and structural conservation, all the studied compounds can fit in the active site of MERS CoV polymerase with a minimum average binding affinity of -7.13 Kcal/mol reported for ribavirin. The four suggested compounds have similar average binding affinities compared to the physiological parent nucleotide (GTP). So, it can compete for MERS CoV polymerase active site with GTP and hence inhibit the polymerase function of MERS CoV.  Besides, these four compounds have slightly lower (better) average binding affinities compared to the antiviral drugs ribavirin and sofosbuvir. Both sofosbuvir and IDX-184 have comparable average binding energies of their physiological molecules from which it was developed (UTP and GTP, respectively).
In order to understand how the interactions are established between the compounds and the MERS HCoV RdRp, we randomly selected two sets of complexes between the compounds and MERS HCoV RdRp and performed in depth analysis of the docking complexes. Tables 1 and 2 show the number of H-bonds formed between the different ligands and the protein active site pocket after performing the docking at 64 ns and 154 ns conformations of MERS CoV polymerase. Also, the amino acids involved in H-bond formation are listed in the tables. The conformations are selected randomly to represent the protein at two different dynamics states, after the equilibration period. As we can conclude from the tables, the number of H-bonds formed is different (from 2 to 7). Additionally, the active site amino acids (D255 and D256) are found to be involved in the H-bonding interaction or the metal (two Mg þ2 ) interactions in all the tested ligands. This is the same mode of interaction reported in other viral polymerases like Hepatitis C Virus and Zika virus Elfiky & Elshemey, 2018;Elfiky & Ismail, 2018J acome et al., 2015;Mayhoub, 2012). Other residues are lining the active site pocket also involved in H-bonding and metal interactions such as D113 and K301 in both conformations of the polymerase (64 and 154 ns conformations). On the other hand, R48 and K302 are slightly involved in the interactions in 64 ns and 154 ns conformations of MERS CoV polymerase, respectively. Metal interactions are reported in all the docking experiments (see Tables 1 and 2). These interactions facilitate the binding of the ligands to the polymerase active site as reported before for other polymerases Elfiky & Elshemey, 2018). Figure 4 shows the docking complexes formed with the protein at 64 ns conformation for GTP and the four suggested guanosine inhibitors. The binding mode of the physiological molecule (GTP) is almost the same as that of the guanosine derivatives. The active site residues are the main site for H- bond formation with the triphosphate groups of the ligands, while the orientation of the guanosine moieties shows a more complicated pattern. This is principally due to the modifications that performed to these guanosine derivatives at 2 0 position of the ribose ring. Compound #1 has (2-hydroxyphenyl)oxidanyl, compound #2 has (3,5-dihydroxyphenyl)oxidanyl, compound #3 has (3-hydroxyphenyl)oxidanyl, while compound #4 has (3-sulfanylphenyl)oxidanyl (Elfiky, 2017). These added groups made new sites for the interaction with active site cavity lining residues (see Tables 1 and 2).
In addition to the H-bonding and metal interactions, water plays an essential role in polymerase function (Bellissent-Funel et al., 2016). Many water molecules complete the coordination of the magnesium ions present in the vicinity of the active site. This facilitates the binding of the ligand to MERS CoV polymerase active site.

Conclusion
As polymerase structure and function is highly conserved, it is possible to target MERS CoV polymerase with anti-virals developed for other viral polymerases. Due to the momentum of research on HCV in the last two decades, a lot of small molecules inhibitors have emerged against the viral polymerase. In this mingled molecular modeling, docking and dynamics study, we demonstrate the ability of some anti-HCV drugs to bind to and consequently inhibit MERS CoV polymerase function. Besides, four suggested guanosine  derivatives are introduced as plausible powerful MERS CoV polymerase blockers compared to ribavirin.