Computational analysis of binding between benzamide-based derivatives and Abl wt and T315I mutant kinases

Abstract

With the aim of discovering new mutation-resistant Abl inhibitors, especially the most difficult to overcome T315I mutant, the interactions between Abl kinases and a series of benzamide-based derivatives were studied using a combined method of a three-dimensional quantitative structure–activity relationship (3D-QSAR), molecular docking and molecular dynamics (MD) simulation. The results show that optimum CoMSIA (comparative molecular similarity index analysis) models have satisfactory internal and external predicted capacity, and the differences of structural features between Abl wt and T315I mutant inhibitors can be well pinpointed by the CoMSIA field plots. The detailed binding process and the in depth comparison of the binding modes of compounds with different activities against both kinases were validated by MD simulation. The binding free energies coincided well with the experimental bioactivities. The energy decomposition demonstrated that the van del Waals interaction is the major driving force for binding, and the hydrogen bond interactions with Met318 and Glu286 are also significant for the Abl potency increase. Finally, the most crucial residues impacting the strong interactions are also identified. These results can provide helpful reference for finding novel potential inhibitors.

---

Introduction

The Bcr-Abl is a non-receptor tyrosine kinase encoded by the Philadelphia chromosome, resulting from a reciprocal translocation between chromosome 9 and 22.¹ Bcr-Abl is involved in various cellular transductions and the progression of hematological malignancies. Many reports have shown that the over expression or deregulation of Bcr-Abl kinase leads to the hyperproliferation of leukemic cells and consequent pathology of the disease. It is the causative agent of chronic myeloid leukemia (CML) and emerging in more than 90% CML patient.² Therefore, Bcr-Abl has been regarded as an important target for the treatment of CML and promotes the finding of many effective inhibitors. After 2001, the first-generation Bcr-Abl kinase inhibitor imatinib has approved for the great clinical use and prominent efficacy for the chronic phase CML patients.³ However, the resistance to the imatinib had been discovered eventually with many patients because of the point mutations in the Bcr-Abl kinase domain.⁴ So far, more than 100 point mutations have been reported. Therefore, the development of new drugs capable of circumventing the imatinib-resistance is the major challenge in CML therapy. Research efforts have been devoted to discover the second-generation Abl inhibitors, including bafetinib,⁵ bosutinib,⁶ nilotinib,⁷ and dasatinib.⁸ But they all show low efficiency against the T315I mutation, which is the most difficult to overcome. The T315I mutation occurs with the reason that the gatekeeper Thr315 is replaced by Ile315, leading to the steric hindrance for imatinib to access the active domain, and the loss of a critical hydrogen bond between inhibitors and Thr315. T315I has been considered as the major frontier in the CML therapy. Delightedly, Bcr-Abl inhibitors GNF-7,⁹ MK-0457,¹⁰ PHA-739358,¹¹ AT9283,¹² and XL-228,¹³ etc., were reported to strongly inhibit not only Bcr-Abl T315I and other mutants with high bioactivity but also the wt kinase.

Recently, Li and co-workers¹⁴ had synthesized a novel series of 3-(1H-1,2,3-triazol-1-yl)benzamide derivatives and discovered that they had remarkable Abl wt and T315I inhibitory activities as well as antiproliferative activity toward some human cancer cell lines, revealing that these derivatives are good potential of developing as a new class of dual Abl wt and T315I inhibitors for the CML therapy. Today, some progresses in experimental researches have been made, but the relevant theoretical studies still be lacked, especially the binding mechanisms of these compounds toward the Abl wt and T315I kinases.

To promote the progress of the drug discovery, it is an important work to reveal the quantitative structural–activity relationship and the ligand–receptor mechanisms¹⁵ by computational studies. A 3D-QSAR CoMSIA method is performed to obtain the key structural factors influencing the inhibitory activity.^16–18 Docking study is also used to predict the optimized conformation of a ligand at the binding site of a receptor.¹⁹ Besides, molecular dynamics simulation is a useful approach to analyse the conformational changes of molecules and further notarize the structural features of ligand–receptor interactions.^20–22 Therefore, a combined method of 3D-QSAR, molecular docking, and MD simulation study can offer the detailed structural information of ligand-Abl wt and ligand-T315I interactions and more proposals for the development of novel inhibitors with improved activity.

For our work, we selected a novel series of 3-(1H-1,2,3-triazol-1-yl)benzamide derivatives as the dual Abl wt and T315I mutant inhibitors to carry out a computational study by using docking, 3D-QSAR and MD simulation methods. All the inhibitors docked into the active pocket of Abl, and the optimized orientation of each compound was selected. The best CoMSIA models were built and the key structural features affecting the bioactivities were verified. Moreover, the key structural features and the detailed interactions for six representative compounds interacting with both Abl kinases were further identified by MD simulations and free energy calculations. We hope that this work can offer useful information for the development of novel inhibitors toward wide-range of Bcr-Abl kinase, including the gatekeeper T315I mutant.^23,24

Materials and methods

Dataset

37 benzamide derivatives (Table S1†) exhibiting potent inhibitory bioactivities against the wt kinase and T315I mutant of Abl were gained from the literature,¹⁵ their general structural formulae were shown in Fig. 1. These derivatives were partitioned into a training set (27 compounds) and a test set (10 compounds) randomly considering the structural diversity and wide range of bioactivities in the data set. The IC₅₀ values were converted into pIC₅₀, and used as dependent variables in the 3D-QSAR analyses.

Fig. 1 General structural formula and numbering of benzamide-based derivatives (a) and template molecule (b, compound 24).

The structures of all compounds were established using the Sybyl 6.9 sketcher, and geometrically minimized by the Tripos force field with MMFF94 charges following a convergence criterion of 0.001 kcal (mol⁻¹ Å⁻¹). The optimal conformation was obtained for molecular docking.

Molecular docking

