Steered Molecular Dynamics Simulations Reveal the Likelier Dissociation Pathway of Imatinib from Its Targeting Kinases c-Kit and Abl

Abstract

Development of small molecular kinase inhibitors has recently been the central focus in drug discovery. And type II kinase inhibitors that target inactive conformation of kinases have attracted particular attention since their potency and selectivity are thought to be easier to achieve compared with their counterpart type I inhibitors that target active conformation of kinases. Although mechanisms underlying the interactions between type II inhibitors and their targeting kinases have been widely studied, there are still some challenging problems, for example, how type II inhibitors associate with or dissociate from their targeting kinases. In this investigation, steered molecular dynamics simulations have been carried out to explore the possible dissociation pathways of typical type II inhibitor imatinib from its targeting protein kinases c-Kit and Abl. The simulation results indicate that the most favorable pathway for imatinib dissociation corresponds to the ATP-channel rather than the relatively wider allosteric-pocket-channel, which is mainly due to the different van der Waals interaction that the ligand suffers during dissociation. Nevertheless, the direct reason comes from the fact that the residues composing the ATP-channel are more flexible than that forming the allosteric-pocket-channel. The present investigation suggests that a bulky hydrophobic head is unfavorable, but a large polar tail is allowed for a potent type II inhibitor. The information obtained here can be used to direct the discovery of type II kinase inhibitors.

About the Authors:

Li-Jun Yang

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Jun Zou

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Huan-Zhang Xie

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Lin-Li Li

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Yu-Quan Wei

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Sheng-Yong Yang

* E-mail: [email protected]

Affiliation: State Key Laboratory of Biotherapy and Cancer Center, West China Medical School, Sichuan University, West China Hospital, Chengdu, People's Republic of China

Introduction

Protein kinases are enzymes essential for cell signal transduction, which regulate a variety of physiological processes including metabolic, cell cycle, apoptosis and cell differentiation [1]–[3]. Dysregulation of protein kinases might lead to some pathological changes, for example, cancer, diabetes, and various autoimmune diseases [4], [5]. Thus protein kinases have been thought as central targets for drug discovery. In the past decade, extensive efforts have been made to develop protein kinase inhibitors as potential drugs against a wide range of diseases [6]–[13]. And it is believed that understanding of issues related to the protein kinase structures, mechanisms underlying enzyme activation and the kinase-inhibitor interaction could benefit the discovery of novel kinase inhibitors.

All protein kinases share a common catalytic domain, which contains two subdomains: the N-terminal lobe and the C-terminal lobe [14]. The two lobes are connected through a flexible chain (hinge region). The natural substrate ATP is bound in the cleft between the two lobes (the ATP binding pocket). The active loop (A-loop), which belongs to the C-terminal lobe but locates outside of the ATP-binding pocket, directly regulates the enzyme activation through its conformational changes. Majority of small molecule kinase inhibitors reversibly occupy the ATP binding pocket, which means that they are ATP-competitive inhibitors. The ATP-competitive inhibitors can be further classified into two categories, type I and type II [15], [16]. Type I inhibitors target the active form of the kinases, in which the A-loop adopts an extended conformation. Such conformational arrangement of A-loop completely exposes the ATP-binding pocket, hence facilitating the entry/exit of ATP or type I inhibitors (this entry/exit pathway will be called as traditional ATP-channel hereafter, see Figure 1A). Type II inhibitors target the inactive form of kinases and bind to an extended ATP-binding site, in contrast to type I inhibitors. In the inactive form, the A-loop crimples outside of the ATP-binding pocket. This conformation of A-loop shrinks the original entry/exit gate, which hinders the access of ATP and protein substrates to the kinase catalytic site. Another concomitant conformational change is the flip of DFG-motif that locates in the beginning of A-loop, which opens a new hydrophobic pocket (usually called allosteric pocket) in the back of the protein [17] (see Figure 1B). Type II inhibitors often occupy both the original ATP-binding pocket and the allosteric pocket. It appears that there are two possible pathways for the entry/exit of type II inhibitors: one is the traditional ATP-channel and the other one is the allosteric-pocket-channel. Now, a question arises that which one is preferred. X-ray crystal structures of kinase-inhibitor complexes show that the allosteric-pocket-channel might be favored since this channel is relatively wider than the ATP-channel [18], [19]. This hypothesis, however, is inconsistent with the fact that many receptor tyrosine kinases have a juxtamembrane region (JMR), which resides close to the gate of allosteric-pocket-channel in the inactive form of kinases. Even so, our previous study on the JMR dynamics did not deny the allosteric-pocket-channel of type II inhibitors [20]. In order to clarify this mechanism, we need to use molecular dynamics simulations, which is mainly due to the fact that the dissociation of ligands from targeting proteins is governed by the dynamic behavior of ligand-protein complexes that is difficult to handle experimentally.

[Figure omitted. See PDF.]

Figure 1. Typical three-dimensional structures of protein kinases shown in Cα ribbon fashion.

(A) is for active conformation, and (B) for inactive conformation. Key structural components of the protein are color coded: A-loop in red, helix αC in purple, others in gray. Type I (for the active conformation) and type II (for the inactive conformation) kinase inhibitors are schematically shown in green wire mesh.

https://doi.org/10.1371/journal.pone.0008470.g001

