Directly Binding Rather than Induced-Fit ... - ACS Publications

Article pubs.acs.org/JCTC

Directly Binding Rather than Induced-Fit Dominated Binding Affinity Difference in (S)- and (R)‑Crizotinib Bound MTH1 Huiyong Sun,† Pengcheng Chen,†,∥ Dan Li,† Youyong Li,§ and Tingjun Hou*,†,‡ †

College of Pharmaceutical Sciences, ‡State Key Lab of CAD&CG, and ∥College of Life Sciences, Zhejiang University, Hangzhou, Zhejiang 310058, China § Institute of Functional Nano & Soft Materials (FUNSOM), Soochow University, Suzhou, Jiangsu 215123, China S Supporting Information *

ABSTRACT: As one of the most successful anticancer drugs, crizotinib is found to be efficient in the suppression of MTH1, a new therapeutic target for RAS-dependent cancers. Deep analysis shows that stereospecificity is prevalent in the binding of crizotinib to MTH1, where the target is more preferred to bind with the (S)-enantiomer of crizotinib. Surprisingly, very similar binding modes were found for the two enantiomers (Huber et al. Nature 2014, 508, 222−227), which puzzled us to ask a question as to why such a subtle structural variation could lead to so large of a binding affinity difference. Thereafter, by using advanced all-atom molecular dynamics simulations, we characterized the free energy surfaces of the binding/unbinding processes of the (S) and (R)-crizotinib enantiomers to/from MTH1. Interestingly, we found that rather than the induced-fit process, which is prevalent in drug selectivity and specificity (Wilson et al. Science 2015, 347, 882−886), the directly binding process has dominated impact on the binding affinity difference of the enantiomers, implying a common mechanism of stereoselectivity of enantiomers.

1. INTRODUCTION

2. RESULTS AND DISCUSSION 2.1. Key Features Reproduction of the apo and holo MTH1. In this study, the Amber99SB-ILDN force field was employed for the MTH1 system simulations due to its good performance in protein conformational sampling.6 It was found that the binding loop (residues 23−30) in the apo-state MTH1 (PDB code: 3ZR1)7 prefers to stay in a more closed state (orange, green, and cyan cartoon models were colored for the apo, (S)-crizotinib-bound, and (R)-crizotinib-bound MTH1, respectively, in Figure 1 top panel). To test whether the Amber99SB-ILDN force field can discriminate the two states of the binding loop (holo-state and apo-state), conventional MD simulation was first performed to generate the distribution of the opening degree of the binding loop in the holo-state MTH1 and the apo-state MTH1. Herein, four systems were performed with the conventional MD simulation, namely, (S)-holo-MTH1, (S)-apo-MTH1, (R)holo-MTH1, and (R)-apo-MTH, where the apo-state MTH1s were constructed by eliminating the cocrystallized ligands in the corresponding holo-state MTH1s. As shown in Figure S1A− S1D, the systems are all stable within the 40 ns MD simulations with two stable hydrogen bonds between MTH1 and crizotinib (H-bonds between D119 and N2 atom of crizotinib and D120 and N1 atom of crizotinib), RMSD = ∼1 Å for the protein, and RMSD = ∼2 Å for the ligand (Figure S2), indicating that it is feasible to analyze the distribution of the binding loop based on

Crizotinib may be one of the most successful anticancer drugs in the past decade, which exhibits a high response ratio to the suppression of several associated drug targets, such as c-MET,1 EML4-ALK,2 CD74-ROS1,3 and the newly found MTH1.4 Compared with the (R)-enantiomer of crizotinib, which is preferred to bind with the tyrosine kinases (such as c-MET, EML4-ALK, and CD74-ROS1), the (S)-enantiomer is more favorable in binding with MTH1.5 Although similar binding poses were revealed by the crystal structures of (R) and (S)crizotinib-MTH1 complexes (PDB codes: 4C9W and 4C9X), the structure of (R)-crizotinib-MTH1 seems less stable because part of the ligand (pyridine ring) and the side chain of a nearby residue (F27) were unresolved.5 However, careful observation of the crystal structure clearly disclosed that an additional hydrogen bond exists between (R)-crizotinib and residue N33 (Figure S1B in the Supporting Information), implying that the (R)-enantiomer may also form a stable interaction with MTH1. Meanwhile, several residues close to (S)-crizotinib (such as F72, M81, and F139) were crystallized with multiple conformations, meaning that the binding of the (S)-enantiomer to MTH1 is not too stable either. Thereafter an intrinsic question is proposed: why so similar binding modes of the two enantiomers result in so different binding affinities. Hence, by using advanced all-atom molecular dynamics (MD) simulations, we revealed the binding/unbinding processes of the two crizotinib enantiomers to/from MTH1 in detail. © 2016 American Chemical Society

Received: October 13, 2015 Published: January 14, 2016 851

DOI: 10.1021/acs.jctc.5b00973 J. Chem. Theory Comput. 2016, 12, 851−860

Article

Journal of Chemical Theory and Computation