To obtain the appropriate binding conformations and orientations of the studied compounds linking with the Abl wt and T315I kinases, molecular docking was applied using the surflex-dock module in Sybyl 6.9 software.²⁵ The crystal structures of Abl wt and Abl T315I kinases (PDB entries: 3CS9 and 3IK3) for docking were gained from the Protein Data Bank. Before the docking process, the water molecules were eliminated, and the wt ligands were extracted. Then the polar hydrogen atoms were assigned and Kollman all atom charges were added to both proteins. The protomol was produced based on the ligand mode, and the protomol_threshold and protomol_bloat values were set to 0.5 and 0, respectively. By default, 20 binding conformations were generated for each ligand and the structure with the highest score was selected for the following studies.

Molecular modeling and alignment

The CoMSIA studies were carried out with SYBYL 6.9 software package. All parameters adopted to CoMSIA were set at the default values unless they are explained.

To develop the best CoMSIA models, two alignment methods were used.^26–28 Firstly, the ligand-based alignment was thought, the most active compound 24 was acted as a template to align the rest of compounds. The second method was receptor-based alignment, the optimal conformations generated from molecular docking were aligned automatically for 3D-QSAR analysis.

Generation of 3D-QSAR model

The CoMSIA descriptors were derived by using a 3D cube lattice with grid spacing of 2 Å beyond the aligned molecules in all directions. The steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields were calculated using an sp³ carbon probe atom of a +1.0 charge, a van der Waals radius of 1.0 Å, hydrophobicity of +1.0, and H-bond donor and acceptor properties of +1. A Gaussian function was used for evaluating the mutual distance between the probe atom and each molecule atom. In addition, the value of attenuation factor α was set to 0.3.^29,30

3D-QSAR models calculation and validation

The partial least-square (PLS) statistical method was exploited to establish the correlation between the CoMSIA fields and the pIC₅₀ values. Running in PLS, the leave-one-out (LOO) cross-validation method was used to compute the highest cross-validation correlation coefficient (q²) and the optimum number of components N. The non-cross-validation was also analysed, and the conventional correlation coefficient R², standard error of estimates (SEE), and F value were yielded. Finally, the bootstrapping analyses³¹ for 100 runs were computed to further evaluate the robustness and stability of the obtained models.

To appraise the predicted ability of the derived models, the pIC₅₀ values of 10 test set compounds belong to the external test set were predicted. The predictive correlation coefficient (R_pred²) was computed by following formula:

where PRESS stands for the sum of squared deviations between the actual and predictive activities of the test set compounds, and SD represents the sum of squared deviation between the experimental activities of the test set compounds and the mean activity of the training set compounds.

Molecular dynamics (MD) simulations

MD simulations were carried out with AMBER 9 software.^32,33 The docked complexes of 3CS9 and 3IK3 with compounds 12, 14, 22, 24, 33 and 36, respectively, were regarded as the starting structures. The electrostatic potentials (ESP) of each compound was computed at the B3LYP/6-31G(d) level in the Gaussian 09 program, and the partial atomic charges for these compounds were generated using the RESP protocol. The AMBER FF03 force field was exploited to describe the proteins and the general AMBER force field (gaff) was used for the ligands.³⁴ The counter ions of Na⁺ were used to neutralize the complexes, and then the whole system was immersed in TIP3P water models³⁵ extending a margin distance of 12 Å.

Prior to MD simulations, two-stages energy minimizations were took to relax the possible stress. In first stage, all parts of the systems was restrained with the constant of 2.0 kcal (mol⁻¹ Å⁻²) and the water molecules and Na⁺ ions were minimized by 2000 steps of steepest descent and 3000 steps of conjugated gradient. In second stage, the whole system was performed a complete minimization with 10 000 steps (5000 steps of steepest descent and 5000 steps of conjugate gradient minimization). Then, the relaxed systems were gradually heated from 0 to 300 K within 200 ps at a constant volume, and equilibrated at 300 K and 1 atm for 500 ps. Finally, 12.0 ns production MD simulation was carried out in a NPT ensemble in which the temperature was targeted at 300 K and the pressure was 1.0 atm. In the course of simulations, we selected the Particle Mesh Ewald (PME) method³⁶ to fix long-range electrostatic interactions with non-bonded cutoff of 8 Å. The SHAKE program³⁷ was employed to constrain all bonds including hydrogen atoms, and the time step was defined as 2 fs. We recorded the coordinate trajectories every 1 ps and confirmed the stability of the complexes by the root mean square deviations (RMSDs).

Binding free energy calculations

To prove the binding stability of the above complexes, the MM-PBSA program in AMBER 9 was used to calculate the binding free energy (see eqn (2)).³⁸

ΔG_bind = G_complex − (G_receptor + G_ligand) = ΔE_MM + ΔG_sol − TΔS = ΔE_vdw + ΔE_ele + ΔG_ele,sol + ΔG_nonpol,sol − TΔS (2)

where ΔE_MM stands for the gas-phase interaction energy, which is the sum of van der Waals energy ΔE_vdw and electrostatic energy ΔE_ele. ΔG_sol represents the solvation free energy determined by the polor solvation free energy ΔG_ele,sol and the nonpolar solvation free energy ΔG_nonpol,sol. The ΔG_ele,sol and ΔG_nonpol,sol were assessed by the Generalized Born (GB) approximation model and the solvent accessible surface area (SASA) that was computed with the Molsurf module in AMBER 9. Due to the expensive computational demand and lack of satisfactory result, the entropy contribution (TΔS) was ignored in this study. Finally, the binding free energies were decomposed to each residue by MM-GBSA method.

Results and discussion

Docking results

The Surflex-dock was applied on this data set, and all 37 studied compounds were docked into the binding pockets of Abl_wt and Abl_T315I, respectively. The reliability of the docking project was first validated by redocking two original ligands NIL and 0LI from their corresponding binding pockets of Abl_wt (PDB entry: 3CS9) and Abl_T315I (PDB entry: 3IK3). As a result, the conformations of docked NIL (or 0LI) and crystallographic NIL (or 0LI) were almost located in the same position, and their root-mean-square deviation (RMSD) values of 0.75 Å and 0.91 Å were obtained, demonstrating that this docking program was reliable. The binding mode of the most active compound 24 with two kinases were selected for more detailed analysis and provided in Fig. 2.

