Introduction

Chronic myeloid leukemia (CML) arises due to the chromosomal aberration in which reciprocal translocation of the Abelson gene on chromosome 9 to break-point cluster region gene on chromosome 22 results in the creation of fusion oncogene, bcr-abl (9; 22) (q34; q11)1. The product of bcr-abl oncogene, BCR-ABL protein is a constitutively active tyrosine kinase that drives the disease CML through phosphorylation of many downstream effector molecules, including Grb2, RAK, ROS, PI3K, JNK, STAT5, AKT and Myc which consequently promoting cell proliferation2,3. Imatinib is the first FDA approved tyrosine kinase inhibitor (TKI) that has shown potent inhibitory effect on the progression of CML4. Despite imatinib's clinical success, resistance due to mutations and side effects are still a limitation of this drug5. To overcome resistance, second generation TKI's inhibitors, nilotinib6 and dasatinib7 were developed. However, neither compound effectively inhibits T315I mutant BCR-ABL, which constitutes 20% of all BCR-ABL mutations8,9,10.

Ponatinib is the only available drug that is designed to overcome T315I gatekeeper mutation11 and is efficient in inhibiting the mutant BCR-ABL12. Ponatinib exhibits triple carbon-carbon (ethynyl linkage) bond between the methyl phenyl and purine groups13,14. It holds the Isoleucine side chain without steric interference and without any loss of hydrogen bond (H_bond)11. It also showed similar binding pattern (DFG-out) as imatinib and nilotinib which interacts with Met318, Asp381 and the side chain of Glu286. In addition, the drug makes van der Waals contacts with Tyr253 and Phe382 because of squeezed conformation of the P-loop and DFG-out mode of activation loop respectively14. Despite potent inhibition capability, ponatinib shows severe side-effects like blood clotting in cardiac valves, chambers and cerebral vessels consequently leading to adverse conditions like myocardial infarction, cardiac stroke and cerebral infraction15. Therefore, a broad spectrum drug capable of inhibiting both wild-type and mutant BCR-ABL with fewer side effects is currently in demand. The availability of 3D-structure of the target protein and the structural details of ponatinib and T315I mutant BCR-ABL protein complex renders an opportunity to identify the most active drug candidate that can efficiently block the catalytic activity of BCR-ABL.

Using a virtual screening approach, we have screened and identified potent drug-like compounds, efficient against both wild-type and mutant BCR-ABL, from a large library of small molecules. Drug molecules exhibiting docking score higher than ponatinib were manually explored in molecular visualization tools in order to analyze their binding patterns. Highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) were calculated using DFT analysis to assess their chemical reactivity. Following which, protein-ligand complexes were incorporated in GROMACS for MD simulation to evaluate their structural stability throughout the trajectory period of simulation. Binding free energy for these complexes were then investigated to cross-check their binding affinity. Results of these studies uncovered seven lead molecules designated as DB07107, DB06977, ST013616, DB04200, ST007180 ST019342 and DB01172. Among these, DB07107, DB06977, ST013616, DB04200 and ST007180 were found to be more effective in blocking drug-resistant T315I mutant than the wild-type BCR-ABL. Interestingly ST019342 and DB01172 were effective only in mutant BCR-ABL.

Methods

Protein selection and preparation

The crystallographic co-ordinates for wild-type BCR-ABL (PDB ID: 3OXZ)16 and mutant T315I (PDB ID: 3QRJ)17 were retrieved from the Protein Data Bank (PDB). Prior to docking, protein structures were prepared by removing water molecules using Schrödinger software. Following which, bond orders were assigned and hydrogen atoms were added to the crystal structures. Finally, a restrained minimization of the protein structure was performed using the default constraint of 0.30Å RMSD and the OPLS-2005 force field18.

Ligand preparation and grid generation

