Repurposing of drugs against methyltransferase as potential Zika virus therapies

In recent years, the outbreak of infectious disease caused by Zika Virus (ZIKV) has posed a major threat to global public health, calling for the development of therapeutics to treat ZIKV disease. Several possible druggable targets involved in virus replication have been identified. In search of additional potential inhibitors, we screened 2895 FDA-approved compounds using Non-Structural Protein 5 (NS5) as a target utilizing virtual screening of in-silco methods. The top 28 compounds with the threshold of binding energy −7.2 kcal/mol value were selected and were cross-docked on the three-dimensional structure of NS5 using AutoDock Tools. Of the 2895 compounds screened, five compounds (Ceforanide, Squanavir, Amcinonide, Cefpiramide, and Olmesartan_Medoxomil) ranked highest based on filtering of having the least negative interactions with the NS5 and were selected for Molecular Dynamic Simulations (MDS) studies. Various parameters such as RMSD, RMSF, Rg, SASA, PCA and binding free energy were calculated to validate the binding of compounds to the target, ZIKV-NS5. The binding free energy was found to be −114.53, −182.01, −168.19, −91.16, −122.56, and −150.65 kJ mol−1 for NS5-SFG, NS5-Ceforanide, NS5-Squanavir, NS5-Amcinonide, NS5-Cefpiramide, and NS5-Ol_Me complexes respectively. The binding energy calculations suggested Cefpiramide and Olmesartan_Medoxomil (Ol_Me) as the most stable compounds for binding to NS5, indicating a strong rationale for their use as lead compounds for development of ZIKV inhibitors. As these drugs have been evaluated on pharmacokinetics and pharmacodynamics parameters only, in vitro and in vivo testing and their impact on Zika viral cell culture may suggest their clinical trials on ZIKV patients.