simulation results from FWTM (based on umbrella sampling, US; and adaptive biasing force, ABF). It is well-known that, for the systems of drug-target interactions, the slowest degree of freedom usually incorporates drug’s binding/unbinding process and large conformational change of the target, such as the induced-fit phenomenon. Hence, the free energy surfaces were designed with two collective variables (CVs): (1) the center of mass (CoM) distance between the heavy atoms of the crizotinib enantiomers and the active pocket of MTH1 (residues within 5 Å of the ligand in the holo-state MTH1) and (2) the loop-sheet distance between the binding loop (residues 23−30) and the β-sheet plane (formed by the Cα atoms of the residues I70, F72, and M81), which represents the opening degree of the binding pocket (also measured by the CoM distance). As shown in Figures 2D and 3D, the black dotted lines denote the favorable binding pathways of the crizotinib enantiomers binding to MTH1 (a 2-D description can be found in Figures 2F and 3F). Consistent with the low opening degree of the binding loop in the apo crystal structure, the loop stays more favorably in the closed state when crizotinib dissociates from the binding pocket of MTH1 (pocket opening degree = 13−14 Å at 15−23 Å of the reaction coordinate, RC; the structural description can be found in Figures 2C and 3C). Besides, the frequently getting into and out of the binding pocket for the crizotinib enantiomers (Figures 2E and 3E; the transition state can be found in Figures 2B and 3B) means that the FESs are convergent for analyses (which has also been validated by the convergent PMF curves in Figure S3A and S3B). Thereafter, the favorable binding pathway can ideally be divided into two independent phases as proposed by Wilson recently:12 namely, (1) the directly binding process (the red dotted lines in Figures 2D and 3D), where the protein target keeps in its closed apo-state and is unfavorable for drug binding (Figures 2A′ and 3A′), and (2) the induced-fit process (the green dotted lines in Figures 2D and 3D), where the binding loop goes from the closed apo-state to the opened holo-state to fit the drug binding (Figures 2A and 3A). By using the convergent FESs, we calculated the binding free energies of the two crizotinib enantiomers along the three pathways (the favorable binding pathway, ΔGMeta‑depth; the directly binding pathway, ΔG1; and the induced-fit pathway, ΔG2). As illustrated in Table 1, the absolute binding free energy of (S)-crizotinib based on FWTM (ΔGbind‑Meta) is much larger than that of (R)-crizotinib (−9.98 kcal/mol versus −5.88 kcal/ mol, ΔΔG = 4.10 kcal/mol), which is reasonably consistent with the experimental data (Table 1).5 It may be a concern that the calculated absolute binding free energies based on FWTM are a bit too lower compared with the experimental data, especially for the system of (R)crizotinib-MTH1 (ΔΔGcal‑exp = 2.74 kcal/mol). This is because the correction of the PMF depth (ΔGMeta‑depth) of 1.38 kcal/ mol as using a restrained funnel. However, as we know, even without any additional restraints in a system, the small ligand will not sample all the reaction-space within limited simulation time. Therefore, it may be too large to minus 1.38 kcal/mol (corresponding to 7.5 Å restrained funnel, eq 2) for the metadynamics based simulation. On the contrary, the PMF depth of FWTM (ΔGMeta‑depth) is closer to the experimental data (with the ΔΔGcal‑exp of −1.03 and 1.36 kcal/mol for (S) and (R)-crizotinib-MTH1, respectively), which is more suitable for the binding pathways discussion. Nevertheless, the calculated binding free energies based on FWTM are all

Figure 1. Validation of the employed force field. The crystallized (S)crizotinib-bound, (R)-crizotinib-bound, and free binding loops in MTH1 are shown in green, cyan, and orange cartoon models, respectively, in the top panel. The distribution of the opening degree of the binding loops is illustrated in the bottom panel, where the free MTH1s were simulated by eliminating the crystallized ligands. It can be found that the same ranking of the most populated opening degree of the binding loops is shown between those in the crystal structures.

the MD trajectories. The reason why we used two apo-state MTH1s derived from holo-state structures for the assessment is that we would use the holo-state proteins to mimic the apo-state proteins in the following enhanced sampling simulations. Therefore, it is necessary to detect whether the employed force field could reproduce the key features of the apo-state protein from a holo one. As shown in Figure 1 bottom panel, the ranking of the most populated opening degree of the binding pocket is consistent with the loop location in the crystal structures ((S)-holo > (R)-holo > apo), implying that the force field is reliable for the studied systems. Thereafter, the holo-state MTH1s (including (S)-crizotinib-MTH1 and (R)crizorinib-MTH1) were employed as the initial structures for the enhanced sampling simulations (including metadynamics and umbrella sampling simulations). 2.2. Absolute Binding Free Energies Calculated by Funnel-Based Well-Tempered Metadynamics and Multidimensional-Constrained Umbrella Sampling. Well-tempered metadynamics has been found effective in characterizing free energy surface (FES) in high dimensional reaction space.8,9 To accelerate the sampling process in associated reaction space, here, well-tempered metadynamics in conjunction with the funnel-restrained algorithm (funnel-based well-tempered metadynamics, FWTM)10 was used to characterize the FES of the binding/unbinding processes of the two crizotinib enantiomers to/from MTH1. Moreover, to guarantee the reliability of the predictions, absolute binding free energy calculation based on Woo and Roux’s scheme11 was also employed to validate the 852


Article


Figure 2. Free energy surface of (S)-crizotinib binding/unbinding from the active pocket of MTH1 (panel D). The favorable binding pathway, the directly binding pathway, and the induced-fit pathway are shown in black, red, and green dotted lines, respectively, in panel D. The corresponding 1D free energy profiles are illustrated in panels F, H, and G, respectively. The time evolution of the (S)-crizotinib binding/unbinding process is illustrated in panel E. Structural descriptions of the favorable binding state (corresponding to the large green point or point A in panels D−H), directly binding state (corresponding to the large red point or point A′ in panels D, G, and H), transition state (corresponding to point B in panels D−F), and free-moving state (corresponding to point C in panels D−F) of the drug are illustrated in panels A−D, respectively. The contours in panel D are spaced at an interval of 1.0 kcal/mol.

2.3. Directly Binding Effect Dominated Binding Affinity Difference upon Crizotinib Enantiomers Binding to MTH1. Recently, Wilson and colleagues12 found that drug selectivity (the same drug targeting homologous targets) is dominated by the induced-fit process during drugs binding, where the directly binding free energies of gleevec (a stereospecific drug targeting Abl) to the Abl and Src kinases are almost the same while the induced-fit free energies, on the contrary, are quite different, which can be explained by the fact that gleevec cannot induce an appropriate P-loop conformation to form favorable interaction with the Src kinase. Interestingly, here we found that, in addition to the enantiomers binding to the same target, the difference of the binding free energies is dominated by the directly binding process rather than the induced-fit process (Table 1). As mentioned above, the favorable binding pathway in FWTM can be divided into two stages (the directly binding stage and the induced-fit stage), where we can imagine a scenario that the ligand at first binds to the binding pocket

reasonably consistent with the experimental data with the binding affinity of (R)-crizotinib much lower than that of (S)crizotinib in MTH1. To further check whether the calculated binding free energies are stable or method dependent, Woo and Roux’s absolute binding free energy calculation methodology was also used to validate the results predicted by FWTM. The conformations involved in the favorable binding pathway in FWTM were collected for the US simulations.13,14 Here, 12 ns US simulation for each window (0.5 Å/window) was employed to guarantee the convergence of the separation PMFs (Figure S3C and S3D, a detailed energetic decomposition can be found in Table S1 and Figure S4). As shown in Table 1, the calculated binding free energies based on Woo and Roux’s scheme (ΔGbind‑US) are very close to the PMF depths calculated by FWTM (−10.77 kcal/mol versus −11.36 kcal/mol and −7.73 kcal/mol versus −7.26 kcal/mol for (S) and (R)-crizotinib, respectively), meaning that the predicted results are reliable and methodindependent. 853