A total of 36,481 small molecules were retrieved in Structure Data Format (SDF) from various small molecule libraries which includes Ligand.info: Small-molecule Meta-Database19,20 (29,090), Drug Bank21,22 (6,825) and PubChem23 database (566). These small molecules were prepared using the LigPrep wizard of Schrödinger by assigning the bond orders and bond angles and then subjected to minimization using OPLS_2005 force field18. For accurate enumeration of ligand protonation states in biological condition we used Epik24,25. Then, Grid box of size 24 × 24 × 24 A3 was generated keeping ATP binding site of BCR-ABL (Glu286, Met318, Ile360, Ala380, Asp381) as centroid using the “Receptor Grid Generation” of Schrödinger Glide26.

Preparation of reference compounds

Reference compounds of FDA approved drugs such as ponatinib, bosutinib, bafetinib, dasatinib, nilotinib and imatinib were docked in the ATP binding site of both wild-type and mutant BCR-ABL. Ligands were prepared using the LigPrep module of Maestro followed by XP docking. The best docked complex score was set as the cutoff value for screening potential inhibitors of wild-type and mutant BCR-ABL.

Virtual screening

High throughput virtual screening (HTVS) was performed against the BCR-ABL using the Schrödinger software. Here, we used the virtual screening workflow of Glide Maestro for the docking. It performed the docking mainly in three phases namely, HTVS, SP (standard-precision) and XP (extra-precision). Using HTVS, we screened the large library of drug like molecules and reduced the intermediate conformations throughout the HTVS docking process. Successful compounds (10%) from HTVS were further assessed in SP docking for reliable docking of the screened compounds with high accuracy. To eliminate false-positive results, the best 10% of thriving compounds from SP docking were further incorporated for XP docking mode using advanced scoring.

Density functional theory analysis

Electronic effects of drug like molecules plays an important role in the pharmacological effects27. Therefore, drug candidates from the best binding poses were exported to Maestro 9.3 of Schrödinger 2012 version and geometry was optimized in the Jaguar panel using Becke's three-parameter exchange potential28 and Lee-Yang-Parr correlation functional (B3LYP) theory29,30 with 6-31G* basis set. Subsequently, surfaces (molecular orbital, density, potential) and atomic electrostatics potential charges (EPS) were monitored to compute the HOMO and LUMO. HOMO energy proposes the region of the small molecules, which can donate electron during the complex formation, while LUMO energy signifies the capacity of the molecule to accept the electrons from the partner protein. The difference in HOMO and LUMO energy, known as HOMO-LUMO gap energy, indicates the electronic excitation energy31 that is necessary to compute the molecular reactivity and stability of the compounds32.

Molecular dynamics simulations for protein ligand complexes

Two sets of Molecular Dynamics Simulations were performed using Gromacs 4.5.533. In the first set, we evaluated the stability of the mutant type of BCR-ABL and selected the best seven drug candidates from the docking study. In the second set, we evaluated drug candidates with wild-type BCR-ABL. We used MD for performing protein and ligand complexes as described34. The topology file for the selected small molecules were generated using the automated topology builder (ATB)35 in the framework of GROMOS 53A6 force-field36. The protein-ligand complexes were then solvated with TIP3P explicit water molecules and placed in the center of a cubic box of size 24 × 24 × 24 A3. Minimum 1.0 Å distance was maintained between the protein and the edge of the simulation box so that protein can fully immerse with water and rotate freely. Then, Particle Mesh Ewald (PME) method37 was used for the electrostatic energy calculation. It permits the use of the Ewald summation at a computational cost comparable to that of a simple truncation method of 10 Å or less and the linear constraint solver (LINCS)38 algorithm was used for covalent bond constraints. Before minimization, the system was neutralized by adding 8 Na+ ions. The steepest descent approach (1000 ps) was used for each protein-ligand complex for energy minimization. Further NVT was performed for 100 ps to equilibrate the system with protein and ligand for constant volume, pressure (1 atm) and temperature (300 K). The final MD run was set to 10000 ps for each protein-ligand complex and trajectories were saved for further analysis using Xmgrace and UCFC Chimera software39.

Rescoring of BCR-ABL and drug candidate complexes using interaction energy and MM-GBSA approach

Interaction energy and Gibbs free energy were calculated using Gromacs and Schrödinger software. Interaction energy for BCR-ABL and drug complexes was calculated by estimating the short range Lennard-Jones and short range Coulomb energy using the g_energy analysis tool of Gromacs software.

