J Mol Model (2014) 20:2399 DOI 10.1007/s00894-014-2399-x

ORIGINAL PAPER

Dual-inhibitors of STAT5 and STAT3: studies from molecular docking and molecular dynamics simulations Shengjuan Shao & Rilei Yu & Yanqing Yu & Yanni Li

Received: 11 January 2014 / Accepted: 23 July 2014 # Springer-Verlag Berlin Heidelberg 2014

Abstract Although molecularly targeted therapy with imatinib has improved treatments of chronic myeloid leukemia (CML), clinical resistance gradually develops in patients with accelerated or blast phase CML. The inability of imatinib to cure CML suggests that inactivation of BCR-ABL kinase activity alone is not sufficient to control the disease. Aberrant STAT signaling and constitutive STAT5 or STAT3 activation are frequently found in both acute and chronic leukemia. Constitutive activation of STAT5 and STAT3 are associated with imatinib resistance on leukemia cells. Development of drugs targeting SH2 domains of STAT5 and STAT3 provides a novel strategy for the treatment of the imatinib-resistant CML. Here, molecular docking and molecular dynamics simulations were used to investigate the interactions of the drugs targeting STAT3 and STAT5 receptors at molecular level. The calculated binding free energies are consistent with the ranking of the experimental affinities and our simulations also explained their differences in binding energy. Then virtual screening based on molecular docking and molecular dynamics was applied to screen a set of ~1500 compounds for dual inhibitors of the SH2 domains of STAT5 and STAT3. Three top score

S. Shao Department of Chemistry and Chemical Engineering, Taiyuan Institute of Technology, Taiyuan 030008, China R. Yu Structural Bioinformatics and Computational Biochemistry, Department of Biochemistry, University of Oxford, OX1 3QU Oxford, United Kingdom Y. Yu : Y. Li (*) Key Laboratory of Systems Bioengineering, Ministry of Education, Department of Pharmaceutical Engineering, School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China e-mail: [email protected]

compounds obtained in virtual screening were compound 660, 304, and 561. Results show that the three predicted dualinhibitors are well fitted within the two binding domains and are predicted to present improved STAT5 and STAT3 SH2 inhibitory activity. Keywords Dual-inhibitor . MM-GBSA . Molecular dynamics . STAT5/STAT3 . Virtual screening

Introduction Chronic myeloid leukemia (CML) is a clonal hematopoietic stem cell disorder that is characterized by the presence of the Philadelphia (Ph+) chromosome and its resulting oncogenic fusion protein BCR-ABL. The BCR-ABL protein shares constitutively active tyrosine kinase activity and plays a central role in the pathogenesis of CML [1, 2]. The constitutive tyrosine kinase activity of BCR-ABL causes activation of a variety of intracellular signaling pathways, such as Ras/ mitogen-activated protein kinase (MAPK), NF-kB, phosphatidylinositol 3-kinase/Akt, Src family kinases (SFK), and STAT signaling cascades [3–5], leading to alterations in the proliferative, adhesive, and survival properties of CML cells. Targeting the tyrosine kinase activity of BCR-ABL is an attractive therapeutic strategy in CML. Imatinib (IM), as the first molecularly targeted therapy drug, specifically inhibits the BCR–ABL tyrosine kinase, and has been used for the firstline treatment of newly diagnosed CML with its excellent hematologic and cytogenetic responses, particularly in patients with chronic phase CML [6]. Nevertheless, resistance to the drug gradually develops in patients with accelerated or blast phase CML, which significantly affects the long-term therapeutic effects of imatinib [7–9]. The most common mechanisms of resistance to imatinib include BCR-ABL kinase domain mutations, BCR-ABL amplification and over-

2399, Page 2 of 17

expression, activation of signal transduction pathways like MAPK, PI3K/Akt, SFK and STATs, P-glycoprotein (P-gp) over-expression, and CML stem cell quiescence [4, 10, 11]. Thus inhibition of STATs is one of the effective ways to decrease the imatinib resistance. Aberrant STAT signaling and constitutive STAT5 or STAT3 activation have been frequently found in both acute and chronic leukemias [12]. Recently, more and more assays suggested that the constitutive activation of STAT5 and STAT3 was associated with imatinib resistance on leukemias cells [13–15]. Wang et al. [13] identified that GM-CSF elicited IM drug resistance via a BCR-ABL-independent activation of JAK2/STAT5 in LAMA84 cell clones with high p-STAT5 levels, and GM-CSF-mediated JAK2/STAT5 pathway activation just circumvented the need for BCR-ABL signaling to maintain survival or proliferation of CML cells. Meanwhile, Warsch et al. [14] suggested that the increased expression of STAT5 itself, independent of JAK2 also sufficed to mediate imatinib resistance. They showed that STAT5 over-expression could lead to a imatinib resistant phenotype. These findings were in accordance with the work of Hoelbl et al. [16]. They proposed that STAT5 activation through direct phosphorylation by BCR-ABL or indirectly through phosphorylation by HCK or JAK2 was associated with IM drug resistance, and STAT5 deletion induced G0/G1 cell cycle arrest and apoptosis of imatinib-sensitive and imatinib-resistant stable leukemic cells. Interestingly, STAT3 was considered as an anti-cancer target protein because of its increased activation in a broad spectrum of cancer cell lines and human tumors [17, 18]. Bewry [15] and his colleagues first investigated the influences of STAT3 on imatinib resistance in CML. The data showed that K562 cells within the bone marrow microenvironment could cause Bcr-Abl-independent resistance to imatinib with increased pTyrSTAT3 and increased levels of the STAT3 target genes Bcl-xl, Mcl-1, and survivin. Furthermore, reducing STAT3 levels with siRNA could reverse the imatinib resistance. Thus it is likely that STAT3 and STAT5 can both be related to IM drug resistance by regulating the expression of critical target genes independence of BCR-ABL. These results suggest that STAT3 and STAT5 could serve as attractive targets to overcome imatinib resistance. As members of the family of STATs, STAT3, and STAT5 are structurally related especially in SH2 domain which is involved in receptor recruitment and dimerization (Fig. 1) [19]. Drugs targeting to the SH2 domain could prevent the tyrosine phosphorylation of STAT monomer, dimerization of phosphorylated STAT monomers, nuclear translocation, and gene expression, and consequently induce the arrest of cell cycle and apoptosis. In fact, tumors dependent on the constitutive activation of one of the STATs may also develop resistance through activation of other STAT family members [20]. As previously described, within the bone marrow

J Mol Model (2014) 20:2399