Fig. 2 The docking modes of compound 24 (yellow) in the binding sites of Abl_wt (a) and Abl_T315I (b). Hydrogen bonds are depicted as red dotted lines.

As displayed in Fig. 2, compound 24 was found to target the ATP binding site of DFG-out form of Abl_wt and Abl_T315I with similar binding modes. The pyrazolo[3,4-b]pyridinyl moiety at C₃ was situated at the hydrophobic pocket created by Phe317, Met318, Thr319 and Leu370/Thr253 for Abl_wt and Abl_T315I, and formed two hydrogen bonds with the NH and the C=O of Met318. The substituent R₁ was penetrated into the solvent accessible region at the entrance of the binding pocket surrounded by Lys285, Glu286, Ile360 and His361. And its terminal N_4′ atom could easily combine H⁺ ion in aqueous solvent and generate a positively charged anime group, implying a favorable hydrogen bond with the carbonyl oxygen atom of Lys285. This is compatible with the previous cocrystal structures studies on Abl inhibitors, verifying that the terminal amine group was significant for the Abl inhibitory activity.^39–41 Meanwhile, the methyl of substituent R₃ lies close to the side chain of Ile313 and Ile314, showing the favorable hydrophobic interaction. Many reports have demonstrated that the methyl group at C₇ was crucial for stabilizing conformation with Abl.^42,43 Finally, the trifluomethylphenyl group is deeply buried into a small hydrophobic pocket formed by Asp381 and Phe382 of DFG-motif. In addition, the linker amide, i.e., the O₁₃ and N₁₄H of compound 24 formed another two hydrogen bonds with the amino moiety of Asp381 and carbonyl oxygen of Glu286.

Through the above analyses, we can find that the binding of 24 to both Abl-wt and Abl-T315I mutant are very similar, but some differences also existed. First, the 1,2,3-triazol group sided close to Thr315 of Abl-wt and formed favorable electrostatic interaction, but no such interaction with Ile315 of Abl-T315I. Nevertherless, the hydrophobic interactions got more enhance in Abl-T315I-inhibitor than in Abl-wt-inhibitor, due to the mutation of Thr315 to Ile315. Moreover, compound 24 could form all four hydrogen bonds with two receptors, with the approximate distance of 1.8 Å, 3.1 Å, 3.0 Å, 2.0 Å for Abl_wt and 1.9 Å, 3.1 Å, 2.9 Å, 2.4 Å for Abl_T315I, indicating that the hydrogen bond interactions got stronger in Abl-wt-24 complex than in Abl-T315I-24 complex. Obviously, the enhancing of hydrogen bond and electrostatic interactions in Abl-wt-24 complex counteracting the reduction of hydrophobic interaction was an occasion for the higher activity of 24 for the wild type than T315I mutant.

3D-QSAR analyses

The CoMSIA analyses were used to generate the 3D-QSAR models of Abl_wt and Abl_T315I kinases, and the detailed statistical parameters were shown in Table 1.

Table 1 Statistical analysis results of the CoMSIA models^a

Statistical parameters	R^2	N	q^2	SEE	F	Rbs^2	SDbs	Rpred^2
a N is the optimal number of components, q^2 is the square of LOO cross-validation coefficient, R^2 is the square of non-cross-validation coefficient, SEE is the standard error of estimation, F is the F-test value, Rbs^2 is the mean R^2 of bootstrapping analysis (100 runs) SDbs is the mean standard deviation by bootstrapping analysis.
Abl_wt	0.998	6	0.853	0.043	2216.90	0.999	0.001	0.702
Abl_T315I	0.993	6	0.682	0.087	487.429	0.996	0.003	0.716

---

The Abl_wt CoMSIA (SEHA) model displayed the optimal results with high R² (0.998), F (2216.90), q² (0.853), and small SEE (0.043), while the corresponding values for the Abl-T315I CoMSIA (SHAD) model were 0.993, 487.429, 0.682 and 0.087, revealing that these models had reliable internal predictivity. To verify the statistical validity and robustness of the constructed models, bootstrapping analysis with 100 runs were performed with R_bs² of 0.993 and 0.996, and SD_bs of 0.003 and 0.003 for these two models, respectively. Moreover, the predictive correlation coefficients R_pred² values for two models were 0.702 and 0.716, demonstrating that these models exhibited satisfactory predictive ability. The steric and hydrophobic fields make the largest contributions (0.294 and 0.342 for Abl_wt model, 0.314 and 0.287 for Abl_T315I model, respectively) to the activity, suggesting their crucial roles in stabilizing the Abl-inhibitor complexes. Table S2† displayed the predicted pIC₅₀ values and the residual ones of compounds for two CoMSIA models. The data observed a good agreement between the experimental and predicted values, where most points were distributed around the regression line (Fig. 3), displaying that the CoMSIA models had good quality.

Fig. 3 The scatter plots of predicted activities vs. actual ones using the training set (black squares) and test set (red triangles) for CoMSIA models of Abl_wt (a) and Abl_T315I (b).

CoMSIA contour map analysis

Abl_wt series. To visualize the results of the best CoMSIA models, we regarded the most potent compound 24 as the template structure to construct the 3D color-coded contour maps (Fig. 4).

Fig. 4 CoMSIA contour maps of the highly active compound 24. Steric contour maps for Abl_wt (a) and Abl_T315I (b). Green contours indicate regions where bulky groups increase activity, yellow contours indicate regions where bulky groups decrease activity. Hydrophobic contour maps for Abl_wt (e) and Abl_T315I (f). Yellow contours indicate regions where hydrophobic substituents enhance activity; white contours indicate regions, where hydrophobic substituents decrease activity. Hydrogen-bond acceptor contour maps for Abl_wt (g), Abl_T315I (h). Magenta contours represent the H-bond acceptor favored regions, red contours indicate the H-bond acceptor disfavored regions. Electrostatic contour map for Abl_wt (c). Blue contours indicate regions where positive charges improve activity; red contours indicate regions where negative charges improve activity. Hydrogen-bond donor contour map for Abl_T315I (d). Cyan contours represent the H-bond donor favored regions, purple contours indicate the H-bond donor disfavored regions.

