Using docking and alchemical free energy approach to determine the binding mechanism of eEF2K inhibitors and prioritizing the compound synthesis

A-484954 is a known eEF2K inhibitor with submicromolar IC₅₀ potency. However, the binding mechanism and the crystal structure of the kinase remains unknown. Here, we employ a homology eEF2K model, docking and alchemical free energy simulations to probe the binding mechanism of eEF2K, and in turn, guide the optimization of potential lead compounds. The inhibitor was docked into the ATP-binding site of a homology model first. Three different binding poses, hypothesis 1, 2, and 3, were obtained and subsequently applied to molecular dynamics (MD) based alchemical free energy simulations. The calculated relative binding free energy of the analogs of A-484954 using the binding pose of hypothesis 1 show a good correlation with the experimental IC₅₀ values, yielding an r² coefficient of 0.96 after removing an outlier (compound 5). Calculations using another two poses show little correlation with experimental data, (r² of less than 0.5 with or without removing any outliers). Based on hypothesis 1, the calculated relative free energy suggests that bigger cyclic groups, at R1 e.g., cyclobutyl and cyclopentyl promote more favorable binding than smaller groups, such as cyclopropyl and hydrogen. Moreover, this study also demonstrates the ability of the alchemical free energy approach in combination with docking and homology modeling to prioritize compound synthesis. This can be an effective means of facilitating structure-based drug design when crystal structures are not available.

Introduction

Computer-based virtual screening (VS) approaches, including docking, pharmacophore, and similarity searching have been proposed and applied at the hit-identification stage of the costly drug discovery process (Singh et al., 2003; Ravindranathan et al., 2010; Kaoud et al., 2012; Wang et al., 2012; Hu et al., 2013; Rea et al., 2013; Sahner et al., 2013; Teli and Rajanikant, 2013; De Luca et al., 2014; Wang et al., 2014; Zhang et al., 2014). These approaches are only effective statistically, e.g., they only provide a meaningful prediction of active compounds in a certain percentile of a large number of samples, rather than information on specific compounds. Such inability hinders their applications in drug discovery processes where relative binding affinities of certain compounds are of interest, for example, in lead-optimization, when tens to hundreds of analogs need to be examined and prioritized for synthesis. Therefore, computational approaches that can correctly predict relative binding affinities are of crucial importance.

Relative free energy of binding as a computable property is a good measure of the binding affinity between two molecules. There is a history of applying physics-based force fields (Jorgensen et al., 1983; Kaminski et al., 2001; Ponder and Case, 2003; Wang et al., 2004; Jiao et al., 2008; Vanommeslaeghe et al., 2010; Ren et al., 2011; Shi et al., 2012; Zhang et al., 2012) or combined quantum mechanics/molecular mechanics (QM/MM) (Raha and Merz, 2005; Wang and Bryce, 2009; Hayik et al., 2010; Min et al., 2010; Wang and Bryce, 2011) methods to predict protein-ligand binding. Specifically, the accuracy of the alchemical free energy approaches in combination with well-tuned force fields has been demonstrated to have a reasonable agreement with experiment results. For example, the typical root mean square differences (RMSD) in hydration free energy of small organic molecules between experiment and prediction are around 1.0 kcal/mol using the commonly used fixed charge force fields, e.g., AMBER/GAFF, CHARMM, and OPLS-AA (Shirts et al., 2003; Shirts and Pande, 2005b; Mobley et al., 2009). While a RMSD of 0.7 kcal/mol has been reported on a smaller scale study of small organic molecules using the multipole-based polarizable AMOEBA force field (Ren et al., 2011). Compared to calculations of hydration free energy, absolute protein-ligand binding free energy calculation is more challenging as the introduction of the protein adds significantly more degrees of freedom to the system (Deng and Roux, 2006; Jayachandran et al., 2006; Mobley et al., 2007; Boyce et al., 2009). Fortunately, the absolute binding free energy is not required to predict the relative binding affinity of two compounds. Instead, the change in binding free energy (ΔΔG) between two compounds is often needed. By constructing an arbitrary thermodynamic path between the two end states (or ligands), the free energy change between them can be obtained by integrating the energy derivative along the path in the so-call thermodynamic integration method (TI) (Kollman, 1993; Shirts and Pande, 2005a). In addition to the TI method which is used in this study, the free energy perturbation (FEP) approach or Bennett acceptance ratio (BAR) are also commonly used in alchemical free energy calculations (Bennett, 1976; Kollman, 1993). In some earlier retrospective studies, the error of the calculated binding free energy is shown to be 1–2 kcal/mol (Jayachandran et al., 2006; Jiao et al., 2008; Michel and Essex, 2008; Boyce et al., 2009; Jiao et al., 2009; Rocklin et al., 2013). Prospective studies, e.g., calculation-driven inhibitor design, have also occasionally been reported in the literature (Jorgensen et al., 2006; Kim et al., 2006; Boyce et al., 2009; Bollini et al., 2011). Building on these recent advances in the theoretical development of the alchemical free energy calculation, we seek here to incorporate this approach to discovery potential inhibitors of eukaryotic elongation factor 2 kinase (eEF2K).