Article


Figure 3. Free energy surface of (R)-crizotinib binding/unbinding from the active pocket of MTH1 (panel D). The favorable binding pathway, the directly binding pathway, and the induced-fit pathway are shown in black, red, and green dotted lines, respectively, in panel D. The corresponding 1D free energy profiles are illustrated in panels F, H, and G, respectively. The time evolution of the (R)-crizotinib binding/unbinding process is illustrated in panel E. Structural descriptions of the favorable binding state (corresponding to the large green point or point A in panels D−H), directly binding state (corresponding to the large red point or point A′ in panels D, G, and H), transition state (corresponding to point B in panels D−F), and free-moving state (corresponding to point C in panels D−F) of the drug are illustrated in panels A−D, respectively. The contours in panel D are spaced at an interval of 1.0 kcal/mol.

(−5.82 kcal/mol), the (R)-crizotinib is very unfavorable in binding with the apo-state MTH1 (−1.25 kcal/mol). Structural analysis gives details upon the directly binding state of the two enantiomers. As has been mentioned above, there should be no conformational change of the target when binding with a ligand in its apo-state. As such, the apo-state conformation of MTH1 was considered as the conformation for directly binding of the crizotinib enantiomers. To analyze the interactions between the crizotinib enantiomers and MTH1 upon the directly binding state, we mapped the crystallized bound-state crizotinib enantiomers into the crystallized apostate MTH1 by overlapping the holo-state and apo-state MTH 1s. As shown in Figure 4C1 and 4C2, besides the binding loop, two residues, N33 and W117, are also involved in large conformational change during the ligands binding. As mentioned above, the only difference of the crizotinib enantiomers is the position of the methyl group (red circle in Figure 4A1 and 4A2). Deep analysis shows that the methyl

without affecting the conformation of the binding pocket of the target (ΔG1, Figures 2A′ and 2H, and Figures 3A′ and 3H) and then induces the conformational change of the binding pocket to fit the ligand binding (ΔG2, Figures 2A and 2G and Figures 3A and 3G). Because the favorable binding free energy in FWTM is just composed of the free energy difference of the two stages (ΔGMeta‑depth = ΔG1+ΔG2), we analyzed the free energy difference between the two crizotinib enantiomers on the level of the two binding pathways. As shown in Table 1, very similar binding free energies were generated at the stage of induced-fit for the two enantiomers binding (−5.54 kcal/mol and −6.01 kcal/mol for (S) and (R)-crizotinib, respectively), meaning that there is little difference of the induced-fit phenomenon when the enantiomers bind to the same target, which is typically different from drug selectivity12 and specificity15 (dominated by induced-fit phenomenon). On the contrary, large binding free energy difference was shown of the directly binding stage, where, compared with the (S)-crizotinib 854


Article


rotates outward during the induced-fit process (the green sticks model in Figure 4A1). Whereas, the methyl group in (R)crizotinib seems very unfavorable in binding with the apo-state MTH1 because the methyl group locates too close to two residues (N33 and W117, Figure 4C2). Although a new hydrogen bond is formed between the ND2 atom of residue N33 and the O1 atom of (R)-crizotinib due to the rotation of the methyl group in (R)-crizotinib (Figure S1B, the green dotted line, which is not found in (S)-enantiomer-MTH1, Figure S1A and S1E), two relatively stable conformations of N33 were still found in the crystallized holo-state MTH1 (the green and cyan sticks model in Figure 4A2), meaning that residue N33 is unfavorable in sitting where it is in the apo-state MTH1 either. Besides, due to the too crowded surroundings, the newly formed hydrogen bond is even forced to increase from 1.65 Å (in the directly binding state) to 1.88 Å (in the favorable-binding state) to reduce the steric hindrance resulting from the unfavorable position of the (R)-methyl (Figure 4A2). Moreover, the (R)-methyl is very close to W117 as well (Figure 4B2), which leads to the decrease of rotation of W117 from 23° (in (S)-crizotinib-MTH1, Figure 4B1) to 16° (in (R)crizotinib-MTH1, Figure 4B2) to relax the too crowded surroundings around the (R)-methyl group. Therefore, the directly binding process is very unfavorable for the (R)crizotinib binding. Although it may be argued that the much closed binding loop in (R)-crizotinib-MTH1 may lead to the decreased induced-fit energy compared with the (S)-complex, the newly formed hydrogen bond between N33 and (R)crizotinib may facilitate the increased induced-fit energy to the system as well, which leads to the comparative induced-fit free energies for the two enantiomers (Table 1).

Table 1. Free Energy Decomposition of the Absolute Binding Free Energy (kcal/mol) free energy decomposition

(S)-crizotinib

(R)-crizotinib

Funnel Based Well-Tempered Metadynamics ΔG1b −5.82 ± 0.14f −1.25 ± 0.74 ΔG2c −5.54 ± 0.18 −6.01 ± 0.68 ΔGMeta‑depthd −11.36 ± 0.11 −7.26 ± 0.32 ΔGbind‑Metae −9.98 ± 0.11 −5.88 ± 0.32 Woo and Roux’s Scheme Based on Umbrella Sampling ΔGbind‑USg −10.77 ± 1.09h −7.73 ± 1.39 Experimental Data ΔGexp‑1i −9.66 −8.02 ΔGexp‑2j −10.33 −8.62

ΔΔGa 4.57 −0.47 4.10 4.10 3.04 1.64 1.71

a Free energy difference calculated by ΔΔG = ΔG(R)‑crizotinib − ΔG(S)‑crizotinib. bDirectly binding free energy, which corresponds to free energy difference between the apo-state and the unfavorable binding state of crizotinib. cInduced-fit binding free energy, which corresponds to the energy difference between the unfavorable binding state and favorable binding state of crizotinib. dTotal binding free energy between the apo-state and the favorable binding state of crizotinib, where ΔGMeta‑depth = ΔG1 + ΔG2. eAbsolute binding free 1 2 ο energy calculated by ΔGbind‑Meta = ΔGMeta‑depth − β ln(πRcyl C ) (β =

1 ), kBT

where the radius was set to 7.5 Å for free rotation of crizotinib in