The CoMSIA steric contour maps regarding Abl-wt are displayed in Fig. 4a, c, e and g. There are a big green, a red and a while contours enclosing the terminal N_4′′ atom of ring-F, demonstrating that a bulky negative and hydrophilic group in this position is favorable. This is in agreement with the fact that compounds 20 and 21, with a terminal (S)- or (R)-dimethylamino group as well as 27 display higher activities than compound 19. The fact that compounds 25 and 26 (having 4-methyl-piperazinyl and 4-dimethylamino group at this area, respectively) are more potent than parent compound 22 is another example. The docking studies suggested that the 4-methyl-1,4-diazepanyl moiety pointed outside the solvent region, and hence, the bulky substituent R₁ with electronegative N atom should possess good physicochemical features to strengthen the binding to plasma. And we can easily find a small (see in Fig. 4a) and a big (see in Fig. 4e) yellow contours located adjacent to substituent R₂, which is plugged by the side chains of Val289 and Met290, meaning that introducing a relatively small and hydrophobic group to this position is beneficial to bioactivity. This can offer an explanation of the fact that the bioactivity of compound 29 with 1-H-1,2,4-triazol as substituent R₂ is higher than that of 30 with 3-methyl-1H-1,2,4-triazol as R₂. R₃ has a small (Fig. 4a) and a big yellow (Fig. 4e) and a blue contours inclose to it, which is blocked by the side chains of Ile313, Ile314 and Ile315. It reveals that including small and hydrophobic alkyl as R₃ can enhance the activity, which can be proved that compounds 17, 35, 36 and 37 have their sequence of activity of 17 > 35 > 36 > 37, bearing the corresponding R₃ of methyl, ethyl, cyclopropyl and isopropyl, respectively. In addition, a large yellow, two green (Fig. 4a), two red and two white (Fig. 4c) contours embed in the substituent Het, showing that negatively charged and medium-size groups at this position are beneficial, in docking, Het is surrounding by Glu316, Phe317, Met318 and Leu370, and too overlarge groups may bring the steric hindrance. Meanwhile, Fig. 4g shows that a magenta and a red contours are found near the H atom connecting with N_7′ and N_3′ atoms, suggesting that hydrogen bond acceptor and donor groups in these two site, respectively, will improve the activity. Consistent with our docking results, the N_7′H group and N_3′ atom can form two hydrogen bonds with Met318. Compounds 13 lacking H atom linking to N_7′ as a hydrogen bond acceptor, as well as compounds 15 and 16 lacking N_3′ as a hydrogen bond donor, showed a respectively 38-, 503-, and 30-fold potency loss than 17. That can illustrate why the potencies of compounds 6, 7, 9, 10 and 11 markedly decreased with the increasing size of substituent Het, because of the large groups affecting the potential hydrogen bond interactions with Met318. Finally, a small red, a white and a magenta contours are found in the vicinity of the O₁₃ position, indicating that the electronegative atom at this site is favored. It is displayed that the O₁₃ atom can form a strong hydrogen bond with Asp381 in docking.

Abl_T315I series. The CoMSIA (SHAD) contour maps for Abl-T315I are shown in Fig. 4b, d, f and h, in which the electrostatic field of Abl-wt is replaced by hydrogen bond donor field. A careful comparison of the contours of two models reveals that the steric, hydrophobic and hydrogen bond acceptor contours are very similar each other, but they also have some differences. For example, as for R₁, the additional sterically unfavorable and hydrophobic favored contours appear round R₁, apart from the big green and white contours of Abl-wt. This suggests that the medium-size groups with certain hydrophibic-lipophilic atoms as R₁ are welcome for activity better. Another additional magenta contour locates at the N_4′′ position of the red contour seen for Abl-wt electrostatic field, demonstrating the hydrogen bond acceptor at N_4′′ position has a significant impact on the Abl-T315I activities. Obviously, the terminal of R₁ is exposed to the solvent-accessible area, the N_4′′ atom can easily associates H⁺ ion in aqueous solvent and bears a positively charged amino group, and involving in a favored hydrogen bond interaction with the carbony atom of Lys285. This can explain well that compound 24 with 4-methyl-1,4-diazepanyl as R₁ exhibits higher potency than compounds 17, 25 and 26 with 4-methylpiperazinyl, 4-(1-(4-methyl)piperazin)yl piperidinyl and 4-dimethylamino-1-piperidinyl as R₁, respectively. The fact that the activities of compounds 19, 22 and 23 against Abl-T315I, which do not have the terminal N_4′′ atom, were obviously decreased also proves this conclusion. So this may be a key function to impact the kinase selectivity for Abl-wt and Abl-T315I.

In the Abl-T315I plot, the contours surrounding substituents Het and R₃ were observed to be almost the same with those in the Abl-wt model, but their sizes are reduced. One notable difference is the green contour near pyridyl ring in Fig. 4a, replacing by the yellow one close to the pyrazolo ring (Fig. 4b), and warning against keeping over large groups. This is in accord with the Abl-wt inhibition and can once again explain why the potencies of compounds 5–11 decreased markedly with the increase size of the N-substituted groups. Meanwhile, the docking studies displayed that the pyrazolo ring of Het may form another π–π stacking with Tyr253, except for the π–π stacking with Phe317, which is the same as that with Abl-wt receptor. This indicated that the binding of pyrazolo ring to Abl-T315I receptor was stronger than that of other groups. This is supported by the fact that compounds 14, 17 and 18, with 1-H-pyrrolo[2,3-b]pyridinyl, 1-H-pyrazolo[3,4-b]phenyl and 1-H-pyrazolo[3,4-b]pyridinyl group as Het, exhibited 20-, 22-, and 10-times increase in Abl-T315I potencies than compound 6 with 2-N-methyl amino pyrimidinyl as Het. But they showed similar Abl-wt activities to compound 6.