Eukaryotic elongation factor 2 kinase is believed to be a regulator of protein synthesis by phosphorylating the eukaryotic elongation factor 2 (eEF2) protein, which promotes the translocation of the ribosome along mRNA. Upon phosphorylation of eEF2 by the eEF2K, translation elongation is impeded (Carlberg et al., 1990; Kruiswijk et al., 2012). A recent study suggests a regulatory connection between eEF2K and nutrition deprivation. Cells undergoing nutrition deprivation survive by blocking the energy-demanding translation elongation process induced by eEF2K. However, the tumor cells can also exploit this pathway to survive from the metabolic stress (Leprivier et al., 2013). This makes eEF2K a potential target of therapeutic interest in drug discovery. Recently, a selective eEF2K inhibitor A-484954 (we refer to it as compound 3 in this study) with a submicromolar IC₅₀ value was reported (Chen et al., 2011). However, the binding mechanism of this compound with eEF2K remains unclear, in part due to the lack of a crystal structure of the kinase.

In this study, we explored the possibility of using in silico approaches to design inhibitors for novel targets that do not have crystal structures. Based on a homology model of eEF2K that we built earlier (Devkota et al., 2014), three hypothetical binding poses of A-484954 were first generated from docking. The relative binding free energies of seven novel analogs of A-484954 were calculated for each hypothetical pose using alchemical free energy approach. The predictions were subsequently compared and validated with the experiment IC₅₀ values we reported earlier (Edupuganti et al., 2014) although docking and alchemical free energy calculations were performed before the actual chemical and biochemical experiments. The computational results were utilized to prioritize the synthesis of the analog compounds in lead-optimization and provide a better understanding of the molecular interaction between eEF2K and the analogs. Based on the correlation between the calculation and experimental data, the most plausible binding mechanism of the compounds was also discussed.

Method

Structure Preparation and Docking

As no X-ray crystal structure for eEF2K is in the public domain, a homology model has been built in our group (Devkota et al., 2014) using the crystal structures of the alpha-kinase domain of myosin heavy chain kinase A (MHCKA, PDB ID: 3LKM) (Ye et al., 2010) and transient receptor potential (TRP) channels (ChaK) (PDB ID: 1IA9) (Yamaguchi et al., 2001). Based on this 3D model structure, compounds were docked into the ATP binding site of eEF2K using the ChemPLP (Korb et al., 2009) and Goldscore (Jones et al., 1995, 1997) scoring functions in the GOLD5.1 software package.

Free Energy Approach

To evaluate the change in the binding free energy between two analog compounds, a two-step free energy calculation scheme was applied. As shown in Figure 1, the change in the binding free energy between compounds A and B can be calculated either by ΔG₃ − ΔG₂ or by ΔG₄ − ΔG₁. Usually, the later formula is chosen since it is computationally more feasible and practical. To obtain ΔG₄ and ΔG₁, two independent perturbations between compounds A and B need to be conducted in the solvent environment (water in this case) with the protein present and absent respectively. The free energy change ΔG of each perturbation can be evaluated by an alchemical approach, such as the thermodynamic integration (TI) approach used in this study (Kollman, 1993).

FIGURE 1

Figure 1. Illustration of the thermodynamic cycle of the perturbation between the analog compounds A and B. Core represents the common part between two analogs. The two panels above represent the ligand in the unbound state in water, while the two panels below represent the ligand in the bound state in water.

For each pair of compounds, an identical protein and octahedral water box were constructed using the leap module in the AMBER12 software package (Case et al., 2012). A buffering region of 10 Å is used to solvate the protein-ligand complex and the ligand in the water box. This results a system of ~30,500 atoms for each protein-ligand complexes. The parameters for protein and water are taken from the ff99SB force field (Hornak et al., 2006) and the TIP3P water model (Jorgensen et al., 1983) respectively. The ligand parameters are obtained from GAFF (Wang et al., 2004) with the charges fitted from HF/6-31G* calculations. All the simulations were started with a quick minimization to remove the close contacts in the structure, followed by a 50 ps NVT simulation to heat the system up to 300 K and another 50 ps NPT simulation to equilibrate the density of the system, both with a time step of 1 fs. Production NVT simulations of 2–4 ns are then conducted for data collection with a time step of 2 fs. Periodic boundary condition and particle mesh Ewald were used to capture long-range effects. The thermodynamic integration along with a softcore potential implementation (Steinbrecher et al., 2011) in AMBER12 was applied to estimate the free energy. Each perturbation used 11 windows with λ values of 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 0.99, where electrostatic and van der Waals interactions were perturbed simultaneously. This saves considerable simulation time than perturbing electrostatic and van der Waals interaction separately. All the molecular dynamics (MD) simulations were performed using the AMBER12 software package (Case et al., 2012). Generally, a good convergence in the thermodynamic integration of the ligands in water can be obtained within 1 ns; in contrast, 2–3 ns are normally required for perturbations with the presence of the kinase using the current setting (see data in Supplementary Material). As a result, by using 11 nodes each with two, eight-core Xeon E5-2680 (Sandy Bridge) processors running at 2.7 GHz on Stampede supercomputer in Texas Advanced Computing Center, it takes about 1 day to obtain the free energy change between two ligands.

Preparation of the Compounds