the bulk. Therefore, the analytical correction corresponds to 1.38 kcal/ mol because of the restrained funnel. fDeviation estimated based on the favorable binding pathway from 15 to 23 Å of the RC based on the last 10 ns metadynamics simulation (with each nanosecond simulation as a sample). gAbsolute binding free energy calculated by ΔGbind‑US = bulk

− β ln(KbindC ο) (Kbind = S*I*e−β[Gc 1

+Gobulk−Gtθ−Gϕt −G0α−G0β−G0γ −Gcsite]

). The

h

energetic components can be found in Table S1. Deviation estimated based on the last 9 ns US simulation from 15 to 20 Å of the RC (with each nanosecond simulation as a sample). iExperimental binding free energy derived from ref 5 at 283 K (∼15 °C). jExperimental binding free energy estimated by ΔGexp = −RT ln KKd based on Kd values at 310 K from ref 5.

3. CONCLUDING REMARKS By using advanced all-atom free energy calculation approaches, we uncovered the mechanism of binding specificity of crizotinib enantiomers to MTH1. We found that different from drug selectivity and specificity, which are typically dominated by an induced-fit process,12,15 the binding affinities difference of the enantiomers may be primarily determined by the directly binding process.

group in (S)-crizotinib locates only close to one residue (N33) in its directly binding state (the orange sticks model of N33, Figure 4C1). To remove the local steric hindrance, the residue

Figure 4. Structural description of the induced-fit mechanism in the two systems (crystal structures are used for the observation). The (S)-crizotinib and (R)-crizotinib systems are shown in panels C1 and C2, respectively, where the orange and green models (sticks and cartoon models) represent the apo-state and the holo-state (induced-fit state, including the cyan stick model of residue N33) of MTH1, respectively. The clear models of the markedly influenced residues (N33 and W117) are shown in panels A1 and B1 for the (S)-enantiomer and A2 and B2 for the (R)-enantiomer, respectively. 855


Article


was not protonated in the apo MTH1), respectively. Cubic TIP3P21 water boxes were added to the systems with the water molecules 10 Å out of the solutes in each direction. 4.2. Conventional MD Simulation. Before the enhanced sampling simulations, conventional MD simulations were performed to assess the performance of the Amber99SBILDN force field on the investigated systems. The cutoff value was set to 10 Å to deal with the short-range electrostatic and van der Waals interactions, and the particle mesh Ewald (PME)22 algorithm was used to handle the long-range electrostatic interactions. All covalent bonds involving hydrogen atoms were restrained with the SHAKE algorithm.23 Langevin dynamics was used to control the temperature of the system, where the Langevin damping coefficient (gamma) was set to 1/ps. Before MD simulation, a four step minimization strategy was used to optimize the systems. At first, all the heavy atoms were restrained with only hydrogen atoms capable of free moving (5000 steps). Second, heavy atoms incorporated in the solvent (including water molecules and ions) were set free as well (5000 steps). Third, heavy atoms in the side chain of protein residues were relaxed (5000 steps) in addition. Finally, the whole system was minimized for 20,000 steps without any restraints. In the stage of MD simulation, the systems were at first gradually heated from 0 to 310 K within 1 ns simulation time in an NVT ensemble (with 5 kcal/mol·Å2 restraints on the heavy atoms of the protein backbone). Then, the systems were gradually relaxed with the restraints decreased from 5 to 0 kcal/ mol·Å2 within 0.2 ns simulation time in an NPT ensemble (P = 1 atm and T = 310 K), where the Poisson Piston algorithm was used for the pressure control.24 At last, a 40 ns production run was performed for each system for the key state analysis, where the collective interval was set to 5 ps (namely, 200 snapshot/ ns), and a total of 8000 snapshots were collected for each system. The time step was set to 2 fs. All the simulations were performed with NAMD version 2.9.25 4.3. Metadynamics Based Absolute Binding Free Energy Calculation. With the rapid free energy surface (FES) characterization capability, metadynamics has become one of the most successful techniques in today’s enhanced sampling simulations. 9,26 In conjunction with different algorithms, standard metadynamics has evolved into various on-the-fly versions of advanced metadynamics, such as welltempered metadynamics,8 funnel metadynamics,10 reconnaissance metadynamics,27 bias-exchange metadynamics,28 and parallel tempering metadynamics,29 etc. Compared with the standard version of metadynamics, the newly developed methodologies are all generated with the same features, namely, accelerated convergence rate of the FES. Therefore, nowadays we see a wide use of metadynamics based techniques in the scopes of protein−protein interaction,30,31 drug-target interaction,32,33 and protein conformational transition or activation,34,35 etc. In the spirit of metadynamics, it adds a history-dependent Gaussian repulsive potential on the original potential of the system as shown in eq 1. Here, the well-tempered algorithm is used as an example

During the enantiomers binding, similar induced-fit free energies of the two crizotinib enantiomers (ΔG2, Table 1) could be observed, which is consistent with the structural features of the crystallized holo-state MTH1s that the similar binding pose of the crizotinib enantiomers and similar induced conformation of the binding loop of MTH1 resulted in at the drugs binding. On the contrary, due to the too crowded surroundings of the apo-state binding pocket of MTH1, the (R)-enantiomer of crizotinib seems very unfavorable in sitting with two residues (N33 and W117) too close to the methyl group, where the two residues both touch or insert into the van der Waals surface of (R)-crizotinib (pink surface in Figure 4A2 and 4B2), meaning that it is indeed very unfavorable for (R)crizotinib binding with the apo-state MTH1 (namely, much lower directly binding free energy resulted compared with that of (S)-crizotinib, ΔG1 in Table 1). Taken together, considering the same binding modes and similar induced conformation of the binding pocket usually occur in the enantiomers binding, the difference of the directly binding free energies may have a dominated effect on the binding specificity of the enantiomers.