Here, Eint stands for interaction energy, ELJ stands for short ranges Lennard-Jones and ECoul denotes short ranges for coulomb energy. It estimates the stability of proteins and drug candidate complexes. The binding free energy was calculated based on the following equation40.

Here, ΔEMM stands for the energy difference between ligands in complex and unliganded receptor.

ΔGsolv is the difference in the GSA solvation energy of the complex and the sum of solvation energy for the unbound ligand protein complex. ΔGSolv is the difference in surface area energy for the complex and the sum of surface area energy for the protein and ligand.

Results and Discussion

Docking of reference compounds for HTVS

To set a cutoff value for docking studies, we docked all FDA approved drug such as ponatinib, bosutinib, bafetinib, dasatinib, nilotinib and imatinib using XP docking. Results showed that ponatinib has the highest binding affinity towards the T315I mutant BCR-ABL with a docking score of −11.050 kcal/mol while bosutinib, bafetinib, dasatinib, nilotinib and imatinib showed −5.513 kcal/mol, −4.689 kcal/mol, −4.702 kcal/mol, −3.772 kcal/mol, −3.593 kcal/mol and −3.78 kcal/mol, respectively (Table 1). Hence, for screening and identifying better lead molecules −11.00 kcal/mol was set as a cutoff value for XP docking.

Table 1 Glide XP results for the existing drug molecules with mutant and wild-type of BCR-ABL (T315I), by Schrödinger 9.3
Full size table

Binding pose analysis of the identified compounds

Results of binding pose analysis of seven lead molecules with better binding affinity and higher binding free energy than the reference compounds. Four of which, designated DB07107, DB06977, DB04200 and DB0117, were from DrugBank and ST007180, ST013616 and ST019342 were from the ligand.info database. The 2D conformations and drug details are given in (Figure 1).

Figure 1
figure 1

2-D Structures of the finally selected seven lead molecules have been shown in the figure.

(1) DB07107, (2) DB06977, (3) ST013616, (4) DB04200, (5) ST007180, (6) ST019342 and (7) DB01172.

Full size image

Interestingly, DB07107 (C23H22N4O) from DrugBank showed the highest binding energy with XP score of −14.045 kcal/mol (Figure 2a). To get an insight into their interacting pattern, we used UCFC Chimera molecular visualization tool and Glide for generating 2D interaction plots. Docking pose analysis revealed four hydrogen bonds (H_bonds) interactions with the ATP binding site residues of BCR-ABL. Here, we observed single H_bond with each Glu286 and Met318 residues with bond length of 2.92Å and 2.97Å while two H_bond formations were observed with Asp381 with the bond length of 3.32Å and 2.72Å.

Figure 2
figure 2

Binding poses of the DB07107, DB06977, ST013616, DB04200, ST007180 and ST019342 lead molecules.

The proposed binding mode of the lead molecules has been shown in the stick format. Residues involved in Hydrogen bonding have been labeled with the Hydrogen bond in dotted red lines and bond length have been shown in Angstrom.

Full size image

Other potent drugs, DB06977 (C23H21N5O), ST013616 (C14H12N4O2), DB04200 (C20H22O6), ST007180 (C18H17FN2OS), ST019342 (C15H11ClO4) (Figure 2b,2c,2d,2e,2f) and DB01172 (C18H36N4O11) (Figure 3) showed high binding affinity with XP scores of −13.163 kcal/mol, −12.065 kcal/mol, −12.041 kcal/mol, −11.555 kcal/mol, −11.033 kcal/mol and −8.433 kcal/mol, respectively (Table 2). Here, it is noticeable that the key residues of Glu286, Meth318 and Asp381 were involved in the H_bond interactions as has been observed with the marketed drug, ponatinib. Although DB01172 exhibits XP docking scored below the cuttoff value, we consider it as a potent inhibitor for BCR-ABL because it showed the highest number of H_bond with the hot spot residues of BCR-ABL. Moreover, it uncovered one side chain hydrogen bond with Glu282, Glu286, Asp381, Lys285 and one backbone hydrogen bond each with Asp381, Ile360. Therefore, we considered DB01172 for our further analysis.