To probe the structure-activity-relationship of the eEF2K inhibitor, a number of analog compounds were designed by replacing different chemical groups at the R1, R2, and R3 sites (Table 1). The analogs were designed based on the predicted importance of the chemical groups in binding from docking and chemical intuition, as well as the ease of chemical synthesis. As a result, three sites of A-484954 (compound 3) were modified with different substituents. The R1 (cyclopropyl) and R2 (ethyl) moieties were substituted by both smaller and bigger hydrophobic groups, including methyl, ethyl, propyl, cyclopropyl, cyclobutyl, and cyclopentyl. At the R3 site, the CONH₂ group in compound 3, which forms a key hydrogen-bond network with the hinge residues of eEF2K, was replaced by a less polar CSNH₂ group. Details about the synthesis and biochemical experiments can be found in reference (Edupuganti et al., 2014).

TABLE 1

Table 1. The analogs tested in free energy calculations.

Results and Discussion

Docking Results and the Hypothetical Binding Poses

Two poses were obtained from two independent docking runs (Figures 2A,B). In the first run, the ChemPLP (Korb et al., 2009) function was used as the primary scoring function with the Goldscore (Jones et al., 1995, 1997) function used for rescoring. By ranking with the ChemPLP score, 7 out of the top 10 poses from 20 generic algorithm optimizations agree with hypothesis 1 (Figure 2A), including the top 3 poses. The ChemPLP score for the 7 poses ranges from 68.6 to 67.3, while their corresponding Goldscore ranges from 46.9 to 37.7. In the second run, the two scoring functions were switched, i.e., the Goldscore was used as the primary scoring function, while ChemPLP used to rescore the poses. By ranking the poses with Goldscore, all the top 10 poses out of 20 generic algorithm optimizations agree with hypothesis 2 (Figure 2B) with Goldscores ranging from 53.6 to 52.8, while the ChemPLP scores range from 52.2 to 43.9 for the same poses. It should be noted that the scores between different scoring functions cannot be compared directly. Only the relative scores of the same scoring function can be compared, ideally the higher the better. The results indicate a clear difference in preference between the two scoring functions, i.e., ChemPLP favors hypothesis 1, while Goldscore favors hypothesis 2.

FIGURE 2

Figure 2. A-484954 (compound 3) and its hypothetical binding poses to eEF2K homology model structure; (A) is hypothesis 1; (B) is hypothesis 2; (C) is hypothesis 3; and (D) is the structure of A-484954 shown in green sticks. The hinge region of eEF2K is shown in purple, the glycine-rich loop is shown in blue, and the pocket deep inside the ATP adenine binding site is shown in gray surface.

In both hypothesis 1 and 2, three residues, K170, I232, and Y236, play a key role in the binding. Y236 is predicted to hold the ligand via pi-pi stacking with the pyridopyrimidine ring of compound 3 in both binding poses. K170 and I232 all form hydrogen bonds with the ligand, however, differently in the two poses. In hypothesis 1, the carbonyl group oxygen (O₂) and amine nitrogen (N₃) of the ligand form hydrogen bonds with the side chain of K170 and the backbone oxygen atom of I232 respectively. The cyclopropyl moiety of compound 3 fits into the deep cavity of the ATP-binding site of the kinase (shown in gray surface in Figure 2), while the ethyl group goes underneath the glycine-rich loop. In contrast, the hydrogen-bonding network changes to O of compound 3 with K170 and O₁ with the backbone nitrogen of I232 in hypothesis 2. The ethyl group still goes underneath the glycine-rich loop, however, the cyclopropyl group is pointing outwards of the pocket. In addition, residue G234 also contributes to the hydrogen bond network with the ligand in hypothesis 1.

Based on the visual inspection of hypotheses 1 and 2, we then deduced hypothesis 3 (Figure 2C), which was obtained by manually flipping the ligand around O₂ in hypothesis 1. This promotes the O and N₄ of the ligand forming the hydrogen bond network with the backbone nitrogen and oxygen of I232 respectively. Another difference between hypothesis 3 and 1 is that the position of the ethyl and cyclopropyl group is switched. The ethyl group flipped into the cavity inside the binding pocket and the cyclopropyl went to underneath the glycine rich loop.

To examine the stability of these hydrogen bonds in each poses, all three structures were minimized and relaxed by means of MD simulations for 10 ns. The measurement of the hydrogen bond distance suggests the hydrogen bonds predicted in docking are reasonably stable in general (Table 2). The average distances between the heavier atoms of the hydrogen bonds are generally around 3 Å. However, in hypothesis 3, the manually docked pose, the hydrogen bond between O₂ of the ligand and K170 diminished quickly during the MD simulation, resulting an average distance of 6.5 Å (Table 2).

TABLE 2

Table 2. Average hydrogen bond distances (in Å) between the kinase and compound 3 from 10 ns molecular dynamics simulations (measured between non-hydrogen atoms).

Determine the Binding Mode

To determine the binding pose of compound 3, the binding free energy change from compound 3 to its analogs are calculated and compared with the experiment IC₅₀ values (Table 3). As we have three hypothetical poses, the calculation was repeated for each of them. As free energy methods have been shown to have good accuracy in predicting the relative binding free energy, the hypothetical binding pose which is closer to the true binding mode should give the best correlation between the predicted binding free energy and the experimentally measured binding affinity.

TABLE 3

Table 3. Calculated binding free energy change (ΔΔG, kcal/mol) of the analogs against A-484954 and their experimental IC₅₀ values (μM).

It should be noted that the binding free energy calculated here is the relative binding free energy, i.e., the cost in free energy by substituting compound 3 with another compound. A positive value suggests a cost in energy, thus the substituted compound (i.e., 3) should be a stronger binder than the other. Conversely, a negative value means the substituent binds better. For example, in Table 3, from compounds 3 to 1, the predicted change in binding free energy is 2.0 kcal/mol based on hypothesis 1, thus compound 3 is predicted to have a higher binding affinity than compound 1.

