In Silico Studies of Potential Selective Inhibitors of Thymidylate Kinase from Variola virus

Continuing the work developed by our research group, in the present manuscript, we performed a theoretical study of 10 new structures derived from the antivirals cidofovir and ribavirin, as inhibitor prototypes for the enzyme thymidylate kinase from Variola virus (VarTMPK). The proposed structures were subjected to docking calculations, molecular dynamics simulations, and free energy calculations, using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method, inside the active sites of VarTMPK and human TMPK (HssTMPK). The docking and molecular dynamic studies pointed to structures 2, 3, 4, 6, and 9 as more selective towards VarTMPK. In addition, the free energy data calculated through the MM-PBSA method, corroborated these results. This suggests that these compounds are potential selective inhibitors of VarTMPK and, thus, can be considered as template molecules to be synthesized and experimentally evaluated against smallpox.


Introduction
Despite the declaration of the World Health Organization (WHO) that the Variola virus was eradicated from the world by the 1980s, smallpox is still a matter of concern. Moreover, many studies on the development of drugs against this disease have been reported in the literature [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. This is due to the fact that Variola virus strains may have been stored in clandestine sites around the world, and their potential use for bioterrorist purposes cannot be ignored [17,18]. Furthermore, the recreation of this kind of virus has proven to be easy to carry out, considering the current technologies [19]. As viruses can survive long periods in nature, there is no guarantee that smallpox will never return as a natural pandemic [20,21]. It is also important to mention that since the 1980s the public vaccination campaigns against smallpox do not exist anymore [4,7,9], a fact that makes all of the world's population under 40 years of age particularly vulnerable.
Currently, there is only one drug approved by the Food and Drug Administration (FDA) of the United States for the treatment of smallpox. This drug, named tecovirimat ( Figure 1), acts as an inhibitor of the protein complex needed for the survival of extracellular viruses, and thereby prevents its spread throughout the body from the infected cells [1,2,10]. However, despite the fact that it was approved by the FDA, tecovirimat is still in phase II of clinical trials and has not been tested in humans yet, due to the lack of patients with smallpox. It has been effective in laboratory tests to protect animals against monkeypox and rabbitpox as well as presented low toxicity to humans [1,2,10,22,23]. The risk of smallpox resurgence either, as a natural pandemic or consequence of a terrorist attack, combined with the huge number of non-immunized people in the world, and the scarcity of drugs to combat this disease, highlight the importance of the search for new drugs against it.
It is known from the literature that the Variola virus belongs to the genus Orthopoxvirus, which, similar to all poxviruses, are capable of encoding their own thymidine and thymidylate kinases (TK and TMPK) [24]. Caillat et al. [24] have detailed the sequence alignment of TMPKs from different Orthopoxviruses (vaccinia virus, smallpox, cowpox, camelpox, and monkeypox) and showed that those enzymes present highly similar amino acid sequencing, which are practically identical, except for a few residues. The majority of the different residues among the viral enzymes belong to the protein loops [24]. Therefore, the crystallographic structure of the Vaccinia virus TMPK (VaccTMPK) can be used as a template for modeling the structure of the Variola virus TMPK (VarTMPK), which is not available yet in the protein data bank (PDB) (https://www.rcsb.org/). Moreover, the study by Caillat et al. [24] included data on practical experiments using compounds such as brivudine monophosphate to test their inhibition of VaccTMPK and human TMPK (HssTMPK). HssTMPK and VaccTMPK are only 42% similar in their amino acid sequence alignments and the dimer interface arrangements of both enzymes are different, as well as their active site geometry [24][25][26]. Due to these differences, as stated by the authors, VaccTMPK is able to accommodate more voluminous compounds, such as brivudine monophosphate, which stabilizes the enzyme and is phosphorylated. The specificity of the bond with brivudine monophosphate shows that selective antipox agents can be developed based on this finding. Therefore, these antivirals can also be investigated for their use in the treatment of TMPK-related diseases from other viruses of the Orthopoxvirus genus [24].
According to the literature, TMPK is responsible for the synthesis of thymidine 5 -triphosphate (TTP) based on the phosphorylation of thymidine 5 -monophosphate (TMP) [27,28]. TTP participates in DNA synthesis (as a building block) and its levels are controlled during the different phases of the cell cycle. Since TMPK is important for the biosynthesis of TTP, the enzyme is considered a molecular target for the development of antiviral drugs against many diseases, and the interruption of the metabolism of TTP can be used to stop the development of said illnesses. In addition to its direct effect on the metabolism of TTP, TMPK also interacts directly with the DNA synthesis process, since it plays an important role in the activation of DNA precursor nucleoside analogues. As a result, its inhibition directly affects the synthesis of genetic material, thus provoking the deactivation (death) of the virus [25,29,30].
Previously, we proposed VarTMPK as a potential target to the design of potential selective inhibitors and studied the behavior of known antivirals inside it, thus pointing to the relevant residues to be targeted in the drug design [7]. Next, we proposed the structures of 10 potential VarTMPK inhibitors based on theoretical studies of another series of antivirals [4,8,24,26]. The best ranked compounds after docking and MD simulations studies were, then, selected to design a new series of compounds based on alterations of their structural features, having in mind the synthetic viability and selectivity towards VarTMPK [9]. A few compounds of this new series presented promising theoretical results [9].
In the present work, we performed the theoretical study of 10 new compounds whose design was based on structural modifications of the antivirals cidofovir and ribavirin (Figure 1), seeking the application of these new compounds against the Variola virus. Cidofovir and ribavirin were chosen to serve as templates for this study, since they have already been extensively studied and validated for diseases caused by Orthopoxviruses. Therefore, they have much data available in specialized databanks [31,32].
Maintaining the main idea of our previous studies [4,7,9], the objective of this work was to minimize the structural complexity of the template compounds, in order to make the synthesis of these structures more feasible and less cumbersome.
To analyze the binding modes and selectivity of the proposed compounds inside VarTMPK and HssTMPK, the molecular docking method was used [33,34]. In order to assess their dynamical behavior and corroborate the docking achievements, molecular dynamics (MD) simulations rounds were performed on the best poses obtained from docking. Finally, free energy calculations, applying the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method, were performed on the most promising compounds to verify the effectiveness, reliability, and selectivity of their binding inside each enzyme [35][36][37][38][39][40][41][42][43][44][45]. Table 1 lists the active site residues of VarTMPK and HssTMPK. It is possible to see that seven of these 14 residues are different. This 50% difference highlights the possibility of designing selective inhibitors for VarTMPK. Similar to our previous works [4,9], the ionization states of the compounds used in this study corresponded to the predominant microspecies at physiological conditions (pH = 7.4). According to the chemicalize server (www.chemicalize.com), our compounds except compound 1, are 100% in the neutral form at pH 7.4. The dominant microspecies for compound 1 predicted by chemicalize (www.chemicalize.com) is the one negatively charged in the phosphate group, with 74% of prevalence at pH 7.4. Moreover, according to the chemicalize server (www.chemicalize.com), all of our compounds meet the Lipinski's criterion of drug likeness [46], as shown in Table 2. The docking studies were meant to identify the compounds with better interactions and higher selectivity towards VarTMPK. The best poses for each compound were selected in accordance with the best (more negative) interaction energies inside VarTMPK and HssTMPK (intermolecular and hydrogen bond) and the best superposition onto TDP. As shown in Figure 2, the poses selected in the docking studies present an optimal superposition onto TDP for both VarTMPK and HssTMPK. These poses were analyzed according to: Their interaction energy values with the enzymes (E interaction ) and the cofactor (E cofactor ); the H-bond energies between the compounds and the enzymes (H bond ); and, finally, the amino acids involved in the interactions with the compounds in both enzymes. The results obtained for those parameters are summarized in Table 3. Figure 3a-j shows the interactions observed between the proposed compounds inside VarTMPK and HssTMPK. The active site residues of both enzymes are shown in red and the red spheres represent water molecules.

Molecular Dynamics Simulations
The best poses of compounds 2, 3, 4, 5, 6, 8, and 9 were selected for additional rounds of MD simulations, since they were pointed as more selective towards VarTMPK by the docking studies. The results are shown in the plots of root-mean-square deviation (RMSD) and number of H-bonds formed during the MD simulations, for their complexes inside VarTMPK and HssTMPK, shown in Figures 4 and 5, respectively.
In order to validate our MD protocol, we performed three rounds of MD simulations for the complexes of compound 2 with both enzymes. The obtained plots of RMSD are shown in Figures S1-S6 of the Supplementary Material.

MM-PBSA Calculations
According to the literature, one of the limitations of docking studies lies in determining the affinities of ligand-protein complexes using calculations based solely on the poses generated for these complexes. In this regard, the MM-PBSA method tends to improve the results related to the binding energy, which is due to the fact that MM-PBSA provides more accurate results when the affinities between the ligand and protein are determined. Moreover, this is achieved by calculating the free binding energy associated with the formation of ligand-protein complexes [47] and decomposing it in the contributing components, as proposed before by Bren et al. [44,45]. Table 4 shows the average binding energy values calculated for compounds 2, 3, 4, 6, and 9 complexed with VarTMPK and HssTMPK, while Figure 6 illustrates the favorable and unfavorable energy contributions of compounds 3 and 6 inside VarTMPK.

Discussion
As can be seen, Table 2 shows the toxicity and mouse carcinogenicity values of the investigated structures and the antivirals used as precursors. In short, with respect to the proposed compounds, it is important to point out that the acute toxicity values in algae calculated for all the structures were lower than the value found for cidofovir (1.19), which is not toxic. Furthermore, except for compounds 1 and 10, all the structures had lower toxicity values in comparison with those of ribavirin (0.54).
According to the literature, the cavity of HssTMPK acquired through the software Molegro virtual docker (MVD) ® , presents a volume of 90.112 Å 3 , which is greater than the value of 76.288 Å 3 observed for the cavity of VarTMPK. It is important to point out that the HssTMPK cavity is more outspread and narrower. Therefore, it is expected to have issues in the entrance and/or permanence of larger inhibitors in the binding site [7].
The results in Table 3 show that all the compounds are capable of binding to the active sites of both enzymes. This is reflected by the negative values of E interaction observed for all of them. Comparing the docking results obtained, we can state that compounds 2, 3, 4, 5, 6, 8, and 9 form more stable complexes with VarTMPK. This enables us to infer that they might show higher affinity for the active site of this enzyme, since the systems formed among these compounds and the viral enzyme presented lower values of E interaction , E cofactor , and H-bond energy. Regarding the ∆E interaction values, determined based on the difference between E interaction inside VarTMPK (lower value) and HssTMPK, the compounds listed above ranked: 2 (−62.85 kcal.mol −1 ) < 5 (−55.40 kcal.mol −1 ) < 9 (−45.90 kcal.mol −1 ) < 4 (−38.57 kcal.mol −1 ) < 3 (−35.41 kcal.mol −1 ) < 6 (−27.85 kcal.mol −1 ) < 8 (−11.32 kcal.mol −1 ).
The best pose of compound 1 interacts with three active site residues of VarTMPK and two active site residues of HssTMPK, as shown in Figure 3a. Among these residues, Arg93 from VarTMPK and Arg97 from HssTMPK are equivalent between both enzymes.
Regarding compound 2 (Figure 3b), we observed interactions with three residues from VarTMPK, of which two belong to the active site, while for HssTMPK, an interaction with only one residue was observed, which also belongs to the active site.
Compounds 3 ( Figure 3c) and 4 ( Figure 3d) interact with the same number of residues inside both enzymes. However, for VarTMPK, two interacting residues belong to the active site versus only one for HssTMPK. Residues Arg93 from VarTMPK and Arg97 from HssTMPK are the only equivalent interacting residues observed. Figure 3e shows that compound 5 interacts with two active site residues from VarTMPK and one from HssTMPK. No matching interacting residues were observed.
According to Figure 3f, compound 6 interacts solely with one active site residue from VarTMPK and has no interaction with the HssTMPK active site residues. On the other hand, for compound 7, Figure 3g shows interactions with two active site residues from VarTMPK and one from HssTMPK. Arg93 from VarTMPK and Arg97 from HssTMPK were the only equivalent interacting residues observed.
Regarding the complexes formed with compound 8 (shown in Figure 3h), we can see interactions with two active site residues from VarTMPK and one from HssTMPK. In this case, Asp13 from VarTMPK and Asp15 from HssTMPK are the equivalent interacting residues observed. Figure 3i shows that compound 9 interacts with three active site residues from VarTMPK, and two from HssTMPK. The equivalent interacting residues observed are: Asp13 (VarTMPK) and Asp15 (HssTMPK), and Arg72 (VarTMPK) and Arg76 (HssTMPK).
Finally, Figure 3j shows that compound 10 interacts with four active site residues from VarTMPK, and only one from HssTMPK. The equivalent residues observed are Arg93 from VarTMPK and Arg97 from HssTMPK.
The RMSD plots for the multiple MD simulations, performed for the systems enzyme/compound 2 ( Figures S1-S6), showed the same behavior for all dynamics, with no variation above 0.1 nm. This validates our MD simulations protocol for the systems under study.
Comparing the RMSD plots in Figure 4 for both enzymes, it is possible to see that for all the systems, the RMSD values never passed 0.4 nm for the enzyme (black lines) and 0.2 nm for the ligands (red lines) during the simulated time. However, for most of the systems, both the enzyme and compounds presented more fluctuations inside HssTMPK. Inside VarTMPK, this fluctuation never passed 0. 05 nm for compounds 2, 3, 4, 6, and 9.
Compounds 5 and 8 were the only ones showing a similar or more instable behavior inside VarTMPK compared to HssTMPK (see Figure 4g,h,k,l). This does not corroborate with the docking results and suggests that these compounds are not selective towards VarTMPK.
In addition, we analyzed the number of H-bonds formed during the MD simulations for the complexes inside VarTMPK and HssTMPK ( Figure 5). The H-bonds profiles of the complexes with compounds 2 (Figure 5a,b), 3 (Figure 5c,d), 4 (Figure 5e,f), 6 (Figure 5i,j), and 9 (Figure 5m,n) show that these compounds were capable of keeping at least three H-bonds with VarTMPK during the whole simulated time, presenting a profile better or similar to HssTMPK. Conversely, the H-bond graphs of the complexes with compounds 5 (Figure 5g,h) and 8 (Figure 5k,l) show a more unstable behavior with the prevalence of no more than two H-bonds during the whole simulation, presenting a worse (compound 5) or similar (compound 8) profile compared to HssTMPK.
The results regarding the formation of H-bonds can be correlated to the RMSD observed during the MD simulations, and point to compounds 2, 3, 4, 6, and 9 as more selective towards VarTMPK. For this reason, these compounds were selected for the determination of their binding energies inside VarTMPK and HssTMPK based on the MM-PBSA calculations.
Regarding the MM-PBSA calculations, the results in Table 4 show that for all compounds, the mean binding energy values of the complexes formed with VarTMPK were lower than the values determined for the complexes formed with HssTMPK. These data confirm the results obtained by the docking and MD simulations, which suggested the selectivity of these compounds towards VarTMPK.
Moreover, the results in Table 4 show that compounds 3 and 6 stand out among the others in terms of binding energy with VarTMPK, which are the most stable ones. Furthermore, the difference related to HssTMPK was higher than 90 KJ/mol in both cases, suggesting a high selectivity towards VarTMPK.
A possible explanation for the high free energy values of the complexes between VarTMPK and compounds 3 and 6, is the presence of the trifluoroacetate group (-COOCF 3 ) attached to the tetrahydrofuran ring. The presence of three fluorine atoms-a chemical element with high electronegativity and with isolated and free electron pairs on its structure-likely favors the selective interactions inside VarTMPK. Therefore, these interactions can help in maintaining the complexes formed between VarTMPK and compounds 3 and 6-thus contributing to the higher stability of these systems.

Compounds Studied
The 2D structures of the compounds studied in this work are shown in Figure 7. As mentioned above, they are all derivatives from the antivirals cidofovir and ribavirin, designed to better explore the potential selective interactions with residues of the active site of VarTMPK.
The ionization states of each compound in Figure 7 at physiological pH (7.4)-as well as the drug likeness criteria established according to the Lipinski's rule of five [48]-were calculated with the aid of the Chemicalize databank (https://chemicalize.com) [46]. Their toxicity and carcinogenicity in mice were calculated with the online program PreADMET (https://preadmet.bmdrc.kr/adme-prediction/).
The tridimensional structures of these compounds, and the calculation of their atomic charges, were done through the Recife model 1 (RM1) semiempirical model [49,50], using the Spartan 08 ® Suite software.

Applied Protocols
In previous works [4,7,9], the tridimensional structures of VarTMPK and HssTMPK that were used for our studies include, respectively, the homology model complexed with TDP and the Mg 2+ cofactor, constructed by Guimarães et al. [7], and the crystallographic structure of HssTMPK complexed with TDP and Mg 2+ , available in the protein data bank (PDB) under the code 1E2G [51]. According to Guimarães et al. [7], the Ramachandran plot of the homology model showed that 99.5% of the residues were situated in the most favored regions.
The docking protocol used, as well as its validation by re-docking studies, was the same as employed by Guimarães et al. [7]. The MolDock algorithm [52] from the Molegro virtual docker (MVD ® ) was used to perform the energy calculations of the ligands inside VarTMPK and HssTMPK. The binding site was restricted to spheres with radii of 6 and 10 Å, respectively, around TDP, and all the residues inside these spheres were set to be flexible. The coordinates were centered on x = 8.95, y = 22.41, and z = 0.69 for VarTMPK, and x = 13.92, y = 75.19, and z = 25.05 for HssTMPK [7]. The best poses for each ligand in the viral and human enzymes were selected for further MD simulations.
The protocol comprising the application of MD simulations rounds to the best poses selected from the docking studies was based on our most recent work [9]. Each ligand was parameterized in order to be recognized by the OPLS-AA force field [53]. The enzyme/ligand complexes were put inside cubic boxes of 450 nm 3 containing around 13,000 Tip4P type water molecules under periodic boundary conditions (PBCs), using the GROMACS 5.1.4 software [54]. The energy minimization steps used were: (1) Steepest descent with position restraint (PR) of the ligand; (2) steepest descent without PR; (3) conjugate gradients; and, lastly, (4) quasi Newton Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [55], with a minimal energy of 1 kcal.mol −1 . Thereafter, the systems were submitted to temperature (NVT) and pressure (NPT) balancing phases in order to attain equilibrium. The equilibrated systems were, then, submitted to 500 ps of MD at 310 K with a PR for the entire system, except for water molecules, in order to ensure the equilibrium of solvent molecules, and finally, 100,000 ps of free MD simulations at 310 K without PR, with 2 fs of integration time, and a cutting radius of 10 Å for long-distance interactions [9]. All the Glu and Asp residues were assigned with negative charges, and the Lys and Arg residues were assigned with positive charges. As mentioned above, the visual molecular dynamics (VMD) software [56] was used to analyze the MD results of the systems. In order to validate our MD simulations protocol, three MD simulations were performed for the systems VarTMPK/compound 2 and HssTMPK/compound 2.
The MM-PBSA method was employed to predict the free binding energies of the ligands inside the VarTMPK and HssTMPK. This study allows us to infer whether the binding process is spontaneous (freeing energy) or nonspontaneous (requiring energy), thus enabling us to point to the prototype that is more interesting and promising regarding the selective inhibition of VarTMPK. The determination of the free energy of formation of the complexes, in association with the MD simulations, take into consideration three energetic terms in the calculation of binding energy: (1) Changes in the potential energy of the system in a vacuum; (2) polar and apolar solvation of the different species; and (3) the entropy related to the formation of the complexes during the gaseous phase [9,36,37,44,45]. For all the enzyme/ligand complexes, MM-PBSA calculations were performed using the g_mmpbsa tool [37] from the GROMACS package. In order to consider non-correlated frames, the structures for the free energy calculations were obtained at 500 ps each after the stabilization of the systems [9].
The parameters used for the docking, MD simulations, and MM-PBSA calculations performed in this work are summarized in Table 5.

Conclusions
According to the docking studies, it was found that only compounds 1, 7, and 10 did not demonstrate selectivity towards VarTMPK. The MD simulations studies corroborated the docking results for compounds 2, 3, 4, 6, and 9, implying that the complexes formed between each of them and VarTMPK are more stable when compared to the same complexes with HssTMPK. This was also observed in the further MM-PBSA calculations of the binding energies that showed lower binding values for these compounds inside VarTMPK than inside HssTMPK. Moreover, the MM-PBSA studies pointed to compounds 3 and 6 as the most promising selective inhibitors of VarTMPK with the higher difference in binding energies. This is probably due to the presence of the-COOCF 3 group on these compounds that could be favoring their stabilization inside VarTMPK.
In conclusion, our results suggest that five of the 10 proposed compounds, derived from the antivirals cidofovir and ribavirin, can be excellent alternatives for the development of new drugs against smallpox. Since TMPK is highly conserved amongst other Orthopoxviruses, similar to the vaccinia and monkeypox viruses, the same theoretical approach used here can be done with the TMPKs of those viruses as a proof of concept. Furthermore, the experimental evaluation of these compounds should be first performed against those viruses as a primarily model to check their potential in vitro. This would help in confirming the activity of those compounds before moving to the more complicated evaluation against the Variola virus, which will demand both an authorized laboratory and a stricter safety protocol.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.3 390/ph14101027/s1, Figure S1: RMSD of the first MD simulation for the system HssTMPK/Compound 2, Figure S2: RMSD of the second MD simulation for the system HssTMPK/Compound 2, Figure  S3: RMSD of the third MD simulation for the system HssTMPK/Compound 2, Figure S4: RMSD of the first MD simulation for the system VarTMPK/Compound 2, Figure S5: RMSD of the second MD simulation for the system VarTMPK/Compound 2, Figure S6: RMSD of the third MD simulation for the system VarTMPK/Compound 2.