In addition, the big purple and cyan contours near N_7′H group and N_3′ atom in Fig. 4d, revealed that hydrogen bond donor substituents are unfavorable and favorable for the activity at these positions, respectively. This observation is in conformity with the forgoing analyses of hydrogen bond acceptor contours (Fig. 4g and h).

MD simulations, MM/GBSA calculations and free energy decomposition

To further verify the binding interactions of the inhibitors and two receptors, 12 ns MD simulations of 12, 14, 22, 24, 33 and 36 with Abl-wt and Abl-T315I complex, respectively, were performed. Six representative compounds was falled into three pairs (12–14, 22–24, 33–36), with the aim of expounding the effects of the substituents R₁, R₃ and Het on inhibitory activity.

To clarity the dynamic stabilities of these systems, the root-mean-square displacements (RMSDs) of the backbone atoms against the starting structures were analyzed. As displayed in Fig. 5, the average RMSD values of Abl-wt/inhibitor complexes (1.45–1.84 Å) showed much lower than those of Abl-T315I/inhibitor complexes (2.25–2.71 Å), and the RMSD fluctuations are very small, revealing that the compounds might make stronger interactions with the Abl-wt kinase in the ATP-binding site than Abl-T315I mutant kinase. It is in agreement with the fact that compounds 12, 14, 22, 24, 33 and 36 have higher inhibitory potency against the wild type than the T315I mutant. In addition, we can obviously discover that compound 12 have the highest RMSD values for two receptors. This may be an explanation for compound 12 with the lowest inhibitory potency against wild type and T315I mutant of Abl.

Fig. 5 The RMSDs of the backbone atoms of the wt–inhibitor complexes (12, 14, 22, 24, 33) (a) and the T315I–inhibitor complexes (12, 14, 22, 24, 33) (b).

Fig. 6 showed the analyses of root-mean-square fluctuation (RMSF) via the residue number of four systems (wt-12, wt-14, T315I-12 and T315I-14). As shown in Fig. 6a, it is obviously seen that residues Ala269, Val289, Met290, Ile293, Val299 and Leu370 in wt-12 and wt-14 have steady quality, because these six residues can make strong van der Waals contacts with two compounds (Fig. 7). Nevertheless, residues Met318, Phe359 and Glu292 located in the active site exhibited larger vibrations for wt-14 than wt-12, especially for Met318, because Met318 can form a stable hydrogen bond with compound 14 but not with 12. In Fig. 6b, we can also note that the residues Ala269, Val289, Val299, Ile360, Leu370 and Ala380 are rather immobile, due to their strong hydrophobic interaction with both systems, moreover, Ile360 can interacts compound 14 with a weak hydrogen bond. However, some apparent differences can also be discovered that residues Lys245, Glu292 and Ala407 for T315I-14 are more flexible than those for T315I-12. Therefore, we can infer that compound 14 may form stronger interactions with two receptors than compound 12, which is in agreement with the experimental activities.

Fig. 6 The RMSFs of backbone atoms versus residue number of the wt-12 and wt-14 complexes (a), and T315I-12 and T315I-14 complexes (b).

Fig. 7 The comparison of the averaged structures for the wt-12 and wt-14 (a) and T315I-12 and T315I-14 (b) complexes (carbon atoms are with cyan and green, respectively); energy difference of each residue to the binding of wt-12 and wt-14 (c) and T315I-12 and T315I-14 (d); the contributions of the individual energy terms for the key residues ((e) and (f), the red and green columns represent ΔG_vdw and ΔG_nonpol,sol, respectively, and the blue and cyan columns represent ΔG_ele and ΔG_ele,sol, respectively).

Besides, we can further analyse the binding conformation by utilizing the hydrogen bond interactions from MD simulations for the above twelve systems (Table 2). From Table 2, we can observe that for each inhibitor pair (12–14, 22–24, 33–36), compounds 14, 24, 33 with higher bioactivity can form more hydrogen bonds with two receptors than 12, 22 and 36 with lower bioactivity. Moreover, some hydrogen bonds existing in the docking model, such as the Met318 NH⋯N₃ in wt-24 and Met318 C=O⋯N₇ hydrogen bond in T315I-24 system, have disappeared, due to the slight movement of the conformations of compound 24.

Table 2 Hydrogen bond analysis from the results of MD simulations

Complex	Donor	Acceptor	Distance/Å	Angle/(°)	Occupancy/%
wt-12	Ligand@O13	Asp381@N-H	2.893	172.23	89.1
wt-14	Ligand@O13	Asp381@N-H	2.778	177.76	99.5
	MET318@O	Ligand@N7′H	2.911	156.36	76.2
	Ligand@N3′	MET318@NH	3.022	166.12	65.3
T315I-12	Ligand@O13	Asp381@N-H	2.733	170.15	92.4
T315I-14	Ligand@O13	Asp381@N-H	2.647	178.90	99.0
	MET318@O	Ligand@N7′H	2.989	173.77	81.3
	Ligand@N3′	MET318@NH	3.242	147.03	52.1
wt-22	Ligand@O13	Asp381@N-H	2.931	160.21	87.6
	MET318@O	Ligand@N7′H	3.378	169.23	67.2
wt-24	Ligand@O13	Asp381@N-H	2.831	164.09	99.6
	MET318@O	Ligand@N7′H	3.133	163.44	78.5
T315I-22	Ligand@O13	Asp381@N-H	3.091	177.04	97.6
	Ligand@N3′	MET318@NH	3.462	140.79	55.8
T315I-24	Ligand@O13	Asp381@N-H	2.862	175.51	99.3
	Ligand@N3′	MET318@NH	3.226	161.20	69.1