To determine the binding mode of compound 3, each of the hypothetical binding poses need to be justified. In hypothesis 2 (Figure 2B), for example, the CONH₂ group forms a hydrogen bond with the hinge residue of eEF2K, which implies that the oxygen atom may be important in binding. Substituting CONH₂ (compound 3) with CSNH₂ (compound 5) yield a 6.9 kcal/mol energy cost based on the free energy calculation (Table 3). This is in accordance with the trend of measured IC50 values, which increases from 0.42 μM to >25 μM for compound 5. However, when the cyclopropyl group in compound 3 is substituted by methyl (compound 8), an increase in affinity would be expected as the cyclopropyl is exposure to solvent in hypothesis 2, substituting it with a smaller hydrophobic group might be more favorable. Not surprisingly, the calculated free energy change based on this pose suggests this is a favorable substitution, given a free energy change of -0.6 kcal/mol. However, the measured IC₅₀ value indicates this is an unfavorable change, where the IC₅₀ of compound 8 is at least 50-fold weaker than compound 3. Nonetheless, such results demonstrate a good prediction of the free energy change, but possibly indicate a bad binding pose predicted in hypothesis 2.

Similarly, when the ethyl group in compound 3 was substituted by the methyl group (compound 2), the free energy calculated based on hypothesis 3 suggested it is a favorable change (-0.5 kcal/mol), while the calculation based on hypothesis 1 suggested an unfavorable change (1.5 kcal/mol). Comparing the IC₅₀ value of 6.1 μM of compound 2 with 0.42 μM of compound 3, it is clear that the prediction based on hypothesis 1 is more accurate. To make a conclusion, a systematic comparison of the correlation between the calculated relative free energy and the experimental IC₅₀ values are conducted.

To facilitate the comparison, relative ln(IC₅₀) values were calculated based on the original experiment IC₅₀ data. For each compound, the natural logarithms were first obtained. Then the relative values were calculated against compound 3, to be consistent with the free energy calculations. As the IC₅₀ values for compounds 5 and 8 have been measured only to be >25 μM, they were all considered to have an IC₅₀ value of 25 μM to simplify the comparison. For compounds 6 and 7, due to the difficulties in synthesis, they were not included in the comparison. This results six data points for each of the hypothetical binding poses.

As shown in Figure 3, the calculated relative binding free energy was plotted against the experimental IC₅₀ values. A good correlation between the calculation and experiment is observed for hypothesis 1 in general (Figure 3). In contrast, neither hypothesis 2 nor 3 shows strong correlation between the predictions and the actual experiment measurements. If all six data points were considered, the r² correlation coefficient for hypothesis 1, 2, and 3 are 0.37, 0.08, and 0.31 respectively, while the Kendall's τ rank correlation coefficient are of 0.55, 0.14, and 0.55. All are not meaningful. However, if one outlier (compound 5) is removed in hypothesis 1, the best r² coefficient becomes 0.96. In contrast, the best r² coefficient of hypothesis 2 and 3 after removing any outliers are 0.42 (removing compound 8) and 0.49 (removing compound 8; removing compound 5 yields a r² of 0.13) respectively. The Kendall's τ rank correlation coefficient becomes 0.97, 0.69, and 0.55 (0.41) by removing the same compounds. With such a good correlation for hypothesis 1, it is likely that this pose is closer to the true binding mode of compound 3 in eEF2K.

FIGURE 3

Figure 3. Plots of the calculated relative binding free energy and the experimental relative ln(IC₅₀) values, i.e., compounds X—compounds 3 (X = 1, 2, 4, 5, 8), use each of the three hypothetical binding poses (dots for hypothesis 1, asterisks for hypothesis 2 and triangles for hypothesis 3). After removing one outlier (either compound 5 or 8), the best r² correlation coefficient for hypothesis 1, 2, and 3 are 0.96 (removing compound 5), 0.42 (removing compound 8), and 0.49 (removing compound 8; removing compound 5 yields a r² of 0.13), respectively, while the Kendall's τ rank correlation are 0.97, 0.69 and 0.55 (0.41), respectively, by removing the same compound in r² calculation.

Ranking and Prioritizing the Compounds

It is shown that hypothesis 1 may be the closest to the true binding mode of compound 3 with eEF2K in previous section. The relative potency of five (out of six) compounds is correctly predicted. The correlation between the prediction and experimental measurement of the five compounds is rather good, with an r² correlation coefficient of 0.96 and a Kendall's τ rank correlation of 0.97 (Figure 3). The calculation correctly predicted the structure-activity-relationship of the R1 and R2 sites of compound 3 (Tables 1, 3). At the R1 site, the calculation suggested that bigger hydrophobic groups, for example cyclobutyl and cyclopentyl tested in calculation, are generally more potent than smaller groups in the order of cyclopentyl (compound 7) > cyclobutyl (compound 6) > cyclopropyl (compound 3) > methyl (compound 8). The change in binding free energy from compound 3 to compounds 7, 6, and 8 are −2.6, −1.1, and 3.1 kcal/mol respectively. A negative change in the binding free energy suggests a favorable change, i.e., a more potent compound, while a positive change suggests an unfavorable change, i.e., a less potent compound. The prediction for the change from compound 3 to compound 8 has been verified by experimental IC₅₀ values of 0.42 and >25 μM for the two compounds respectively. However, for compound 6 and compound 7, experimental results are not obtained due to the difficulties in synthesis of the two compounds. Another compound, which uses ethyl to substitute cyclopropyl at the R1 site, has an IC₅₀ value between 1 ~ 2 μM (Edupuganti et al., 2014). This is in accordance with our predicted structure-activity-relationship.