In this account, steered molecular dynamics (SMD) simulations [21]–[23] will be employed to explore the possible dissociation channel for type II kinase inhibitors from inactive form of kinases. Despite several methods of molecular dynamics simulations that are currently available for handling this type of rare events (see the section of Models and Methods), SMD has been seen as the most suitable method to address the binding/unbinding issues of biomolecules [24]–[30]. (A brief description about the SMD algorithm, as well as its advantages and disadvantages compared with other approaches, will be given in the section of Models and Methods). Here, we chose imatinib (see Figure 2A) as our model of type II kinase inhibitors; imatinib is the first type II protein kinase inhibitor on the market and has become a reference compound for the kinase community. The cellular targets of imatinib include the Abl/Bcr-Abl protein kinase that is involved in the development of chronic myeloid leukaemia (CML) [31], the c-Kit kinase that is implicated in gastrointestinal stromal tumours (GIST) [32], [33]. Thanks to the advance of structural biology, the crystal structures of c-Kit and Abl in complex with imatinib have been available [18], [19] (see Figure 2B and 2C). Another important reason why we chose both c-Kit and Abl protein kinases as our models is that c-Kit and Abl represent two distinct regulation modes of catalytic activation. Abl is a typical cytoplastic kinase, which activation depends directly on the phosphoralation of TYR412 on A-loop [34]. Contrastively, c-Kit is a typical receptor tyrosine kinase, which activation depends on not only the phosphoralation of the A-loop but also the JMR (TYR1-ASP33) [20], [35]. Since JMR locates around the gate of allosteric-pocket-channel, it seems that the JMR could have some impact on the entry/exit of ligand. In a word, we expected that SMD simulations could provide us with more insights into the entry/exit mechanism of type II inhibitors from the inactive form of protein kinases. We also believe that a detailed understanding of the entry/exit dynamics will ultimately help with the design of new type II kinase inhibitors.

[Figure omitted. See PDF.]

Figure 2. Chemical structure of imatinib and its binding modes to c-Kit and Abl.

(A) Chemical structure of imatinib. (B) 3D structure of imatinib-c-Kit complex (PDB entry: 1T46). (C) 3D structure of imatinib-Abl complex (PDB entry: 2HYY). In (B) and (C), protein structures are in Cα ribbon fashion and imatinib in green CPK fashion. Key structural components of the protein are color coded: A-loop in red, helix αC in purple, JMR in blue, others in gray.

https://doi.org/10.1371/journal.pone.0008470.g002

Methods

Computational Models

The initial structures for imatinib-c-Kit and imatinib-Abl complexes were prepared based on their corresponding X-ray crystal structures (PDB entry is 1T46 for imatinib-c-Kit and 2HYY for imatinib-Abl). The unresolved residues in imatinib-c-Kit, including TYR547-ASN564, ILE690-THR694, SER753-GLU761, HIS934-ILE935, were constructed by using the homology modeling module Modeler within InsightII (Accelrys, San Diego, CA), in which the reference protein structure was taken from the PDB code 1T45.

Steered Molecular Dynamics Simulations

SMD is an extended MD simulation method that mimics the basic idea of atomic force microscopy (AFM). SMD allows exploring the binding and unbinding properties of a variety of biomolecules as well as their responses to external mechanical manipulations at the atomic level [36], [37]. It is particularly suitable for the investigation of dissociation of a ligand from its binding protein. In such case, a ligand is forcedly pulled out of the binding site of a given protein along the route that is specified in advance. The pulling force F is applied according to the following equation:(1)where is the force constant, is the pulling velocity, is the pulling direction normal, and are the position of ligand at time t and initial time. Briefly, if ligand moves forward along the leaving pathway, the increases and the force decreases; if the pulling force can not displace the ligand, the force will increase because of the increase of t. So, it is possible to draw a profile to express the force required to dissociate the ligand as the function of simulation time, and hence determine the difficulty to dissociate ligand along the different pathway.

Actually several methods that utilize the principles of molecular dynamics to investigate this type of rare processes have been developed. In these methods, some external perturbations are applied to the system in order to accelerate the dynamic process and guide the system toward a target state. Such methods typically include, for examples, targeted molecular dynamics (TMD) [38], [39] and biased molecular dynamics (BMD) [40], [41], in addition to SMD used here. The TMD methodology imposes a time-dependent holonomic or harmonic constraint on the RMSD (root mean square deviation) in order to drive the system to a target structure. BMD provides the least perturbation compared with SMD and TMD, in which the system feels no external force if it moves toward the target and the biasing potential is nonzero only if the system moves away from the target.

In contrast to TMD and BMD, SMD closely resembles the AFM experiment. In addition, the exact target structure is not necessary to know in advance for SMD simulations. These make the SMD very suitable to be used for exploring the dynamics behaviors of binding/unbinding processes of biomolecules. Although SMD provides an approach to investigate the binding/unbinding events, yet the pulling direction of the force in SMD is chosen randomly or by guessing on the basis of structural information, which makes that the force applied to the ligand in such chosen directions may not move it along a favorable pathway. In order to overcome this shortcoming, a generally used strategy is to perform a series of SMD simulations with different directions of pulling force. Dynamic properties of the system can be obtained by averaging those of all the individual simulations. This strategy has also been used in this investigation, which will be mentioned again later.

Computational Details