Materials and methods
Structure and ligand retrieval. The protein 3D structure of the ZIKV, methyltransferase or NS5 (PDB ID: 5MRK, X-ray, 1.9 Å) 27 was retrieved from the Protein Data Bank (PDB). It is co-complexed with an inhibitor, Sinefungin (SFG) 27 . We compared the ZIKV, methyltransferase (PDB ID: 5MRK) with other reported ZIKV methyltransferase structures and found that all the structures share similar folds. We observed that binding site residue orientation is also similar for all the structures. The RMSD between the 5MRK and 5GOZ, 5M5B, 5WZ1, 5NUJ, and 5NJV was 0.479, 0.362, 0.396, 0.349 and 0.390 Å respectively (Fig. 2). SFG is isolated from the Streptomyces species and it is a natural nucleoside related to S-adenosylmethionine. It has antiviral, antifungal, and antiparasitic activity. It is a well-known antifungal antibiotic and competes with S-adenosyl-1-methionine (SAM), the natural substrate of MTases 28 . SFG has shown intramolecular interactions with ZIKV methyltransferase 27 ; therefore, we utilized SFG as a control compound for comparing binding of the predicted ligands. FDA-approved compounds (n = 2895) were retrieved from the ZINC database (https:// zinc. docki ng. org/) in .mol2 file format. ZINC is a large database and comprises millions of compounds in different subsets.  www.nature.com/scientificreports/ Protein and ligand preparation. Heteroatoms such as Cl -, SFG, and water, etc. were removed from the NS5 protein structure using Chimera 1.13.2 and then the apo-protein structure was employed for energy minimization using Chimera 1.13.2 29 . We used Amber ff99SB force field 30 and the 100 steepest descent steps were run for removing the steric clashes in the protein. The protein structures was converted into .pdbqt file format by using Autodock 31 prior to the virtual screening. All the desired hydrogen atoms were added in the structure and the Kollmann charges were set. The ligands from ZINC database in .mol2 format contains all the hydrogens with proper geometry; therefore, they were directly used for the .pdbqt conversion instead of energy minimization. The ligands were taken and converted using prepare_ligand4.py python script provided in the AutoDock Tools. This automated script adds all the charges, removes the steric clashes, and converts the .mol2 file into .pdbqt file format. Then the protein and ligands (.pdbqt format) were used for the virtual screening through Autodock Vina v.1.2.0 32 using the in-house developed pipeline (https:// github. com/ shukl arohi t815/ pyVSv ina).
High-throughput screening. High-throughput virtual screening (HTVS), a widely used technique to identify the potential inhibitors from a large dataset 33 , avoids the extensive effort of synthesizing individual molecules and the respective cost to identify the lead molecule. Due to invention of computational power, the HTVS technique is able to screen compound datasets within hours. We prepared a grid based on the co-crystallized ligand SFG. The SFG showed the interactions with Ser56, Gly86, Trp87, Lys105, His110, Glu111, Asp131, Val132, Asp146, and Ile147 in the co-crystallized structure 27 . We selected the same binding residues and created a grid (Center_X = 18, Center_Y = 6, Center_Z = 5 and Size_X = 20, Size_Y = 20 and Size_Z = 20) and utilized the same grid for the three docking software analyses. We next screened the FDA-approved compound dataset against the NS5 enzyme.

Molecular docking analysis.
Cross-docking from multiple software is a widely used approach to remove false positive binders 34,35 . We selected 28 compounds from the virtual screening analysis based on their binding properties with the human receptors. These compounds were employed for cross-docking analysis using Autodock Tools (ADT). The Autodock uses the Lamarckian genetic algorithm (LGA) for producing the docking pose. The same grid box, which was set for virtual screening, was utilized for the docking. The grid box was set based on SFG binding residues to these coordinates (Center_X = 18, Center_Y = 6, Center_Z = 5 and Size_X = 20, Size_Y = 20 and Size_Z = 20). In the Autodock tools, the binding affinity was evaluated in the following two steps: first the energy was calculated in an unbound state followed by energy calculation of the NS5-ligand complexes. The difference between the first and second steps was then evaluated. The sum of the freedom of torsional degree was used for the calculation of conformational entropy. The free energy of binding was calculated as: where the protein and ligand are referred to as P and L, pairwise evaluation is denoted by V, and ΔS conf denotes the loss of conformational entropy during binding. One hundred binding poses for each ligand were generated by using the Lamarckian genetic algorithm. The complexes were selected through the binding mode and binding affinity analysis. Finally, the top five compounds were selected based on their interaction and binding energy and the control ligand in complex form were employed for MDS.

Molecular dynamics simulation.
MDS is a widely used technique for predicting protein-ligand complex stability. We performed 150 ns MDS for evaluating the stability of the protein-ligand complexes. We created seven systems (one for apo-NS5, one for NS5-SFG, and the other 5 for the selected hits) and employed them for MDS studies. All the systems were placed in a cubic box and water molecules were filled using the gmx solvate tool. Ions were added for the neutralization of the systems. The neutralized systems were employed for energy minimization using the steepest descent algorithm for 5000 steps for removing the steric clashes and irregular geometry. After that Constant Number of particles, Volume, and Temperature (NVT) and Constant Number of particles, Pressure, and Temperature (NPT) simulation of 1 ns each was performed for maintaining the temperature, pressure, and volume of the systems. The van der Waals and short-range electrostatic interactions applied a cutoff of 1.0 nm. The Particle Mesh Ewald method was used for treating the large range electrostatic interaction 36 . The LINCS algorithm was used for constraining all the bonds 37 . Finally, all the systems were employed for 150 ns MDS studies and the coordinates were saved for every 0.2 fs. Analyses of RMSD, RMSF, Rg, SASA, intramolecular interaction, PCA, and MM-PBSA were performed with gmx rms, gmx rmsf, gmx gyration, and gmx sasa, etc. tools respectively as earlier described 38,39 . The trajectory was visually analyzed using Chimera 1.13.2 29 .
Principal component analysis. The principal component analysis (PCA) is a widely used method to identify the correlated motions which are induced by ligand binding 40 . The conformational ensemble of MD simulation was used to conduct PCA analysis. A covariance matrix was constructed on the basis of Cα atom displacement 41,42 . The inbuilt tool of Gromacs gmx covar and gmx anaeig was used for the eigenvector and eigenvalue generation. The covariance matrix (C ij ) is defined from the following equation: where, the r i and r j are the mass-weighted Cartesian coordinates of the i th and j th Cα atoms.
(1) www.nature.com/scientificreports/ The eigenvalues and their respective eigenvectors were obtained by the diagonalization of the covariance matrix. Finally, the plots were drawn by the Origin 6.0 software. We also took the first two Principal Components (PCs) and plotted them against each other to generate the 2D projection plot. This plot was used for the threedimension free energy landscape analysis as described earlier 43 . Binding free energy calculations. The binding free energy for all the complexes was calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) tool 44 . The final 500 snapshots generated from each trajectory was used for the analysis. The binding free energy was calculated using this equation: The binding affinity ΔG mm (molecular mechanics energy) was estimated considering the electrostatic and van der waals interactions. The polar and non-polar contributions estimate the ΔG sol (solvation free energy). The solvent accessible surface area (SASA) describes the nonpolar solvation free energy. The entropy (TΔS) was excluded due to its high computational cost.