wt-33	Ligand@O13	Asp381@N-H	2.776	162.68	99.7
	MET318@O	Ligand@N7′H	3.190	145.03	62.9
	Glu286@O	Ligand@N14H	3.397	131.44	53.2
wt-36	Ligand@O13	Asp381@N-H	3.073	157.09	96.9
T315I-33	Ligand@O13	Asp381@N-H	2.882	172.16	97.4
	MET318@O	Ligand@N7′H	3.142	162.76	71.6
	Glu286@O	Ligand@N14H	2.782	169.04	87.2
T315I-36	Ligand@O13	Asp381@N-H	3.062	161.03	86.3

---

The binding free energies of the twelve complexes were computed using MM-PBSA method and shown in Table 3. It is observed that their binding affinities are in line with the experimental bioactivities. To gain further insight into the ligand–receptor interactions, we use the free energy decomposition to establish the key residue–ligand interaction spectra. Accordingly, three compound pairs (12–14, 22–24, 33–36) for wild type and T315I mutant of ABL were systematically compared, depending on the binding free energies and the key residue–inhibitor interaction spectra.

Table 3 Binding free energy (kcal mol⁻¹) of the twelve systems

System	Polar contributions		Nopolar contributions		ΔGbind
	ΔEele	ΔGele,sol	ΔEvdw	ΔGnonpol,sol
wt-12	−14.24	38.79	−68.80	−9.86	−54.11
wt-14	−27.38	40.80	−72.68	−6.26	−65.52
wt-22	−29.79	44.53	−74.27	−5.89	−65.41
wt-24	−23.95	41.15	−80.48	−6.46	−69.74
wt-33	−25.47	37.98	−73.22	−6.10	−66.81
wt-36	−24.80	44.05	−76.37	−6.32	−63.44
T315I-12	−13.61	37.86	−70.11	−6.51	−52.37
T315I-14	−28.74	39.88	−70.20	−6.08	−65.13
T315I-22	−13.46	33.80	−71.82	−6.36	−57.84
T315I-24	−31.20	46.03	−75.16	−8.62	−68.95
T315I-33	−26.29	45.42	−76.68	−6.31	−63.86
T315I-36	−17.06	44.05	−74.54	−6.27	−53.82

---

Comparing 12 and 14, it shows that their substituents Het are different. When the 4-pyrimidinepiperidinyl group in 12 was replaced by the pyrrolo[2,3-b]pyridinyl group, the resulting compound 14 has dramatic better potency than 12, indicating that a hydrogen bond donor group at N₇′-position can cause significant potency increase. This is also compatible with the results of binding free energy that the ΔG_bind values for wt-14 (−65.52 kcal mol⁻¹) and T315I-14 (−65.13 kcal mol⁻¹) are significantly higher than wt-12 (−54.11 kcal mol⁻¹) and T315I-12 (−56.37 kcal mol⁻¹), which demonstrate that 14 binds more tightly with two receptors than 12. In Fig. 7a and b, compound 14 can form two hydrogen bonds Met318 for wild type and T315I mutant, respectively, and those for wt-12 and T315I-12 are disappeared, which can be proved by the binding free energy calculations, because the polar contributions (ΔG_ele + ΔG_ele,sol) of wt-14 (13.42 kcal mol⁻¹) and T315I-14 (11.14 kcal mol⁻¹) are apparently stronger than those of wt-12 (24.55 kcal mol⁻¹) and T315I-14 (24.25 kcal mol⁻¹). Meanwhile, the energy differences of key residues to wt-14 and wt-12 systems are shown in Fig. 7e, where the polar interactions of Met318, Phe317 and Lys271 with 14 are largely higher than those of 12, but Tyr253 and Ile293 locate near compound 12 and make stronger non-polar interactions with 12 than with 14. However, the strengthening hydrophobic interactions of compound 12 with Tyr253 and Ile293 are not sufficient to compensate the weakening interactions with other residues. As for the T315I mutant, we can observe that the polar residues Met318 and Phe317 and the non-polar residue Ile360 contribute obviously more to 14 than 12 in Fig. 7f, but Leu248 connects tighter with 12 than 14 by the non-polar interaction, due to its smaller distance to 12. Nevertheless, Table 3 shows that the non-polar contribution (ΔG_vdw + ΔG_nonpol,sol) of 14 (−76.28 kcal mol⁻¹) is approximatively equal to that of 12 (−76.62 kcal mol⁻¹), which means that the polar interaction plays the key role in the potency difference between T315I-12 and T315I-14. From the above, we can conclude that the polar interactions are a key role in the stronger binding affinities of 14 over 12 in both wild type and T315I mutant of Abl. We also reveal that introducing a hydrogen bond acceptor and/or donor atom at 3′- and/or 7′-position of molecule should be prioritize in designing the substituent Het.

The differences between compounds 22 and 24 are R₁ substituent, and the piperidine ring in 22 is replaced by 4-methyl-1,4-diazepane ring in 24. From Table 3, it is obvious to find that the non-polar energy makes prominent contribution to the binding, where the value (ΔG_vdw + ΔG_nonpol,sol) of wt-24 is −6.78 kcal mol⁻¹ over wt-22, and that of T315I-24 is −5.6 kcal mol⁻¹ over T315I-22. For wt-22 and wt-24, Fig. 8c and e show that three important residues Ile293, Ala380 and Phe382 exhibit large difference further demonstrated the is major contribution of the non-polar terms in two compounds. With regard to T315I-22 and 24 systems (Fig. 8d and f), Ile293, Phe359 and Ile360 apparently connect compound 24 with stronger interactions than 22, and the non-polar terms are also the chief contributions. Here Ile293 interacting directly with R₁ is the key residue to influence the different bioactivity of compound 24 and 22 in complex with wild type and T315I mutant of Abl. Meanwhile, as shown in Fig. 8f, we can also observe that the polar term of lle360 for T315I-24 is distinctly higher than that of T315I-22, because the C=O motif of Ile360 can form a weak hydrogen bond of 3.5 Å with N_4′′ of compound 24. This is in line with the energy individual terms in Table 3, which indicates that the polar contribution of T315I-24 (14.83 kcal mol⁻¹) is stronger than that of T315I-22 (20.34 kcal mol⁻¹). Therefore, introducing the relatively large and hydrophobic R₁ substituent possessing the terminal anime group may be favorable for the potency of both kinases.