At the R2 site, the calculated free energy suggested that the ethyl group of compound 3 is the optimal group. Replacing it with bigger or smaller hydrocarbon groups all lowers the potency of the compound. The calculated binding free energy change from substituting ethyl with hydrogen (compound 1), methyl (compound 2) and propyl (compound 4) are 2.0, 1.5, and 0.7 kcal/mol, respectively. This is in an excellent agreement with the experimental IC₅₀ values of 0.42, 6.6, 6.1, and 0.93 μM for compounds 3, 1, 2, and 4.

An earlier study (Lockman et al., 2010) also reported a compound (Figure 4B), which has similar scaffold as compound 3, favors relatively bigger groups (e.g., furan and isobutyl) than smaller groups (e.g., hydrogen atom and methyl) at R1 in general. After docking this compound to eEF2K in reference to the binding poses of compound 3 in hypothesis 1 and 2, two binding poses are obtained. We then calculated the relative binding free energy upon the substitution of the hydrogen atom with a furan group at R1 of the compound (Figure 4). The calculated changes in binding free energy are of −0.8 ± 0.4 and −2.6 ± 0.3 kcal/mol for the two poses, respectively. Given the experimental IC₅₀ values of 2.5 and 0.64 μM of the compounds with hydrogen and furan at R1, respectively, the calculated results based on hypothesis 1 (ΔΔG = −0.8 ± 0.4 kcal/mol) is clearly better correlated with the 4-fold increase in binding affinity measured in experiment.

FIGURE 4

Figure 4. Two eEF2K inhibitors, (A) A-484954 and (B) compound 2 in the reference (Lockman et al., 2010), (C) in complex with eEF2K with the conformations of hypothesis 1.

Problem in Free Energy Calculation between Compounds 3 and 5

As discussed in previous sections, the calculated change in the binding free energy from compound 3 to 5 in hypothesis 1 was the only outlier in comparison with the experimental IC₅₀ values (Figure 3). However, the independent MD simulations of compounds 3 and 5 clearly indicates that the later lost one key hydrogen bond between the N₃ atom of the compound and the Gly234 residue in the hinge region of the kinase (Figure 2 and Table 2). A short hydrogen bond distance of 2.9 Å between Gly234 and compound 3 suggests a strong interaction between the kinase and the ligand. However, after the CONH₂ group is substituted with CSNH₂, this hydrogen bond is completely lost, given an average distance of 5.2 Å in the MD simulation of compound 5 (Table 2). Quantum mechanical calculations (HF/6-31G*) for compounds 3 and 5 suggest that the substitution of CONH₂ with CSNH₂ is associated with significant change on the charge distribution of the atoms nearby. As a result, the restrained electrostatic potential (RESP) method fitted partial charges of the neighbor NH₂ group, as well as the carbon atoms, are less polar in compound 5 comparing to the charges in compound 3. The calculated electrostatic interaction energy between compound 5 and eEF2K is about 9.1 kcal/mol less favorable than that of compound 3, while the contribution from vdWs interaction of compounds 5 and 3 are about the same, of −40.5 and −40.8 kcal/mol, respectively (averaged energy from MD simulations). This shift in energy is in accordance with the observed shift in hydrogen bond distances. These evidences all suggest that compound 5 may be a weaker binder than 3. The failure of the free energy calculation between compounds 3 and 5 might be a combination of many factors, e.g., the force field and the free energy method. More studies are needed before a conclusion can be drawn. This, however, is out the scope of our study.

Conclusion

In this study, the binding mechanism of a known eEF2K inhibitor (Chen et al., 2011) (compound 3) has been studied using docking and alchemical free energy approaches in conjunction with experimental measurements. The inhibitor was firstly docked into the ATP-binding site of a homology model (Devkota et al., 2014) we built earlier. Then, three different binding poses, hypothesis 1, 2, and 3, were used to predict the structure-activity-relationship (SAR) of the inhibitor. The calculated relative binding free energy of analogs of compound 3 using the pose in hypothesis 1 shows a good correlation with the experimental IC₅₀ values, giving an r² coefficient of 0.96 after removing a suspicious outlier (compound 5). Calculations using another two poses merely show correlations with experiment, given the r² coefficients of <0.5 either or not remove any single outliers (Figure 3). Based on hypothesis 1, further free energy calculations suggest bigger cyclic groups, e.g., cyclobutyl and cyclopentyl, at R1 might be more favorable in binding than smaller groups, e.g., cyclopropyl and hydrogen. This is in accordance with previous study of an eEF2K inhibitor which has similar core scaffold as compound 3 (Lockman et al., 2010). At R2, the ethyl group might be close to the optimal size, as either longer or shorter side chains do not increase the potency of the compound. Moreover, this study also demonstrates the ability of the alchemical free energy approach in prioritizing the compound synthesis in combination with docking and homology modeling when experimental protein structures are not available, suggesting an alternative way to reduce the dependence of crystal structures in structure-based drug design, as well as the possibility to reduce the total number of compounds that need to be synthesized and tested.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