microenvironment, STAT3 can compensate for the loss of STAT5 and cause drug resistance in K562 cells. In all, STAT3 and STAT5 share high sequence similarity (59 %) of SH2 domains [21], overlapping regulatory pathways, and even functionally complemented, which forms the molecular basis for the discovery of new small-molecule compounds that can target the two proteins simultaneously. Since STAT3 was regarded as a candidate target for cancer therapy, many efforts have been geared toward design of its inhibitors. These efforts include computational threedimensional database screening [22–24], experimental screening [25, 26] and structure-based de novo design [27–29]. Earlier, the blockers of STAT3 SH2 domain mainly included some small peptides and peptidomimetics. Turkson et al. [30] identified the phosphotyrosyl peptide (Pro-pTyr-Leu-Lys-ThrLys) that could block STAT3-mediated DNA binding activity, gene regulation, and cell transformation. This peptide inhibitor provides the conceptual basis for using of the peptide as a lead for novel peptide-mimetic drug design and opens the door for the discovery of novel small-molecule inhibitors that can mimic its function. Recently, a few non-peptide STAT3 SH2 inhibitors have been developed to inhibit STAT3 dimerization, including Stattic [21], STA-21 [23], S3I-201 [24], FLLL32 [28], and LLL12 [29]. Several new inhibitors of JAK2, the upstream kinase of STAT3, such as curcumin [31], AG490 [32], WP-1066 [33], and WP-1034 [34] have also been reported. These molecules are more selective for STAT3 than for its homolog STAT5, of which, the apparent IC50 values for Stattic and S3I-201 in binding to STAT3 SH2 domain are 5.1 μM [21] at 37 °C and 86 μM [24] at 25 °C. In contrast, the IC50 values in binding to STAT5 SH2 domain are more than 100 μM. So far, there are no nonpeptidic molecules that explicitly inhibit the function of the STAT5 SH2 domain except for SC-355979 with an IC50 of 47±17 μM [35]. Müller et al. [35] identified this compound along with a series of chromone derivatives, such as Acyl hydrazone 2 and Acyl hydrazone 3 [35], that could inhibit the function of the STAT5 SH2 domain when they screened diverse chemical libraries by using a homogeneous assay based on fluorescence polarization (Fig. 2). Their results showed that some compounds also cause comparable inhibition of STAT3 to STAT5. Thus the development of novel compounds that are STAT5 specific or optimized for dual-inhibition of STAT3 and STAT5 is a necessary step to block the STAT signaling and to overcome the resistance of imatinib in CML cells. In the present study, we investigated the interactions between three known ligands and STAT3 and STAT5 proteins using molecular docking, molecular dynamics simulations and binding free energy calculations. Results from our simulations were consistent with experimental studies. In addition, we conducted a virtual screening against the STAT3 and STAT5 SH2 domains and discovered several novel dualinhibitors of STAT3 and STAT5, which were predicted to

J Mol Model (2014) 20:2399

Page 3 of 17, 2399

Fig. 1 SH2 domains of STAT3 and STAT5. The SH2 domains are presented in cartoon (a, c) and surface (b, d) respectively. SH2 domain is characterized by a central βB, βC, βD anti-parallel β-sheet that is

flanked by adjacent α-helices. (e) A superimposed structure comparison between STAT5 and STAT3 SH2 domains

possess better affinities than the known STATs inhibitors. These novel inhibitors can be tested in future experimental studies. It is expected that our results will eventually be used in the design of more potent dual-inhibitors of STAT3 and STAT5 and help overcome the resistance of imatinib clinically.

Materials and methods

Fig. 2 Chemical structures of the three positive controls [24, 35]

Preparations of ligand database and receptors Due to the structural similarity of the SH2 domains within the two proteins, a STAT5-inhibitor can be a derivative structure

2399, Page 4 of 17

of the known STAT3-inhibitors, and vice versa. We created a library of ~1500 derivative structures similar to STAT3 inhibitor S3I201 [24], STAT5 inhibitor SC-355979 [35], JAK2 inhibitor curcumin [31] and a novel JAK-STAT inhibitor WP-1034 [34] by searching the PubChem database with a similarity threshold of 90 %. Compounds similar to the query structure are measured using the Tanimoto score [36]. Various predefined thresholds between 99–60 % are allowed. Using the similarity search based on the known active inhibitors can provide a reasonable number of compounds with the ones that are not similar to the original query structure returned and yet yield higher hit rates [37]. Our aim is to find small molecules that can be docked into the SH2 domains of STAT3 and STAT5, and that potentially disrupt the STAT-STAT dimerization. We retrieved residues 580–681 from the PDB entry 1BG1 [38] and 589–682 from the crystal structure (PDB ID:1Y1U) [39] as the STAT3 and STAT5 SH2 domain, respectively. In addition, S3I201 [24] and two analogs of SC355979 [35], Acyl hydrazone 2 and Acyl hydrazone 3 [35] that have known binding affinities with both STAT3 and STAT5 were also used as the positive controls in this study (Fig. 2). Ligand screening Ligand screening on the pTyr-binding sites within the SH2 domain of STAT3 and STAT5 was performed using AUTODOCK, version 4.2 [40]. Hydrogen atoms were added to the receptors and ligands and partial atomic charges were then assigned using the Gasteiger–Marsili [41] method. For STAT3, a docking grid map with 70×70×70 points and grid spacing of 0.375 Å was set up by centering near the centroid of the pTyr705 of the other monomer. Rotatable bonds of each ligand were automatically assigned using AUTOTORS utility. Docking was performed using the Lamarckian Genetic Algorithm (LGA) method with an initial population of 300 random individuals, a maximum number of 25,000,000 energy evaluations, ten trials, 27,000 maximum generations, a mutation rate of 0.02, a crossover rate of 0.80, and the requirement that only one individual can survive into the next generation. With the structural similarity of the two SH2 domains, the above parameters were also applied to STAT5 with only slight alterations, when its structure was converted into the coordinate system of the STAT3. At the end of the docking, the docked conformations were clustered into subgroups within the RMSD values of 2.0 Å. The complex structures were ranked based on the lowest energy of the largest cluster. In all, the full set of compounds was first screened against the STAT5 SH2 domain, and the top 200 STAT5-hits from this screening result were then filtered for compounds that can bind to STAT3 SH2 domain with high affinity. The ranking of compounds involved two successive iterations. The AUTODOCK method was used to place the compounds within the binding sites of

J Mol Model (2014) 20:2399

STAT3 and STAT5 and to search for their minimal energy conformations. Then, the top ten hits obtained from the AUTODOCK screening process were rescored using MMGBSA method [42] after molecular dynamics simulations. Molecular dynamics simulations Molecular dynamics (MD) simulations were carried out on the protein-ligand complex and protein in the apo state using Amber 10 package [43]. FF99SB and GAFF force field were employed for the receptors and ligands respectively [44]. Each protein-ligand complex was solvated into an 8 Å thick octahedral TIP3P water box and neutralized with Na+ or Cl− according to the net charge of the complex. Then it was minimized to remove unfavorable van der Waals interactions. The sander module of Amber 10 was used. After 2000 cycles of constrained minimization including 1000 steps of steepest descent method and 1000 steps of conjugate gradient method, the whole system was then relaxed with 1000 steps of steepest descent method and 2000 steps of conjugate gradient method. The cutoff of the non-bonded interactions was set to 10 Å for the energy minimization process. The system was then heated up from 0 to 300 K over 60 ps via the Langevin temperature regulation scheme with constant volume restraints of 10.0 kcal mol −1 Å −2 weight. Following heating, the system was gradually equilibrated using a small constant pressure restraints of 1.0 kcal mol −1 Å −2 for 100 ps. Finally, the whole system was relaxed and a 5 ns long molecular dynamics simulation was carried out. Coordinate trajectories were saved every 2 ps throughout the production runs, and 1000 conformations from the last 2 ns MD trajectory were extracted for the further binding free energy calculations. The simulations used the SHAKE [45] algorithm to constrain the covalent bonds involving hydrogen atoms and the particle mesh Ewald (PME) [46] method to calculate the long-range electrostatic interactions with a time step of 2 fs. MM-GBSA calculation Binding energies were calculated for the snapshots collected from the equilibrated trajectories of the simulation systems via MM-GBSA method [42]. The values of the binding free energy (ΔG) of each inhibitor were calculated according to the equation: ΔGbind ¼ Gcom − Grec − Glig ;