Fig. 8 The comparison of the averaged structures for the wt-22 and wt-24 (a) and T315I-22 and T315I-24 (b) complexes (carbon atoms are with cyan and green, respectively); energy difference of each residue to the binding of wt-22 and wt-24 (c) and T315I-22 and T315I-24 (d); the contributions of the individual energy terms for the key residues ((e) and (f), the red and green columns represent ΔG_vdw and ΔG_nonpol,sol, respectively, and the blue and cyan columns represent ΔG_ele + ΔG_ele,sol, respectively).

At the R₃ position, the Cl atom in 33 is replaced with the cyclopropyl group in 36, the resulting compound 36 displayed markedly potency loss in 33 with higher inhibition against the wild type and T315I mutant of ABL than 33, which coincided with the results of binding free energy that the ΔG_bind values of wt-33 (−66.81 kcal mol⁻¹) and T315I-33 (−63.26 kcal mol⁻¹) are substantial higher than those of wt-36 (−63.44 kcal mol⁻¹) and T315I-36 (−53.82 kcal mol⁻¹). For two wild type systems, the polar terms (ΔG_ele + ΔG_ele,sol) of wt-33 (12.51 kcal mol⁻¹) are much higher than that of wt-36 (19.25 kcal mol⁻¹) in Table 3. As shown in Fig. 9e, the residues with mainly responsible for the differences of energy are Lys271, Glu286 and Thr315, which are all polar residues, especially for Glu286, whose energy difference between 33 and 36 is up to −2.66 kcal mol⁻¹. Compared with compound 36 in Fig. 9a, compound 33 can form an additional hydrogen bond with Glu286, resulting in the much higher polar interaction over compound 36. Moreover, the NH group of Thr315 and C=O group of Lys271 can directly interact with the strong electronegative chlorine atom in R₃ position. Compound 33 forms stronger polar interactions with these residues than compound 36. As for two T315I mutant systems, we can also discover that the polar terms of T315I-33 (19.13 kcal mol⁻¹) display visibly better than that of T315I-36 (26.99 kcal mol⁻¹) (Table 3). As it can be seen in Fig. 9f, the important residues with favorable large energy difference between 33 and 36 are Ala269, Glu286 and Met318. Obviously, the energy difference of these residues is primarily contributed by the polar energy terms. Compared with compound 36 in Fig. 9b, it is noted that compound 33 can form two additional hydrogen bonds with Glu286 and Met318, suggesting that 33 have higher polar interaction with T315I mutant of Abl than 36. Furthermore, the triazol ring in 33 may form a π–π stacking with Phe382, leading to the higher non-polar interaction of Phe382 with compound 33 than 36, because of the different conformations between 33 and 36. From Fig. 9b, we can note that the triazol ring of 36 is rotated by about 90° compared with 33, due to the blocking of Ile315 to the cyclopropyl ring of 36, resulting in the weak interaction with Phe382. In summary, selecting small and electronegative atom or group as R₃ may enhance the bioactivities in the wild type and T315I mutant of Abl. We should take more considerations in the interactions with Glu286 when designing novel substituent R₃.

Fig. 9 The comparison of the averaged structures for the wt-33 and wt-36 (a) and T315I-33 and T315I-36 (b) complexes (carbon atoms are with cyan and green, respectively); energy difference of each residue to the binding of wt-33 and wt-36 (c) and T315I-33 and T315I-36 (d); the contributions of the individual energy terms for the key residues ((e) and (f) the red and green columns represent ΔG_vdw and ΔG_nonpol,sol, respectively, and the blue and cyan columns represent ΔG_ele + ΔG_ele,sol, respectively).

Conclusions

In this work, we analyzed the interactions between the wild type and T315I mutant of Abl and a series of benzamide-based derivatives by utilizing an integrated computational method. The generated 3D-QSAR models revealed the important structural features affecting the inhibitory activities against both kinases. The docking studies shown that inhibitors bound to two receptors with the similar modes, and its reliability was confirmed by MD simulation. Decomposition analysis of the free energy further demonstrated that the hydrogen bond interactions with Met318 and Glu286 were the key element for distinguishing the different binding and activity of compounds. The comparison between the interactions of three pairs of inhibitors (12–14, 22–24, 33–36) with both kinases might shed light on the design of novel inhibitors of wild type and T315I mutant of Abl.

Acknowledgements

The authors gratefully acknowledge supports of this research by the Science and Technology planning Project of Guangzhou (No. 2013J4100071). We also heartily thank the computation environment support by the Department of Biochemistry, College of Life Sciences, SunYat-Sen University for SYBYL 6.9.

Notes and references