PR thanks Welch Foundation (F-1691), NIH (R01GM106137) and XSEDE (TG-MCB100057); KD. thanks NIH (GM059802 and CA167505), Welch Foundation (F-1390) and CPRIT (RP110532). The authors also thank the high performance computing support provided by Texas Advanced Computing Center.

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fmolb.2015.00009/abstract

Detailed free energy results is provided in the Supporting Material, including the perturbation from bounded (protein-ligand complex in water) and unbounded (ligand in water) states, as well as the convergence of the simulations.

References

Bennett, C. H. (1976). Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 22, 245–268. doi: 10.1016/0021-9991(76)90078-4

CrossRef Full Text | Google Scholar

Bollini, M., Domaoal, R. A., Thakur, V. V., Gallardo-Macias, R., Spasov, K. A., Anderson, K. S., et al. (2011). Computationally-guided optimization of a docking hit to yield catechol diethers as potent anti-HIV agents. J. Med. Chem. 54, 8582–8591. doi: 10.1021/jm201134m

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Boyce, S. E., Mobley, D. L., Rocklin, G. J., Graves, A. P., Dill, K. A., and Shoichet, B. K. (2009). Predicting ligand binding affinity with alchemical free energy methods in a polar model binding site. J. Mol. Biol. 394, 747–763. doi: 10.1016/j.jmb.2009.09.049

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Carlberg, U., Nilsson, A., and Nygård, O. (1990). Functional properties of phosphorylated elongation factor 2. Eur. J. Biochem. 191, 639–645. doi: 10.1111/j.1432-1033.1990.tb19169.x

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Case, D. A., Darden, T. A., Cheatham, I. T. E., Simmerling, C. L., Wang, J., and Kollman, P. A. (2012). AMBER12. San Francisco, CA: University of California.

Chen, Z., Gopalakrishnan, S. M., Bui, M.-H., Soni, N. B., Warrior, U., Johnson, E. F., et al. (2011). 1-Benzyl-3-cetyl-2-methylimidazolium iodide (NH125) induces phosphorylation of eukaryotic Elongation Factor-2 (eEF2): a cautionary note on the anticancer mechanism of an eEF2 kinase inhibitor. J. Biol. Chem. 286, 43951–43958. doi: 10.1074/jbc.M111.301291

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

De Luca, L., Ferro, S., Damiano, F. M., Supuran, C. T., Vullo, D., Chimirri, A., et al. (2014). Structure-based screening for the discovery of new carbonic anhydrase VII inhibitors. Eur. J. Med. Chem. 71, 105–111. doi: 10.1016/j.ejmech.2013.10.071

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Deng, Y., and Roux, B. (2006). Calculation of standard binding free energies: aromatic molecules in the T4 Lysozyme L99A mutant. J. Chem. Theory Comput. 2, 1255–1273. doi: 10.1021/ct060037v

CrossRef Full Text | Google Scholar

Devkota, A. K., Edupuganti, R., Yan, C., Shi, Y., Jose, J., Wang, Q., et al. (2014). Reversible Covalent Inhibition of eEF-2K by Carbonitriles. Chembiochem. 15, 2435–2442. doi: 10.1002/cbic.201402321

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Edupuganti, R., Wang, Q., Tavares, C. D. J., Chitjian, C. A., Bachman, J. L., Ren, P., et al. (2014). Synthesis and biological evaluation of pyrido[2,3-d]pyrimidine-2,4-dione derivatives as eEF-2K inhibitors. Bioorg. Med. Chem. 22, 4910–4916. doi: 10.1016/j.bmc.2014.06.050

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hayik, S. A., Dunbrack, R., and Merz, K. M. (2010). Mixed quantum mechanics/molecular mechanics scoring function to predict protein-ligand binding affinity. J. Chem. Theory Comput. 6, 3079–3091. doi: 10.1021/ct100315g

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hu, X., Compton, J. R., Abdulhameed, M. D. M., Marchand, C. L., Robertson, K. L., Leary, D. H., et al. (2013). 3-substituted indole inhibitors against Francisella tularensis FabI identified by structure-based virtual screening. J. Med. Chem. 56, 5275–5287. doi: 10.1021/jm4001242

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jayachandran, G., Shirts, M. R., Park, S., and Pande, V. S. (2006). Parallelized-over-parts computation of absolute binding free energy with docking and molecular dynamics. J. Chem. Phys. 125, 084901–084912. doi: 10.1063/1.2221680

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jiao, D., Golubkov, P. A., Darden, T. A., and Ren, P. (2008). Calculation of protein–ligand binding free energy by using a polarizable potential. Proc. Natl. Acad. Sci. U.S.A. 105, 6290–6295. doi: 10.1073/pnas.0711686105

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jiao, D., Zhang, J., Duke, R. E., Li, G., Schnieders, M. J., and Ren, P. (2009). Trypsin-ligand binding free energies from explicit and implicit solvent simulations with polarizable potential. J. Comput. Chem. 30, 1701–1711. doi: 10.1002/jcc.21268

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jones, G., Willett, P., and Glen, R. C. (1995). Molecular recognition of receptor sites using a genetic algorithm with a description of desolvation. J. Mol. Biol. 245, 43–53. doi: 10.1016/S0022-2836(95)80037-9

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jones, G., Willett, P., Glen, R. C., Leach, A. R., and Taylor, R. (1997). Development and validation of a genetic algorithm for flexible docking. J. Mol. Biol. 267, 727–748. doi: 10.1006/jmbi.1996.0897

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869