Table 2 Glide XP results for the seven lead molecules with the mutant (T315I) and wild-type of BCR-ABL, by Schrödinger 9.3
Full size table
Figure 3
figure 3

Binding poses of the DB01172 lead molecule.

The proposed binding mode of the lead molecules has been shown in the stick format. Residues involved in Hydrogen bonding have been labeled with the Hydrogen bond in dotted red lines and bond length have been shown in Angstrom.

Full size image

Assessment of identified drugs in Wild-type of BCR-ABL

To determine the efficacy of the selected drug in binding to and inhibiting wild-type of BCR-ABL we performed docking studies as described in the methods section. To compare their efficacy and binding pattern with our screened compounds we again docked FDA approved drugs with wild-type of BCR-ABL.

As shown in table 2, DB07107 and DB06977 showed better docking scores than existing drugs while others showed results comparable to those of marketed drugs like ponatinib, bosutinib, bafetinib, dasatinib, nilotinib and imatinib. In the case of DB06977, although one H_bond is missing (Met318) with the wild-type of BCR-ABL complex compared to T315I mutant type, it still exhibits a favorable binding score of −10.94 kcal/mol. While, in the case of ST013616, one extra H_bond interaction is found to interact with Lys271. However, it did not show H-bond interaction with Asp381. Most of the selected drugs in these studies showed slightly less efficacy against wild-type compared to mutant BCR-ABL but better when compared to ponatinib. These investigations lead us to conclude that DB07107, DB06977, ST013616, DB04200 and ST007180 are more effective against the mutant than the wild-type BCR-ABL. ST019342 and DB01172 are effective only against mutant BCR-ABL.

DFT analysis

HOMO-LUMO plays an important role in stabilizing the interactions between drug and receptor protein32. Hence, the orbital energy of both HOMO and LUMO and the gap between them was calculated to estimate the chemical reactivity of the selected compounds using DFT (Table 3). HOMO-LUMO plots were generated to analyze the atomic contribution for these orbitals. The plots of HOMO and LUMO show the positive electron density in red color and negative electron density in blue (Figures 4 and 5). All the selected drug candidates showed minimal HOMO-LUMO gap with the average energy difference of 0.13 eV, signifying molecular reactivity. The lowest energy gap was exhibited by DB06977 with the HOMO-LUMO gap energy of 0.09911 eV and the highest energy gap was observed in DB04200 with the gap energy of 0.18442 eV (Table 3). HOMO energy was higher for DB07107 (Figure 4) and DB01172 (Figure 5), when compared to LUMO energy indicating an ability to donate the electrons rather than accept electrons with their partner receptor-binding site region.

Table 3 Frontier orbital energies of best seven lead compounds
Full size table
Figure 4
figure 4

Plots of highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of DB07107, DB06977, ST013616 and DB04200.

The positive electron density has been shown in red color while negative have been shown in blue.

Full size image
Figure 5
figure 5

Plots of highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) ST007180, ST013942 and DB01172.

The positive electron density has been shown in red color while negative have been shown in blue.

Full size image

Molecular dynamics simulations

To compare the structural behavior and flexibility of the wild-type and mutant BCR-ABL, all lead compounds were incorporated in Gromacs4.5.5 and MD was performed for 10 ns of each complexes. Root mean square deviations (RMSD) of the wild-type and mutant BCR-ABL were calculated against their initial structure in the protein-ligand complexes and graphs were generated to compare the flexibility of the backbone of the proteins using the Xmgrace software. Throughout the simulation period, no significant fluctuations were observed in the backbone of the wild-type and mutant T315I BCR-ABL, implying that the binding of drug candidates at the active site of the proteins is not only stable and strong but also does not disturb the protein backbone stability (Figure 6).

Figure 6
figure 6

Time dependence of root mean square deviations (RMSDs) of the Backbone of mutant (T3I51)and wild type of BCR-ABL have been shown in figure against the initial structure during 10,000 ps molecular dynamics (MD) simulation.

Full size image

To ensure the binding stability of the drug candidates in the active site of proteins, ligand positional RMSD of each lead molecule were generated and analyzed as described study41. ST019342 showed more and continues fluctuations in the noticeable window size of 0.2–2.5 nm (Figure 7a) for wild-type and 0.2–4 nm for the mutant BCR-ABL complex (Figure 7b). Co-ordinates of ST019342-BCR-ABL (wild-type and mutant) complex were downloaded from the trajectory in the interval of 1000 ps and investigated in PyMOL for protein-ligand interactions. Our investigation uncovered that ST019342 has an unstable binding affinity towards BCR-ABL. It exhibits weak hydrogen bond interactions with the receptor binding site, which leads to the inefficient inhibition of BCR-ABL. Upon investigation of downloaded timeframe of MD, we observed that after 2 to 5 ns ST019342 was not bound to the binding pocket. It suggests its inability to inhibit the protein target efficiently while other drug candidates have shown stable and strong binding affinity.

Figure 7
figure 7

Backbone RMSD values of drug candidates from both the types (wild and mutant) of protein ligand complexes were generated against the initial structures of protein-ligand complexes during 10,000 ps of molecular dynamics (MD) simulation period.

Graphs were plotted using Xmgrace, a 3-D plotting tool.

Full size image

In order to calculate the residual mobility of each lead molecules in BCR-ABL-ligand complexes (wild-type and mutant), Root Mean Square Fluctuation (RMSF) was calculated in each complexes and the graph was plotted against the residue number based on the trajectory period of MD simulation. None of the drug candidates showed noticeable changes in their residual level as the general profile of residual fluctuation of wild-type and T315I is minimal in each complex without any abnormal fluctuation (Figure 8).

Figure 8
figure 8

Root mean square fluctuations (RMSF) from the initial structures of protein-ligand complexes during the trajectory period of simulation.

Full size image

Hydrogen bond analysis

To determine the stability of hydrogen bonds with the ATP binding site of protein, MD analysis of BCR-ABL and the selected drug candidate's complex stability were monitored during the trajectory period. Hydrogen bond profiles between the selected drugs and BCR-ABL (T315I and wild-type) were calculated using the g_hbond utility of GROMACS. This analysis revealed that DB01172 comprises 6–7 (highest) average H_bonds during the simulation period sharing four H_bonds with Glu282, Lys285, Glu286 and Ile360 and two H_bonds with Asp381 (Figures 9 and 10) in T315I (mutant) while Wild-type showed poor H_bond interactions with 1–2 H_bonds on average during the trajectory period. The same pattern was also observed in the case of ST007180 and ST019342. ST007180 and ST019342 showed less number of hydrogen bond plots in the mutant BCR-ABL compared to wild type. To gain insight into the hydrogen bond instability, we downloaded the drug-receptor complexes in the interval of one ns for DB01172, ST007180 and ST019342. Following which, the structures were superimposed to monitor the binding pose of the drug candidates throughout the trajectories. Our investigation revealed that DB01172 is moving away from the ATP binding site of BCR-ABL in both mutant and wild-type towards the north-west (130°) and east direction (180°) (Figure 11). Compound ST007180 showed satisfactory results, while ST019342 showed slight movement in the binding site.

Figure 9
figure 9

Total number of inter-molecular H_bond interactions between the drug compounds (DB07107, DB06977, ST013616 and DB04200) and BCR-ABL (mutant and wild).

Full size image
Figure 10
figure 10

Total number of inter-molecular H_bond interactions between BCR-ABL (mutant and wild) and ST007180, ST019342 and DB01172.

Full size image
Figure 11
figure 11

Spnapshot of the superimposed structures of mutant and wild type of DB01172 and ST007180 and ST019342.

Stuctures were downloaded from the trajectory file in the interval of 1 ns.

Full size image