4. METHODOLOGIES IN DETAIL 4.1. Structure Preparation. To accurately reproduce the key states of the studied systems, crystal structures of MTH1 complexed with (S)-crizotinib (PDB code 4C9X, resolution 1.20 Å)5 and (R)-crizotinib (PDB code 4C9W, resolution 1.65 Å)5 were employed as the initial structures for the holo-state MD simulations. Due to the imperfect crystallization of (R)crizotinib-MTH1, where the side chain of F27 (close to the ligand) is missing, we modeled the side chain of F27 by considering the position of this residue in (S)-crizotinibMTH1. As shown in Figure S5, due to the more closed state of the binding loop in the (R)-crizotinib bound MTH1, it is very unfavorable for F27 (cyan sticks model) staying where it is in the (S)-crizotinib-MTH1 (the orange sticks model). Thereby, the side chain of F27 can only be rotated outside of the binding pocket to reduce the local steric hindrance (the green sticks model). PROPKA (version 3.1)16 was used to determine the protonation state of the protein residues (such as histidines and cysteines). Specially, residue D119 was also protonated to keep stable the hydrogen bond between D119 and the crizotinib enantiomers (Figure S1C). To obtain the parameters for the small molecules simulations, we first optimized the ligand structures and calculated the electrostatic potentials at the Hartree−Fock 631G* level of theory using the Gaussian 09 program.17 Afterward, the restrained electrostatic potential technique (RESP)18 in the Ambertools 1.5 package19 was employed to fit the atomic partial charges of the optimized structures. The Amber99SB-ILDN force field6 and the General AMBER Force Field (GAFF)20 were employed for the protein and ligand simulations, respectively (the parameter files for (S) and (R)crizotinib can be found in the Supporting Information). To test whether the Amber99SB-ILDN force field can discriminate different conformational states of the studied system (especially the conformation of the binding loop), the crystallized ligands ((S) and (R)-crizotinib) were eliminated from the corresponding holo-state MTH1 to construct the apo-state initial structures. The crystallized apo MTH1 (PDB code 3ZR1, resolution 1.90 Å)7 was employed to give a comparison of the simulated result. Counterions of 6 Na+ and 5 Na+ were added to neutralize the holo-state and apo-state systems (residue D119

⎛ (s − s(t ′))2 ⎞ ωeV (s , t ′)/ ΔT exp⎜− ⎟ 2σ 2 ⎝ ⎠ t ′= τ ,2τ ,... t

V (s , t ) =

∑

(t ′ < t ) (1)

where V(s, t) denotes the history-dependent biasing potential, and t is the cumulative time. In metadynamics simulation, a 856


Article

Journal of Chemical Theory and Computation Gaussian potential with the hill height ωeV(s,t′)/ΔT will be added to the concurrent position s(t′) of the biased molecule in each time interval τ, where ΔT is a bias-factor to adjust the change rate of the hill height in the well-tempered algorithm. In welltempered metadynamics, the hill height will be gradually scaled down with the exponential of V(s, t′)/ΔT, which is superior to the standard version of metadynamics that uses an invariable hill height (ω) for the simulation (and typically needs a very long simulation time for the FES convergence). Therefore, at the beginning of the well-tempered simulation, one can use a large hill height for the deep potential-well filling (i.e. 1 kcal/ mol versus standard version of 0.01 kcal/mol, that is 100 times larger than the standard version of metadynamics), and after sufficient simulation time, the hill height will go close to ∼0 kcal/mol, namely, convergence of the FES. As mentioned above, in the system of drug-target interaction, the slowest degree of freedom usually incorporates the drug’s binding/unbinding process and large conformational change of the target, such as the induced-fit phenomenon. Therefore, herein, the free energy surfaces were designed with two collective variables (CVs): (1) the center of mass (CoM) distance between the heavy atoms of crizotinib enantiomers and the active pocket of MTH1 (residues within 5 Å of the ligand in the holo-state MTH1) and (2) the loop-sheet distance between the binding loop (residues 23−30) and the β-sheet plane made from the Cα atoms in residues I70, F72, and M81, which represents the opening degree of the binding pocket and was also measured by the CoM distance. To attenuate the random influence of MD simulation, here, the crystal structure of (S)- and (R)-crizotinib-MTH1 complexes were employed as the initial structures for the metadynamics simulation. Before metadynamics simulation, we at first detected the largest pocket direction by using Caver 2.0,36 which was shown to be effective in drug-dissociation based enhanced sampling simulations in our previous work.37−39 Afterward, the complexes were rotated with the largest pocket direction toward the X-axis, and a rectangular water box was added to each drug-target complex with the water molecules 30 Å out of the solutes in the X-axis and 10 Å in the Y- and Z-axes. Third, the systems were optimized with all the heavy atoms in the drug-target complexes restrained in 100 kcal/mol·Å2. Finally, the solvent molecules (including water molecules, ions, and hydrogen atoms in the solutes) were equilibrated for 1 ns in an NPT ensemble (P = 1 atm and T = 310 K; the heavy atoms in the drug-target complexes were still restrained). The equilibrated structures were used for the following metadynamics simulations. In metadynamics simulations, to prevent drifting of the systems, only heavy atoms incorporated in the β-sheet which have been used for the pocket opening degree detection (residues I70, F72, and M81) were restrained with 5 kcal/mol· Å2.40 The initial hill height (ω), the deposition rate (τ), and the bias-factor (ΔT) were set to 1 kcal/mol, 1 kcal/mol·ps, and 3100, respectively. As mentioned above, two CVs were used for the FES construction: in the first CV, namely the drug-target distance, the width of the Gaussian hill (σ) was set to 0.4 Å (0.05 × 8 in the NAMD algorithm), and 500 bins were divided from 0 to 25 Å to detect the free energy change along the drugtarget distance; and in the second CV, known as the pocket opening degree, the same width of the Gaussian hill was set with the first CV (0.4 Å), while 400 bins ranging from 10 to 30 Å were used for the detection of free energy upon the pocket opening degree.

Due to the fact that the biased molecule may absorb on the reaction unrelated region of the target in the distance-restrained metadynamics simulation, here, we merged the funnelrestraining algorithm10 into the well-tempered metadynamics, making it funnel-based well-tempered metadynamics (FWTM).14 In the spirit of funnel based metadynamics simulation, it actually adds a harmonic restrained wall around the selected CV, which means that once the biased molecule goes out of the reaction associated sampling space (namely, the restrained funnel), a restrained energy will be added to force it back into the reaction associated sampling space. Herein, a cylinder-shaped restrained funnel was constructed along the first CV of the system (drug-target distance). The radius of the restrained funnel was set to 7.5 Å to provide enough sampling space for the rotation of the crizotinib enantiomers (∼10 Å in diameter). The elastic constant of the restrained funnel was set to 100 kcal/mol·Å2. Due to the fact that the use of cylindrical restraints may reduce the sampling space of the biased molecule in the bulk, it is necessary to correct the absolute binding free energy according to eq 2 ΔG bind‐Meta = ΔG Meta‐depth −

⎛ 1 1 ⎞ 2 ο ln(πR cyl C )⎜ β = ⎟ β kBT ⎠ ⎝ (2)

where ΔGbind‑Meta denotes the absolute binding free energy of the FWTM; ΔGMeta‑depth is the PMF depth of the crizotinib enantiomers going from bulk to the active site of MTH1; πR2cyl is the surface of the cylinder funnel; C° is the standard concentration (1/1661 Å3), and kB is Boltzmann constant. The detailed descriptions of the algorithm can be found in ref 10. Due to the constraint of the funnel, there is 1.38 kcal/mol analytical correction of the PMF depth (ΔGMeta‑depth), which corresponds to the 7.5 Å radius of the restrained funnel (at 310 K). To determine the binding stages of crizotinib enantiomers binding, the favorable binding pathway (black dotted lines in Figures 2D and 3D, and a 2-D description could be found in Figures 2F and 3F) was constructed along the X-axis (drugtarget distance) by connecting bins with the lowest ensemble energies along the Y-axis (pocket opening degree). The directly binding pathway was determined by averaging bins with the lowest ensemble energies from 15 to 23 Å along the X-axis (red dotted lines in Figures 2D and 3D, and a 2-D description was shown in Figures 2H and 3H), which represents the apo-state MTH1. While the induced-fit pathway (green dotted lines in Figure 2D and 3D, and a 2-D description was shown in Figures 2G and 3G) was determined by extending the stable binding position (the huge green dots in Figures 2D and 3D) along the Y-axis until intersecting with the directly binding pathway (the huge red dots in Figure 2D and 3D). 4.4. Umbrella Sampling Based Absolute Binding Free Energy Calculation. As a classic algorithm, umbrella sampling (US) should be one of the most widely used enhanced sampling methodologies in the scope of free energy calculation.41−44 Through adding biasing potentials to the given reaction coordinate (RC), the system can be driven from one thermodynamic state to another.45 For the issue of sufficient sampling, it is necessary to divide the RC into different pieces, namely windows, during the US simulations. In each window, a biasing potential (usually harmonic potential) is added in the middle of the RC for the purpose of biased 857


Article

Journal of Chemical Theory and Computation sampling. In the case of harmonic potential, the biasing potential can be added according to the following equation

uibias =

1 ki(r − riref )2 2

(3)

where ki is the biasing force constant in window i. r is the current position of a given reference state; and rref I represents the reference state of window i, which can denote the conformation (RMSD), angle (α, β, γ, θ, and Φ), and position (r) of the biased molecule. After sufficient sampling, the distribution of the biased samples can be unbiased with the reweighted methods, such as weighted histogram analysis method (WHAM).46,47 Finally, the normal-distributed samples in each window can be integrated into a potential of mean force curve (PMF curve) to obtain the total free energy change along the RC. Besides US technology, another newly developed enhanced sampling methodology, adaptive biasing force (ABF),48,49 is also widely used in the field of drug-target interactions due to the superiority of fewer simulation windows and priori parameters needed to be divided and tested.50−52 Different from US that uses sample probability to deduce the free energy change along the RC, ABF adds biasing force to the CV directly. With sufficient sampling, the biased molecule can go with no energy barriers along the RC (namely, free-diffusionlike behavior).37,50 Thereafter, the free energy change along the RC can be obtained from the equilibrated biasing force directly. In the scope of drug-target interactions, it may be hard for traditional US/ABF simulations (without using additional restraints) to obtain the absolute binding free energy upon drugs binding as the fact that too many degrees of freedom are usually incorporated in the reaction pathway. Fortunately, Woo and Roux have developed a well-characterized and very effective absolute binding free energy calculation strategy for the enhanced sampling simulations.11 Through restraining the conformation, the orientation, and the translation of the biased molecule, the sampling in the given RC can converge rapidly with the fact that the calculation strategy can constrain the relative external degrees of freedom. 53−55 As a sound methodology, the absolute binding free energy of Woo and Roux’s scheme is derived from the reaction equilibrium constant (Kbind) of the drug-target interaction. By calculating the restraints associated energies upon the conformation (Gbulk c bulk α β γ and Gsite c ), the orientation (Go , Gc , Gc , and Gc ), and the translation (Gθt , and GΦ t ) of the biased molecule, the absolute binding free energy can be recomposed according to the following equations 1 ΔG bind‐US = − ln(KbindC ο) β bulk

Kbind = S*I *e−β[Gc

+ Gobulk − Gtθ − GtΘ − Goα − Goβ − Goγ − Gcsite]

Figure 5. Rotational and translational restraints on crizotinib enantimers. The defined restrained points are shown in green sphere models (P1, P2, and PC) for MTH1 (yellow) and cyan/blue/white sphere models (L1, L2, and LC) for the crizotinib enantiomers (white), respectively. P1 and P2 points are Cα atoms of residues Y136 and R102, respectively, and PC point is a fictitious atom (modeled by overlapping with Cα atom in Y122) sited about 12 Å away from LC along the X-axis. L1, L2, and LC are the atoms of F1, N5, and C13, respectively, in the crizotinib enantiomers. The rotational restraints are defined as α(L1, LC, PC), β(L1, LC, PC, P2), and γ(L2, L1, LC, PC), and the translational restraints are designed with θ(LC, PC, P2), Φ(LC, PC, ⎯⎯⎯⎯→ P2, P1), and r (along the vector of PCLC , which is the dissociation direction of the ligand). The restrained angles are not labeled for the issue of clarity.

molecule) associated with free energy calculations. To give a comparison between the absolute binding free energy calculated by FWTM, we used the same structures of crizotinib-MTH1 complexes in metadynamics simulations as the initial conformations for the Woo and Roux scheme based absolute binding free energy calculations. In the stage of conformational restrained simulations, the conformational change of the crizotinib enantiomers (characterized by RMSD) was employed as the reaction coordinate for the US simulation. The elastic constant was set to 0.01 kcal/mol·Å2. Eleven and 7 windows were used for the apo-state and holostate crizotinib simulations, respectively. A 5 ns US simulation was performed for each window (with the last 3 ns samples used for the conformational PMF construction). The window size was set to 0.5 Å, making the RC range from 0 to 5 Å and 0 to 3 Å for the apo-state and the holo-state crizotinib. In addition to the angular restrained simulations, we orderly performed the ABF simulations for the five angles (α, β, γ, θ, and Φ) with the crizotinib enantiomers restrained in their initial states (the crizotinib enantiomers were restrained with 1 kcal/mol·Å2 at 0 Å of the RMSD and followed by the orderly angular restraints of 0.3 kcal/mol·Å2 in each angle). To improve the computational efficiency, only one window was used for the angular restrained ABF simulations.14 Here, fictitious atoms were employed to construct a generalized coordinate to overcome the shortcomings of the ABF algorithm in the NAMD code.56 A 5 ns ABF simulation was performed for each angle with the bin size of 0.2 Å. As for the separation degree of simulation (r), 41 windows were divided with each window 0.5 Å in length and ranged from 12 to 32 Å of the RC, where the RC was designed ⎯⎯⎯⎯→ along the vector of PCLC , and the PC point is a dummy atom (overlapped with the Cα atom in residue Y122, Figure 5)