ð1Þ

where com, rec, and lig denote the complex, receptor, and ligand, respectively. The free energy of each was estimated as a sum of the four terms: G ¼< E MM > þ < GGB > þ < F SUR > −T < S >;

ð2Þ

where EMM is the molecular mechanics energy expressed as the sum of the internal energy of the molecule plus the

J Mol Model (2014) 20:2399

electrostatics and van der Waals interactions, GGB is the polar contribution to the solvation free energy calculated by the GB module at igb=5 and saltcon=0.15, GSUR is the nonpolar solvation energy calculated by an empirical model with the solvent accessible surface area (SASA) coefficient of 0.0072, T is the absolute temperature, and S is the entropy of the molecule estimated using the normal mode analysis. Considering the high computational cost, the produced 1000 frames of the last 2 ns production were extracted at 20 frames intervals, i.e., a total of 50 frames were utilized for entropy calculations [42]. To investigate the energetic contributions of key residues to the ligand binding, we also performed pairwise per-residue energy decomposition with the GB implicit solvent model via MM-GBSA [47]. Using this model the calculated energies can be further decomposed into backbone, side chain, and total contributions to their decomposition energies.

Results and discussion Constitutive activations of STAT5 and STAT3 were reported to be associated with imatinib resistance on leukemias cells [12]. As a critical step in STAT activation, the dimerization between two STAT monomers represents an attractive target to abolish STATs DNA-binding and transcriptional activity and to inhibit STATs biological functions. The dimerization relies on the reciprocal binding of the SH2 domain of one monomer to the pTyr peptide of the other, i.e., pTyr–SH2 interactions [30, 48]. Consequently, targeting the SH2 domains of STAT5 and STAT3 will be a promising and effective therapeutic approach to overcome the imatinib resistance. Recently, many STAT3 SH2 inhibitors have been discovered via both computational and experimental methods [23–28]. STATs SH2 domains are quite conserved with similar backbone conformation between different members [19] and merely small differences are identified between STAT3 and STAT5 SH2 domains. Thus, it offers promise for the discovery of small-molecule inhibitors which can target the two SH2 domains simultaneously. Results of positive inhibitors S3I201, Acyl hydrazone 2 and Acyl hydrazone 3 were reported to have binding affinities with both STAT3 and STAT5 SH2 domains [24, 35] and thus were selected as the positive controls in the present study. The apparent IC50 values for S3I201, Acyl hydrazone 2 and Acyl hydrazone 3 in bindings to STAT3 are 86 μM [24], 54 μM [35], and 159 μM [35], respectively. Regarding the bindings to STAT5, the IC50 values are 166 μM [24], 53 μM [35], and 79 μM [35], respectively. To make a more intuitive comparison, we converted IC50 values from the fluorescence polarization (FP)

Page 5 of 17, 2399

competition assay to the inhibition constant Ki values and calculated the related ΔG with the formula Ki = IC50/(1+ [L*]/Kd) and ΔG = 1.36logKi (kcal mol-1) at 300 K [49]. The Ki value describes the strength of the protein-ligand interactions and it is equivalently used with the dissociation constant Kd and the association constant Ka. However, if Ki is not available, IC50 value can be used as replacement of the Ki, given that both values run approximately parallel and the latter is easier to determine and is also well suited for the characterization of a ligand in comparison to other structures [50]. Table 1 lists the estimated binding energies for the three compounds compared to the experimental values. The predicted values of the binding energies obtained by docking method were in better agreement with the experimental data comparing to values given by MM-GBSA. Interestingly, MM-GBSA approach seems highly sensitive and the results are generally evident to distinguish compounds with lower activity. S3I201 was regarded as a selective chemical probe inhibitor of STAT3 and a weaker binder to STAT5 [24] with corresponding IC50 data of 86 and 166 μM. As shown in Table 1, S3I201 was predicted to bind STAT5 with a highly positive binding energy 4.83 kcal mol−1 given by MM-GBSA. Acyl hydrazone 3 was also tested as a weak binder to STAT3 with the positive value of 3.86 kcal mol−1. Besides, MMGBSA was more robust than molecular docking to rank of the binding affinities of these positive controls. As shown in Table 1, the ranking of the binding affinities to STAT5 given by IC50 data is Acyl hydrazone 2>Acyl hydrazone 3>S3I201, which is consistent with theoretical ranking given by MMGBSA but not by molecular docking. Additionally, molecular docking also failed to rank the binding affinities of these positive controls to STAT3, while MM-GBSA method gave correct ranking (Table 1). These results are not a surprise. Typically, docking-score function which is experienced based can fast yield good predictions for the binding affinities of small-molecule ligands to their protein targets with high resolution. However, it lacks the accuracy to rank the activity of the ligands with different structures in comparison with the more rigorous MMPBSA method which is physical-function based and more time consuming. To take advantage of the different scoring functions, we combined the molecular docking and MM-GBSA methods for database screening in the following studies [51, 52]. Dynamics of the positive controls It has been established that MD simulations can be used as a post-docking tool for the refinement of the final docked complexes. Although this procedure is computationally expensive, it is able to refine the protein-ligand interactions by taking the flexibility of the receptor and ligand into consideration. Moreover, MD simulations can be used to evaluate the stability of the complex, as a direct measure of the consistency of binding,

2399, Page 6 of 17

J Mol Model (2014) 20:2399

Table 1 Relative ranking of positive controls using the two scoring methods compared to experimental dataa STAT5 ranking

STAT3 ranking

Ligand

Autodock

MM-GBSA

IC50

Ki b

ΔG’c

Autodock

MM-GBSA

IC50

Ki b

ΔG’c

Acyl hydrazone 2 Acyl hydrazone 3 S3I201

−6.23 −6.00 −6.42

−6.53±0.94 −4.64±1.08 4.83±0.97

53d 79d 166f

27 39 /

−6.22 −5.99 /

−6.49 −6.02 −7.44

−5.18±1.00 3.86±1.56 −4.05±0.38

54e 159e 86f

27 75 /

−6.21 −5.61 /

a

All energies are in kcal mol−1

b

The Ki values are calculated in the formula Ki = IC50/(1+ [L*]/Kd)

c

The experimental values are calculated in the formula ΔG=1.36logKi (kcal mol−1 ) at 300 K

d

values are from the FP competition assay with native pTyr peptide (Kd =105 nM and [L*]=10 nM)

e

values are from the FP competition assay with native pTyr peptide (Kd =150 nM and [L*]=10 nM)

f

IC50 data of cellular level measurement

since improperly docked structures are expected to produce unstable trajectories. It is essential to make sure that equilibrium has successfully been achieved before moving on to consider the structures and detailed interactions. The potential energy plots for the three positive controls with the two targets help to see if the systems have reached an equilibrium (see Fig. 3). From the plot, we can see that all of the potential energy terms initially increased, corresponding to the heating from 0 to 300 K, and then leveled off and plateaued when we switched to smaller restraints in the constant pressure stage. Subsequently, the potential energy decreased and remained constant indicating that the systems were totally relaxed and reached an equilibrium. It can be seen that the potential energy was fluctuating around a constant mean value, and the plots for the STAT5 simulations were very similar. A detailed interactions analysis of the three positive controls with the two targets provides us insights into the design of future STAT5/STAT3 SH2 dual-inhibitors. Here, we focused on several aspects involved in the protein-ligand binding, such as structure based dynamics, H-bond analysis, and decompositions of the binding energy.