Other complexes, DB07107, DB06977, ST013616 and DB04200 exposed an average number of hydrogen bonds 2–3 (Glu286, Met318 and Asp381), 3 (Glu286, Met318 and Asp381), 3–4 (Glu286, Asp381, Glu316, Met318), 2–3 (Glu286 and Met318), 1 (322) and 1–2 (Met318 and Asn322), respectively in both wild-type and mutant BCR-ABL. The average number of H_bonds throughout the MD period signifies their continuous contribution in H_bond interactions, which provides stability to the complexes in holding it to the ATP binding position. This suggests that the functionality and ability of these compounds can efficiently inhibit BCR-ABL.

Interaction energy and binding free energy

To estimate the binding association between wild-type and mutant of BCR-ABL, two different approaches were carried out. Firstly, crude interaction energy was estimated based on the short-range energy of the system. Next, to determine the binding affinity between BCR-ABL and drug complexes we scored the MM-GBSA energy of the complexes. Interaction energy data suggests that all the selected drug compounds have higher interaction energy of 203444.3 kcal/mol, −203449.57 kcal/mol, −292066.20 kcal/mol, −310571.94 kcal/mol, −295967.25 kcal/mol, −292055.68 kcal/mol, −292134.56 kcal/mol in DB07107, DB06977, ST013616, DB04200, ST007180, ST019342 and DB01172 respectively, for the mutant than wild-type which showed −232088.91 kcal/mol, −232068.12 kcal/mol, −232062.86 kcal/mol, −232033.94 kcal/mol, −232030.35 kcal/mol, −232091.54 kcal/mol and −232062.86 kcal/mol for DB07107, DB06977, ST013616, DB04200, ST007180, ST019342 and DB01172 respectively (Table 4).

Table 4 Interactions energy and binding free energy of mutant and wild-type of BCR-ABL-ligand complexes
Full size table

Re-scoring of complexes using Prime MM-GBSA results indicated that all the selected drug candidates displayed higher binding free energy (Table 4). Moreover, all of them displayed results ideal for being a drug candidate in the treatment of CML bearing wild-type or mutant BCR-ABL. Our studies show that these drugs can efficiently inhibit BCR-ABL by binding to DFG out conformation of BCR-ABL protein and can obstruct its catalytic activity. Hence, we recommend DB07107, DB06977, ST013616, DB04200 and ST007180 to be experimentally tested and further authenticated. Since ST019342 and DB01172 are not recommended in the broad-spectrum purpose due to their unbalanced backbone conformation in wild-type during the MD study, we recommend for further experimentation and use against T315I. The entire work flow chart of this study is shown in Figure 12.

Figure 12
figure 12

Graphical abstract.

Full size image

Conclusions

Using a virtual screening approach, we performed a molecular docking study for both wild-type and mutant (T315I) of BCR-ABL. The study yielded seven lead molecules which showed promising results against both wild-type and T315I mutant BCR-ABL. The bioactivity of these selected lead compounds were determined using frontier orbital study and the results showed that the selected drug candidates are chemically reactive. The stability of the protein-ligand complexes were investigated through 10 ns of MD simulation of each complex for both types. MD simulation studies illustrated the dynamic behavior of protein-ligand complexes. Results showed that the protein backbone and ligand backbone of DB07107, DB06977, ST013616, DB04200, ST007180 and DB01172 compounds are stable throughout the simulation period without any significant fluctuation within this period. In contrast, backbone of ST019342 showed anomalous fluctuations throughout the simulation period while DB01172 showed instability in RMSD graph with the wild-type BCR-ABL. Hydrogen bond analysis revealed that DB01172 has six hydrogens on average, while in wild-type it has 2–3 on average, implying low efficacy towards the wild-type of target. Therefore, in terms of broad spectrum therapy for CML ST019342 and DB01172 drugs are not recommended while DB07107, DB06977, ST013616, DB04200 and ST007180 showed remarkable results during the trajectory period. These drugs promise a new gateway for the further development of anti CML therapeutics targeting BCR-ABL. In summary, DB07107, DB06977, ST013616, DB04200 and ST007180 are more effective in inhibiting mutant than wild-type of BCR-ABL. ST019342 and DB01172 are effective only in mutant BCR-ABL.