All the simulations, including conventional molecular dynamics simulations (CMD) which purpose was to make the model systems to reach their thermodynamic equilibration states, and SMD, were carried out using the programs NAMD 2.6 [42] with constant molecular number, pressure, and temperature (NPT) and periodic boundary conditions. The CHARMM 27 force field [43], [44] was employed for the protein and the gaff force field [45] for imatinib. The model systems were solvated in a bath of TIP3P water molecules [46] with the buffering distance set to 15 Å which is large enough for the conformational changes. By assuming normal charge states of ionizable groups corresponding to pH 7, sodium (Na+) and chloride (Cl−) counter-ions at physiological concentration of 0.15 mol·L−1 were added in the box in random positions to ensure the global charge neutrality. All Na+ and Cl− ions were placed more than 8 Å away from any protein atoms and from each other.

Constant temperature (T = 310 K) and constant pressure (P = 1 bar) were maintained using Langevin piston coupling algorithm [47]. The integration time step of the simulations was set to 2.0 fs, the SHAKE algorithm was used to constrain the lengths of all chemical bonds involving hydrogen atoms at their equilibrium values [48], and the water geometry was restrained rigid by using the SETTLE algorithm [49]. Non-bonded van der Waals interactions were treated by using a switching function at 10 Å and reaching zero at a distance of 12 Å. Long-range electrostatic forces were handled by using the particle-mesh Ewald (PME) algorithm [50], which is an efficient full electrostatics method for use with periodic boundary conditions.

The systems were firstly minimized using steepest-descents algorithm. Afterward, the systems were heated from 0 K to 310 K in two steps: the first step was a 100ps heating process with the protein and ligand fixed; the second step was another 100 ps heating process with the whole system relaxed. Thirdly, a 200 ps CMD simulation was performed with the protein and ligand fixed, followed by a 4 ns CMD simulation with the whole system relaxed. The SMD simulations were then performed, in which initial structures were the snapshots randomly taken from the equilibration states in the CMD simulations. The acting point of pulling point for the ATP-channel is the mass center of pyridine group of imatinib, and that for the allosteric-pocket-channel is the mass center of piperazinyl group of imatinib. The harmonic constraints were applied to the Cα atoms of helix in the Abl and c-Kit and the exponent to be used in the harmonic constraint energy functions was 2. The force constant and velocity used in the SMD simulations are 4 kcal mol−1 Å−2 (277.9 pN Å−1) and 0.00005 Å timestep−1 (0.025Å ps−1), respectively.

By the way, it is necessary to point out that the affinity of the ligand for the two pockets can not be directly given by the force measurements in the current SMD simulations, although SMD provides a possibility to calculate the binding affinity by, for example, applying the Jarzynski equality [36]. This is mainly due to that (1) the system concerned here experiences a large conformational change, and (2) the SMD simulations here are still very limited. To obtain the ligand binding affinity, a much larger set of simulations are required. Fortunately, the calculation of ligand binding affinity is not necessory for the objective in this investigation that is to explore the possible dissociation pathways of typical type II inhibitor imatinib from its targeting protein kinases c-Kit and Abl.

Results and Discussion

Dissociation of Imatinib from c-Kit

Before performing SMD simulations to pull imatinib out of c-Kit, 4 ns CMD simulation was first carried out which purpose was to make the whole system equilibrated at the designated temperature. Figure 3 shows the variations of RMSD values for the Cα atoms of JMR (TYR1-ASP33), A-loop (CYS205-LEU227) and the remaining regions (RR, that are all the residues of c-Kit except those in JMR and in A-loop), relative to the crystal structure during the MD simulations. Obviously, in the beginning, the RMSD values for the JMR have a significant change, which is understandable since JMR is a loop and completely exposed to the solvent. In contrast to JMR, a small fluctuation of RMSD values is observed for A-loop and RR. After 3 ns, the RMSD values for the whole system fluctuate steadily, implying that the whole system reaches its equilibration state.

[Figure omitted. See PDF.]

Figure 3. RMSD (root mean square deviation) profiles of key motifs in imatinib-c-Kit complex during conventional molecular dynamics simulation.

The reference structure from which the RMSD was calculated is the crystal structure 1T46.

https://doi.org/10.1371/journal.pone.0008470.g003

After that, ten snapshots were randomly chosen from the equilibration state (3ns–4ns) as the initial structures for the subsequent SMD simulations, with five snapshots for each pathway (ATP-channel or allosteric-pocket-channel). Before executing the actual SMD simulations, all the possible directions for the pulling force were identified based on the relative orientation of ligand and its binding site. Then five independent SMD simulations for each pathway were performed, in which the direction of pulling force was randomly chosen in the range of possible directions just identified.. The simulation time for each SMD simulation is 1 ns.

Figure 4A presents the profiles of pulling force during the SMD simulations for the ATP-channel. Clearly, increased force appears when the ligand imatinib begins to move out of the binding site, irrespective of initial structure and the direction of pulling force, which imply that the ligand encounters energy barriers. The maximum pulling forces for the ATP-channel are between 737–1080 pN, and averaged at 816 pN (see Figure 4C). Finally the pulling force becomes to zero, indicating that the ligand completely dissociated from c-Kit. The pulling force profiles for dissociation along the allosteric-pocket-channel are shown in Figure 4B. Very similar to the ATP-channel, the forces increase in the beginning of ligand unbinding, and reaches a maximum around 0.4 ns, followed by returning to zero.

[Figure omitted. See PDF.]

Figure 4. Force profiles of imatinib dissociation from imatinib-c-Kit complex.

(A) Force profile of imatinib unbinding along traditional ATP-channel. (B) Force profile of imatinib unbinding along allosteric-pocket-channel. (C) Averaged force profiles of imatinib unbinding along ATP-channel (red) and allosteric-pocket-channel (black).