Plots of the RMSDs for the backbone atoms from the initial coordinates of the STAT3 and STAT5 SH2 domains for the 5 ns simulation, both in the apo and holo states, illustrated the inherent stability of the complex and the best fit interactions (see Fig. 4). The RMSDs of all the systems gradually became stable after 1 ns, indicating that they reached equilibrium. For STAT3 SH2 domain, all the holo-systems induced fluctuations equivalent to the apo-STAT3. For the apo-STAT5 SH2 domain simulations, the protein backbone RMSD fluctuated around a large value of 3.75 Å that may be related to the low resolution of 3.21 Å for the crystal structure of STAT5a [39]. However, the RMSDs of its holo-systems fluctuated around lower values of less than 2.5 Å, indicating mutually stabilizing effects induced by protein–ligand interactions. In addition, the RMSFs of the protein backbone atoms for residues constructing the hot spots within STAT3 and STAT5 SH2 domains also illustrated the inherent stability of the complex during the molecular dynamic simulations (see Fig. 5a, b). As shown in the Fig. 5a, there were three highly fluctuating regions within STAT3 SH2 domain, i.e., the binding hot spots pY-X (Δ592–605), pY+0 (591, Δ609–620), and pY+1 (Δ626– 639) site, respectively, which was in line with the previous

Fig. 3 The potential energies of STAT3 and STAT5 SH2 domains in complex with Acyl hydrazone 2 (red), Acyl hydrazone 3 (blue), and S3I201 (green) versus simulation time

J Mol Model (2014) 20:2399

Page 7 of 17, 2399

Fig. 4 a The RMSDs of STAT3 SH2 domain in the apo state (black), and in complexes with Acyl hydrazone 2 (red), Acyl hydrazone 3 (blue), S3I201 (green). b The RMSDs of STAT5 SH2 domain in the apo state (black), and in complexes with Acyl hydrazone 2(red), Acyl hydrazone 3(blue), S3I201(green). It could be seen that the RMSDs of the three

systems gradually became stable after 1 ns, indicating well equilibrated models were reached, which also meant that it was reasonable for us to analyze the detailed interactions between the protein and the ligands further

modeling study of Park et al. [49]. The so-called pY-X, pY+0, and pY+1 sites are the pTyr 705 binding pocket, Leu706 binding pocket and residues at pY-X position binding pocket, relatively [30, 49]. The hydrophobic side pocket pY-X was unique to STAT3 and was drug targetable, holding all three inhibitors tightly. The RMSFs of this pocket in the holosystems were much lower than that of the apo states. The STAT5 SH2 domain possessed two sub-sites pY+0 (Δ618– 628) and pY+1 (Δ632–645) that were also characterized by noticeable fluctuations shown in Fig. 5b, which correlated well with previous studies: reference [19, 53]. The RMSFs of residues in holo-systems were lower than that in aposystems, suggesting tight binding interactions between ligands and these proteins. For the pY+0 pocket, the region within STAT5 SH2 domain was more stabilized, and the conformation fluctuations induced by bindings were almost unnoticeable, indeed, this pocket was much shallower, and less flexible than that of STAT3. Larger fluctuations were found in the regions of pY+1 pocket, indicating the bindings to pY+1 subsite were overall loose, which confirmed the findings of Park et al. [49] that the pY+1 pocket was the most dynamic and difficult to target.

SH2 domain and maintained a stable value of~2 during the overall dynamic process. Analysis of H-bond donor-acceptor pairs provided us a detailed understanding of protein-ligand interactions. As we can see in Table 2, there were five remarkably strong hydrogen bonds between S3I201 and the polar pY+0 pocket of STAT3 SH2 domain that consists of Lys591, Arg609, and Ser611, which was consistent with the researches of Siddiquee et al. [24]. The negatively charged salicylic acid group of S3I201 played an important role in forming hydrogen bonds with STAT3 and can be regarded as a pTyr-mimetic [54, 55]. For STAT5 SH2 domain, Acyl hydrazone 2 formed most of hydrogen bonds with residues surrounding the pY+0 pocket (Ser620, Asp621, and Thr628). It suggests that the dominant residues in the pY+0 pocket are polar and frequently form Hbonds with the polar functional groups of Acyl hydrazone 2. In contrast, most of the residues in pY+1 and pY-X pockets are hydrophobic [49] and they form a smaller number of Hbonds with these compounds.

Hydrogen bonds analysis Hydrogen bonds are important indicators of the electrostatic interactions between the protein and ligands, and play significant roles in the protein–ligand recognition process [28]. The hydrogen bonds density was calculated along the 2 ns trajectory (see Fig. 6). The plots gave us an overview of the general evolutions of hydrogen bond density between proteins and ligands. It can be easily seen that S3I201 formed remarkable hydrogen bonds with STAT3 SH2 domain, implying strong electrostatic interactions with the target. For Acyl hydrazone 2 and 3, strong hydrogen bonds were identified with STAT5

MMGBSA decompositions of the binding energy To evaluate the energetic profiles of interactions between proteins and the three inhibitors, energy decompositions were performed using MM-GBSA approach [56, 57]. The binding free energy includes enthalpy term (ΔH) that is the sum of van der Waals (vdW), electrostatic (EEL) as well as the solvation interactions and entropic component (TΔS) (see Table 3). Typically, a negative TΔS value is expected for a biological complex. This symbolizes a decrease in available microstates as a result of the ligand being trapped and having limited mobility. However, such a negative value also indicates an entropically unfavorable protein-ligand complex. Since the binding affinity is determined by the Gibbs energy (Ki =exp (-ΔG/RT)) and ΔG is given by ΔG= ΔH – TΔS, it is

2399, Page 8 of 17

Fig. 5 The root mean square fluctuations of the binding hot spots of STAT3 (588–640) and STAT5 (615–645) SH2 domains in complex with the top three hits, Acyl hydrazone 2 (red), Acyl hydrazone 3 (blue), and S3I201 (green), compared with that in the apo state (black). a apo-STAT3

J Mol Model (2014) 20:2399

and holo-STAT3; b apo-STAT5 and holo-STAT5; c, d The binding hot spots of STAT3 and STAT5 SH2 domains with pY+0, pY+1, and pY-X designations

Fig. 6 Evolution of hydrogen bond density for the protein-ligand complex. a holo-STAT3, b holo-STAT5. Ligands: Acyl hydrazone2(red), Acyl hydrazone3(blue), S3I201(green)

J Mol Model (2014) 20:2399

Page 9 of 17, 2399

Table 2 The occupancies of hydrogen bonds formed between the ligands and receptors during the simulation time STAT3

STAT5

Ligand

Residue: atom ID

Group: atom ID

Occupancy/%

Residue: atom ID

Group: atom ID

Occupancy/%

Acyl hydrazone2

SER636: N / / SER636: N / ARG609: N LYS591: NZ ARG609: N SER611: OG SER636: N

: N9 / / : N8 / : O7 : O7 : O6 : O6 : O4

35.5 / / 65.1 / 83.4 64.0 49.2 40.10 91.7

SER620: OG ASP621: N THR628: OG1 ASP621: N ASN642: OD1 ASN642: O ASN642: OD1 / / /

:N8 :O23 :N9: :N6 :N8 :N1 : N1 / / /

79.6 77.3 41.6 75.2 95.6 90.0 77.0 / / /

Acyl hydrazone3 S3I201

