The early stage of chronic myeloid leukemia is triggered by the tyrosine kinase Bcr-Abl. Imatinib mesylate, a selective inhibitor of Bcr-Abl, has been successful in chronic myeloid leukemia clinical trials, but short-lived remissions are usually observed in blast crisis patients. Sequencing of the BCR-ABL gene in relapsed patients revealed a set of mutants that mediate drug resistance. Previously reported work postulated that the missense T315I mutation both alters the three-dimensional structure of the protein binding site, thus decreasing the protein sensitivity for the drug, and does not feature a fundamental hydrogen bond that is critical for binding with imatinib. These speculations, however, were not supported by investigations at the molecular modeling level. Here, we present the results obtained from the application of molecular dynamics simulations to the study of the interactions between T315I Bcr-Abl and imatinib. For the first time, we show that, with respect to the wild-type system, the absence of the supposedly critical H-bond is not the only cause for the failure of receptor inhibition by imatinib, but also a plethora of other protein/drug interactions are drastically and unfavorably changed in the mutant protein.

Chronic myeloid leukemia (CML) is a clonal disease involving the pluripotent hematopoietic stem cell compartment and is associated with the Philadelphia chromosome, a reciprocal translocation between chromosome 9 and 22 (14). This translocation links the c-Abl tyrosine kinase oncogene on chromosome 9 to the 5′ half of the BCR gene on chromosome 22 and originates the fusion gene BCR-ABL (4, 5). The fusion gene produces a chimeric 8.5 kb transcript that codes for the p210BCR-ABL protein (6), which possesses constitutive tyrosine kinase activity and is the pathogenic agent of CML. From the therapeutic perspective, as CML arises from a single genetic lesion, a research goal has been to develop kinase-specific inhibitors. The development of imatinib (also known as imatinib mesylate, Glivec, or Gleevec) is a landmark achievement in this respect as it has shown great promise in the chronic phase (7, 8), and some expectations in the accelerated and blastic phase of CML, as well as in BCR-ABL–expressing acute lymphoblastic leukemias (9). Unfortunately, however, most responding blast-phase patients relapse despite continued chemotherapy (10) and, as frequently happens with cancer chemotherapeutics, resistance to imatinib has been reported in both Bcr-Abl–expressing cell lines and in patients with CML (1114).

Recently, point mutations within the DNA sequence coding for the ATP-binding pocket of the BCR-ABL gene were identified in cells from some patients with CML who had imatinib refractory disease or who had a relapse during the treatment (1518). In particular, a point mutation resulting in a threonine-to-isoleucine change at amino acid position 315 (T315I) has been described in detail (14). This mutation, when engineered into wild-type (WT) p210BCR-ABL and transiently infected into 293T cells or Ba/F3 cells, interfered with the inhibition of Bcr-Abl kinase activity in cells exposed to imatinib (14, 18). Further, Corbin et al. (15) showed that the T315I-mutated Abl kinase domain exhibited no significant inhibition at imatinib concentrations 200-fold higher than the IC50 value of the WT kinase; it also showed a 2-fold increase in its ATP affinity relative to the WT protein. However, the extent to which this mutation contributes to imatinib resistance in vivo remains a matter of lively literature dispute (14, 17, 1922).

On the basis of these evidences, both Schindler et al. (23) and Corbin et al. (15) postulated that the above-mentioned point mutation T315I in the tyrosine kinase domain of Bcr-Abl alters the three-dimensional structure of the ATP pocket, thus decreasing the sensitivity of Bcr-Abl to imatinib, and does not feature a fundamental hydrogen bond that is critical for binding with the drug. These speculations, however, were not supported by further investigations at a molecular modeling level. Accordingly, we present the results obtained from a detailed molecular simulation study of both WT and T315I-mutated Abl kinase domain in complex with imatinib, aiming at offering, for the first time, both a qualitative and a quantitative picture of the molecular mechanism of failure of the tyrosine kinase inhibitor binding to the mutated protein.

The original three-dimensional structure of the kinase domain of Abl in complex with imatinib (chain A, Protein Data Bank code 1IEP; ref. 24), including the 178 crystallographic water molecules, was used as a starting point. All missing hydrogen atoms were added to the protein backbone and side chains with the parse module of AMBER 7 (25). All ionizable residues were considered in the standard ionization state at neutral pH. The geometry of added hydrogen was refined for 200 steps (steepest descent) in vacuo. The parm94 all-atom force field by Cornell et al. (26) was applied using the sander module of AMBER 7. Eventual missing force field parameters for inhibitors were generated as follows. AM1 (27) geometry optimization of the structure was followed by restricted Hartree-Fock/6-31G* single-point calculation to obtain the electrostatic potentials. Next, the restrained potential method (28) was used for charge fitting. The missing bonds, angle torsions, or van der Waals parameters not included in the parm94 force field were transferred from the newly developed parm99 parameter set (29) or from a general AMBER force field (30).

Mutation T315I was introduced into the WT crystal structure of Abl/imatinib complex using the Biopolymer module of Insight II (v. 2001, Accerlys, Inc., San Diego, CA) by swapping the mutant residue into the specific site (31). As this mutant has several possible rotamers, each rotamer structure was generated and minimized separately. A systematic conformational search for 100 conformations was carried out for the mutant residue. For each conformation, a steric bump check was executed first to avoid serious sterical clash between the mutant residue and any other residue in the molecule. If any atom of the built-in residue was closer than 1 Å to any atom of any other residue in that molecule, that conformation of the built-in residue was discarded. Each surviving mutant complex was then minimized using the sander module of AMBER 7 with steepest descent algorithm for the first 100 steps, and then subsequently with a conjugate gradient approach. A cutoff of 10 Å and a distance-dependent dielectric were applied. The minimization protocol was as follows: First, the mutated residue was minimized for a maximum of 500 steps, whereas the rest of the system was held fixed; second, only the hydrogen atoms were minimized for a maximum of 500 steps; finally, all atoms were minimized for 1,000 steps. Six rotamers were obtained from the above-described procedure, all compatible with the local backbone conformation. The energy difference among these conformers was small and confined to the range 0.28 to 0.82 kcal/mol. The structure whose rotamer led to the lowest final energy value was selected for the subsequent molecular dynamics experiments. To set up each complex for the simulation, we adopted the following ansatz: first, to let the protein relax in an aqueous environment, each complex was immersed in a 75 Å radius sphere of TIP3P water molecules (32). The resulting system was minimized with a gradual decrease in the position restraints of the protein atoms. At the end of the relaxation process, all water molecules beyond the first hydration shell (i.e., at a distance >3.5 Å from any protein atom) were removed. Finally, to achieve electroneutrality, a suitable number of counterions were added in the positions of largest electrostatic potential as determined by the module cion within AMBER 7. To reduce computational time to reasonable limits, all protein residues with any atom closer than 30 Å from the center of mass of considered mutation were chosen to be flexible in the dynamic simulations. Subsequently, a spherical TIP3P water cap of radius equal to 30 Å was centered on imatinib in the corresponding complex, including the hydrating water molecules within the sphere resulting from the previous step. After energy minimization of the new water cap for 1,500 steps, keeping the protein, the inhibitor, and the preexisting waters rigid, followed by a molecular dynamics equilibration of the entire water sphere with fixed solute for 20 picoseconds, further unfavorable interactions within the structures were relieved by progressively smaller positional restraints on the solute [from 25 to 0 kcal/(mol Å2)] for a total of 4,000 steps. Each system was gradually heated to 27°C in three intervals, allowing a 5-picosecond interval per each 100°C, and then equilibrated for 50 picoseconds at 27°C, followed by 400 picoseconds of data collection runs, necessary for the estimation of the free energy of binding (vide infra). After the first 20 picoseconds of molecular dynamics equilibration, additional TIP3P water molecules were added to the 30 Å water cap to compensate for molecules that were able to diffuse into gaps of the enzyme. The molecular dynamics simulations were done at a constant temperature (T = 27°C) using the Berendsen coupling algorithm (33) with separate coupling of the solute and solvent to the heat, an integration time step of 2 femtoseconds, and the applications of the SHAKE algorithm (34) to constrain all bonds to their equilibrium values, thus removing high-frequency vibrations. All intrasolute pairwise interactions (i.e., the bond, angle, dihedral, van der Waals, and electrostatics terms of the molecular mechanics force field) were accumulated without a distance-based cutoff. Long-range nonbonded interactions were truncated by using a dual cutoff of 12 and 30 Å, respectively, where energies and forces due to interactions between 12 and 30 Å were updated every 20 time step. The same frequency of update was used for the nonbonded list. For the calculation of the binding free energy between kinase domains of Abl and imatinib in water, a total of 400 snapshots were saved during the molecular dynamics data collection period described above, one snapshot per each 1 picosecond of molecular dynamics simulation.

All energetics analyses were done for only a single molecular dynamics trajectory of each Abl/imatinib complex considered, with unbound protein and substrate snapshots taken from the snapshots of that trajectory. According to the so-called molecular mechanics/Poisson-Boltzmann surface area method (35), the binding free energy between imatinib and Abl can be calculated as:

\[\mathrm{{\Delta}}\mathit{G}_{\mathrm{bind}}\ =\ \mathrm{{\Delta}}\mathit{E}_{\mathrm{MM}}\ +\ \mathrm{{\Delta}}\mathit{G}_{\mathrm{solv}}\ {-}\ \mathit{T}\mathrm{{\Delta}}\mathit{S}\]

where:

\[\mathrm{{\Delta}}\mathit{G}_{\mathrm{solv}}\ =\ \mathrm{{\Delta}}\mathit{G}_{\mathrm{PB}}\ +\ \mathrm{{\Delta}}\mathit{G}_{\mathrm{NP}}\]

ΔEMM denotes the sum of molecular mechanics energies, and can be further split into contributions from electrostatic (EEL) and van der Waals (EvdW) energies:

\[\mathrm{{\Delta}}\mathit{E}_{\mathrm{MM}}\ =\ \mathrm{{\Delta}}\mathit{E}_{\mathrm{EL}}\ +\ \mathrm{{\Delta}}\mathit{E}_{\mathrm{vdW}}\]

The terms in Eq. C were calculated by using the carnal and anal modules of AMBER 7. To examine the individual contribution of the protein residues to ΔEMM, the collected frames were further processed, and the interaction energy between imatinib and the binding site was decomposed on a residue basis again using the anal module of AMBER 7.

ΔGsolv represents the solvation free energy. The polar solvation process (i.e., ΔGPB in Eq. B) is equivalent to the transfer of a protein from one medium with dielectric constant equal to that of the interior of the protein to another medium with dielectric constant equal to that of the exterior of the protein. This term yields a free energy because it corresponds to the work done to reversibly charge the solute, and it is a polarization free energy because the work goes to the polarization of the solvent. The nonpolar solvation contribution (i.e., ΔGNP in Eq. B) includes cavity creation in water and van der Waals interactions between the modeled nonpolar protein and water molecules. This term can be imagined as transferring a nonpolar molecule with the shape of the protein from vacuum to water. The polar component of ΔGsolv was evaluated with the Poisson-Boltzmann approach (36). This procedure involves using a continuum solvent model, which represents the solute as a low dielectric medium (i.e., of dielectric constant ε = 1) with embedded charges and the solvent as a high dielectric medium (ε = 80) with no salt. All atomic charges were taken from the Cornell et al. (26) force field and from our ab initio calculations (see above), because these are consistent with the molecular mechanics energy calculations. However, as suggested by Chong et al. (37), the atomic radii were taken from the parse parameter set instead of the parm94 force field set because of the small size of hydrogens in the latter. The dielectric boundary is the contact surface between the radii of the solute and the radius (1.4 Å) of a water molecule. The numerical solution of the linearized Poisson-Boltzmann equations were solved on a cubic lattice by using the iterative finite-difference method implemented in the DelPhi software package (38). The grid size used was 0.5 Å. Potentials at the boundaries of the finite-difference lattice were set to the sum of the Debye-Hückel potentials. The nonpolar contribution to the solvation energy was calculated as ΔGNP = γ (solvent-accessible surface area) + β, where γ = 0.00542 kcal/Å2, β = 0.92 kcal/mol, and solvent-accessible surface area is the solvent-accessible surface estimated with the MSMS program (39).

The normal-mode analysis approach was followed to estimate the last parameter (i.e., the change in solute entropy upon association −TΔS; ref. 40). In the first step of this calculation, an 8 Å sphere around imatinib was cut out from a molecular dynamics snapshot for each inhibitor-protein complex. This value was shown to be large enough to yield converged mean changes in solute entropy. On the basis of the size-reduced snapshots of the complex, we generated structures of the uncomplexed reactants by removing the atoms of the protein and ligand, respectively. Each of those structures was minimized using a distance-dependent dielectric constant, ε = 4r, to account for solvent screening; its entropy was calculated using classic statistical formulas and normal-mode analysis. To minimize the effects due to different conformations adopted by individual snapshots, we averaged the estimation of entropy over 10 snapshots.

In our comparison of binding free energies of WT and mutant protein/ligand complexes, it is important to ensure that each of these systems is thoroughly equilibrated ahead of data collection. Starting from the solvated and equilibrated systems, the evolution of the dynamics of each simulated complex is consistent with a stabilized system with minimal fluctuations in total energy being exhibited as a function of time (data not shown). For instance, in the case of the WT protein complex, the overall protein structure remains well defined, with an average root-mean-square deviation from the crystal structure of ∼1.5 Å. By comparison, the root-mean-square deviation for all the heavy atoms is predictably larger, ∼2.1 to 2.3 Å. These root-mean-square deviation values are typical of the deviations observed between crystal and solution conformation in careful molecular dynamics simulations. It is worth noting that no restraints were required to keep the imatinib molecule properly in place within its binding site. This implies that the force field and water-type molecules adopted here are appropriate to model the attraction between the two moieties.

From the standpoint of the energetics of interactions of this system, Table 1 reports the binding free energy, all energy terms contributing to ΔGbind, and the corresponding ΔGbind value calculated for both WT (see Fig. 1A) and T315I (Fig. 1C) Abl/imatinib complex. Table 2 lists the calculated individual contributions of imatinib-contacting protein residues to EMM. Clearly, the native protein binds tightly and efficiently to the inhibitor as expected. Among all van der Waals interactions between residues contacting imatinib and belonging to the ATP-binding loop (or P-loop), only Y253 yields a considerable favorable dispersive term (ΔEvdW,Y253 = −2.28 kcal/mol). Moreover, this residue also contributes to the electrostatics, with a notable value of ΔEEL,Y253 = −1.12 kcal/mol (Table 2). The analysis of the molecular dynamics trajectory reveals that the average dynamic distance between the C-51 atom of L248 and the C3 atom of the imatinib pyridine ring is equal to 3.61 Å, that between the C-3 atom of A269 and N-3 of the pyrimidine ring is 3.50 Å, and that between the C-41 atom of V256 and C-6 of the same cyclic moiety is 4.06 Å (see Table 3). Accordingly, the side chains of L248, V256, and A269 indeed encase the pyridine ring system of the inhibitor with hydrophobic interactions, but these forces are weak due to the relative distance of these residues from the inhibitor (ΔEvdW,L248 = −0.73, ΔEvdW,V256 = −1.41, and ΔEvdW,A269 = −1.09 kcal/mol, respectively; Tables 2 and 3). In the case of residue Y253, an important π-π stabilizing interaction occurs between the ring system of imatinib and the phenyl side chain of Y253, resulting from the edge-on geometry of these moieties. Aromatic interactions consist of van der Waals, hydrophobic, and electrostatic forces. The electronic nature of the π-π interactions, however, indeed favor the stacking of aromatic rings either by parallel-displaced (off-center) or edge-on (T stacking) geometries, whereas the face-to-face geometry is unfavorable (particularly in environments where there is a low effective dielectric constant), because the dominant interaction is π-electron repulsion. Moreover, the water-mediated hydrogen bond between the -OH group of Y253 and the side chain of N322 is still observed in our simulated system, with an average dynamic length of 2.80 Å on both sides. This feature plays an important role in imatinib binding as it contributes to keep in place the P-loop in an induced-fit mechanism that increases surface complementarity with the drug.

Table 1.

Calculated free energy of binding ΔGbind (kcal/mol) and energy contributions to ΔGbind (kcal/mol) of WT and T315 mutant Abl kinase domains and imatinib

ContributionWTT315I
ΔEvdW −71.35 ± 0.61 −67.80 ± 0.53 
ΔEEL −57.44 ± 0.11 −55.12 ± 0.20 
ΔEMM −128.79 ± 0.11 −122.92 ± 0.18 
ΔGPB 108.00 ± 1.20 105.18 ± 1.21 
ΔGNP −6.70 ± 0.03 −6.18 ± 0.01 
TΔS 17.10 17.37 
ΔGbind −10.39 (−10.37) −6.55 (<−7.23) 
ContributionWTT315I
ΔEvdW −71.35 ± 0.61 −67.80 ± 0.53 
ΔEEL −57.44 ± 0.11 −55.12 ± 0.20 
ΔEMM −128.79 ± 0.11 −122.92 ± 0.18 
ΔGPB 108.00 ± 1.20 105.18 ± 1.21 
ΔGNP −6.70 ± 0.03 −6.18 ± 0.01 
TΔS 17.10 17.37 
ΔGbind −10.39 (−10.37) −6.55 (<−7.23) 

NOTE: The experimental ΔGbind value (reported in parenthesis for comparison) was obtained from the corresponding experimental IC50 value (15) using the following relationship (41): ΔGbind = RTlnKdiss = RTln(IC50 + 0.50Cenz) ≅ RTlnIC50, where R is ideal gas constant, T is the temperature in Kelvin, and Cenz is the concentration of the enzyme, which is a small number after equilibration and can be omitted in most cases.

Figure 1.

A, ribbon representation of the structure of the WT Abl kinase domain in complex with imatinib. The inhibitor and its van der Waals surface are colored khaki. B, details of residues contacting imatinib: L248, Y253, V256 (of the P-loop), and A269 are colored pink; D381, F382, and G383 of the DFG motif of the A-loop are colored light green; E286 is colored dark green; N322 is colored purple; and T315 is colored red. C, comparison between WT (gold) and T315 (brown) mutant Abl kinase domains. D, details of residues contacting imatinib: L248, Y253, V256 (of the P-loop), and A269 are colored pink; D381, F382, and G383 of the DFG motif of the A-loop are colored light green; E286 is colored dark green; N322 is colored purple; and I315 is colored red.

Figure 1.

A, ribbon representation of the structure of the WT Abl kinase domain in complex with imatinib. The inhibitor and its van der Waals surface are colored khaki. B, details of residues contacting imatinib: L248, Y253, V256 (of the P-loop), and A269 are colored pink; D381, F382, and G383 of the DFG motif of the A-loop are colored light green; E286 is colored dark green; N322 is colored purple; and T315 is colored red. C, comparison between WT (gold) and T315 (brown) mutant Abl kinase domains. D, details of residues contacting imatinib: L248, Y253, V256 (of the P-loop), and A269 are colored pink; D381, F382, and G383 of the DFG motif of the A-loop are colored light green; E286 is colored dark green; N322 is colored purple; and I315 is colored red.

Close modal
Table 2.

Residue-based decomposition of the noncovalent interaction energies (kcal/mol) between WT and T315I mutant Abl kinase domains and imatinib

ResidueΔEvdW WTΔEEL WTΔEvdW T315IΔEEL T315I
L248 −0.73 +0.36 −1.12 +0.39 
Y253 −2.28 −1.12 −0.66 −1.23 
V256 −1.41 −0.15 −1.11 −0.30 
A269 −1.09 −0.13 −1.01 −0.26 
K271 −1.70 +1.87 +2.96 +1.50 
E286 −2.81 −7.92 −0.65 −7.33 
V289 −1.54 +0.49 −1.64 +0.28 
M290 +0.33 −0.15 +0.95 −0.17 
V299 −1.77 −0.85 −1.85 −0.51 
I313 −1.90 −0.12 −2.26 −0.01 
T315/I315 −1.31 −2.77 −1.14 −0.57 
F317 −1.27 −0.51 −2.51 −0.91 
M318 +0.94 −1.17 −1.06 −2.44 
G321 −0.99 +1.53 −0.80 +0.92 
I360 −2.23 +1.09 −2.15 +0.80 
H361 −1.70 +0.65 −3.77 +1.02 
R362 −1.95 −0.92 −0.10 +1.10 
L370 −2.42 −0.59 −1.98 −0.660 
A380 +1.84 −1.55 +4.07 −1.76 
D381 −4.83 −8.79 −1.89 +0.20 
F382 −3.42 −1.87 +2.10 −0.87 
ΔEtot −32.24 −22.62 −15.62 −10.81 
ResidueΔEvdW WTΔEEL WTΔEvdW T315IΔEEL T315I
L248 −0.73 +0.36 −1.12 +0.39 
Y253 −2.28 −1.12 −0.66 −1.23 
V256 −1.41 −0.15 −1.11 −0.30 
A269 −1.09 −0.13 −1.01 −0.26 
K271 −1.70 +1.87 +2.96 +1.50 
E286 −2.81 −7.92 −0.65 −7.33 
V289 −1.54 +0.49 −1.64 +0.28 
M290 +0.33 −0.15 +0.95 −0.17 
V299 −1.77 −0.85 −1.85 −0.51 
I313 −1.90 −0.12 −2.26 −0.01 
T315/I315 −1.31 −2.77 −1.14 −0.57 
F317 −1.27 −0.51 −2.51 −0.91 
M318 +0.94 −1.17 −1.06 −2.44 
G321 −0.99 +1.53 −0.80 +0.92 
I360 −2.23 +1.09 −2.15 +0.80 
H361 −1.70 +0.65 −3.77 +1.02 
R362 −1.95 −0.92 −0.10 +1.10 
L370 −2.42 −0.59 −1.98 −0.660 
A380 +1.84 −1.55 +4.07 −1.76 
D381 −4.83 −8.79 −1.89 +0.20 
F382 −3.42 −1.87 +2.10 −0.87 
ΔEtot −32.24 −22.62 −15.62 −10.81 
Table 3.

Average dynamic distances (Å) between selected Abl kinase domain residues and imatinib

Average dynamic distanceWTT315I
C-51 of L248–C-3 atom of imatinib pyridine ring 3.61 3.75 
C-3 of A269–N-3 of imatinib pyrimidine ring 3.50 3.25 
C-41 of V256–C-6 of imatinib pyrimidine ring 4.06 3.82 
C-3 of F382–C-6 of imatinib pyrimidine ring 2.96 3.38 
C-3 of K271–C-7 of imatinib toluene ring 3.66 3.22 
N-δ of R362–C-7 of imatinib piperazinyl ring 4.18 3.33 
Average dynamic distanceWTT315I
C-51 of L248–C-3 atom of imatinib pyridine ring 3.61 3.75 
C-3 of A269–N-3 of imatinib pyrimidine ring 3.50 3.25 
C-41 of V256–C-6 of imatinib pyrimidine ring 4.06 3.82 
C-3 of F382–C-6 of imatinib pyrimidine ring 2.96 3.38 
C-3 of K271–C-7 of imatinib toluene ring 3.66 3.22 
N-δ of R362–C-7 of imatinib piperazinyl ring 4.18 3.33 

Another important steric role is played by residue F382, belonging to the activation loop (or A-loop), which give rise to a strong off-center π-π stacking interaction with the pyrimidine ring of imatinib (average dynamic distances between the two ring centers = 3.79 Å), aptly testified by the two corresponding energy terms (ΔEvdW,F382 = −3.42 kcal/mol, ΔEEL,F382 = −1.87 kcal/mol). These interactions contribute to maintain this residue, belonging to the well-conserved DFG motif of the A-loop, in the proper position of pointing toward the ATP docking site, peculiar to Abl and crucial for effective binding of the inhibitor (Fig. 1B).

Hydrogen bonding has been claimed as the main force in imatinib binding, and the absence of some of these nonbonded interactions have been pointed out as the major cause for imatinib failure. The analysis of the molecular trajectory for the WT Bcr-Abl protein/imatinib complex reveals that an important H-bond takes place between the carboxylate oxygens of residue E286 and the amide hydrogen of the inhibitor, and is characterized by an average dynamic length of 3.28 Å. At the same time, the two hydrogens on the adjacent phenyl ring contribute to further stabilize the negative charge of this residue, being their average dynamic distance values of 2.96 and 3.45 Å, respectively. Energetically, this reflects into two of the highest favorable binding energy components for the interaction of this residue with imatinib: ΔEvdW,E286 = −2.81 kcal/mol and ΔEEL,E286 = −7.92 kcal/mol (Table 2).

The residue that exerts the strongest attraction toward imatinib is D381, which combined an optimal sterical position to the formation of a persistent hydrogen bond with the carbonyl oxygen of the inhibitor (average dynamic length = 2.62 Å); the intermolecular dispersive and electrostatic interactions are accordingly strong (ΔEvdW,D381 = −4.83 kcal/mol, ΔEEL,D381 = −8.79 kcal/mol). The high value of the intermolecular electrostatics is also accounted for by the vicinity of the negatively charged carboxylate group of D381 and some of the hydrogen atoms of the piperazinyl ring of imatinib, which bear a weak positive partial charge.

Residue 315 deserves a special comment because this position is found mutated in most CML patients. Indeed, the resulting average interaction energy between this amino acid and imatinib is less important than the values discussed above for other residues, being ΔEvdW,T315 = −1.31 kcal/mol and ΔEEL,T315 = −2.77 kcal/mol, respectively. The confined dimensions of the amino acid, with a calculated surface area of 140 Å2 and a molecular volume of 117 Å3, justify the limited value of the van der Waals term, whereas the electrostatic component is in line with the estimated value for a buried intermolecular hydrogen bond in protein/ligand systems. The average dynamic length of the H-bond between the -OH moiety of T315 and the secondary amine group of the inhibitor is 2.83 Å, in agreement with the 2.88 Å length predicted from X-ray structure (24). The same moiety is involved in another alternative hydrogen bond interaction with the backbone carbonyl of E316, characterized by an average dynamic length of 3.11 Å.

Mutating T to I at position 315 in the ATP-binding domain of Abl results in a calculated ΔGbind for imatinib of −3.84 kcal/mol lower than the corresponding WT value (Table 1). The detailed analysis of the free energy of binding components reveals that both the van der Waals and the electrostatic interactions drastically decrease in the case of the isoleucine mutant (Tables 1 and 2). The absence of the hydrogen bond between the -OH side chain of T315 and the amino group of imatinib is correctly reflected in a decrease of the electrostatic stabilizing contribution of 2.20 kcal/mol (Table 2). On the other hand, the loss of the hydrogen bond between T315 and E316 is counterbalanced by the permanence of two further hydrogen bonds between the side-chain carboxylate group of E286 and the backbone -NH group of Q300 and water 99 (average dynamic lengths = 3.01 and 2.73 Å, respectively), and a salt bridge between the same negative moiety and the positively charged side chain nitrogen atom of K378 (average dynamic distance = 2.70 Å). Further, the analysis of the corresponding molecular dynamics trajectory reveals the formation of another hydrogen bond between the negative oxygens of E286 and water 62, its average dynamic length being equal to 3.07 Å.

There are numerous instances in which the structures of threonine to nonpolar valine, leucine, and isoleucine variants of different proteins have been solved; the protein structures seem to be very tolerant to substitution. Valine and threonine are both β-branched side chains and can be considered isosteric. Despite this similarity in shape and atomic volumes, valine seems to have a somewhat larger van der Waals volume than threonine (58.4 versus 48.6 Å3). Thus, the average distance from the neighboring atoms will be longer for the hydroxyl group of threonine than for the methyl group of valine. Because packing interactions are strongly distance dependent, this should lead to a decrease in the strength of packing interactions. However, the polarizability of hydroxyl oxygen is higher than that for carbon, and thus the minimum of the packing interaction (as modeled by the 6–12 potential used here) with the hydroxyl oxygen atom is deeper than for the carbon atom. Thus, one should expect an overall similar strength in packing interaction for substitutions of threonine with valine because of compensation due to increase in distance and polarizability. This could explain why this mutation has not been reported yet in patients. On the other hand, the isoleucine side chain is larger than that of threonine (48.6 versus 73.0 Å3, respectively); from the model, it seems, however, that it can be accommodated without severe van der Waals repulsions (ΔEvdW,I315 = 1.14 kcal/mol). Keeping in mind the foregoing discussion, the only possible reason for the destabilization of the threonine-to-isoleucine change can be the difference in Gibbs free energy of hydration. Indeed, the calculated difference in ΔGsolv of threonine versus isoleucine is equal to 7.19 kcal/mol.

In a sort of a domino effect, however, the substitution of a single residue in this position induces several, sometimes substantial, modifications in the conformation of other residues, both belonging to the binding site and the surrounding areas as a direct consequence of the different positions assumed by imatinib in the binding site (Fig. 1C and D). The residues that register the most consistent loss of favorable interactions with the inhibitor in the presence of the I315 substitutions are K271, E286, R362, and D381. The readjustment of K271 results in an unfavorable sterical accommodation of its side chain toward the methyl-substituted ring of imatinib (ΔEvdW,K271 = +2.96 kcal/mol; Table 3), the electrostatic interactions being almost unaltered (ΔEEL,K271 = +1.50 kcal/mol; Table 2). Although E286 remains one of the strongest imatinib contact points, its slight conformational modification results in a decrease of their mutual dispersive interactions as evidenced by the corresponding value of ΔEvdW,E286 = −0.65 kcal/mol (Table 2). The electrostatic environment described for the WT case remains almost unaltered in the presence of the isoleucine in position 315, as detected both from the visual inspection of the corresponding molecular trajectory and the corresponding Coulombic term (ΔEEL,E286 = − 7.33 kcal/mol). R362 is forced to find a new conformation in which it is not able to optimize the corresponding interactions with imatinib (ΔEvdW,R362 = −0.10 kcal/mol, ΔEEL,R362 = +1.10 kcal/mol). This is basically due to the proximity of the R362 side chain to the piperazinyl ring of the inhibitor, which results in both a sterical and electrostatic repulsion (Table 3). The residue most destabilized in its interaction with imatinib in the presence of the T315I mutation is D381, for which ΔEvdW,D381 = −1.89 kcal/mol and ΔEEL,D381 = +0.20 kcal/mol. The corresponding trajectory shows that there is a different positioning of this residue with respect to imatinib, leading to a less favorable geometry for the hydrogen bonding. Moreover, this has an influence on the conformation of F382 of the DFG motif, which is no longer maintained in the apt position in the mutant trajectory and is reflected in the net loss of a favorable stabilizing interaction (ΔEvdW,F382 = +2.10 kcal/mol, ΔEEL,F382 = −0.87 kcal/mol; Fig. 1D).

The results obtained from the application of the molecular mechanics/Poisson-Boltzmann surface area method to the evaluation of the binding free energy of WT and T315I mutant Abl-core domains with imatinib led us to the conclusion that the analyzed mutation can be claimed to be one of the molecular causes of resistance when identified in patients. Although other authors have used molecular modeling to study the T315I mutation and other known Bcr-Abl mutations from a rigorous but qualitative point of view (42), we showed for the first time that the missing-hydrogen-bond hypothesis is not the leading cause of failure of imatinib binding to the protein; instead, it results from a domino effect induced by the conformational readjustment necessary to accommodate the mutant residue and involves several other important drug contact point.

Undoubtedly, any simulation result strongly relies on the quality of the starting model and of the procedure adopted. Indeed, the fact that the protein may adopt a different, “active” conformation has been invoked to explain resistance in the presence of the I315 mutation. On the other hand, many authoritative papers have presented only speculative explanations just based on the WT structure. Simulating the transition between an inactive (T315) to an active (I315) conformation would be, in our opinion, not only impractical but also somewhat arbitrary. The goal of this paper was to present a computational approach that is able to quantify, and not only to qualify, the interactions between imatinib and Abl core domain in wild type as well as in those mutated at position 315. Further, we wanted to show that other important chemicophysical parameters so often neglected, such as the energy of solvation, can play a major role in the presence of a mutated residue.

Finally, the major purpose of this work, which is part of a larger, ongoing project, is to propose a reliable and relatively fast computational recipe for the classification of different mutations in Abl kinase domain with respect to their affinity toward imatinib. Once we have tested a number of Abl kinase domain mutations, differing by position and affinity toward the inhibitor, and verified the correspondence between experimental and computational evidences, we could claim that this technique could be adopted as a routine analysis and can be flanked to all other techniques already used in cancer research to define a better therapeutic strategy.

Grant support: Italian Association for Cancer Research, Italian Ministry of Health, and Italian Ministry for University and Scientific Research.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1
Fialkow PJ, Jacobson RJ, Papayannopoulou T. Chromic myelocytic leukemia: clonal origin in a stem cell common to granulocyte, erythrocyte, platelet and monocyte/macrophage.
Am J Med
1977
;
63
:
125
–30.
2
Caspersson T, Gahrton G, Lindsten J, Zech L. Identification of the Philadelphia chromosome on a number 22 by quinacrine mustard fluorescence analysis.
Exp Cell Res
1970
;
63
:
238
–40.
3
Rowley JD. A new consistent chromosomal abnormality in chronic myelogenous leukemia identified by crinaquine, fluorescence and Giemsa staining.
Nature
1973
;
243
:
290
–3.
4
Heisterkamp N, Stephenson JR, Groffen J, et al. Localization of the c-abl oncogene adjacent to a translocation break point in chronic myelocytic leukemia.
Nature
1983
;
306
:
239
–42.
5
Groffen J, Stephenson JR, Heisterkamp N, de Klein A, Bartram CR, Grosveld G. Philadelphia chromosomal breakpoints are clustered within a limited region, bcr, on chromosome 22.
Cell
1984
;
36
:
93
–9.
6
Shtivelman E, Lifshitz B, Gale RP, Canaani E. Fused transcript of ABL and BCR genes in chronic myelogeneous leukemia.
Nature
1985
;
315
:
550
–4.
7
Druker BJ, Tamura S, Buchdunger E, et al. Effects of a selective inhibitor of the Abl tyrosine kinase on the growth of Bcr-Abl positive cells.
Nat Med
1996
;
2
:
561
–6.
8
Druker BJ, Talpaz M, Resta DJ, et al. Efficacy and safety of a specific inhibitor of the BCR-ABL tyrosine kinase in chronic myeloid leukemia.
N Engl J Med
2001
;
344
:
1031
–7.
9
Druker BJ, Sawyers CL, Kantarjian H, et al. Activity of a specific inhibitor of the BCR-ABL tyrosine kinase in the blast crisis of chronic myeloid leukemia and acute lymphoblastic leukemia with the Philadelphia chromosome.
N Engl J Med
2001
;
44
:
1038
–42.
10
Savage D, Antman K. Imatinib mesylate: a new oral targeted therapy.
N Engl J Med
2002
;
346
:
683
–93.
11
le Coutre P, Tassi E, Varella-Garcia M, et al. Induction of resistance to the Abelson inhibitor STI571 in human leukemic cells through gene amplification.
Blood
2000
;
95
:
1758
–66.
12
Weisberg E, Griffin JD. Mechanism of resistance to the Abl tyrosine kinase inhibitor STI571 in BCR/ABL-transformed hematopoietic cell lines.
Blood
2000
;
95
:
3498
–505.
13
Mahon FX, Deininger MW, Schultheis B, et al. Selection and characterization of BCR-ABL positive cell lines with differential sensitivity to the tyrosine kinase inhibitor STI571: diverse mechanism of resistance.
Blood
2000
;
96
:
1070
–9.
14
Gorre ME, Mohammed M, Ellwood K, et al. Clinical resistance to STI571 cancer therapy caused by BCR-ABL gene mutation or amplification.
Science
2001
;
293
:
876
–80.
15
Corbin AS, Buchdunger E, Furet P, Druker BJ. Analysis of the structural basis of specificity of the Abl kinase by STI571.
J Biol Chem
2001
;
277
:
32214
–9.
16
Shah NP, Nicoll JM, Nagar B, et al. Multiple BCR-ABL kinase domain mutations confer polyclonal resistance to the tyrosine kinase inhibitor imatinib (STI571) in chronic phase and blast crisis chronic myeloid leukemia.
Cancer Cell
2002
;
2
:
117
–25.
17
Ricci C, Scappini B, Divoky V, et al. Mutation in the ATP-binding pocket of the ABL kinase domain in an STI571-resistant BCR/ABL-positive cell line.
Cancer Res
2002
;
62
:
5995
–8.
18
von Boubnoff N, Schneller F, Peschel C, Duyster J. BCR-ABL gene mutations in relation to clinical resistance of Philadelphia chromosome-positive leukemia to STI571: a prospective study.
Lancet
2002
;
359
:
487
–91.
19
Barthe C, Cony-Makhoul P, Melo JV, Reiffers J, Mahon FX. Roots of clinical resistance to STI-571 cancer therapy.
Science
2001
:
293
:
2163a
.
20
Hochhaus A, Kreil S, Corbin AS, et al. Molecular and chromosomal mechanisms of resistance to imatinib (STI571) therapy.
Leukemia
2002
;
16
:
2190
–6.
21
Gambacorti-Passerini C, Corneo G, D'Incalci M. Roots of clinical resistance to STI-571 cancer therapy.
Science
2001
:
293
:
2163a
.
22
Gorre ME, Shah N, Ellwood K, Nicoll J, Sawyers CL. Roots of clinical resistance to STI-571 cancer therapy.
Science
2001
;
293
:
2163a
.
23
Schindler T, Bornmann W, Pellicena P, Miller WT, Clarkson B, Kuriyan J. Structural mechanism for STI-571 inhibition of Abelson tyrosine kinase.
Science
2000
;
289
:
938
–1942.
24
Nagar B, Bornmann WG, Pellicena P, et al. Crystal structures of the kinase domain of c-Abl in complex with the small molecule inhibitors PD173955 and Imatinib (STI-571).
Cancer Res
2002
;
62
:
4236
–43.
25
Case DA, Pearlman DA, Caldwell JD, et al. AMBER 7. San Francisco (CA): University of California; 2000.
26
Cornell WD, Cieplak P, Bayly CI, et al. A second generation force field for the simulation of proteins and nucleic acids.
J Am Chem Soc
1995
;
117
:
5179
–97.
27
Dewar MJS, Zoebisch EG, Healy EF, Stewart JJP. AM1: a new general purpose quantum mechanical model.
J Am Chem Soc
1985
;
107
:
3902
–9.
28
Bayly CI, Cieplak P, Cornell WD, Kollman PA. A well-behaved electrostatic potential based method using charge restraints for determining atom-centered charges: the RESP model.
J Phys Chem
1993
;
97
:
10269
–80.
29
Wang J, Cieplak P, Kollman PA. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules.
J Comput Chem
2001
;
21
:
1049
–74.
30
Wang J, Wang W, Kollman PA, Case DA. Antechamber, an accessory software package for molecular mechanical calculations.
J Comput Chem
2005
;
25
:
1157
–74.
31
Reyes CM, Kollman PA. Investigating the binding specificity of U1A-RNA by computational mutagenesis.
J Mol Biol
2000
;
295
:
1
–6.
32
Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water.
J Chem Phys
1983
;
79
:
926
–35.
33
Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath.
J Chem Phys
1984
;
81
:
3684
–90.
34
Ryckaert JP, Ciccotti G, Berendsen HJC. Numerical integration of the Cartesian equations of motion of a system with constraints; molecular dynamics of n-alkanes.
J Comp Phys
1977
;
23
:
327
–41.
35
Srinivasan J, Cheathem TE III, Cieplak P, Kollman PA, Case DA. Continuum solvent studies of the stability of DNA, RNA and phosphoramidate-DNA helices.
J Am Chem Soc
1998
;
120
:
9401
–9.
36
Gilson MK, Sharp KA, Honig B. Calculating the electrostatic potential of molecules in solution: method and error assessment.
J Comput Chem
1987
;
9
:
327
–35.
37
Chong LT, Duan Y, Wang L, Massova I, Kollman PA. Molecular dynamics simulation and free energy calculations applied to affinity maturation in antibody 48G7.
Proc Natl Acad Sci U S A
1990
;
96
:
14330
–5.
38
Sitkoff D, Sharp KA, Honig B. Accurate calculation of hydration free energies using macroscopic solvent models.
J Phys Chem
1994
;
98
:
1978
–88.
39
Sanner MF, Olson AJ, Spehner JC. Reduced surface: an efficient way to compute molecular surfaces.
Biopolymers
1996
;
38
:
305
–20.
40
Wilson EB, Decius JC, Cross PC. Molecular vibrations. New York: McGraw-Hill; 1995.
41
Ho W, Kukla MJ, Breslin HJ, et al. Synthesis and anti-HIV-1 activity of 4,5,6,7-tetrahydro-5-methylimidazo-[4,5,1-jk]-[1,4]benzodiazepine-2(1H)-one (TIBO) derivatives. 4.
J Med Chem
1995
;
38
:
794
–802.
42
Gambacorti-Passerini CB, Gunby RH, Piazza R, Galietta A, Rostagno R, Scapozza L. Molecular mechanisms of resistance to imatinib in Philadelphia-chromosome-positive leukaemias.
Lancet Oncol
2003
;
4
:
75
–85.