CrossRef Full Text | Google Scholar

Jorgensen, W. L., Ruiz-Caro, J., Tirado-Rives, J., Basavapathruni, A., Anderson, K. S., and Hamilton, A. D. (2006). Computer-aided design of non-nucleoside inhibitors of HIV-1 reverse transcriptase. Bioorg. Med. Chem. Lett. 16, 663–667. doi: 10.1016/j.bmcl.2005.10.038

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kaminski, G. A., Friesner, R. A., Tirado-Rives, J., and Jorgensen, W. L. (2001). Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides†. J. Phys. Chem. 105, 6474–6487. doi: 10.1021/jp003919d

CrossRef Full Text | Google Scholar

Kaoud, T. S., Yan, C., Mitra, S., Tseng, C.-C., Jose, J., Taliaferro, J. M., et al. (2012). From in silico discovery to intracellular activity: targeting JNK–protein interactions with small molecules. ACS Med. Chem. Lett. 3, 721–725. doi: 10.1021/ml300129b

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kim, J. T., Hamilton, A. D., Bailey, C. M., Domoal, R. A., Wang, L., Anderson, K. S., et al. (2006). FEP-guided selection of bicyclic heterocycles in lead optimization for non-nucleoside inhibitors of HIV-1 reverse transcriptase. J. Am. Chem. Soc. 128, 15372–15373. doi: 10.1021/ja066472g

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kollman, P. (1993). Free energy calculations: applications to chemical and biochemical phenomena. Chem. Rev. 93, 2395–2417. doi: 10.1021/cr00023a004

CrossRef Full Text | Google Scholar

Korb, O., Stützle, T., and Exner, T. E. (2009). Empirical scoring functions for advanced protein-ligand docking with PLANTS. J. Chem. Inf. Model. 49, 84–96. doi: 10.1021/ci800298z

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kruiswijk, F., Yuniati, L., Magliozzi, R., Low, T. Y., Lim, R., Bolder, R., et al. (2012). Coupled activation and degradation of eEF2K regulates protein synthesis in response to genotoxic stress. Sci. Signal. 5, ra40. doi: 10.1126/scisignal.2002718

PubMed Abstract | Full Text | CrossRef Full Text

Leprivier, G., Remke, M., Rotblat, B., Dubuc, A., Mateo, A.-R., Kool, M., et al. (2013). The eEF2 kinase confers resistance to nutrient deprivation by blocking translation elongation. Cell 153, 1064–1079. doi: 10.1016/j.cell.2013.04.055

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Lockman, J. W., Reeder, M. D., Suzuki, K., Ostanin, K., Hoff, R., Bhoite, L., et al. (2010). Inhibition of eEF2-K by thieno[2,3-b]pyridine analogues. Bioorg. Med. Chem. Lett. 20, 2283–2286. doi: 10.1016/j.bmcl.2010.02.005

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Michel, J., and Essex, J. W. (2008). Hit identification and binding mode predictions by rigorous free energy simulations. J. Med. Chem. 51, 6654–6664. doi: 10.1021/jm800524s

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Min, D., Zheng, L., Harris, W., Chen, M., Lv, C., and Yang, W. (2010). Practically efficient QM/MM alchemical free energy simulations: the orthogonal space random walk strategy. J. Chem. Theory Comput. 6, 2253–2266. doi: 10.1021/ct100033s

CrossRef Full Text | Google Scholar