https://doi.org/10.1371/journal.pone.0008470.g004

Visibly, the ligand suffers a smaller force when dissociation along the ATP-channel than through the allosteric-pocket-channel. Quantitatively, the maximum average force along the ATP-channel is 816 pN, which is smaller by 19% than that along the allosteric-pocket-channel (1008 pN). The integrated forces shown in Table 1 reflect even clearer how dissociation is easier along ATP-channel: the total force required for the ligand to be extracted is 332 pN·ns, which is 16% smaller than that in allosteric-pocket-channel (397 pN·ns).

[Figure omitted. See PDF.]

Table 1. Properties (maximum force and integrated force) of the average forces suffered by imatinib during dissociation processes along the ATP-channel and allosteric-pocket-channel.

https://doi.org/10.1371/journal.pone.0008470.t001

What are the determining factors to result in that the ligand dissociation along the allosteric-pocket-channel is difficult than the ATP-channel? Here one may simply think that the JMR, which locates outside of the allosteric-pocket-channel, might be the most important factor. Thus, in order to further explore this matter, another model system, namely, imatinib-Abl, in which Abl is a cytoplastic kinase, was chosen to perform the same SMD simulation. Since no JMR in the Abl kinase, imatinib-Abl is expected to be a good model to judge the role of JMR.

Dissociation of Imatinib from Abl

Again a 4 ns CMD simulation was firstly performed in order to make the imatinib-Abl equilibrated in the solvent at the designated temperature. Figure 5 depicts the variation of the RMSD values relative to the initial structure for the Cα atoms of A-loop (ALA146-PRO168) and the remaining regions (RR, that are all the residues of Abl except those in A-loop). Clearly, after 2.5 ns, the RMSD values for A-loop, RR, and hence the overall system, fluctuate steadily, implying that the system reaches its equilibration state.

[Figure omitted. See PDF.]

Figure 5. RMSD (root mean square deviation) profiles of key motifs in imatinib-Abl complex during conventional molecular dynamics simulation.

The reference structure from which the RMSD was calculated is the crystal structure 2HYY.

https://doi.org/10.1371/journal.pone.0008470.g005

Similarly, we selected randomly ten snapshots in the equilibration state as the initial structures for the SMD simulations. Five of them are for the ATP-channel and the other five for the allosteric-pocket-channel. Figure 6A and 6B show the pulling force profiles for dissociation along the ATP-channel and allosteric-pocket-channel, respectively. For both pathways, the pulling forces increase in the beginning of ligand unbinding process, irrespective of initial structure and the direction of pulling force, indicating again that the ligand encounters energy barriers. Finally the pulling forces return to zero, denoting that the ligand completely dissociated from the Abl.

[Figure omitted. See PDF.]

Figure 6. Force profiles of imatinib dissociation from imatinib-Abl complex.

(A) Force profile of imatinib unbinding along ATP-channel. (B) Force profile of imatinib unbinding along allosteric-pocket-channel. (C) Averaged force profiles of imatinib unbinding along ATP-channel (red) and allosteric-pocket-channel (black).

https://doi.org/10.1371/journal.pone.0008470.g006

Comparisons of the force profiles for both pathways show that the maximum pulling forces for the ATP-channel are obviously smaller than that for the allosteric-pocket-channel. The average forces together with the integrated forces shown in Figure 6C and Table 1 further clearly demonstrate that dissociation of imatinib is easier along ATP-channel than allosteric-pocket-channel: the total force required for the ligand to be extracted is 324 pN·ns, which is 16% smaller than that in allosteric-pocket-channel (384 pN·ns).

Intrinsic Reasons Why the Ligand Dissociation through the ATP-Channel Is Favored

The above simulations clearly indicate that the ligand imatinib dissociates from its binding protein kinases readily through the ATP-channel rather than the allosteric-pocket-channel, irrespective of the receptor, initial structure, and the direction of pulling force. Further, the simulation results also demonstrate that the JMR that resides outside of the allosteric-pocket-channel is not a factor that results in that the allosteric-pocket-channel is not preferred. What are the determining factors leading to that the ligand dissociation through the ATP-channel is favored? In order to address this question, we analyzed the key components of interaction energies between the ligand and protein kinases during SMD simulations, including electrostatic, van de Waals (vdW), and hydrogen-bonding interactions. It was found that the profiles of electrostatic and hydrogen-bonding interaction energy were very similar for the both pathways (data not shown). These are understandable due to the following reasons: (1) the ligand imatinib is a very weak polar molecule, which results in that the ligand suffers very similar electrostatic interaction when dissociating through different pathways; (2) hydrogen bonds are always in the courses of breaking and forming in the ligand unbinding process, and the overall number of hydrogen bonds keep relatively invariable.

For the vdW interaction energy, however, obvious differences could be found for these profiles along the ATP-channel and allosteric-pocket-channel. Figure 7A presents the profiles of average vdW interaction energy between imatinib and c-Kit during dissociation through the two pathways. Obviously, there is some fluctuation for both of the vdW energy profiles during the SMD simulations. The highest vdW energy peak (the energy barrier), however, is on the curve corresponding to the allosteric-pocket-channel. The time point for the highest vdW energy peak is around 0.4 ns, which is consistent with the time when the ligand suffers the maximum pulling force (see Figure 4C). The vdW interaction energy profiles for the imatinib-Abl system during dissociation through ATP-channel and allosteric-pocket-channel are shown in Figure 7B. Again, a fluctuation is observed for both of the vdW energy curves during the SMD simulations. The highest vdW energy peak (the energy barrier) occurs on the profile corresponding to the allosteric-pocket-channel. The time point for the highest vdW energy peak is around 0.3 ns, which is also consistent with the time when the ligand suffers the maximum pulling force (see Figure 6C).

[Figure omitted. See PDF.]

Figure 7. Energy profiles of van der Waals (vdW) interaction between imatinib and its targeting kinases during dissociation.

(A) vdW interaction of imatinib-c-Kit system. (B) vdW interaction of imatinib-Abl system. The vdW profiles shown in solid line are for the ATP-channel and those shown in dotted line are for the allosteric-pocket-channel.

https://doi.org/10.1371/journal.pone.0008470.g007

The preceding energy analysis indicates that the difference of vdW interaction energy is possibly a key factor determining that the ligand dissociation through the ATP-channel is favored. But what are the intrinsic reasons leading to the difference of vdW interaction energies? In order to answer this question, we further analyzed the flexibility of residues forming the ATP-channel and allosteric-pocket-channel. Figure 8A and 8B present the variations of RMSD values for the residues composing the ATP-channel and allosteric-pocket-channel during the SMD simulations. Very clearly, for both of the imatinib-c-Kit and imatinib-Abl model systems, the variations of RMSD values for ATP-channel residues are much larger than that for allosteric-pocket-channel, implying that the residues forming the ATP-channel are more flexible than that forming the allosteric-pocket-channel. These are consistent with B-factor values experimentally measured [see crystal structures of imatinib-c-Kit complex (PDB entry: 1T46) and imatinib-Abl complex (PDB entry: 2HYY)]: the values of B-factor for the residues forming the ATP-channel are generally larger than that for the residues forming the allosteric-pocket-channel. Further, the three-dimensional structures of ATP-channel and allosteric-pocket-channel provide additional evidences: the residues forming the ATP-channel belong to certain β-sheet or loop, nevertheless, in addition to β-sheets and loops, part of residues composing the allosteric-pocket-channel also belong to the α-helix (helix αC, see Figure 1). In general, flexibility of residues in α-helix is smaller than that in β-sheet and loop.

[Figure omitted. See PDF.]

Figure 8. Variations of RMSD (root mean square deviation) values for the residues composing the ATP-channel (solid line) and allosteric-pocket-channel (dotted line) during the steered molecular dynamics simulations.

(A) is for imatinib-c-Kit system, and (B) for imaintib-Abl system.

https://doi.org/10.1371/journal.pone.0008470.g008

Implications to the Development of Type II Kinase Inhibitors

Finally, what implications to the development of type II protein kinase inhibitors could be obtained from the present investigation? As discussed above, the likelier dissociation pathway for the ligand from its targeting kinases corresponds to the ATP-channel. Rationally we can deduce that the ATP-channel would be also preferred for the ligand association. This implies that, in the association process of a type II kinase inhibitors, the hydrophobic head of the inhibitor that is the motif locating in the allosteric pocket, has to first pass through the ATP-channel and then the allosteric-pocket-channel. It is clear that a bulky hydrophobic head in type II kinase inhibitors is unfavorable. On the contrary, a bulky polar tail in type II kinase inhibitors is allowed. These agree very well with chemical structures of most of the known type II kinase inhibitors that they often have a smaller hydrophobic head. Several typical examples are shown in Figure 9.

[Figure omitted. See PDF.]

Figure 9. Chemical structures of selected type II kinase inhibitors.

The substituent in dotted rectangle corresponds to the head of inhibitor and that in solid rectangle corresponds to the tail of inhibitor.

https://doi.org/10.1371/journal.pone.0008470.g009

In conclusion, molecular dynamics simulations reveal that the most favorable pathway for imatinib dissociation corresponds to the ATP-channel rather than the relatively wider allosteric-pocket-channel, which is mainly due to the different van der Waals interaction that the ligand suffers during dissociation. Nevertheless, the direct reason might come from the fact that the residues composing the ATP-channel are more flexible than those forming the allosteric-pocket-channel. This investigation, in addition to clarify the dissociation mechanism of type II kinase inhibitor imatinib from its targeting protein kinases, also suggests us that a bulky hydrophobic head is unfavorable, but a large polar tail is allowed for a potent type II inhibitor. The information obtained here can be used to direct the discovery of type II kinase inhibitors.

Author Contributions

Conceived and designed the experiments: L-JY L-LL Y-QW S-YY. Performed the experiments: L-JY JZ H-ZX. Analyzed the data: L-JY JZ H-ZX L-LL S-YY. Wrote the paper: L-JY S-YY.

Citation: Yang L-J, Zou J, Xie H-Z, Li L-L, Wei Y-Q, Yang S-Y (2009) Steered Molecular Dynamics Simulations Reveal the Likelier Dissociation Pathway of Imatinib from Its Targeting Kinases c-Kit and Abl. PLoS ONE4(12): e8470. https://doi.org/10.1371/journal.pone.0008470

References

1. Manning G, Whyte DB, Martinez R, Hunter T, Sudarsanam S (2002) The protein kinase complement of the human genome. Science 298: 1912–1934.G. ManningDB WhyteR. MartinezT. HunterS. Sudarsanam2002The protein kinase complement of the human genome.Science29819121934

2. Cohen P (2002) Protein kinases - the major drug targets of the twenty-first century? Nat Rev Drug Discov 1: 309–315.P. Cohen2002Protein kinases - the major drug targets of the twenty-first century?Nat Rev Drug Discov1309315

3. Dancey J, Sausville EA (2003) Issues and progress with protein kinase inhibitors for cancer treatment. Nat Rev Drug Discov 2: 296–313.J. DanceyEA Sausville2003Issues and progress with protein kinase inhibitors for cancer treatment.Nat Rev Drug Discov2296313

4. Noble MEM, Endicott JA, Johnson LN (2004) Protein kinase inhibitors: Insights into drug design from structure. Science 303: 1800–1805.MEM NobleJA EndicottLN Johnson2004Protein kinase inhibitors: Insights into drug design from structure.Science30318001805

5. Thaimattam R, Banerjee R, Miglani R, Iqbal J (2007) Protein kinase inhibitors: Structural insights into selectivity. Curr Pharm Des 13: 2751–2765.R. ThaimattamR. BanerjeeR. MiglaniJ. Iqbal2007Protein kinase inhibitors: Structural insights into selectivity.Curr Pharm Des1327512765

6. Sawyers C (2004) Targeted cancer therapy. Nature 432: 294–297.C. Sawyers2004Targeted cancer therapy.Nature432294297

7. Dowell J, Minna JD, Kirkpatrick P (2005) Erlotinib hydrochloride. Nat Rev Drug Discov 4: 13–14.J. DowellJD MinnaP. Kirkpatrick2005Erlotinib hydrochloride.Nat Rev Drug Discov41314

8. Arora A, Scholar EM (2005) Role of tyrosine kinase inhibitors in cancer therapy. J Pharmacol Exp Ther 315: 971–979.A. AroraEM Scholar2005Role of tyrosine kinase inhibitors in cancer therapy.J Pharmacol Exp Ther315971979

9. Biscardi JS, Tice DA, Parsons SJ (1999) c-Src, receptor tyrosine kinases, and human cancer. Adv Cancer Res 76: 61–199.JS BiscardiDA TiceSJ Parsons1999c-Src, receptor tyrosine kinases, and human cancer.Adv Cancer Res7661199

10. Wan PTC, Garnett MJ, Roe SM, Lee S, Niculescu-Duvaz D, et al. (2004) Mechanism of activation of the RAF-ERK signaling pathway by oncogenic mutations of B-RAF. Cell 116: 855–867.PTC WanMJ GarnettSM RoeS. LeeD. Niculescu-Duvaz2004Mechanism of activation of the RAF-ERK signaling pathway by oncogenic mutations of B-RAF.Cell116855867

11. Murphy E, Steenbergen C (2005) Inhibition of GSK-3b as a target for cardioprotection: the importance of timing, location, duration and degree of inhibition. Expert Opin Ther Targets 9: 447–456.E. MurphyC. Steenbergen2005Inhibition of GSK-3b as a target for cardioprotection: the importance of timing, location, duration and degree of inhibition.Expert Opin Ther Targets9447456

12. Vlahos CJ, McDowell SA, Clerk A (2003) Kinases as therapeutic targets for heart failure. Nat Rev Drug Discov 2: 99–113.CJ VlahosSA McDowellA. Clerk2003Kinases as therapeutic targets for heart failure.Nat Rev Drug Discov299113

13. Giamas G, Stebbing J, Vorgia CE, Knippschild U (2007) Protein kinases as targets for cancer treatment. Pharmacogenomics 8: 1005–1016.G. GiamasJ. StebbingCE VorgiaU. Knippschild2007Protein kinases as targets for cancer treatment.Pharmacogenomics810051016

14. Morgan H, John K (2002) The conformational plasticity of protein kinase. Cell 109: 275–282.H. MorganK. John2002The conformational plasticity of protein kinase.Cell109275282

15. Liu Y, Gray NS (2006) Rational design of inhibitors that bind to inactive kinase conformations. Nat Chem Biol 2: 358–364.Y. LiuNS Gray2006Rational design of inhibitors that bind to inactive kinase conformations.Nat Chem Biol2358364

16. Okram B, Nagle A, Adrian FJ, Lee C, Ren P, et al. (2006) A general strategy for creating “Inactive-conformation” Abl inhibitors. Chem Biol 13: 779–786.B. OkramA. NagleFJ AdrianC. LeeP. Ren2006A general strategy for creating “Inactive-conformation” Abl inhibitors.Chem Biol13779786

17. Pargellis C, Tong L, Churchill L, Cirillo PF, Gilmore T, et al. (2002) Inhibition of p38 MAP kinase by utilizing a novel allosteric binding site. Nat Struct Mol Biol 9: 268–272.C. PargellisL. TongL. ChurchillPF CirilloT. Gilmore2002Inhibition of p38 MAP kinase by utilizing a novel allosteric binding site.Nat Struct Mol Biol9268272

18. Mol CD, Dougan DR, Schneider TR, Skene RJ, Kraus ML, et al. (2004) Structural basis for the autoinhibition and STI-571 inhibition of c-Kit tyrosine kinase. J Biol Chem 279: 31655–31663.CD MolDR DouganTR SchneiderRJ SkeneML Kraus2004Structural basis for the autoinhibition and STI-571 inhibition of c-Kit tyrosine kinase.J Biol Chem2793165531663

19. Cowan-Jacob SW, Fendrich G, Floersheimer A, Furet P, Liebetanz J, et al. (2007) Structural biology contributions to the discovery of drugs to treat chronic myelogenous leukaemia. Acta Crystallogr D Biol Crystallogr 63: 80–93.SW Cowan-JacobG. FendrichA. FloersheimerP. FuretJ. Liebetanz2007Structural biology contributions to the discovery of drugs to treat chronic myelogenous leukaemia.Acta Crystallogr D Biol Crystallogr638093

20. Zou J, Wang YD, Ma FX, Xiang ML, Shi B, et al. (2008) Detailed conformational dynamics of juxtamembrane region and activation loop in c-Kit kinase activation process. Proteins 72: 323–332.J. ZouYD WangFX MaML XiangB. Shi2008Detailed conformational dynamics of juxtamembrane region and activation loop in c-Kit kinase activation process.Proteins72323332

21. Leech J, Prins J, Hermans J (1996) SMD: Visual steering of molecular dynamics for protein design. IEEE Comput Sci Eng 3: 38–45.J. LeechJ. PrinsJ. Hermans1996SMD: Visual steering of molecular dynamics for protein design.IEEE Comput Sci Eng33845

22. Isralewitz B, Gao M, Schulten K (2001) Steered molecular dynamics and mechanical functions of proteins. Curr Opin Struct Biol 11: 224–230.B. IsralewitzM. GaoK. Schulten2001Steered molecular dynamics and mechanical functions of proteins.Curr Opin Struct Biol11224230

23. Grubmüller H, Heymann B, Tavan P (1996) Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science 271: 997–999.H. GrubmüllerB. HeymannP. Tavan1996Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force.Science271997999

24. Gua Y, Shrivastava IH, Amara SG, Bahar I (2009) Molecular simulations elucidate the substrate translocation pathway in a glutamate transporter. Proc Natl Acad Sci U S A 106: 2589–2594.Y. GuaIH ShrivastavaSG AmaraI. Bahar2009Molecular simulations elucidate the substrate translocation pathway in a glutamate transporter.Proc Natl Acad Sci U S A10625892594

25. Martinez L, Polikarpov I, Skaf MS (2008) Only subtle protein conformational adaptations are required for ligand binding to thyroid hormone receptors: Simulations using a novel multipoint steered molecular dynamics approach. J Phys Chem B 112: 10741–10751.L. MartinezI. PolikarpovMS Skaf2008Only subtle protein conformational adaptations are required for ligand binding to thyroid hormone receptors: Simulations using a novel multipoint steered molecular dynamics approach.J Phys Chem B1121074110751

26. Martinez L, Webb P, Polikarpov I, Skaf MS (2006) Molecular dynamics simulations of ligand dissociation from thyroid hormone receptors: Evidence of the likeliest escape pathway and its implications for the design of novel ligands. J Med Chem 49: 23–26.L. MartinezP. WebbI. PolikarpovMS Skaf2006Molecular dynamics simulations of ligand dissociation from thyroid hormone receptors: Evidence of the likeliest escape pathway and its implications for the design of novel ligands.J Med Chem492326

27. Sonoda MT, Martinez L, Webb P, Skaf MS, Polikarpov I (2008) Ligand dissociation from estrogen receptor Is mediated by receptor dimerization: evidence from molecular dynamics simulations. Mol Endocrinol 22: 1565–1578.MT SonodaL. MartinezP. WebbMS SkafI. Polikarpov2008Ligand dissociation from estrogen receptor Is mediated by receptor dimerization: evidence from molecular dynamics simulations.Mol Endocrinol2215651578

28. Genest D, Garnier N, Arrault A, Marot C, Morin-Allory L, et al. (2008) Ligand-escape pathways from the ligand-binding domain of PPARg receptor as probed by molecular dynamics simulations. Eur Biophys J 37: 369–379.D. GenestN. GarnierA. ArraultC. MarotL. Morin-Allory2008Ligand-escape pathways from the ligand-binding domain of PPARg receptor as probed by molecular dynamics simulations.Eur Biophys J37369379

29. Carlsson P, Burendahl S, Nilsson L (2006) Unbinding of retinoic acid from the retinoic acid receptor by random expulsion molecular dynamics. Biophys J 91: 3151–3161.P. CarlssonS. BurendahlL. Nilsson2006Unbinding of retinoic acid from the retinoic acid receptor by random expulsion molecular dynamics.Biophys J9131513161

30. Huang H, Ozkirimli E, Post CB (2009) Comparison of Three Perturbation Molecular Dynamics Methods for Modeling Conformational Transitions. J Chem Theory Comput 5: 1304–1314.H. HuangE. OzkirimliCB Post2009Comparison of Three Perturbation Molecular Dynamics Methods for Modeling Conformational Transitions.J Chem Theory Comput513041314

31. Piccaluga PP, Martinelli G, Baccarani M (2006) Advances in the treatment for haematological malignancies. Expert Opin Pharmacother 7: 721–732.PP PiccalugaG. MartinelliM. Baccarani2006Advances in the treatment for haematological malignancies.Expert Opin Pharmacother7721732

32. Joensuu H, Roberts PJ, Sarlomo-Rikala M, Andersson LC, Tervahartiala P, et al. (2001) Effect of the tyrosine kinase inhibitor STI571 in a patient with a metastatic gastrointestinal stromal tumor. N Engl J Med 344: 1052–1056.H. JoensuuPJ RobertsM. Sarlomo-RikalaLC AnderssonP. Tervahartiala2001Effect of the tyrosine kinase inhibitor STI571 in a patient with a metastatic gastrointestinal stromal tumor.N Engl J Med34410521056