1. O. Hantschel and G. Superti-Furga, Nat. Rev. Mol. Cell Biol., 2004, 5, 33–44.
2. J. M. Goldman and J. V. Melo, N. Engl. J. Med., 2003, 349, 1451–1464.
3. M. Deininger, E. Buchdunger and B. J. Druker, Blood, 2005, 105, 2640–2653.
4. E. Jabbour and H. Kantarjian, Am. J. Hematol., 2012, 87, 1037–1045.
5. S. Schenone, C. Brullo and M. Botta, Curr. Med. Chem., 2010, 17, 1220–1245.
6. T. O'Hare, C. A. Eide and M. W. Deininger, Expert Opin. Invest. Drugs, 2008, 17, 865–878.
7. H. Kantarjian, F. Giles, L. Wunderle, K. Bhalla, S. O'Brien, B. Wassmann, C. Tanaka, P. Manley, P. Rae and W. Mietlowski, N. Engl. J. Med., 2006, 354, 2542–2551.
8. N. P. Shah, C. Tran, F. Y. Lee, P. Chen, D. Norris and C. L. Sawyers, Science, 2004, 305, 399–401.
9. X. Y. Lu, Q. Cai and K. Ding, Curr. Med. Chem., 2011, 18, 2146–2157.
10. F. J. Giles, J. Cortes, D. Jones, D. Bergstrom, H. Kantarjian and S. T. Freedman, Blood, 2007, 109, 500–502.
11. P. Carpinelli, R. Ceruti, M. L. Giorgini, P. Cappella, L. Gianellini, V. Croci, A. Degrassi, G. Texido and M. Rocchetti, Mol. Cancer Ther., 2007, 6, 3158–3168.
12. S. Howard, V. Berdini, J. A. Boulstridge, M. G. Carr, D. M. Cross, J. Curry and L. A. Devine, J. Med. Chem., 2009, 52, 379–388.
13. J. Cortes, R. Paquette, M. Talpaz, J. Pinilla, E. Asatiani, M. Wetzler, J. H. Lipton, C. Kasap, L. A. Bui, D. O. Clary and N. Shah, Blood, 2008, 112, 3232.
14. Y. Li, M. Shen, Z. Zhang, J. Luo, X. Pan, X. Lu, H. Long, D. Wen, F. Zhang and F. Leng, J. Med. Chem., 2012, 55, 10033–10046.
15. F. Ntie-Kang, J. N. Yong, L. C. O. Owono, W. Sippl and E. Megnassan, Curr. Med. Chem., 2014, 21, 3466–3477.
16. A. Ashek and S. J. Cho, Bioorg. Med. Chem., 2006, 14, 1474–1482.
17. H. Assefa, S. Kamath and J. K. Buolamwini, J. Comput.-Aided Mol. Des., 2003, 17, 475–493.
18. F. Cheng, J. Shen, X. Luo, W. Zhu, J. Gu, R. Ji, H. Jiang and K. Chen, Bioorg. Med. Chem., 2002, 10, 2883–2891.
19. S. Durdagi, T. Mavromoustakos and M. G. Papadopoulos, Bioorg. Med. Chem. Lett., 2008, 18, 6283–6289.
20. S. Kaur, A. V. Shivange and N. Roy, Mol. Diversity, 2010, 14, 169–178.
21. C. King, H. Diaz, D. Barnard, D. Barda, D. Clawson, W. Blosser, K. Cox, S. Guo and M. Marshall, Invest. New Drugs, 2014, 32, 213–226.
22. S. J. Ma, G. H. Zeng, D. Q. Fang, J. P. Wang, W. J. Wu, W. G. Xie, S. P. Tan and K. C. Zheng, Mol. BioSyst., 2015, 11, 394–406.
23. X. Li, L. Ye, X. Wang, X. Wang, H. Liu, X. Qian, Y. Zhu and H. Yu, Sci. Total Environ., 2012, 441, 230–238.
24. Y. Yang, J. Qin, H. Liu and X. Yao, J. Chem. Inf. Model., 2011, 51, 680–692.
25. A. N. Jain, J. Med. Chem., 2003, 46, 499–511.
26. S. J. Cho and A. Tropsha, J. Med. Chem., 1995, 38, 1060–1066.
27. B. Hoffman, S. J. Cho, W. Zheng, S. Wyrick, D. E. Nichols, R. B. Mailman and A. Tropsha, J. Med. Chem., 1999, 42, 3217–3226.
28. A. Tropsha and S. J. Cho, 3D QSAR Drug Des., 1998, 3, 57–69.
29. R. D. Cramer, J. D. Bunce, D. E. Patterson and I. E. Frank, Quant. Struct.-Act. Relat., 1988, 7, 18–25.
30. R. D. Cramer, D. E. Patterson and J. D. Bunce, J. Am. Chem. Soc., 1988, 110, 5959–5967.
31. Y. Yang, Y. Shen, H. Liu and X. Yao, J. Chem. Inf. Model., 2011, 51, 3235–3246.
32. D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang and R. J. Woods, J. Comput. Chem., 2005, 26, 1668–1688.
33. D. A. Case, T. Darden, T. E. Cheatham III, C. Simmerling, J. Wang, R. E. Duke, R. Luo, K. M. Merz, D. A. Pearlman and M. Crowley, Amber 9, University of California, San Francisco, 2006.
34. M. J. Hsieh and R. Luo, Proteins, 2004, 56, 475–486.
35. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
36. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
37. J. P. Ryckaert, G. Ciccotti and H. J. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
38. T. Laitinen, J. A. Kankare and M. Peräkylä, Proteins, 2004, 55, 34–43.
39. T. Horio, T. Hamasaki, T. Inoue, T. Wakayama, S. Itou, H. Naito, T. Asaki, H. Hayase and T. Niwa, Bioorg. Med. Chem. Lett., 2007, 17, 2712–2717.
40. T. Schindler, W. Bornmann, P. Pellicena, W. T. Miller, B. Clarkson and J. Kuriyan, Science, 2000, 289, 1938–1942.
41. T. Zhou, L. Commodore, W. S. Huang, Y. Wang, M. Thomas, J. Keats, Q. Xu, V. M. Rivera, W. C. Shakespeare and T. Clackson, Chem. Biol. Drug Des., 2011, 77, 1–11.
42. W. S. Huang, C. A. Metcalf, R. Sundaramoorthi, Y. Wang, D. Zou, R. M. Thomas, X. Zhu, L. Cai, D. Wen and S. Liu, J. Med. Chem., 2010, 53, 4701–4719.
43. X. Ren, L. Duan, Q. He, Z. Zhang, Y. Zhou, D. Wu, J. Pan, D. Pei and K. Ding, ACS Med. Chem. Lett., 2010, 1, 454–459.

---

Footnote

† Electronic supplementary information (ESI) available: Tables S1 and S2. See DOI: 10.1039/c6ra19494j

---