Mobley, D. L., Bayly, C. I., Cooper, M. D., Shirts, M. R., and Dill, K. A. (2009). Small molecule hydration free energies in explicit solvent: an extensive test of fixed-charge atomistic simulations. J. Chem. Theory Comput. 5, 350–358. doi: 10.1021/ct800409d

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mobley, D. L., Graves, A. P., Chodera, J. D., Mcreynolds, A. C., Shoichet, B. K., and Dill, K. A. (2007). Predicting absolute ligand binding free energies to a simple model site. J. Mol. Biol. 371, 1118–1134. doi: 10.1016/j.jmb.2007.06.002

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ponder, J. W., and Case, D. A. (2003). Force fields for protein simulations. Adv. Protein Chem. 66, 27–85. doi: 10.1016/S0065-3233(03)66002-X

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Raha, K., and Merz, K. M. (2005). Large-scale validation of a quantum mechanics based scoring function: predicting the binding affinity and the binding mode of a diverse set of protein-ligand complexes. J. Med. Chem. 48, 4558–4575. doi: 10.1021/jm048973n

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ravindranathan, K. P., Mandiyan, V., Ekkati, A. R., Bae, J. H., Schlessinger, J., and Jorgensen, W. L. (2010). Discovery of novel fibroblast growth factor receptor 1 kinase inhibitor by structure-based virtual screening. J. Med. Chem. 53, 1662–1672. doi: 10.1021/jm901386e

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Rea, V. E. A., Lavecchia, A., Di Giovanni, C., Rossi, F. W., Gorrasi, A., Pesapane, A., et al. (2013). Discovery of new small molecules targeting the vitronectin binding site of the urokinase receptor that block cancer cell invasion. Mol. Cancer Ther. 12, 1402–1416. doi: 10.1158/1535-7163.MCT-12-1249

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ren, P., Wu, C., and Ponder, J. W. (2011). Polarizable atomic multipole-based molecular mechanics for organic molecules. J. Chem. Theory Comput. 7, 3143–3161. doi: 10.1021/ct200304d

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Rocklin, G. J., Boyce, S. E., Fischer, M., Fish, I., Mobley, D. L., Shoichet, B. K., et al. (2013). Blind prediction of charged ligand binding affinities in a model binding site. J. Mol. Biol. 425, 4569–4583. doi: 10.1016/j.jmb.2013.07.030

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Sahner, J. H., Groh, M., Negri, M., Haupenthal, J., and Hartmann, R. W. (2013). Novel small molecule inhibitors targeting the “switch region” of bacterial RNAP: structure-based optimization of a virtual screening hit. Eur. J. Med. Chem. 65, 223–231. doi: 10.1016/j.ejmech.2013.04.060

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shi, Y., Zhu, C. Z., Martin, S. F., and Ren, P. (2012). Probing the effect of conformational constraint on phosphorylated ligand binding to an SH2 domain using polarizable force field simulations. J. Phys. Chem. 116, 1716–1727. doi: 10.1021/jp210265d

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shirts, M. R., and Pande, V. S. (2005a). Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration. J. Chem. Phys. 122, 144107–144122. doi: 10.1063/1.1873592

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shirts, M. R., and Pande, V. S. (2005b). Solvation free energies of amino acid side chain analogs for common molecular mechanics water models. J. Chem. Phys. 122, 134508–134520. doi: 10.1063/1.1877132

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shirts, M. R., Pitera, J. W., Swope, W. C., and Pande, V. S. (2003). Extremely precise free energy calculations of amino acid side chain analogs: comparison of common molecular mechanics force fields for proteins. J. Chem. Phys. 119, 5740–5761. doi: 10.1063/1.1587119

CrossRef Full Text | Google Scholar

Singh, J., Chuaqui, C. E., Boriack-Sjodin, P. A., Lee, W.-C., Pontz, T., Corbley, M. J., et al. (2003). Successful shape-based virtual screening: the discovery of a potent inhibitor of the type I TGFbeta receptor kinase (TbetaRI). Bioorg. Med. Chem. Lett. 13, 4355–4359. doi: 10.1016/j.bmcl.2003.09.028

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Steinbrecher, T., Joung, I., and Case, D. A. (2011). Soft-core potentials in thermodynamic integration: comparing one- and two-step transformations. J. Comput. Chem. 32, 3253–3263. doi: 10.1002/jcc.21909

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Teli, M. K., and Rajanikant, G. K. (2013). Computational repositioning and experimental validation of approved drugs for HIF-Prolyl hydroxylase inhibition. J. Chem. Inf. Model. 53, 1818–1824. doi: 10.1021/ci400254a

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Vanommeslaeghe, K., Hatcher, E., Acharya, C., Kundu, S., Zhong, S., Shim, J., et al. (2010). CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 31, 671–690. doi: 10.1002/jcc.21367

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, J., Chen, L., Sinha, S. H., Liang, Z., Chai, H., Muniyan, S., et al. (2012). Pharmacophore-based virtual screening and biological evaluation of small molecule inhibitors for protein arginine methylation. J. Med. Chem. 55, 7978–7987. doi: 10.1021/jm300521m

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, Q., and Bryce, R. A. (2011). Accounting for non-optimal interactions in molecular recognition: a study of ion-[small pi] complexes using a QM/MM model with a dipole-polarisable MM region. Phys. Chem. Chem. Phys. 13, 19401–19408. doi: 10.1039/c1cp21944h

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, Q., Park, J., Devkota, A. K., Cho, E. J., Dalby, K. N., and Ren, P. (2014). Identification and validation of novel PERK inhibitors. J. Chem. Inf. Model. 54, 1467–1475. doi: 10.1021/ci500114r

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, Q. T., and Bryce, R. A. (2009). Improved hydrogen bonding at the NDDO-type semiempirical quantum mechanical/molecular mechanical interface. J. Chem. Theory Comput. 5, 2206–2211. doi: 10.1021/ct9002674

CrossRef Full Text | Google Scholar

Yamaguchi, H., Matsushita, M., Nairn, A. C., and Kuriyan, J. (2001). Crystal structure of the atypical protein kinase domain of a TRP channel with phosphotransferase activity. Mol. Cell 7, 1047–1057. doi: 10.1016/S1097-2765(01)00256-8

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ye, Q., Crawley, S. W., Yang, Y., Cote, G. P., and Jia, Z. (2010). Crystal structure of the {alpha}-kinase domain of Dictyostelium myosin heavy chain kinase A. Sci. Signal. 3, ra17. doi: 10.1126/scisignal.2000525

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Zhang, J., Yang, W., Piquemal, J.-P., and Ren, P. (2012). Modeling structural coordination and ligand binding in zinc proteins with a polarizable potential. J. Chem. Theory Comput. 8, 1314–1324. doi: 10.1021/ct200812y

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Zhang, Z., Martiny, V., Lagorce, D., Ikeguchi, Y., Alexov, E., and Miteva, M. A. (2014). Rational design of small-molecule stabilizers of spermine synthase dimer by virtual screening and free energy-based approach. PLoS ONE 9:e110884. doi: 10.1371/journal.pone.0110884

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar