
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
1 Molecular Simulation Engineering Laboratory, Department of Chemical Engineering, University of Trieste, Trieste, Italy and 2 Experimental Molecular Pathology, Department of Pathology, Istituto Nazionale per lo Studio e la Cura dei Tumori, Milan, Italy
Requests for reprints: Sabrina Pricl, Molecular Simulation Engineering Laboratory, Department of Chemical Engineering, University of Trieste, Piazzale Europa 1, 34127 Trieste, Italy. Phone: 39-40-5583750; Fax: 39-40-569823. E-mail: sabrina.pricl{at}dicamp.units.it
| Abstract |
|---|
|
|
|---|
| Introduction |
|---|
|
|
|---|
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.
| Materials and Methods |
|---|
|
|
|---|
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:
![]() | (A) |
![]() | (B) |
EMM denotes the sum of molecular mechanics energies, and can be further split into contributions from electrostatic (EEL) and van der Waals (EvdW) energies:
![]() | (C) |
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.
| Results and Discussion |
|---|
|
|
|---|
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.
|
|
|
|
-
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 612 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).
| Conclusions |
|---|
|
|
|---|
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.
| Footnotes |
|---|
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.
Received 4/ 4/05; revised 5/19/05; accepted 6/14/05.
| References |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
E. J. Kennedy, G. Ghosh, and L. Pillus Identification of Functionally Distinct Regions That Mediate Biological Activity of the Protein Kinase A Homolog Tpk2 J. Biol. Chem., January 11, 2008; 283(2): 1084 - 1093. [Abstract] [Full Text] [PDF] |
||||
![]() |
F. E. Nicolini, S. Hayette, S. Corm, E. Bachy, D. Bories, M. Tulliez, F. Guilhot, L. Legros, F. Maloisel, J.-J. Kiladjian, et al. Clinical outcome of 27 imatinib mesylate-resistant chronic myelogenous leukemia patients harboring a T315I BCR-ABL mutation Haematologica, September 1, 2007; 92(9): 1238 - 1241. [Abstract] [Full Text] [PDF] |
||||
![]() |
G. Martinelli, S. Soverini, I. Iacobucci, and M. Baccarani MK-0457: a light at the end of the tunnel? Blood, January 15, 2007; 109(2): 396 - 397. [Full Text] [PDF] |
||||
![]() |
E. Banfi, G. Scialino, D. Zampieri, M. G. Mamolo, L. Vio, M. Ferrone, M. Fermeglia, M. S. Paneni, and S. Pricl Antifungal and antimycobacterial activity of new imidazole and triazole derivatives. A combined experimental and computational approach J. Antimicrob. Chemother., July 1, 2006; 58(1): 76 - 84. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. Ferrone, F. Perrone, E. Tamborini, M. S. Paneni, M. Fermeglia, S. Suardi, E. Pastore, D. Delia, M. A. Pierotti, S. Pricl, et al. Functional analysis and molecular modeling show a preserved wild-type activity of p53C238Y. Mol. Cancer Ther., June 1, 2006; 5(6): 1467 - 1473. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Cancer Research | Clinical Cancer Research |
| Cancer Epidemiology Biomarkers & Prevention | Molecular Cancer Therapeutics |
| Molecular Cancer Research | Cancer Prevention Research |
| Cancer Prevention Journals Portal | Cancer Reviews Online |
| Annual Meeting Education Book | Meeting Abstracts Online |