Results
Virtual screening. Structure-based virtual screening was conducted for the identification of novel inhibitors against NS5 (PDB ID: 5MRK). We utilized Autodock Vina, Smina, and I-dock for this purpose. From all these tools we observed that ZINC03830383 (Carbinoxamine) was the highest energy compound based on analyses of all three docking software. A binding affinity of −11.2 kcal/mol, −11.42 kcal/mol, and −12.9 kcal/ mol were obtained from Autodock Vina, I-dock, and Smina for ZINC03830383. ZINC08034121 showed that the least binding binding affinity from all the docking software, which was −2.5, −2.39, and −2.6 kcal/mol from Autodock Vina, I-dock, and Smina. The binding affinity of 2895 compounds is shown in Supplementary Table 1. From this screening we selected compounds based on their binding energy as well as their interaction analysis. Further selected compounds were used for the cross-docking analysis using Autodock Tools.

Docking result analysis.
Based on the virtual screening results, we selected the 28 best compounds for the cross-docking studies using Autodock. ZINC3914596 (Saquinavir) showed the highest binding affinity of −9.77 kcal/mol, while ZINC3938482 (Posaconazole) showed the best binding energy of −6.03 kcal/mol. All 28 compounds showed a favorable binding affinity toward NS5, indicating that they can act as inhibitors. From this result we selected five compounds that showed good binding affinity from the three docking software which do not directly target human receptors except Ol_Me; however, they showed binding to the bacterial or viral receptors [45][46][47] . The details of binding energy, hydrogen bond interactions, and other interacting residues are shown in Supplementary Table 2. The ZINC ID, drug name, and binding affinity of selected 5 compounds with the control compound are shown in Table 1.
NS5-SFG interaction. SFG was used as the positive control and demonstrated binding affinities of −7.2, −7.7, −7.76, and −8.1 kcal/mol for Autodock, Autodock Vina, idock, and Smina analyses respectively. The NS5-SFG complex was stabilized by various forces including hydrophobic, hydrogen bond, and electrostatic interactions. SFG formed hydrogen bonds with Ser56, Gly85, Gly86, Trp87, and Glu111 residues in NS5. The Interaction details of SFG are shown in Fig. 3A. The residue details are given in Supplementary Table 2. NS5-Ceforanide interaction. Ceforanide showed binding affinities of −9.54, −9.2, −9.09, and −9.2 kcal/mol from Autodock, Autodock Vina, I-dock, and Smina respectively. The interaction details are shown in Fig. 3B from Autodock. The Autodock complex was stabilized by various interactions including hydrogen bond interactions, hydrophobic interactions, and ionic interactions. The Ceforanide complex was mainly stabilized by hydrogen bonds with Ser56, Lys61, Gly85, Trp87, and Lys182 of NS5. The interaction details are given in Supplementary Table 2.
NS5-Saquinavir interactions. The binding affinities for Saquinavir toward the NS5 receptor were calculated to be −9.77, −9.5, −9.31, and −9.6 kcal/mol based on Autodock, Autodock Vina, I-dock, and Smina respectively. Figure 3C illustrates the brief interaction of ligand with the NS5 protein. The complex was stabilized by various interactions including hydrogen bond interactions, hydrophobic interactions, and ionic interactions. Saquinavir forms hydrogen bonds with the NS5 residues including Lys105, Gly106, Glu111, and Glu149. The interaction details are given in Supplementary Table 2. NS5-Amcinonide interaction. Amcinonide showed binding affinities of −9.52, −10.0, −10.7, and −10.07 kcal/ mol from Autodock, Autodock Vina, I-dock, and Smina respectively. We analyzed the binding of the drug using Autodock (Fig. 3D). The complex was stabilized by hydrogen bonds and hydrophobic interactions. Amcinonide formed hydrogen bonds with NS5 residues including Ser56, Gly86, Gly148, and Glu149 residues of NS5 (Supplementary Table 2).  Table 2. www.nature.com/scientificreports/ predicted complexes for 150 ns with SFG as the control compound along with apo-NS5. Finally, seven systems were created and simulated to achieve the protein-ligand dynamics.     www.nature.com/scientificreports/ Interaction analysis. NS5-ligand interactions are key to stabilization of the complex. Hydrogen bonds are very specific and play a pivotal role in the stability of the protein-ligand complex. We calculated the number of hydrogen bonds for the MD trajectory of NS5-ligand complexes (Fig. 5C). The average number of hydrogen bonds for NS5-SFG, NS5-Ceforanide, NS5-Squanavir, NS5-Amcinonide, NS5-Cefpiramide, and NS5-Ol_Me was 7, 8, 6, 7, 12, and 10 respectively. The results illustrate that most of the complexes showed a higher number of hydrogen bonds compared to the number of hydrogen bonds of the control ligand (SFG) NS5 complex, with the single exception of the NS5-Squanavir complex. The NS5-Cefpiramide and NS5-Ol_Me showed multiple numbers of hydrogen bonds, 10 and 12 respectively, compared to the other selected ligands. The overall result suggests that these two complexes are more stable than other predicted hits.
Principal component analysis. We next predicted the correlated motions of NS5-drug complexes using PCA. All the NS5 bound complexes along with the apo-NS5 were analyzed using gmx covar and gmx anaeig tools. The eigenvalues were obtained by the diagonalization of the covariance matrix of atomic fluctuations. They are plotted in Fig. 6A. The first few eigenvectors play a key role in the overall dynamics of the system. We selected the first 50 eigenvectors which showed the motions of 79.33%, 80.82%, 78.14, 81.40%, 72.44%, 76.53%, and 80.29% for apo-NS5, NS5-SFG, NS5-Ceforanide, NS5-Squanavir, NS5-Amcinonide, NS5-Cefpiramide, and NS5-Ol_Me respectively. Figure 6A clearly shows that apo-NS5 has lower motions while the highest motion was observed in the NS5-Squanavir complex. Compared to the control compound, the predicted ligands showed lesser motion. It is evident that NS5-Ceforanide, NS5-Amcinonide, and NS5-Cefpiramide induced less motion and showed higher stability in the complex. We selected the first two principal components and plotted them against each other (Fig. 6B). Apo-NS5 showed a very dense and stable cluster as compared to the control compound and other ligands. The control complex, NS5-SFG, showed a very large and unstable cluster while the predicted drugs also showed dense and stable clusters. The cluster of the predicted hits, NS5-Amcinonide, NS5-Cefpiramide, and NS5-Ol_Me, were more stable as compared to the other complexes. www.nature.com/scientificreports/ After analysis we selected the first eigenvector to predict the correlated motions on the basis of residues. All the values were calculated based on eigenvalues vs. residues (Fig. 6C) Gibbs free energy landscape analysis. The Gibbs free energy landscape analysis was performed using the first two principal components to get the free energy pattern of the apo-NS5, NS5-SFG, and predicted hits (Fig. 7). The apo-NS5 shows two energy funnels separated by energy barriers. The energy funnels are separated by various small energy minima and occupy the most space; the large area with intense blue color represents the stable cluster. The NS5-SFG also shows an apo-NS5-like pattern but here one energy funnel clearly represents a deep well and the other is bifurcated into four small energy minima, illustrating that the NS5-SFG complex changes through many conformational states to reach a stable state. The NS5-Ceforanide clearly displays two deep wells, with one well occupying area, as seen with the deep blue color. It indicates that this complex has two thermodynamically conformational states. The NS5-Saquinavir shows four clear deep wells which are separated by a high-energy barrier, symbolizing the four conformational states for this complex. The NS5-Amcinonide shows a single energy funnel with a deep blue color. It represents that this complex has only one conformation, indicating that it is the most stable complex as compared to the others. The NS5-Cefpiramide shows five energy funnels and they are very close to each other, forming a stable cluster. They are not separated by any high-energy barrier, again representing that this complex is also stable. The NS5-Ol_Me displays three deep wells or energy funnels, www.nature.com/scientificreports/ two of which show small conformational states with light blue color that are not thermodynamically favorable, while larger area with deep blue color illustrates the stable conformation of this complex. From the overall FEL results, we have concluded that the predicted hits undergo several transition states to reach the stable states and show the NS5-SFG like pattern.
To evaluate the hotspot residues important for ligand stabilization, we calculated the per residue energy contribution (Fig. 8). To reduce the complexity of data for representation, we selected a few key binding site residues for the analysis. We observed that Lys105, Asp131, Asp143, and Asp146 residues play a crucial role in drug binding.