Hydrogen bond criteria are 3.5 Å for donor-acceptor distance and 135° for donor-H-acceptor angle

apparent that a favorable binding can be accomplished either by making ΔH more negative, or making ΔS more positive, and the total enthalpy and entropy changes upon binding will determine the binding affinity [58]. Just as shown in Table 3, Acyl hydrazone 2 and S3I201 showed strong binding behaviors with STAT3 SH2 domain, mainly due to their preferable enthalpy contributions, and the latter also benefited from the minimum entropic penalties. For STAT5 SH2 domain, Acyl hydrazone 2 and 3 attributed their tight binding interactions to the large negative enthalpy contributions. In contrast, Acyl hydrazone 3 and S3I201 showed a weak binding behavior with STAT3 and STAT5 SH2 domain, respectively, as a result of minimum enthalpy contributions. Generally, the favorable EEL term and the unfavorable polar solvation energy component in GB often offset with

each other (see Table 3) [49]. The polar energy contributions ΔGpol, the sum of Coulombic electrostatic interaction EEL and polar solvation energy GB term, often show positive values and are unfavorable for the bindings of these ligands. Here, the only exception was S3I201. Its negative value guaranteed a higher potency of being an effective STAT3specific inhibitor [24]. The more preferable electrostatic interaction was resulted from the polar interactions between the negatively charged salicylic acid group and polar residues in the pY+0 site [55]. For all the systems, the non-polar interactions ΔGnon-pol acted as the main driving force, and the van der Waals interactions made major contributions for the binding. Acyl hydrazone 3 and S3I201 mainly attributed their weak binding affinities with STAT3 and STAT5 SH2 domains, to the minimum non-polar energy contributions, especially the

Table 3 Energy terms contributing to the binding free energy STAT3

STAT5

Contributiona

Acyl hydrazone2

Acyl hydrazone3

S3I201

Acyl hydrazone2

Acyl hydrazone3

S3I201

ΔEVDW ΔGSUR ΔEEEL ΔGGB

−29.84±0.38 −4.14±0.03 −8.75±0.64 21.87±0.60

−16.13±0.56 −2.91±0.07 −6.46±1.05 14.70±0.99

−14.24±0.59 −2.89±0.07 −91.79±2.90 91.09±2.59

−26.68±0.42 −3.36±0.04 −23.21±1.34 30.70±1.02

−24.35±0.46 −3.54±0.04 −24.47±0.67 31.57±0.52

−20.01±0.33 −3.51±0.03 −80.35±2.24 91.63±2.07

ΔGpolb ΔGnon-polc ΔHd TΔS ΔGe

13.13±0.21 −33.98±1.18 −20.85±0.35 −15.67±1.24 −5.18±1.00

8.23±0.27 −19.04±0.84 −10.81±0.43 −14.67±2.61 3.86±1.56

−0.70±2.48 −17.13±0.86 −17.83±1.06 −13.77±0.93 −4.05±0.38

7.49±1.62 −30.04±0.72 −22.55±0.54 −16.02±1.33 −6.53±0.94

7.10±1.15 −27.89±0.75 −20.78±0.39 −16.15±1.62 −4.64±1.08

11.30±2.20 −23.52±0.63 −12.24±0.31 −17.07±1.16 4.83±0.97

a

All energies are in kcal mol−1

b

ΔGpol = ΔEEEL + ΔGGB

c

ΔGnon-pol = ΔEVDW + ΔGSUR

d

ΔH = ΔGpol + ΔGnon-pol

e

ΔG = ΔH – TΔS

2399, Page 10 of 17

weak van der Waals interactions. According to the energy decompositions, we can draw the following conclusions: the differences of the binding free energies of the three positive controls are primarily associated with the enthalpy contributions, especially the non-polar energy contributions and the differences of the entropy penalties are not significant. In order to get a deeper insight into the protein ligand interaction patterns, the binding free energy was also decomposed into specific residue-ligand pairs based on MM-GBSA [59]. Figure 7 shows the energy contribution of each residue to the binding free energy. In the bindings with STAT3, Acyl hydrazone 2 primarily attributed its binding free energies to the polar contributions by Lys591 in the pY+0 pocket and Glu592, Arg595 in the pY-X pocket. In addition, its non-polar contributions by five residues in the pY+1 site, which includes Ile634, Gln635, Ser636, Val637, and Glu638 also partially contributed to the binding free energies. The non-polar binding contributions for Acyl hydrazone 3 and S3I201 are mainly resulted from residues around the pY+1 site, and they formed few contacts with residues in the pY-X site. Interestingly, S3I201 mainly attributes its remarkable contributions to the polar interactions with residues Lys591,

J Mol Model (2014) 20:2399

Arg609, Ser611, Gln635, and Ser636, mainly arising from the pTyr mimetic. Overall, the major contributions for Acyl hydrazone 2 and S3I201 are the non-polar and polar energy contributions, respectively, and the critical residues related are Lys591, Arg609, Ser611, Gln635, and Ser636. Given the fewer non-polar and polar interactions with STAT3 SH2 domain, Acyl hydrazone 3 shows a weak binding behavior with STAT3. In the bindings with STAT5 for the three positive controls, the significant contributions mainly came from key residues Arg618, Ser620, Asp621, Ser622, Thr628 of the pY+0 site and Asn642, Leu643, Lys644, Pro645 of the pY+1 site, which confirmed the work of Page, Brent et al. [53]. For compound Acyl hydrazone 2 and 3, they formed strong polar and nonpolar contacts with both pY+0 and pY+1 sites. In comparison, S3I201 made fewer contributions with residues surrounding the pY+1 site within STAT5 SH2 domain. Due to the lack of strong interactions with residues in the pY+0 site, S3I201 shows a weak binding behavior with STAT5 SH2 domain [24], which may be related with the poor physical occupation of the large volume salicylic acid group for the shallow pY+0 pocket. The loose binding interactions (van der Waals,

Fig. 7 Energetic contributions of key residues within STAT3/STAT5 SH2 domain to the binding affinity of different systems. Ligands: Acyl hydrazone 2 (red), Acyl hydrazone 3 (blue), S3I201 (green)

J Mol Model (2014) 20:2399

Page 11 of 17, 2399

hydrogen bonds, etc.) directly reflects a decrease in the enthalpy and the ‘induced fit’ model with tight binding interactions is preferred. Virtual screening Virtual screening is a well-established technique to discover novel molecular inhibitors that are complementary to a target protein in terms of shape, charge and several additional biophysical or biochemical properties [52]. Here, we applied this

technique to screen STAT5/STAT3 dual-inhibitors. A compound library of ~1500 small molecules were screened against STAT5 SH2 domain via AUTODOCK, and top 200 hits that showed strong affinity for STAT5 SH2 domain were used in a second round of screening against the SH2 domain of STAT3. The MM-GBSA method was applied to predict the absolute binding energies for the top ten hits from the second round of screening. The binding energies were obtained by averaging 20 snapshots evenly extracted from the last 2 ns MD trajectories.

Table 4 STAT5/STAT3 top 35 hits ranked by their binding energies as were calculated using the AUTODOCK method. Compounds in bold were predicted to bind to STAT5 and STAT3 with high affinities STAT5 ranking

STAT3 ranking

Rank

BE (kcal mol−1)

CID

BE (kcal mol−1)

CID

1 2 3 4 5 6 7 8 9

−8.40 −7.93 −7.89 −7.78 −7.49 −7.38 −7.36 −7.28 −7.26