(4) (5)

where ΔGbind‑US is the absolute binding free energy of drugtarget interaction, and S* and I* correspond to the translational angle restraints in bulk and the separation PMF based integral, respectively.11 To calculate the absolute binding free energies of the two crizotinib enantiomers to MTH1, here, we use US for the conformational (RMSD) and the separating (r) restraints associated with free energy calculations and ABF for the angular restraints (including restraints on rotational (α, β, γ) and translational (θ, Φ) degrees of freedom on the biased 858


Article


(4) Gad, H.; Koolmeister, T.; Jemth, A.-S.; Eshtad, S.; Jacques, S. A.; Ström, C. E.; Svensson, L. M.; Schultz, N.; Lundbäck, T.; Einarsdottir, B. O. MTH1 inhibition eradicates cancer by preventing sanitation of the dNTP pool. Nature 2014, 508, 215−221. (5) Huber, K. V.; Salah, E.; Radic, B.; Gridling, M.; Elkins, J. M.; Stukalov, A.; Jemth, A.-S.; Göktürk, C.; Sanjiv, K.; Strömberg, K. Stereospecific targeting of MTH1 by (S)-crizotinib as an anticancer strategy. Nature 2014, 508, 222−227. (6) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (7) Svensson, L. M.; Jemth, A.-S.; Desroses, M.; Loseva, O.; Helleday, T.; Högbom, M.; Stenmark, P. Crystal structure of human MTH1 and the 8-oxo-dGMP product complex. FEBS Lett. 2011, 585, 2617−2621. (8) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. (9) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (10) Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel metadynamics as accurate binding free-energy method. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 6358−6363. (11) Woo, H.-J.; Roux, B. Calculation of absolute protein−ligand binding free energy from computer simulations. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6825−6830. (12) Wilson, C.; Agafonov, R.; Hoemberger, M.; Kutter, S.; Zorba, A.; Halpin, J.; Buosi, V.; Otten, R.; Waterman, D.; Theobald, D. Using ancient protein kinases to unravel a modern cancer drug’s mechanism. Science 2015, 347, 882−886. (13) Zhang, Y.; Voth, G. A. Combined metadynamics and umbrella sampling method for the calculation of ion permeation free energy profiles. J. Chem. Theory Comput. 2011, 7, 2277−2283. (14) Sun, H.; Li, Y.; Tian, S.; Wang, J.; Hou, T. P-loop conformation governed crizotinib resistance in G2032R-mutated ROS1 tyrosine kinase: clues from free energy landscape. PLoS Comput. Biol. 2014, 10, e1003729. (15) Lau, A. Y.; Roux, B. The hidden energetics of ligand binding and activation in a glutamate receptor. Nat. Struct. Mol. Biol. 2011, 18, 283−287. (16) Søndergaard, C. R.; Olsson, M. H.; Rostkowski, M.; Jensen, J. H. Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa values. J. Chem. Theory Comput. 2011, 7, 2284−2295. (17) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K.; Burant, J. Gaussian 03; Gaussian. Inc.: Wallingford, CT, 2004. (18) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269−10280. (19) Case, D.; Darden, T.; Cheatham, T., III; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Walker, R.; Zhang, W.; Merz, K. AMBER 12; University of California: San Francisco, CA, 2012. (20) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (21) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (22) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (23) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341.

placed at approximately 12 Å away from the LC point (C13 atom in the crizotinib enantiomers). The elastic constant in each window was set to 5 kcal/mol·Å2. The initial structure of each window was derived from the corresponding confromation of the favorable pathway in the metadynamics trajectory.13 In the separation US simulation, each window was preformed under the conformational and angular restraints in crizotinib. A 12 ns simulation was preformed for each window, and a total of 492 ns simulation was used for each separation sampling of the crizotinib enantiomers (a detailed simulation protocol can be found in Table S2). Due to the isotropic property of the bulk, it could use direct numerical integation for the orientation associated energies calculation rather than actual MD simulations.11

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00973. Figure S1, H-bond interactions between the two crizotinib enantiomers and MTH1; Figure S2, systems stability across 40 ns conventional MD simulation; Figure S3, convergence of the separation PMFs for Metadynamics (panels A and B) and US (panels C and D) simulations; Figure S4, 1-D free energy profiles of (S)-crizotinib (blue) and (R)-crizotinib (orange); Figure S5, modeled conformation of residue F27 in the (R)crizotinib system; Table S1, free energy decomposition of Woo and Roux’s absolute binding free energy (kcal/ mol); and Table S2, simulation details in this study (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 86-571-88208412. E-mail: [email protected] or [email protected]. Author Contributions

H.S. and P.C. contributed equally to this work. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This study was supported by the National Science Foundation of China (21575128; 81302679) and the National Science Foundation for Postdoctoral Scientists of China (2015M581953).

■

REFERENCES

(1) Cui, J. J.; Tran-Dubé, M.; Shen, H.; Nambu, M.; Kung, P. P.; Pairish, M.; Jia, L.; Meng, J.; Funk, L.; Botrous, I. Structure Based Drug Design of Crizotinib (PF-02341066), a Potent and Selective Dual Inhibitor of Mesenchymal−Epithelial Transition Factor (c-MET) Kinase and Anaplastic Lymphoma Kinase (ALK). J. Med. Chem. 2011, 54, 6342−6363. (2) Kwak, E. L.; Bang, Y. J.; Camidge, D. R.; Shaw, A. T.; Solomon, B.; Maki, R. G.; Ou, S. H. I.; Dezube, B. J.; Jänne, P. A.; Costa, D. B. Anaplastic lymphoma kinase inhibition in non-small-cell lung cancer. N. Engl. J. Med. 2010, 363, 1693−1703. (3) Davies, K. D.; Le, A. T.; Theodoro, M. F.; Skokan, M. C.; Aisner, D. L.; Berge, E. M.; Terracciano, L. M.; Cappuzzo, F.; Incarbone, M.; Roncalli, M. Identifying and targeting ROS1 gene fusions in non− small cell lung cancer. Clin. Cancer Res. 2012, 18, 4570−4579. 859


Article


(45) Kästner, J. Umbrella sampling. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 932−942. (46) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (47) Souaille, M.; Roux, B. t. Extension to the weighted histogram analysis method: combining umbrella sampling with free energy calculations. Comput. Phys. Commun. 2001, 135, 40−57. (48) Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 144120−144132. (49) Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. J. Chem. Theory Comput. 2010, 6, 35−47. (50) Dehez, F.; Pebay-Peyroula, E.; Chipot, C. Binding of ADP in the mitochondrial ADP/ATP carrier is driven by an electrostatic funnel. J. Am. Chem. Soc. 2008, 130, 12725−12733. (51) Neumann, A.; Baginski, M.; Czub, J. How do sterols determine the antifungal activity of amphotericin B? Free energy of binding between the drug and its membrane targets. J. Am. Chem. Soc. 2010, 132, 18266−18272. (52) Wereszczynski, J.; McCammon, J. A. Nucleotide-dependent mechanism of Get3 as elucidated from free energy calculations. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 7759−7764. (53) Faraldo-Gómez, J. D.; Roux, B. On the importance of a funneled energy landscape for the assembly and regulation of multidomain Src tyrosine kinases. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 13643− 13648. (54) Lin, Y.-L.; Meng, Y.; Jiang, W.; Roux, B. Explaining why Gleevec is a specific and potent inhibitor of Abl kinase. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1664−1669. (55) Lin, Y.-L.; Roux, B. Computational Analysis of the Binding Specificity of Gleevec to Abl, c-Kit, Lck, and c-Src Tyrosine Kinases. J. Am. Chem. Soc. 2013, 135, 14741−14753. (56) Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013, 9, 794−802.

(24) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constantpressure molecular-dynamics simulation-the Langevin piston method. J. Chem. Phys. 1995, 103, 4613−4621. (25) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (26) Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601. (27) Tribello, G. A.; Ceriotti, M.; Parrinello, M. A self-learning algorithm for biased molecular dynamics. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 17509−17514. (28) Piana, S.; Laio, A. A bias-exchange approach to protein folding. J. Phys. Chem. B 2007, 111, 4553−4559. (29) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-energy landscape for β hairpin folding from combined parallel tempering and metadynamics. J. Am. Chem. Soc. 2006, 128, 13435−13441. (30) Barducci, A.; Bonomi, M.; Prakash, M. K.; Parrinello, M. Freeenergy landscape of protein oligomerization from atomistic simulations. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, E4708−E4713. (31) Corbi-Verge, C.; Marinelli, F.; Zafra-Ruano, A.; Ruiz-Sanz, J.; Luque, I.; Faraldo-Gómez, J. D. Two-state dynamics of the SH3−SH2 tandem of Abl kinase and the allosteric role of the N-cap. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, E3372−E3380. (32) Li, L.; Martinis, S. A.; Luthey-Schulten, Z. Capture and Quality Control Mechanisms for Adenosine-5′-triphosphate Binding. J. Am. Chem. Soc. 2013, 135, 6047−6055. (33) Limongelli, V.; Marinelli, L.; Cosconati, S.; La Motta, C.; Sartini, S.; Mugnaini, L.; Da Settimo, F.; Novellino, E.; Parrinello, M. Sampling protein motion and solvent effect during ligand binding. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1467−1472. (34) Berteotti, A.; Barducci, A.; Parrinello, M. Effect of urea on the βhairpin conformational ensemble and protein denaturation mechanism. J. Am. Chem. Soc. 2011, 133, 17200−17206. (35) Palazzesi, F.; Barducci, A.; Tollinger, M.; Parrinello, M. The allosteric communication pathways in KIX domain of CBP. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 14237−14242. (36) Petřek, M.; Otyepka, M.; Banás,̌ P.; Košinová, P.; Koča, J.; Damborský, J. CAVER: a new tool to explore routes from protein clefts, pockets and cavities. BMC Bioinf. 2006, 7, 316. (37) Sun, H.; Li, Y.; Li, D.; Hou, T. Insight into Crizotinib Resistance Mechanisms Caused by Three Mutations in ALK Tyrosine Kinase using Free Energy Calculation Approaches. J. Chem. Inf. Model. 2013, 53, 2376−2389. (38) Pan, P.; Li, Y.; Yu, H.; Sun, H.; Hou, T. Molecular Principle of Topotecan Resistance by Topoisomerase I Mutations through Molecular Modeling Approaches. J. Chem. Inf. Model. 2013, 53, 997−1006. (39) Sun, H.; Tian, S.; Zhou, S.; Li, Y.; Li, D.; Xu, L.; Shen, M.; Pan, P.; Hou, T. Revealing the favorable dissociation pathway of type II kinase inhibitors via enhanced sampling simulations and two-end-state calculations. Sci. Rep. 2015, 5, 8457. (40) Zeller, F.; Zacharias, M. Adaptive biasing combined with Hamiltonian replica exchange to improve umbrella sampling free energy simulations. J. Chem. Theory Comput. 2014, 10, 703−710. (41) Bock, L. V.; Blau, C.; Schröder, G. F.; Davydov, I. I.; Fischer, N.; Stark, H.; Rodnina, M. V.; Vaiana, A. C.; Grubmüller, H. Energy barriers and driving forces in tRNA translocation through the ribosome. Nat. Struct. Mol. Biol. 2013, 20, 1390−1396. (42) Fowler, P. W.; Sansom, M. S. The pore of voltage-gated potassium ion channels is strained when closed. Nat. Commun. 2013, 4, 1872. (43) Ostmeyer, J.; Chakrapani, S.; Pan, A. C.; Perozo, E.; Roux, B. Recovery from slow inactivation in K+ channels is controlled by water molecules. Nature 2013, 501, 121−124. (44) Setny, P.; Baron, R.; Kekenes-Huskey, P. M.; McCammon, J. A.; Dzubiella, J. Solvent fluctuations in hydrophobic cavity−ligand binding kinetics. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1197−1202. 860


Directly Binding Rather than Induced-Fit ... - ACS Publications

Recommend Documents