33. van Oosterom AT, Judson I, Verweij J, Stroobants S, di Paola ED, et al. (2001) Safety and efficacy of imatinib (ST1571) in metastatic gastrointestinal stromal tumours: a phase I study. Lancet 358: 1421–1423.AT van OosteromI. JudsonJ. VerweijS. StroobantsED di Paola2001Safety and efficacy of imatinib (ST1571) in metastatic gastrointestinal stromal tumours: a phase I study.Lancet35814211423

34. Nagar B, Hantschel O, Young MA, Scheffzek K, Veach D, et al. (2003) Structural basis for the autoinhibition of c-Abl tyrosine kinase. Cell 112: 859–871.B. NagarO. HantschelMA YoungK. ScheffzekD. Veach2003Structural basis for the autoinhibition of c-Abl tyrosine kinase.Cell112859871

35. Hubbard SR (2004) Juxtamembrane autoinhibition in receptor tyrosine kinases. Nat Rev Mol Cell Biol 5: 464–470.SR Hubbard2004Juxtamembrane autoinhibition in receptor tyrosine kinases.Nat Rev Mol Cell Biol5464470

36. Jarzynski C (1997) Nonequilibrium equality for free energy differences. Phys Rev Lett 78: 2690–2693.C. Jarzynski1997Nonequilibrium equality for free energy differences.Phys Rev Lett7826902693

37. Liphardt J, Dumont S, Smith SB, Tinoco I, Bustamante C (2002) Equilibrium information from nonequilibrium measurements in an experimental test of Jarzynski's equality. Science 296: 1832–1835.J. LiphardtS. DumontSB SmithI. TinocoC. Bustamante2002Equilibrium information from nonequilibrium measurements in an experimental test of Jarzynski's equality.Science29618321835

38. Schlitter J, Engels M, Krüger P, Jacoby E, Wollmer A (1993) Targeted molecular dynamics simulation of conformational change-application to the t↔r transition in insulin. Mol Simul 10: 291–308.J. SchlitterM. EngelsP. KrügerE. JacobyA. Wollmer1993Targeted molecular dynamics simulation of conformational change-application to the t↔r transition in insulin.Mol Simul10291308

39. Schlitter J, Engels M, Krüger P (1994) Targeted molecular dynamics: a new approach for searching pathways of conformational transitions. J Mol Graphics 12: 84–89.J. SchlitterM. EngelsP. Krüger1994Targeted molecular dynamics: a new approach for searching pathways of conformational transitions.J Mol Graphics128489

40. Marchi M, Ballone P (1999) Adiabatic bias molecular dynamics: A method to navigate the conformational space of complex molecular systems. J Chem Phys 110: 3697–3702.M. MarchiP. Ballone1999Adiabatic bias molecular dynamics: A method to navigate the conformational space of complex molecular systems.J Chem Phys11036973702

41. Paci E, Karplus M (1999) Forced unfolding of fibronectin type 3 modules: an analysis by biased molecular dynamics simulations. J Mol Biol 288: 441–459.E. PaciM. Karplus1999Forced unfolding of fibronectin type 3 modules: an analysis by biased molecular dynamics simulations.J Mol Biol288441459

42. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, et al. (2005) Scalable molecular dynamics with NAMD. J Comput Chem 26: 1781–1802.JC PhillipsR. BraunW. WangJ. GumbartE. Tajkhorshid2005Scalable molecular dynamics with NAMD.J Comput Chem2617811802

43. MacKerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, et al. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 102: 3586–3616.AD MacKerellD. BashfordM. BellottRL DunbrackJD Evanseck1998All-atom empirical potential for molecular modeling and dynamics studies of proteins.J Phys Chem B10235863616

44. Beglov D, Roux B (1994) Finite representation of an infinite bulk system: Solvent boundary potential for computer simulations. J Chem Phys 100: 9050–9063.D. BeglovB. Roux1994Finite representation of an infinite bulk system: Solvent boundary potential for computer simulations.J Chem Phys10090509063

45. Wang JM, Wolf RM, Caldwell JW, Kollman PA, Case DA (2004) Development and testing of a general amber force field. J Comput Chem 25: 1157–1174.JM WangRM WolfJW CaldwellPA KollmanDA Case2004Development and testing of a general amber force field.J Comput Chem2511571174

46. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML (1983) Comparison of simple potential functions for simulating liquid water. J Chem Phys 79: 926–935.WL JorgensenJ. ChandrasekharJD MaduraRW ImpeyML Klein1983Comparison of simple potential functions for simulating liquid water.J Chem Phys79926935

47. Feller SE, Zhang YH, Pastor RW, Brooks BR (1995) Constant-pressure molecular-dynamics simulation - the langevin piston method. J Chem Phys 103: 4613–4621.SE FellerYH ZhangRW PastorBR Brooks1995Constant-pressure molecular-dynamics simulation - the langevin piston method.J Chem Phys10346134621

48. J-P R, G C, HJC B (1977) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys 23: 327–341.R. J-PC. GB. HJC1977Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.J Comput Phys23327341

49. S M, PA K (1992) SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models.. J Comput Chem 13: 952–962.M. SK. PA1992SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models..J Comput Chem13952962

50. Darden T, York D, Pedersen L (1993) Particle mesh ewald - an n.log(n) method for ewald sums in large systems. J Chem Phys 98: 10089–10092.T. DardenD. YorkL. Pedersen1993Particle mesh ewald - an n.log(n) method for ewald sums in large systems.J Chem Phys981008910092