Pub#7515545 Pub#1878702 Pub#24233809 Pub#25423195 Pub#3959254 Pub#3871088 Pub#46441391 Pub#4023335 Pub#4926251

561 660 496 33 625 627 420 624 45

1 2 3 4 5 6 7 8 9

−9.28 −9.06 −8.97 −8.84 −8.63 −8.58 −8.55 −8.53 −8.48

Pub#1008300 Pub#25423195 Pub#3903522 Pub#5882222 Pub#16750713 Pub#1878702 Pub#1083660 Pub#25423183 Pub#7515545

668 33 343 120 516 660 667 9 561

10 11 12 13 14 15 16 17 18 19 20

−7.25 −7.23 −7.22 −7.21 −7.20 −7.20 −7.20 −7.19 −7.19 −7.18 −7.18

Pub#23013654 Pub#52221631 Pub#39821522 Pub#25423215 Pub#4926221 Pub#52725077 Pub#2131903 Pub#8199626 Pub#25424243 Pub#4926157 Pub#8544853

35 381 1 271 18 223 326 307 488 610 304

10 11 12 13 14 15 16 17 18 19 20

−8.47 −8.33 −8.29 −8.23 −8.17 −8.13 −8.11 −8.07 −8.06 −8.03 −8.00

Pub#46497477 Pub#25424243 Pub#16632472 Pub#5438977 Pub#24233809 Pub#25423177 Pub#21209573 Pub#17620651 Pub#8199425 Pub#8544853 Pub#23013652

238 488 517 600 496 10 13 512 23 304 11

21 22 23 24 25 26 27

−7.15 −7.14 −7.13 −7.12 −7.11 −7.09 −7.08

Pub#2939668 Pub#4926215 Pub#9479783 Pub#5438977 Pub#1371415 Pub#51781159 Pub#1371416

348 336 524 600 356 386 355

21 22 23 24 25 26 27

−7.99 −7.97 −7.94 −7.86 −7.86 −7.85 −7.85

Pub#26714006 Pub#43911521 Pub#25423215 Pub#9410928 Pub#25423171 Pub#32260135 Pub#32260110

264 8 271 290 99 151 255

28 29 30 31 32 33 34 35

−7.07 −7.06 −7.06 −7.05 −7.05 −7.04 −7.04 −7.03

Pub#46761880 Pub#40691534 Pub#9278507 Pub#8912181 Pub#18283431 Pub#25423183 Pub#45754571 Pub#3104366

415 455 293 537 279 9 424 636

28 29 30 31 32 33 34 35

−7.83 −7.79 −7.74 −7.74 −7.72 −7.71 −7.70 −7.69

Pub#30885190 Pub#54570516 Pub#46761880 Pub#6232638 Pub#9278507 Pub#9303813 Pub#8503186 Pub#4926221

466 361 415 576 293 529 183 18

Rank

2399, Page 12 of 17

J Mol Model (2014) 20:2399

Table 5 The listed compounds were predicted to bind to the STAT5 and STAT3 SH2 domains. The compounds were ranked by their binding energies as were calculated using the MM-GBSA method

J Mol Model (2014) 20:2399

Page 13 of 17, 2399

Table 5 (continued)

Although STAT5 and STAT3 SH2 domains share high structural similarities, STAT3 SH2 domain was predicted to

be more drug targetable than that of STAT5. Table 4 lists the STAT5/STAT3 top 35 hits ranked by their binding energies as

2399, Page 14 of 17

J Mol Model (2014) 20:2399

Table 6 Energy terms contributing to the binding free energy of the three top-hits in complex with STAT3 and STAT5 SH2 domains STAT3

STAT5

Contributiona

660

561

304

660

561

304

ΔEVDW ΔGSUR ΔEEEL ΔGGB ΔGpolb ΔGnon-polc ΔHd TΔS ΔGe

−19.59±0.94 −4.36±0.07 −132.80±5.00 124.50±3.49 −8.30±0.22 −23.95±1.90 −32.25±2.11 −17.91±2.00 −14.33±0.91

−32.55±0.65 −4.50±0.05 −33.67±0.92 34.92±0.64 1.25±0.15 −37.05±0.34 −35.81±0.72 −21.09±2.18 −14.72±0.64

−28.97±0.44 −3.98±0.03 −47.17±2.13 53.66±2.09 6.49±0.15 −32.95±0.34 −26.47±0.41 −21.37±1.63 −5.10±0.45

−26.79±0.85 −4.88±0.04 −189.28±4.1 179.31±3.44 −9.97±0.19 −31.67±1.69 −41.64±1.31 −20.78±1.24 −20.85±0.57

−31.30±0.55 −4.45±0.06 −16.52±0.84 28.80±0.68 12.28±0.13 −35.75±0.34 −23.46±0.57 −20.44±1.03 −3.02±0.35

−35.83±0.98 −5.03±0.08 −15.40±1.85 24.59±1.86 9.19±0.23 −40.86±0.82 −31.67±1.05 −17.56±0.89 −14.10±0.43

a

All energies are in kcal mol−1

b

ΔGpol = ΔEEEL + ΔGGB

c

ΔGnon-pol = ΔEVDW + ΔGSUR

d

ΔH = ΔGpol + ΔGnon-pol e ΔG = ΔH – TΔS

were calculated using the AUTODOCK method. These hits shared the same sites with the pTyr-containing peptide, formed a number of hydrogen bonds with key residues surrounding the pockets and were suggested to function as dualSTAT5/STAT3-inhibitors. We selected ten candidate compounds that potentially bind to STAT5 SH2 domain with predicted affinities comparable to STAT3 SH2 domain to undergo a further screening. Table 5 shows the final ranking of the ten hits after rescoring their interactions using the MMGBSA method. A considerable number of compounds showed positive binding energies after rescoring them using the MM-GBSA method, supporting the capability of this method to discriminate inactive compounds from the AUTODOCK suggested hits. For instance, compound 33 was predicted to bind to STAT5 and STAT3 SH2 domains with a relatively lower binding energy of −7.78 and −9.06 kcal mol−1, respectively, and ranked in the top five using the AUTODOCK method (see Table 4). However, rescored with MM-GBSA method, the binding energies became positive, which were 1.33 and 1.35 kcal mol−1, respectively (see Table 5). This was mainly due to its poor enthalpy contributions and large entropic penalties in the dynamics simulations. The docking model showed a suitable conformation of compound 33’s ethyl phenyl and vinyl phenyl groups in the hydrophobic regions of pY+1 and pY-X pocket, respectively. However, MD simulations showed that the ethyl phenyl group caused significant conformational fluctuations, and the vinyl phenyl group gradually moved out of the pY-X pockets reducing the overall binding affinity. It is well known that the enthalpy associates with any possible conformational change in the protein or ligand upon binding [60] and the entropy change also mainly reflects two contributions: changes in solvation entropy and changes in conformational entropy