Discussion
At present, various promising drug targets for ZIKV have been identified which can inhibit the viral growth. The anti-ZIKV compounds can reduce or halt the progression of ZIKV infection through the direct inhibition of key enzymes responsible for maintaining the biological function of ZIKV. Currently, many anti-ZIKV compounds are available which can inhibit various enzymes such as Temoporfin (NS2B/NS3 protease inhibitor) 48 , Suramin (NS3 helicase inhibitor) 49 , Nanchangmycin (envelope glycoprotein inhibitor) 50 , Sofosbuvir (NS5 RdRp inhibitor) 51 , and SFG (NS5 methyltransferase inhibitor) 13 . Although several approaches to developing new therapies have been tested over many years in the anti-ZIKV drug discovery program, no medications are available in the market against ZIKV. The compounds mentioned above are either cytotoxic or require high doses to inhibit the targets, which is a main reason that they have not been proposed as drugs in the market 52 . However, several targets are available for anti-ZIKV inhibitor identification with ZIKV NS5 as one of the major and most promising targets to date. The ZIKV NS5 crystal structure showed that the SAM binding site domain, called the SAM-dependent domain, could be a promising drug target 53 . Various structures complexed with the inhibitors targeting the SAM-dependent domain are available in the protein data bank. An analogue compound of SAM called SFG has been shown to inhibit the Dengue virus (DENV-2) and West Nile virus (WNV) with values of 0.63 mM and 14 mM IC 50 , respectively 54 . Various structures are available complexed with the SFG. Another SAM analogue called NSC 12155, identified in a virtual screening study, showed in vitro inhibitory activity against the WNV and Yellow fever virus (YFV) NS5 with an IC 50 ranging from 0.50 to 3.00 mM. NSC 12155 showed an EC 50 ranging from 1.00 to 7.00 mM against WNV, DENV-2 and Japanese encephalitis (JEV) in a cell-based assay 55 ; However, NSC 12155 did not exhibit any activity against ZIKV. Theaflavin, a main component of tea, was also tested against NS5 of ZIKV, and showed 10.10 µM IC 50 and 8.19 µM EC 50 values against NS5. The authors also employed a computational study and predicted that theaflavin would show strong binding with D146 residue, which is a key catalytic residue 56 .
Drug repurposing is an effective approach to find potential drugs against many diseases, including ZIKV. It involves identifying drugs previously approved by the FDA for various diseases for use against other diseases for which they were not intended. For example, chloroquine, a widely used antimalarial drug, can significantly inhibit in vitro ZIKV infection with a value of 1-5 µM EC 50 57,58 . Chloroquine alters the low pH-induced conformational changes which are necessary for fusion of the envelope protein with the endosomal membrane 59 ; hence early stage ZIKV replication is blocked by chloroquine 60 . It is generally known that RNA viruses can mutate quickly; therefore, the identification of potential compounds with different mechanisms is required. Due to the importance of ZIKV NS5, we selected this enzyme as a drug target in this study and used computational approaches to find a novel drug through the drug repurposing approach. Computational approaches are practical 61 as they avoid the time-consuming and expensive hurdles posed by ADMET in traditional drug designing 62 . Therefore, in this study, we retrieved the FDA-approved compound and screened it against NS5. We employed SFG as a control compound because it is co-crystallized with the crystal structure and has also shown good inhibitory activity against NS5 in several other studies. Virtual screening revealed various high-energy compounds; however, we selected only the five best antibacterial and antiviral compounds against NS5. The binding free energy and various structural parameters were calculated using the MDS trajectory, which revealed that Cefpiramide and Ol_Me can act as anti-ZIKV compounds with good binding affinity as compared to SFG. The drugs showed binding with key catalytic residue such as D146 which are required to perform the native enzyme activity 27,56 . Our study further revealed that these predicted compounds also showed binding with key catalytic residues, which proves that they can alter enzyme activity and reduce the viral burden. The selected drugs belong to antiviral www.nature.com/scientificreports/ and antibacterial families; therefore, they can be easily repurposed against ZIKV because they cannot inhibit the human targets and cause toxicity. Cefpiramide is a well-known anti-bacterial compound active against Pseudomonas aeruginosa. It belongs to the semi-synthetic, broad-spectrum, beta-lactam, third-generation cephalosporin antibiotic family 63 . It inhibits bacterial cell wall synthesis by inactivating penicillin binding proteins (PBPs) through alteration of the final transpeptidation step, which is necessary for peptidoglycan cross-linking, a major cell wall component. 47,64 . Cefpiramide also showed antiviral activity against SARS-CoV-2. An in-silico study showed that Cefpiramide can inhibit the spike protein of SARS-CoV-2 with good binding affinity 65 . In another study, Cefpiramide was tested against the ACE2 receptor of SARS-CoV-2 and was observed to inhibit the receptor with −9.1 kcal/mol binding affinity 66 . In the current study, Cefpiramide showed both good binding affinity and stability in the MDS analysis, which indicates that it can act as a potential anti-ZIKV compound.
Ol_Me, also known as Benicar, is used to treat the high blood pressure (hypertension). It is an Angiotensin II Receptor Blocker 67 . It is a prodrug which is hydrolysed into Olmesartan during absorption from the gastrointestinal tract 68 . Although it is not known for its antiviral or antibacterial activity, in our study we found Ol_Me to be a potential NS5 inhibitor in the docking, binding free-energy and other MDS analyses. However, this drug has been shown to directly interact with human targets so minor changes would be required in this drug to propose it as an anti-ZIKV compound. In the future, researchers can test the efficacy of this drug against ZIKV NS5.
The overall analysis revealed that drug repurposing is a very powerful approach that can predict which FDAapproved drugs can be directly repurposed for diseases other than those for which they were originally intended.

Conclusions
ZIKV is an arthropod-borne virus (arbovirus) belonging to the family of Flaviviridae and genus Flavivirus. As a single-stranded positive RNA virus, the genome of ZIKV is approximately 10 kb that encodes for three structural proteins and seven non-structural proteins. Inhibiting the proteins that are crucial for viral activity can reduce the growth of ZIKV. In this study, we used the available crystal structure of the NS5 protein and screened FDAapproved compounds to search for a potential inhibitor. A total of 2866 FDA-approved drugs were screened using virtual screening methods and the top five drugs were selected for evaluation of their ability to bind to NS5 using 150 ns MDS studies. Based on various MDS parameters such as RMSD, RMSF, Rg, SASA, PCA, and binding free energy, we found that, out of the five selected hits, Cefpiramide and Ol_Me formed stable interactions with NS5. Since the selected molecules are FDA-approved drugs, they may have an advantage in terms of their previous pharmacodynamics and pharmacokinetics properties. After evaluation, we propose that these may act as lead compounds for the development of potential inhibitors of NS5.

Data availability
All data generated or analyzed during this study are included in this published article (and its Supplementary Information files).