[61]. Moreover, compound 33 has the most torsions among all the selected compounds while the vibrational entropy penalty of ligand shows a linear relationship to the number of torsions it possessed [49]. The entropy penalty of compound 33 is ~ 5 kcal mol−1 larger than other ligands as calculated. We were interested in the top three compounds 660, 304, and 561 as shown in Table 6. The three hits shared the same sites with the pTyr-containing peptide, formed a number of hydrogen bonds with key residues surrounding the pockets and showed much lower binding energies compared with the positive controls. The binding energies for compound 660 with STAT5 and STAT3 SH2 domain were −20.85 and −14.33 kcal mol−1 respectively, much lower than that of the positive controls. For compound 304, it also showed lower binding energies of −14.10 and −5.10 kcal mol−1 compared with the positive controls. While compound 561 was suggested to function as a STAT3 specific inhibitor more than a dua l- in hibitor with a m uch strong er affinity of −14.72 kcal mol−1 to STAT3 but only −3.02 kcal mol−1 to STAT5. Given its possibility of being a specific STAT3 inhibitor, it was taken into our considerations as well. The three top hits showed stronger binding behaviors with SH2 domains of both STAT3 and STAT5, especially for the compound 660. That is mainly attributed to the favorable non-polar contributions compared with the positive controls (see Table 6). It is easily seen that the non-polar contributions were lower than −30 kcal mol−1, except for compound 660. Actually, the phenyl group as the common feature of the three inhibitors Fig. 8 The binding modes of the three top-hits in complex with the„ STAT3 and STAT5 SH2 domain. a 660/STAT3, b 304/STAT3, c 561/ STAT3, a* 660/STAT5, b* 304/STAT5, c* 561/STAT5

J Mol Model (2014) 20:2399

Page 15 of 17, 2399

2399, Page 16 of 17

made key hydrophobic interactions with residues surrounding the hydrophobic pY+1 pocket. Interestingly, compound 660 attributed its strong binding affinity to the polar energy contributions as well as the positive control S3I201. The negative value of the polar energy contribution may guarantee a higher potency of being an effective dual-inhibitor and the more preferable electrostatic interaction probably results from the polar interactions between the negatively charged carboxylate group with the electrostatically driven residues in the pY+0 site. The carboxylic acid group that is complementary to the positively charged residues contributes most to the binding energy and can be regarded as a pTyr-mimetic. The same is true for the S3I201 and its favorable polar interactions is from the negatively charged salicylic acid group with polar residues in the pY+0 site of STAT3 SH2 domain. However, compared with compound 660, the salicylic acid group seemed less active than the carboxylic acid group, which may be related to the effect of intramolecular hydrogen bonding within the salicylic acid group. Tanimoto similarity scores were also calculated for the three hits using OpenBabel, of which, Tanimoto score of compound 561 from 304 was 0.85, and Tanimoto values of compound 304 and 561 from compound 660 were 0.63 and 0.58, respectively. The three hits were structured similarly to a certain extent, while the carboxyl group of compound 660 distinguished it from the other compounds and guaranteed its high possibility of being an effective dual-inhibitor. Figure 8 shows the binding modes of the three inhibitors in complex with STAT3 and STAT5 SH2 domains. In the bindings with STAT3, compound 660 formed strong interactions with the five polar residues in the pY+0 site, which includes Lys591, Arg609, Ser611, Glu612, and Ser613. The negatively charged carboxyl group of compound 660 played an important role in forming hydrogen bonds and can be regarded as a pTyr-mimetic [22]. Due to the strength of such interactions between the oppositely charged ions, a considerable portion of the binding between SH2 domain and compound 660 were just from the pTyr-mimetic. In addition, Ile634, Ser636, Pro639 in the pY+1 pocket and Glu594 in the pY-X pocket also significantly contributed to the binding affinities. Compound 304 formed extensive contacts with residues Ser636, Val637, Glu638, Pro639 of the pY+1 site, whereas it formed few contacts with the key residues in the pY+0 site. The binding contributions for compound 561 were mainly resulted from key residues Arg609, Lys591, Ser613 of the pY+0 site and residues Ile634, Gln635, Ser636, Val637, Glu638, Pro639 of the pY+1 site. Compound 561 and 304 formed fewer hydrogen bonds with pY+0 pocket, because the nitro groups were too weak to act as proton acceptors in hydrogen bonding. In the bindings with STAT5, compound 660 also formed strong contacts with both pY+0 and pY+1 sites, in which, key residues Arg618, Ser620, Ser622, Thr628, Asn642, Lys644

J Mol Model (2014) 20:2399

contributed most of the binding affinity. In comparison, compound 561 and 304 mainly attributed their remarkable contributions to residues surrounding the pY+1 site within STAT5 SH2 domain. Due to the most favorable interactions with the electrostatically driven residues in the pY+0 site, compound 660 is regarded as the representative for the dual-inhibitors against the two targets.

Conclusions Targeting the SH2 domains of STAT5 and STAT3 may provide a novel strategy for treating the imatinib-resistant CML harboring aberrant STAT5 and STAT3 signaling, considering that SH2 domains are critical for the binding of STATs to the tyrosine-phosphorylated receptors and for the dimerization that is necessary for DNA binding and gene expression. Currently, many researchers are striving to uncover STAT3 SH2 small-molecule inhibitors and limited work has been done on targeting STAT5 SH2 domain. This suggests a need to develop novel compounds that are STAT5-specific or optimized for dual-inhibition of STAT3 and STAT5 to fully block the STAT signaling and overcome the imatinib-resistant-CML. Our results confirmed the experimental findings on the binding affinities between inhibitors S3I201, Acyl hydrazone 2 and Acyl hydrazone 3 to STAT5 and STAT3 SH2 domains, suggesting that the simulation method by combination of molecular docking and MD simulations in this study are reasonable and can be used in predicting the dual-inhibitors with good affinity and binding modes. Besides, the detailed interactions of the three positive controls with STAT5 and STAT3 SH2 domains based on structure dynamics, H-bond analysis, and energetic contributions were investigated. These simulations gave us insights into the inhibition mechanism of inhibitors targeting STAT5 and STAT3 SH2 domains. The final top hits (compound 660, 304, and 561) from our screening are well fitted within the two binding domains and are predicted to present improved inhibitory activity on the SH2 domains of STAT5 and STAT3 by comparing their estimated binding energies to the positive controls. Among the three top hits, compound 660 showed a stronger binding behavior with the SH2 domains of both STAT3 and STAT5 and exhibited the most possibility of being a real dualinhibitor against the two targets. That is mainly attributed to the existence of a carboxylic acid group which is complementary to the positively charged residues and contributes to the binding energy. In comparison, the nitro group within the three compounds that was predicted to bind in the pY+0 site turned out not strong enough to hold the interaction with residues over the entire MD simulations, suggesting that a more polar group prone to being proton acceptors will be beneficial for the binding. The phenyl group as the common

J Mol Model (2014) 20:2399

feature of the three inhibitors made key hydrophobic interactions with residues surrounding the hydrophobic pY+1 pocket, and it should be maintained in the future structural modifications. Acknowledgments The authors gratefully acknowledge a grant support from the Scientific Research Foundation for Returned Overseas Chinese Scholars and the State Education Ministry.

References 1. An X, Tiwari AK, Sun Y, Ding PR, Ashby CR, Chen ZS (2010) Leuk Res 34:1255–1268 2. Bauer S, Romvari E (2009) Clin J Oncol Nurs 13:523–534 3. Vaidya S, Ghosh K, Vundinti BR (2011) Eur J Haematol 87:381–393 4. Pene-Dumitrescu T, Smithgall TE (2010) J Biol Chem 285:21446– 21457 5. Noor SM, Bell R, Ward AC (2011) Crit Rev Oncol Hematol 78:33–44 6. O’Brien SG, Guilhot F, Larson RA, Gathmann I, Baccarani M, Cervantes F, Cornelissen JJ, Fischer T, Hochhaus A, Hughes T (2003) N Engl J Med 348:994–1004 7. Hochhaus A (2006) Ann Oncol 17:274–279 8. Mahon FX, Deininger MWN, Schultheis B, Chabrol J, Reiffers J, Goldman JM, Melo JV (2000) Blood 96:1070–1079 9. Gorre ME, Mohammed M, Ellwood K, Hsu N, Paquette R, Rao PN, Sawyers CL (2001) Sci Signal 293:876–880 10. Bixby D, Talpaz M (2010) Leukemia 25:7–22 11. Roychowdhury S, Talpaz M (2011) Blood Rev 25:279–290 12. Benekli M, Baer MR, Baumann H, Wetzler M (2003) Blood 101: 2940–2954 13. Wang Y, Cai D, Brendel C, Barett C, Erben P, Manley PW, Hochhaus A, Neubauer A, Burchert A (2007) Blood 109:2147–2155 14. Warsch W, Kollmann K, Eckelhart E, Fajmann S, Cerny-Reiterer S, Hölbl A, Gleixner KV, Dworzak M, Mayerhofer M, Hoermann G (2011) Blood 117:3409–3420 15. Bewry NN, Nair RR, Emmons MF, Boulware D, Pinilla-Ibarz J, Hazlehurst LA (2008) Mol Cancer Ther 7:3169–3175 16. Hoelbl A, Schuster C, Kovacic B, Zhu B, Wickre M, Hoelzl MA, Fajmann S, Grebien F, Warsch W, Stengl G (2010) EMBO Mol Med 2:98–110 17. Leeman RJ, Lui VWY, Grandis JR (2006) Expert Opin Biol Ther 6: 231–241 18. Darnell JE (2005) Nat Med 11:595–596 19. McMurray JS (2008) Biopolymers 90:69–79 20. Nelson EA, Sharma SV, Settleman J, Frank DA (2011) Oncotarget 2: 518–524 21. Schust J, Sperl B, Hollis A, Mayer TU, Berg T (2006) Chem Biol 13: 1235–1242 22. Hao W, Hu Y, Niu C, Huang X, Chang CPB, Gibbons J, Xu J (2008) Bioorg Med Chem Lett 18:4988–4992 23. Song H, Wang R, Wang S, Lin J (2005) Proc Natl Acad Sci U S A 102:4700–4705 24. Siddiquee K, Zhang S, Guida WC, Blaskovich MA, Greedy B, Lawrence HR, Yip MLR, Jove R, McLaughlin MM, Lawrence NJ (2007) Proc Natl Acad Sci U S A 104:7391–7396 25. Schust J, Berg T (2004) Anal Biochem 330:114–118 26. Shahani VM, Yue P, Haftchenary S, Zhao W, Lukkarila JL, Zhang X, Ball D, Nona C, Gunning PT, Turkson J (2010) ACS Med Chem Lett 2:79–84

Page 17 of 17, 2399 27. Bhasin D, Cisek K, Pandharkar T, Regan N, Li C, Pandit B, Lin J, Li PK (2008) Bioorg Med Chem Lett 18:391–395 28. Lin L, Deangelis S, Foust E, Fuchs J, Li C, Li PK, Schwartz EB, Lesinski GB, Benson D, Lu J (2010) Mol Cancer 9:217 29. Lin L, Hutzen B, Li PK, Ball S, Zuo M, DeAngelis S, Foust E, Sobo M, Friedman L, Bhasin D (2010) Neoplasia 12:39–50 30. Turkson J, Ryan D, Kim JS, Zhang Y, Chen Z, Haura E, Laudano A, Sebti S, Hamilton AD, Jove R (2001) J Biol Chem 276:45443–45455 31. Bharti AC, Donato N, Aggarwal BB (2003) J Immunol 171:3863– 3871 32. Meydan N, Grunberger T, Dadi H, Shahar M, Arpaia E, Lapidot Z, Leeder JS, Freedman M, Cohen A, Gazit A (1996) Nature 379:645–648 33. Iwamaru A, Szymanski S, Iwado E, Aoki H, Yokoyama T, Fokt I, Hess K, Conrad C, Madden T, Sawaya R (2006) Oncogene 26:2435– 2444 34. Faderl S, Ferrajoli A, Harris D, Van Q, Priebe W, Estrov Z (2005) Anticancer Res 25:1841–1850 35. Müller J, Sperl B, Reindl W, Kiessling A, Berg T (2008) Chembiochem 9:723–727 36. Tanimoto TT (1957) Nov 17:1957 37. Zhu T, Cao S, Su P-C, Patel R, Shah D, Chokshi HB, Szukala R, Johnson ME, Hevener KE (2013) J Med Chem 56:6560–6572 38. Becker S, Groner B, Muller CW (1998) Nature 394:145–151 39. Neculai D, Neculai AM, Verrier S, Straub K, Klumpp K, Pfitzner E, Becker S (2005) J Biol Chem 280:40782–40787 40. Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ (1998) J Comput Chem 19:1639–1662 41. Gasteiger J, Marsili M (1980) Tetrahedron 36:3219–3228 42. Hou T, Wang J, Li Y, Wang W (2011) J Chem Inf Model 51:69–82 43. Case DA, Cheatham TE, Darden T, Gohlke H, Luo R, Merz KM, Onufriev A, Simmerling C, Wang B, Woods RJ (2005) J Comput Chem 26:1668–1688 44. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA (2004) J Comput Chem 25:1157–1174 45. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) J Chem Phys 23: 327–341 46. Darden T, York D, Pedersen L (1993) J Chem Phys 98:10089 47. Yang Y, Shen Y, Li S, Jin N, Liu H, Yao X (2012) Mol BioSyst 8: 3049–3060 48. Turkson J, Kim JS, Zhang S, Yuan J, Huang M, Glenn M, Haura E, Sebti S, Hamilton AD, Jove R (2004) Mol Cancer Ther 3:261–269 49. Park IH, Li C (2010) J Mol Recognit 24:254–265 50. Klebe G (2009) The foundations of protein–ligand interaction In: From molecules to medicines. Springer, Heidelberg, pp 79-101 51. Wang J, Morin P, Wang W, Kollman PA (2001) J Am Chem Soc 123: 5221–5230 52. Barakat K, Mane J, Friesen D, Tuszynski J (2010) J Mol Graph Model 28:555–568 53. Page BD, Khoury H, Laister RC, Fletcher S, Vellozo M, Manzoli A, Yue P, Turkson J, Minden MD, Gunning PT (2012) J Med Chem 55: 1047–1055 54. Zhang X, Yue P, Fletcher S, Zhao W, Gunning PT, Turkson J (2010) Biochem Pharmacol 79:1398–1409 55. Fletcher S, Page BD, Zhang X, Yue P, Li ZH, Sharmeen S, Singh J, Zhao W, Schimmer AD, Trudel S (2011) ChemMedChem 6:1459–1470 56. Shi C, Yu R, Shao S, Li Y (2012) J Mol Model 19:871–878 57. Yu R, Craik DJ, Kaas Q (2011) PLoS Comput Biol 7:e1002011 58. Velazquez-Campoy A, Todd MJ, Freire E (2000) Biochemistry 39: 2201–2207 59. Li L, Li Y, Zhang L, Hou T (2012) J Chem Inf Model 52:2715–2729 60. Luque I, Freire E (2002) Proteins Struct Funct Bioinforma 49:181– 190 61. Leavitt S, Freire E (2001) Curr Opin Struct Biol 11:560–566