Molecular Determinants Underlying Binding Specificities of the ABL Kinase Inhibitors: Combining Alanine Scanning of Binding Hot Spots with Network Analysis of Residue Interactions and Coevolution

Quantifying binding specificity and drug resistance of protein kinase inhibitors is of fundamental importance and remains highly challenging due to complex interplay of structural and thermodynamic factors. In this work, molecular simulations and computational alanine scanning are combined with the network-based approaches to characterize molecular determinants underlying binding specificities of the ABL kinase inhibitors. The proposed theoretical framework unveiled a relationship between ligand binding and inhibitor-mediated changes in the residue interaction networks. By using topological parameters, we have described the organization of the residue interaction networks and networks of coevolving residues in the ABL kinase structures. This analysis has shown that functionally critical regulatory residues can simultaneously embody strong coevolutionary signal and high network centrality with a propensity to be energetic hot spots for drug binding. We have found that selective (Nilotinib) and promiscuous (Bosutinib, Dasatinib) kinase inhibitors can use their energetic hot spots to differentially modulate stability of the residue interaction networks, thus inhibiting or promoting conformational equilibrium between inactive and active states. According to our results, Nilotinib binding may induce a significant network-bridging effect and enhance centrality of the hot spot residues that stabilize structural environment favored by the specific kinase form. In contrast, Bosutinib and Dasatinib can incur modest changes in the residue interaction network in which ligand binding is primarily coupled only with the identity of the gate-keeper residue. These factors may promote structural adaptability of the active kinase states in binding with these promiscuous inhibitors. Our results have related ligand-induced changes in the residue interaction networks with drug resistance effects, showing that network robustness may be compromised by targeted mutations of key mediating residues. This study has outlined mechanisms by which inhibitor binding could modulate resilience and efficiency of allosteric interactions in the kinase structures, while preserving structural topology required for catalytic activity and regulation.

modulate resilience and efficiency of allosteric interactions in the kinase structures, while preserving structural topology required for catalytic activity and regulation.

Introduction
Protein kinases act as dynamic molecular switches in cellular signaling and their functional activity is essential for the integrity and viability of signaling pathways involved in cell cycle control, organism development, and stress response [1][2][3][4][5][6][7][8][9][10][11][12]. The human protein kinases represent one of the largest protein families that orchestrate functional processes in cellular networks and comprise an important class of therapeutic targets, owing to the presence of a highly conserved ATP binding pocket that can be exploited by small molecule inhibitors [13][14][15][16][17]. Due to evolutionary conservation of the ATP binding site and structural similarity of the protein kinase folds, most ATP-competitive kinase inhibitors can promiscuously inhibit multiple kinases. Understanding of the molecular determinants underlying binding specificities of the kinase inhibitors and the development of selective and multi-target kinase drugs with a desirable activity profile are of fundamental and practical importance and remain to be highly challenging. The continuously growing body of structural and functional studies has revealed that protein kinase activity and binding can be regulated via a dynamic equilibrium between distinct functional states: active, inactive, and Src-like inactive conformations [18][19][20][21][22][23][24]. A diverse repertoire of crystallographic conformations has also indicated that molecular switching mechanism of protein kinases may not necessarily imply an on-off binary operation (from inactive to active), but could rather represent a continuous dynamic process in which kinases may adopt a wide spectrum of inactive-like and active-like conformations exhibiting a range of activity levels. Conformational transitions between kinase states are orchestrated by three conserved structural motifs in the catalytic domain: the αC-helix, the DFG-Asp motif (DFG-Asp in, active; DFG-Asp out, inactive), and the activation loop (A-loop open, active; A-loop closed, inactive). The conserved His-Arg-Asp (HRD) motif in the catalytic loop and the DFG motif are coupled with the αC-helix to form conserved intramolecular networks termed regulatory spine (R-spine) and catalytic spine (C-spine) whose assembly and stabilization are intimately linked with the conformational transformations and kinase activation [25,26].
The equilibrium between functional kinase states can be modulated and often redistributed by activation mutations, posttranslational modifications, protein interactions, and binding of small molecule inhibitors. On the basis of the molecular mechanism of action, one can distinguish three major classes of kinase inhibitors (types 1, 2 and 3) [14][15][16][17]. Type 1 inhibitors target the catalytically competent active (DFG-in) conformation of the kinase domain, while type 2 inhibitors recognize the inactive DFG-out kinase conformation. It has been long assumed that type 1 inhibitors are less specific than type 2 inhibitors because the active conformation is very similar in most kinases. Type 3 allosteric inhibitors do not compete with ATP and tend to be more selective than types 1 and 2 inhibitors by binding to the regulatory sites outside of the ATP binding site. The discovery of Imatinib, a highly selective type 2 inhibitor targeting ABL, KIT, and PDGFR kinases by recognizing a specific inactive conformation, launched a new era in the development of the tyrosine kinase inhibitors and galvanized a stream of structure-functional studies of ABL kinase regulation [27][28][29][30]. Nilotinib is a second generation type 2 inhibitor and a 20-fold more potent ABL inhibitor than Imatinib, but both Imatinib and Nilotinib share high specificity towards ABL kinase and bind to a specific inactive conformation (DFGout/αC-helix-in, A-loop closed) [31,32]. Ponatinib is a third-generation type 2 inhibitor that circumvents drug resistance against the T315I mutation and is more potent and less selective inhibitor than both Imatinib and Nilotinib [33,34]. Despite structural similarities of the ABL kinase complexes with the type 2 inhibitors, Ponatinib inhibits activity of multiple kinases including members of the VEGFR, PDGFR, FGFR, EPH receptor and SRC families of kinases, KIT, RET, TIE2, and FLT3. Dasatinib is a type 1 inhibitor that targets the catalytically competent active ABL conformation [23] and has a broad spectrum of activity against the SRC, CSK, TEC, and EPH families of tyrosine kinases. Large-scale profiling studies have revealed a significant promiscuity of Dasatinib that can target over 30 different tyrosine and Ser/Thr kinases [35][36][37][38][39]. Crystal structures of Dasatinib complexes with ABL [23], SRC [40], LYN [41], EPHA4 [42], BMX [43], and BTK kinases [44] have supported the notion that functional promiscuity of this inhibitor may be associated with conformational tolerance to multiple kinase forms. Structural plasticity of protein kinase conformations in complexes with Dasatinib ranged from fully active DFG-in states (in ABL, SRC, LYN, EPHA4, and BTK kinases) to an inactive DFG-out conformation of BMX [43], and an intermediate DFG-out form of BTK kinase [44]. NMR studies of Dasatinib binding with ABL kinase [45] have reassured structural preferences of the inhibitor to bind an active DFG-in form, but this conformational state may not be mandatory for a productive inhibitor binding. Recent studies have revealed that Dasatinib binding is phosphorylation stateindependent and can be compatible with the DFG-out inactive ABL conformation, as the thermodynamic preferences for association with the active kinase state appeared to be only marginal [46]. Bosutinib is another type 1 inhibitor with a broad spectrum of activity against SRC and Ser/ Thr kinases [47]. The binding modes of Bosutinib and Dasatinib in complexes with ABL kinase are similar, but the kinase domain can adopt a specific conformation in response to Bosutinib binding with an inactive DFG-out conformation and an open, active conformation of the A-loop [48]. Structural and spectroscopic analysis showed that binding affinities of Bosutinib with the phosphorylated and unphosphorylated ABL kinase were virtually indistinguishable, suggesting that this promiscuous inhibitor is also conformationally permissive and can bind to both DFG-in and DFG-out kinase conformations [48]. Recently, the crystal structure of Bosutinib complex with SRC kinase revealed an active DFG-in protein conformation and a similar inhibitor binding mode to that observed in the ABL complex [49].
Comprehensive analyses of kinase inhibitor selectivity probed the interactions of many known kinase inhibitors against panels of kinase targets representing more than 80% of the human protein kinome [50][51][52]. A complex interaction pattern has unveiled a wide spectrum of promiscuity in which off-target interactions could often occur with a different kinase subfamily than the subfamily of the intended kinase target. In the kinase inhibitors developed against tyrosine kinases a spectrum of binding promiscuity was less broad, with~24% of offtargets representing Ser/Thr kinases [52]. A mass spectrometry-based proteomic approach identified the direct targets and downstream signaling effect of Imatinib, Nilotinib, Dasatinib, and Bosutinib in epidermoid carcinoma cells [53]. This study confirmed similarities in the binding promiscuity of Dasatinib and Bosutinib that displayed a broad target spectrum (19 and 18 tyrosine kinases, respectively) including members of the SRC kinase family (e.g., YES, LCK, LYN, and FRK) and members of the ephrin receptor family (EPHA2, EPHB2, EPHB3, and EPHB4 kinases). A systematic analysis of type 2 inhibitors has recently shown that high selectivity and strong sensitivity towards the unphosphorylated form of ABL kinase is not a universal characteristic of type 2 inhibitors that bind to the DFG-out conformation [54]. A survey of type 2 kinase inhibitors has concluded that many kinases may be inherently predisposed for binding with this class of small molecules since the DFG-out conformation may be thermodynamically native state for some kinases [55]. Although type 2 inhibitors are generally more selective than type 1 inhibitors, selectivity may be also achieved with a type 1 binding mode and is not guaranteed with a type 2 binding mode.
Computational studies have investigated molecular mechanisms of protein kinases and structural effects of drug resistance mutations [70][71][72][73][74][75][76][77][78][79][80][81]. Molecular dynamics (MD) simulations and free energy calculations have shown that conformational selection and the stability difference between the inactive DFG-out and active DFG-in conformations may be the primary determinant underlying Imatinib specificity [76,77]. Integration of modeling approaches and isothermal titration calorimetry experiments have quantified conformational dynamics and stability of the ABL conformations that could define the mechanism of Imatinib selectivity [78]. The absolute binding free energies of Imatinib binding with ABL and SRC kinases were computed in a large-scale simulation study [79] and scrutinized at atomic details the proposed mechanisms of Imatinib selectivity, suggesting that conformational selection mechanism and subtle differences in the interactions with the binding site residues may be collectively responsible for binding specificity. MD simulations have also explored thermodynamic preferences of Dasatinib for binding to the active ABL state, suggesting that favorable inhibitor interactions with the assembled R-spine may contribute to the stability of the active complex [82]. Mechanisms of Dasatinib and Ponatinib resistance to clinically relevant ABL polymutants were recently investigated using a combination of simulations and biochemical assays [83], showing that Ponatinib can maintain binding affinity not only for most single ABL mutants but also for some of the polymutants detected in the studied patient pool. Recently published MD simulations of Ponatinib binding with a large panel of Imatinib-resistant ABL mutants have reached similar conclusions [84]. The thermodynamic factors underlying kinase inhibitor specificities and drug resistant profiles can be evaluated by combining biophysical simulations and free energy simulations with alanine scanning mutagenesis. Computational alanine scanning approach pioneered by late Peter Kollman and colleagues [85,86] can estimate the energetic contribution of each residue to the total binding energy through systematic alanine modifications by employing the molecular mechanics (MM) force field [87] combined with the generalized Born and solvent accessible surface area (GB/SA) solvation model [88,89]. This protocol assumes that mutations would not cause significant conformational changes and global rearrangements in the residue interactions.
Complex changes in the kinase-inhibitor interactions obtained from X-ray crystallography, NMR studies, and large-scale computer simulations often reflect not only local variations at the binding interfaces, but also subtle yet important changes in the residue interaction networks. A graph-based representation of protein structures yields a convenient description of residue interaction networks [90][91][92][93], providing a robust framework for understanding allosteric communications in protein systems. Structure-based network models of protein structures often employ various measures of node centrality (degree, closeness, and betweenness) to characterize local and global connectivity in the residue interaction networks. The network centrality parameters have been robust predictors of functional residues that regulate protein-protein interactions [94,95], ligand binding [96,97], enzyme catalysis [98] and allosteric signaling [99]. Functional sites that mediate stability of the residue interaction networks are also involved in the networks of coevolving protein residues. Statistical coupling analysis (SCA), mutual information (MI) model and related covariance-based approaches have employed sequencebased analysis of residue coevolution in homologous families to show that functional residues are connected via strong coevolutionary relationships [100][101][102][103][104][105][106][107]. Coevolution of protein residues can reflect a coordinated involvement of these sites in mediating residue-residue contacts [108] thus promoting protein folding [109], facilitating protein recognition and allosteric signaling in multi-protein complexes [110], controlling an enzymatic activity and progression of disease-associated phenotypes and drug resistant variants [111,112]. Coevolving residues tend to be spatially coupled and correspond to functionally important sites exhibiting correlated and compensatory mutations in homologous proteins [113]. These residues could also form networks with connections corresponding to coevolutionary interaction strengths between nodes, where the underlying small-world topology of such networks is similar to the structurebased residue interaction networks [114][115][116]. Furthermore, coevolving residues can assemble into structurally stable and quasi-independent modules of physically interacting residues termed 'protein sectors' [103,104]. Integrated analysis of residue coevolution networks and residue interaction networks can capture the emergence of independent modules of functional residues that cooperatively mediate structural stability and conformational transitions required for diverse functions in the dynamic protein environment [117][118][119][120][121][122][123].
In the current work, we present a theoretical framework for rationalizing binding specificities of the ABL kinase inhibitors by dissecting a relationship between ligand binding and residue interaction networks in the kinase structures. The primary hypothesis addressed in this study is that kinase inhibitors could effectively intervene into organization of the residue interaction networks through binding hot spots that also function as global mediating sites in protein kinase structures. We characterize the residue interaction networks and networks of coevolving residues in the kinase structures, revealing a selected group of high centrality residues that may be involved in kinase activity, regulation and binding. MD simulations and MM-GBSA alanine scanning of binding interactions are combined with the structure-based network analysis to probe ligand-induced changes in the residue interaction networks. This analysis demonstrates that binding specificity and drug resistance effects may be adequately described by the unique networking signatures of the kinase complexes. We show that binding sensitivities of the kinase inhibitors may be exposed through energetic coupling with the hot spot residues that regulate conformational equilibrium and can modulate stability of the kinase states. Our study offers a systems-based perspective on drug design by investigating how selective and promiscuous inhibitor binding can be interrelated with the efficiency and robustness of the residue interaction networks in the kinase structures.

Conformational Dynamics of the Protein Kinase Complexes
Using MD simulations, we first characterized conformational dynamics of the ABL kinase complexes with the type 2 inhibitors (Nilotinib, Ponatinib) and type 1 inhibitor Bosutinib. Although Nilotinib and Ponatinib share a similar binding mode targeting a specific inactive ABL conformation (DFG-out/αC-helix-in, A-loop closed) (Fig 1), these drugs have a drastically different kinase specificity profiles and binding sensitivities towards Imatinib-resistant ABL mutants. In the ABL kinase complex with the type 1 inhibitor Bosutinib, the kinase domain revealed an uncharacteristic inactive position of the DFG motif (DFG-out/αC-helix-in) but retained an open active conformation of the A-loop that is similar to that observed in the active kinases (Fig 2A). However, in the Src kinase complex with Bosutinib, the kinase domain adopts an active DFG-in/αC-helix-in conformation and an open A-loop conformation that are structurally similar to Dasatinib complexes with ABL and SRC kinases (Fig 2B). In the course of simulations, we tracked conformational variations of the αC-helix and the A-loop regions, particularly focusing on structural stability and ligand-induced changes in the R-spine and Cspine intramolecular networks. The R-spine in ABL consists of M290 from the C-terminal end Conformational dynamics profiles were computed by projecting MD trajectories onto the space of three lowest frequency modes. The color gradient from blue to red indicates the decreasing structural rigidity (or increasing conformational mobility) of the protein residues and refers to an average value over the backbone atoms in each residue. The R-spine residues are annotated in spheres and colored according to their degree of structural stability. A partially disjointed architecture of the R-spine is characteristic of the inactive ABL conformation in the crystal structures. The inhibitors are shown in sticks and atom-based color-coded. The inhibitor binding modes and binding site residues are shown for Nilotinib (C) and Ponatinib (D). The highlighted αC-helix position (αC-in) and DFG-out conformation are characteristic of the inactive ABL conformation in the complexes with type 2 inhibitors. of the αC-helix, L301 from the β4-strand, F382 of the DFG motif in the beginning of the Aloop, H361 of the HRD motif in the catalytic loop, and D421 of the αF-helix (Figs 1 and 2). The C-spine is comprised of hydrophobic residues (V256, A269, L323, C369, L370, V317, L438, and I432) that connect the kinase lobes anchoring catalytically important sites to the hydrophobic C-terminus of the αF-helix (L438, I432).
To characterize functional motions of the kinase complexes we determined conformational mobility profiles using principal component analysis (PCA) [124]. In the current study, we performed PCA of protein conformational dynamics based on the backbone heavy atoms (N, Cα, Conformational dynamics profiles were computed by projecting MD trajectories onto the space of three lowest frequency modes. The color gradient from blue to red indicates the decreasing structural rigidity (or increasing conformational mobility) of the protein residues and refers to an average value over the backbone atoms in each residue. The R-spine residues are annotated in spheres and colored according to their degree of structural stability. A fully assembled architecture of the R-spine is characteristic of the active kinase conformations in the crystal structures. The inhibitors are shown in sticks and atom-based color-coded. A close up of the inhibitor binding mode and interacting residues is shown for ABL kinase (C) and SRC kinase (D). The different orientation of the regulatory kinase motifs in otherwise similar complexes with Bosutinib is highlighted: DFG-out/αC-helix-in conformation in ABL kinase (C) and DFG-in/αC-helix-in conformation in SRC kinase (D).
Cβ, C, O) and the Cα atoms only, both producing very similar profiles. The conformational mobility profiles were computed along the three low frequency modes revealing conservation of the kinase dynamics as well as structural stability of the catalytic core and the R-spine residues. Although the R-spine was partially disassembled in the ABL complexes with Nilotinib ( Fig 1A), Ponatinib ( Fig 1B) and Bosutinib (Fig 2A), such structural arrangement of the intramolecular network remained stable during simulations of the kinase complexes. This indicated that both types of kinase inhibitors may be compatible with the DFG-out conformation and a disjointed structure of the R-spine. Furthermore, similar functional dynamics of Bosutinib complexes with ABL (DFG-out conformation) and SRC kinases (DFG-in conformation) supported the notion that this type 1 inhibitor may be dynamically adaptable and relatively insensitive to the R-spine architecture (Fig 2). The conformational fluctuations based on the computed B-factors of the backbone residues were generally similar in all complexes (Fig 3). While the P-loop (residues 248-255) displayed small fluctuations in the ABL complexes with Nilotinib and Ponatinib (Fig 3A and 3B), the greater conformational flexibility of this motif was observed in the Bosutinib complexes (Fig 3C and 3D). A partially increased mobility of the P-loop in the ABL-Bosutinib complex may reflect coupling between the P-loop (Y253) and DFG motif (F382) that maintains the DFG-out conformation (Fig 3C). In the SRC-Bosutinib complex, the conformational flexibility of the P-loop (residues 273-279) further increased, while the αC-helix (residues 303-317), and A-loop (residues 403-430) displayed similar mobility to that seen in the ABL-Bosutinib complex (Fig 3D).
We also simulated conformational dynamics of Dasatinib binding using crystal structures of Dasatinib complexes with ABL [23], SRC [40], LYN [41], EPHA4 [42], BMX [43], and BTK [44], and P38 kinases. The projection of MD trajectories onto the essential space of low frequency modes revealed characteristic similarities in the distribution of structurally rigid and conformationally flexible regions of the kinase core (S1 Fig). The computed fluctuations of the backbone residues indicated that the binding site residues interacting with Dasatinib were stable (S2 Fig). Conformational dynamics of Dasatinib binding with kinase targets was consistent with the HX-MS experiments of the inhibitor binding with ABL kinase that revealed only small dynamic changes localized near the ATP binding site, while other regions in the kinase domain remain mostly unaffected [125]. Structural alignment of the ensemble-average positions for the binding site residues revealed their considerable structural conservation, which was particularly pronounced for T315, F317, M318, and L370 (S3 Fig). Structural rigidity of these residues may be contrasted with some positional variability of the DFG motif. In Dasatinib complexes with the DFG-out conformations of BMX and BTK kinases the αC-helix moved away from the catalytically competent position and the characteristic salt bridge was broken (S4 Fig). Nonetheless, the binding mode of Dasatinib was virtually unchanged in the respective crystal structures and remained stable during simulations.
Computational Alanine Scanning of the Binding Site Residues: Energetic Hot Spots of the Type 1 and Type 2 Kinase Inhibitors MD simulations were employed to perform binding free energy calculations followed by a systematic alanine scanning of the binding site residues in the studied complexes. The central objective of this analysis was to determine key functional residues and identify energetic hot spots that may contribute to the binding preferences of the inhibitors. We first evaluated the equilibrium distributions of the intermolecular contacts in simulated complexes using the Ligand Protein Contact (LPC) program [126]. Structural similarity of Nilotinib and Ponatinib binding modes produced nearly identical profiles of the inhibitor-kinase contacts (Fig 4A and 4B). Both Nilotinib and Ponatinib could form a number of stable intermolecular contacts with the P-loop residue Y253, catalytic residues K271, E286, the R-spine residue M290, the gate-keeper residue T315, the catalytic HRD triad (H361, R362), and the DFG motif (D381, F382). In the ABL-Nilotinib complex, the P-loop folds over the adenine pocket in the kinked conformation, and a group of hydrophobic residues (L248, Y253, F317, L370) encloses the pyridine and pyrimidine groups of the inhibitor. Binding free energies and alanine scanning confirmed that specific Nilotinib binding may be primarily associated with the favorable interactions made with Y253, M290, T315, F317, F359, L370, and F382 ( Fig 5A). Structural coupling between these functional sites from the P-loop, αC-β4-loop, A-loop, HRD and DFG motifs is important to ensure thermodynamic stability of the inactive ABL conformation structure and thus may contribute to the binding specificity of Nilotinib. It is worth stressing that some of the energetic hot spots anchoring specific Nilotinib binding included sites of drug resistant mutations Y253, T315, and F317. In patients who have developed Nilotinib resistance, frequently detected mutations included Q252H, Y253H, E255K/V, V299L, T315I/A, F317L, and F359C/V [56][57][58].
An important feature of Ponatinib binding is a network of hydrophobic contacts formed by the trifluoromethyl substituent on the phenyl ring with the I293, L298, V299, M290, V359, V379, A380, H361, and L354 residues ( Fig 4B). The energetic analysis suggested that multiple contacts formed by Ponatinib could provide a more balanced distribution of binding interactions since the individual alanine substitutions resulted in a moderate loss of binding energies ( Fig  5B). Consistent with the experiments [33,34], we determined that mutations of Ponatinib-interacting residues may produce only a moderate reduction in binding affinity. In particular, a relatively minor change in the binding energetics could be seen for T315A, which is consistent with the fact that Ponatinib is capable of escaping mutations at the critical gate-keeper position. At the same time, a few residues including E286, M290, L370, F317, D381, and F382 caused a . The intermolecular contacts were evaluated by applying the LPC program to the MD-based ensemble of structures and using a LPC-based classification scheme that includes 8 classes of atom types to define interactions: hydrophilic, hydrophobic, aromatic, acceptor, donor, neutral, neutral-donor, neutral-acceptor. The computed intermolecular contacts accounted for protein kinase residues and ligand atoms that interact through hydrogen bonds, hydrophobic contacts, aromatic-aromatic, and hydrophilic-hydrophobic interactions. larger change in the binding interactions, suggesting that these sites may act as the energetic anchors of Ponatinib binding to ABL kinase. These residues are the integral elements of the Rspine (M290, F382) and C-spine (L370) subnetworks. Accordingly, the energetic hot spots of these type 2 inhibitors corresponded to the key functional residues that define structural architecture of the specific ABL conformation. Nonetheless, binding sensitivities of Ponatinib may not be immediately apparent from alanine scanning and analysis of local binding interactions. Based on this evidence, we suggested that Ponatinib binding may induce global changes in the residue interaction networks.
For type 1 inhibitor Bosutinib, the total number of the intermolecular contacts with ABL was relatively moderate (Fig 4C). The favorable inhibitor interactions were primarily made with L248, C-spine residues (V256, A269, L370), and the hinge residues T315, M318. Interestingly, this type 1 inhibitor maintained only a limited number of contacts with the DFG-out . Computational alanine scanning evaluated the effect of mutations in the active site residues on binding affinity using MD trajectories of the wild type (WT) complexes and MM-GBSA calculations. The protocol involved a systematic modification of the inhibitor-interacting residues to alanine by eliminating side-chain atoms beyond C β , and measuring the effect of each mutation on binding affinity. motif and the A-loop during simulations. Despite a different orientation of the DFG motif (DFG-in) in the SRC complex, the distribution of the intermolecular contacts was mostly unchanged, with the most favorable interactions formed by the hinge residues (T338, Y340) and C-spine residues V281, A293, L393 ( Fig 4D). Hence, binding preferences of Bosutinib may be relatively insensitive to the DFG conformation and therefore better tolerate conformational fluctuations of the kinase domain. Computational alanine scanning confirmed that Bosutinib binding may be primarily determined by contact residues L248, T315, F317, and L370 (Fig 5C  and 5D). Importantly, these residues serve as energetic hot spots for both type 1 and type 2 inhibitors. These findings were also consistent with the recent experiments in which L248R and T315V showed resistance to both type 1 (Imatinib, Nilotinib, Ponatinib) and type 2 inhibitors (Bosutinib, Dasatinib) [64]. In addition, F317R mutant was moderately resistant to Imatinib and Nilotinib, but highly resistant to Dasatinib, Bosutinib and Ponatinib.
The distribution of Dasatinib-kinase contacts was fairly consistent across kinase targets, dominated by the favorable binding interactions with the hinge residues. The primary residues involved in stable ABL interactions with Dasatinib included L248, K271 (catalytic salt bridge), T315 (gate-keeper), M318 (hinge), L370 (C-spine) ( Fig 6A). In general, alanine scanning of the binding site residues across the Dasatinib complexes resulted in similarly small changes of the binding energies (Fig 6), which is markedly different from selective Nilotinib binding. The noticeable changes were only observed upon substitutions of the gate-keeper residue (T315 in ABL) that is conserved in all crystal structures and may function as a primary hot spot of Dasatinib binding. In addition, mutations of L248, K271, F317, M318, and L370 residues appeared to detrimentally affect binding energetics of Dasatinib to ABL kinase. Some of these residues (K271, M318, and L370) are important for the kinase activity, while mutations in other residues (L248, T315, and F317) may be associated with Dasatinib resistance [59][60][61]64]. Importantly, alanine mutations of the corresponding residues in other kinases led to similar appreciable reductions in the binding affinities (Fig 6). We found that Dasatinib binding to kinase targets may be primarily associated with the following energetic hot spots: I627, T669, Y470, L753 (EPHA4); L253, T319, F321, L374 (LYN); L423, T489, Y491, L543 (BMX); L408, T474, Y476, L528 (BTK). In the course of simulations, Dasatinib binding pose remained stable in all kinase complexes, confirming that the inhibitor binding energetics may be relatively insensitive to variations of the αC-helix, A-loop, the DFG motif and the R-spine residues (S4 Fig). These factors may promote conformational tolerance and encourage structural adaptability of the active kinase states in formation complexes with Bosutinib and Dasatinib.

The Residue Interaction Networks of the Kinase Structures and Binding Sensitivities of the ABL Kinase Inhibitors
Our results indicated that computational alanine scanning experiments may primarily reflect the effect of local binding interactions. However, functional role of the binding hot spots may be also linked with their global centrality in the residue interaction networks. Mutations of these residues may simultaneously alter multiple interactions, causing reorganizations in the global network and triggering shifts in the conformational equilibrium of kinase states. To probe binding sensitivities of the kinase inhibitors, we combined computational alanine scanning with the structure-based network analysis of the kinase-inhibitor complexes. We evaluated ligand-induced changes in the residue interaction networks to determine a mechanism by which selective and promiscuous inhibitors can use energetic hot spots to modulate stability of the residue interaction networks. Residue betweenness which is a global centrality measure was used as a stability probe in evaluating residue interaction networks in kinase structures. The betweenness of a node is defined as the number of shortest paths that pass through that node in the network, representing a global measure of the node contribution to the communication within the network. The residue centrality can characterize and differentiate highly connected residues that mediate stable interaction networks and allosteric communications in protein structures [94][95][96][97][98][99]. We also probed local residue environment using an energetics-based evaluation of relative solvent accessibility (RSA) [127,128]. A residue-specific RSA measure is defined as the ratio of the observed solvent-accessible surface area for a residue to the expected unfolded state value for the respective amino acid type [129].
We first examined the residue interaction networks of the ABL complexes with Nilotinib, Ponatinib and Bosutinib. For this analysis, we constructed joint distributions of the global residue betweenness and local residue flexibility measured by the computed B-factors (Fig 7A- 7C), and residue-based RSA values (Fig 7D-7F). It could be seen that locally stable residues may display a broad range of centrality values (0.4-0.12), suggesting that only a small fraction of stable residues in the kinase core may serve as global hot spots mediating stability and allosteric interactions in the ABL kinase complexes. For the ABL-Nilotinib complex, we detected a small "isolated" cluster of high centrality residues that seemed to be separated from other residues ( Fig 7D). This indicated that selective Nilotinib binding may be coupled with a small number of mediating residues. In some contrast, the residue interaction network of the ABL-Ponatinib complex revealed a denser distribution (Fig 7E). More specifically, Ponatinib binding may render a network where a significant number of residues (both buried and solvent-exposed) exhibited moderate-to-high centrality. This may render a more balanced connectivity of binding site residues that can reduce the inhibitor dependence on interactions with a few specific residues and thus promote a greater binding promiscuity. Although the crystal structures of Ponatinib and Bosutinib complexes with ABL kinase are radically different, the residue centrality distributions of these complexes were quite similar (Fig 7E and 7F). This indicated that common functional traits that dictate similar binding promiscuity of these inhibitors may be connected with the global networking properties.
To directly evaluate ligand-induced effects on the residue interaction networks, we examined the residue centrality profiles of the kinase structures in their apo and inhibitor-bound forms. For the ABL-Nilotinib complex, the binding site and R-spine residues corresponded to local maxima of the distribution, suggesting that these functional residues tend to have distinct networking characteristics (Fig 8A). We found that Nilotinib binding induced a sharply increased centrality for a subset of functional residues: Y253 (P-loop), T315 (gate-keeper), D363 (HRD motif), L370 (C-spine), R367, Y393 (phosphorylation site in the A-loop), and F382 (DFG motif) (Fig 8B). Some of these residues (R367, Y393) are not involved in direct contacts with the inhibitor, but are important for regulation and stability of the inactive kinase conformation. Importantly, several globally connected residues (Y253, T315, and F382) are also among energetic hot spots of Nilotinib binding (Y253, M290, T315, F317, F359, L370, and F382). In this mechanism, through direct contacts with the globally connected energetic hot spots, Nilotinib binding may affect allosteric interactions in the kinase domain and modulate stability of the residue interaction network in the inactive structure. According to our model, the enhanced centrality of these residues would strengthen cooperative interactions between the kinked P-loop (Y253), catalytic core (D393) and the unphosphorylated A-loop (Y393), that are central in promoting stabilization of the specific inactive state [54]. Selective Nilotinib binding could also enhance the stability and mediating capabilities of the R-spine (H361, F382) and C-spine residues (C369, L370, V371), which may additionally optimize allosteric communications between the N-terminal and C-terminal lobes in the inactive conformation. Overall, a pronounced network-bridging effect exerted by this selective inhibitor could stabilize structural environment favored by the inactive kinase form and thus shift thermodynamic preferences towards the specific kinase conformation.
We used the experimental data [67][68][69] that quantified changes in the inhibition constants of kinase drugs (Imatinib, Nilotinib, Ponatinib, Dasatinib, Bosutinib) against ABL kinase upon drug resistant mutations (S5 Fig) in our comparison with the computational results. In this analysis, we employed the nomenclature of resistant, moderately resistant and sensitive (nonresistant) mutations based on the experimental data and proposed classification [56]. A rigorous analysis of drug resistant effects would require conducting independent MD simulations for all studied mutants and their complexes. In the context of this study, we simplified the fundamental problem by suggesting that direct mapping of drug resistant sites onto the centrality profiles may help to differentiate the severity of mutations and provide a robust metric for assessing molecular source of resistance.
By mapping Nilotinib resistant sites onto the centrality profile, we found that this network parameter can distinguish between positions of highly resistant mutations (Y253H, T315I) and less sensitive or moderately resistant positions (Q252H, E255K/V, F359C/V) [56] (Fig 8C). In fact, highly resistant sites of Nilotinib binding corresponded to global mediating residues whose exceptionally high centrality is induced by the inhibitor binding. As a result, targeted mutations of these high centrality sites could disrupt allosteric coupling between functional regions, leading to the weakening and fragmentation of the residue interaction network. A strong network dependency on high centrality residues in the specific complex may explain vulnerability of Nilotinib binding to their targeted mutations and the emergence of drug resistance. It is worth noting that mutational mapping of the residue-based RSA profile (Fig 8D) was unable to distinguish between resistant and sensitive residues. In general, our results suggested that the severity of drug resistance mutations may be associated with the global mediating role of targeted residues.
Despite similar crystal structures and nearly identical inhibitor binding modes, the residue interaction networks in the ABL complexes with Nilotinib and Ponatinib complexes are different. The Ponatinib-induced changes in the residue centrality profile manifested in the increased betweenness for a significant number of kinase residues (Fig 9A). Most notably, the majority of the Ponatinib-interacting residues acquired higher centrality upon binding, leading to a broadly distributed and more uniform allocation of global mediating sites. The networkbridging effect of Ponatinib resulted in the increased centrality of many functional residues including Y253, K271, E286, M290, L298, V299, T315, F317, M318, D381, and F382 ( Fig 9B).
Accordingly, Ponatinib-induced reorganization of the interaction network could produce clusters of highly connected hot spots. Several distinct peaks of the centrality distribution (Y253, T315, and F317) also corresponded to the ABL residues associated with moderate Ponatinib resistance ( Fig 9C). Other peaks represented functionally important residues interacting with Ponatinib that can preserve tumorigenic potential of the oncogenic protein [65]. In particular, mutations of E286, M318, and D381 can abolish kinase activity and compromise protein stability. Hence, the network organization of ABL-Ponatinib complex revealed the increased number and the broader connectivity of high centrality residues. These network signatures may reduce Ponatinib dependence on interactions with specific hot spot residues and facilitate greater promiscuity in binding with other kinase targets. To illustrate how binding of type 2 inhibitors is coupled with the residue interaction networks, we mapped residues corresponding to the peaks of the centrality profile onto the inactive ABL structure (Fig 10). The θ-like shape of the interaction network provided two lines of communication that connected the ATP binding site and the substrate binding regions ( Fig  10A). Strikingly, the C-spine and R-spine represent subnetworks of the global interaction network of high centrality sites. Indeed, one side of the network could be traced from the active site to the αD-helix through the C-spine residues to the integral αF-helix and the substrate binding region. Another line of communication linked the ATP binding site with the αC-helix and via the R-spine residues to the substrate region. The θ-shaped residue network is rested on allosteric coupling between two lines of communications that is provided by structurally rigid and strategically positioned H361, D363, and Y393 residues. High centrality residues that mediate structural stability in the kinase structures often correspond to functionally important sites involved in the kinase activity, regulation and binding. We observed that selective Nilotinib binding tends to maximize complementarity and allosteric coupling with the global interaction network (Fig 10B). Nilotinib binding is coupled with the global interaction network via direct contacts with the energetic hot spots T315, M318, L370, D381, F382, and H361 residues. The interaction network of these high centrality sites may also act as a direct and optimal path for transmitting allosteric signals between functional kinase regions. By targeting the energetic hot spots and increasing centrality of these global mediating residues, Nilotinib binding may optimize the residue interaction networks and induce further stabilization of the inactive ABL conformation. Remarkably, relatively minor structural deviations in the binding mode of Ponatinib may trigger global changes in the ligand coupling with the residue interaction network that alter the binding specificity profile. Our analysis detected a weakened coupling of Ponatinib to key mediating residues T315 and F382 that control access to the hinge region and conformational preferences of the DFG regulatory motif ( Fig  10B). Collectively, these subtle changes may alleviate inhibitor dependency on the specific kinase conformation of ABL kinase and switch inhibitor preferences towards more promiscuous binding.
We also examined ligand-induced changes in the centrality profiles for the kinase complexes with type 1 inhibitors. The betweenness profile of the Bosutinib-ABL complex revealed the increased centrality of the αC-β4-loop/αC-helix region and hinge residues that correspond to the major distribution peaks (Fig 11A). The network bridging effect appeared to be significant only in these sites, which is consistent with the notion that structural stability of the αC-β4-loop/αChelix (E286, V299, M290) and the hinge region (T315, F317) are vital for the active ABL conformation. The close-up of the centrality profile further highlighted the network-bridging effect in E286, V289, M290, V299, T315, and F359 residues, indicating that the assembly of the R-spine and active position of the αC-helix residues are central in mediating stability of the Bosutinib complex ( Fig 11B). These residues also included Bosutinib resistant sites (V299, T315, F317, and F359), showing that residue centrality could differentiate between inhibitor-resistant and inhibitor-sensitive residues (Fig 11C). Hence, ligand-induced network bridging effects may be linked with binding sensitivities of Bosutinib towards mutational variants of ABL kinase.
Using MD simulations of Dasatinib complexes with different kinases, we aggregated the results of network analysis and analyzed coupling between residue betweenness and local structural parameters (Fig 12A and 12B). While low centrality residues are typically flexible and solvent-exposed, we detected a dense area of the medium-betweenness residues that have a similar degree of local mobility as the high centrality sites. The centrality distribution density revealed the characteristic small-world topology, in which the high centrality residues in Dasatinib complexes formed a relatively small population (Fig 12C). The number of highly connected central nodes seemed to exponentially decay as the residue centrality increased. In agreement with network-based studies of protein structures [93], the distribution of global mediating residues may obey the Poisson model, as the residue nodes that significantly deviated from the average degree were rare. We also mapped out the centrality distribution of the Dasatinib-interacting residues that is shallow and shifted towards the low and medium betweenness values (Fig 12D). It implied that the network-bridging effect of Dasatinib binding was relatively minor since most of the binding site residues displayed an average betweenness similar to the rest of the kinase molecule. The observed effect was uniform across all Dasatinib-kinase complexes (S6 Fig). Notably, the only energetic hot spots with high centrality were the conserved residues in the hinge region (T315 and M318 in ABL) and C-spine (L370 in ABL). Hence, binding sensitivities of Dasatinib may be primarily controlled by the identity and structural conservation of the energetic hot spots in the hinge region (gate-keeper T315 and M318), and could be also coupled to the C-spine network (via interactions with L370). The results are consistent with recent observations pointing to strong evolutionary conservation of the gate-keeper residue in the kinases that bind Dasatinib, since among 94 kinases with the conserved threonine at the gatekeeper position~60 kinases could be targeted by Dasatinib [39].
We also mapped residues corresponding to the centrality peaks in the complexes with type 1 inhibitors onto the active ABL conformation (Fig 13). The key functional residues occupy strategic central positions in this network including catalytic residues (K271 and E286), T315, M318 (hinge region), HRD and DFG motifs as well as the substrate binding motif 405-WTAPE-409. The θ-like shape of the residue interaction network appeared to be a conserved feature of the inactive and active kinase structures. In both cases, the C-spine and R-spine residues contribute to the global interaction network. Noteworthy is coupling of Dasatinib binding with the subnetwork that links the binding site via the C-spine to the αF-helix, and the substrate binding region (Fig 13). At the same time, Dasatinib binding may be more tolerant to conformational variations in the other communication line that connects the active site with the R-spine and DFG residues. Hence, binding of type 1 inhibitors can incur relatively moderate network changes and allow for a greater conformational tolerance of the active kinase states.

The Mutual Information Network of Coevolving Residues and Binding Sensitivities of the Kinase Inhibitors
In this section we investigated how binding preferences of the kinase inhibitors may be linked with the networks of coevolving residues in the kinase structures. Prediction of strongly coevolving residue pairs can help to clarify molecular mechanisms underlying kinase function and binding. We characterize coevolved residue positions that are associated with conformational stability, catalytic and ligand binding function. The underlying premise of this analysis is that strongly coevolving kinase residues are often located near active sites and coevolve with many other residues, thus forming interaction clusters that modulate stability of the kinase conformations. We constructed mutual information networks of coevolving residues in protein tyrosine kinases and mapped residues with the high coevolutionary signal onto ABL kinase structures. In this analysis, we compared the organization and topology of the coevolutionary residue network and the residue interaction network. According to our hypothesis, a subset of kinase residues exhibiting strong coevolutionary signal may be also associated with high centrality and correspond to energetic hot spots of kinase inhibitors. By using mutual information (MI) approach [130] we evaluated coevolutionary relationships between position pairs in the tyrosine kinase family. This provided a sufficient accuracy of the sequence analysis and identified positional correlations that guided subsequent structural mapping onto the kinase conformations and analysis of functionally important residues. We computed residue-based MI score that characterizes the extent of mutual information shared by a given residue with other protein residues. Based on MI scores and using MISTIC server [130], a network of coevolutionary residues was constructed where nodes are residues and links between nodes represent a coevolutionary signal between residue pairs. To analyze the network of coevolving residues, we explored the following parameters: the Kullback-Leibler conservation score (Fig 14A), the number of coevolving residue interactions per residue (Fig 14B), a cumulative mutual information score (cMI) for each residue (Fig 14C), and the proximity mutual information score (pMI), which measures the mutual information shared by the residues in the proximity of a given residue within a certain distance threshold (Fig 14D).
The conservation score ( Fig 14A) and cMI profile (Fig 14C) were fairly similar, both revealing that highly coevolving residues were primarily localized in the catalytic loop (residues 357- The residue-based pMI score. The coevolutionary parameters were calculated using MISTIC approach [130]. The coevolutionary parameters were mapped onto ABL kinase and residue annotation is consistent with residue numbering in ABL complexes. The coevolutionary parameters are highlighted for the R-spine residues (green circles), C-spine residues (red squares), and Imatinib-resistant mutations (maroon diamonds).
doi:10.1371/journal.pone.0130203.g014 372 in ABL) and substrate binding site (residues 399-416 in ABL). This loop includes the critical 361-HRD-363 motif that is involved in kinase activity and regulation. The catalytic loop sequence 363-DLAARN-368 is also highly conserved throughout the majority of protein tyrosine kinases and includes a critically important for kinase activity R367 [131,132]. A strong correspondence between coevolving residue positions and conserved residues in the catalytic and substrate binding sites is consistent with the previous studies [112][113][114][115][116][117][118][119][120][121][122][123]. The pMI profile evaluated the propensity to form modules of coevolving residues, which revealed that dense clusters tend to be formed in the proximity of the catalytic and substrate binding sites (Fig 14D). Hence, residues that are functionally critical for kinase activity coevolve and coevolutionary relationships often include pairs of residues from allosterically communicating active sites. Moreover, catalytic residues and their structural neighbors tend to coevolve with each other and form independent structural modules stabilized by the physical interactions [120,121]. To facilitate functional analysis of coevolving residues, we also mapped positions of the R-spine and C-spine residues onto MI profiles (Fig 14). The network connectivity of the coevolving residue pairs is illustrated by a sequential circular representation of the MI scores (Fig 15). The connections in the center of the circle link pairs of positions with significant MI score. We also mapped highly coevolving groups of residues onto the inactive and active kinase structures to explore a relationship between organization of these networks and inhibitor binding (Fig 15).
The rigid R-spine residues H361 (HRD), D382 (DFG), and D421 (αF-helix) represented conserved and highly coevolving positions in all four distributions (Fig 14). Similarly, highly coevolving C-spine residues are either proximal to the catalytic region (L370, V371) or belong to the highly conserved αF-helix (L428, I432). The high conservation of the αF-helix stems from its central role in anchoring the R-and C-Spines subnetworks with the rest of the kinase core. Notably, the decisive role of the R-spine residues was especially striking in the conservation score (Fig 14A) and pMI distributions (Fig 14D), showing that modules of highly coevolving residues can link the highly conserved and rigid catalytic residues with the more flexible residues (such as R-spine M290 and L301 from the αC-helix) that facilitate conformational changes in the kinase domain. The pMI profile revealed clusters of coevolving residues in the ATP binding site, catalytic loop, the αF-helix, and the substrate binding site that have distinguishable network signatures and may have evolved independently (Fig 14D). This finding is consistent with the notion that most of the coevolving residues are spatially proximal in the protein structures and tend to form relatively independent protein sectors of functional residues [103,104,117,118]. However, the dense network connectivity demonstrated that coevolved clusters could emerge at different distance regimes, including coupling between distant functional modules such as the ATP binding site (hinge residues 313-322), αC-helix (residues 279-292), and A-loop (residues 378-415) (Fig 15). The appearance of coevolutionary clusters that include mobile residues in the αC-helix and A-loop is likely associated with dynamic requirements for global rearrangements. In particular, we observed that modules of coevolving residues can connect the ATP binding site residues with the αC-β4-loop/αC-helix hinge region and DFG motif of the A-loop. The formation and disruption of residue interaction networks near the hinge site is a significant driver of conformational transitions between the kinase states.
The network of coevolving residues presented a characteristic θ-shaped topology, when projected onto the kinase structures, with two routes linking the ATP binding sites and the substrate binding via the R-spine and C-spine (Fig 15). We observed that high centrality and highly coevolving kinase residues often coincide and correspond to the conserved functional sites mediating stability and allosteric communications in the kinase domain. These residues included the gate-keeper T315 from the ATP site, HRD motif of the catalytic loop, and substrate binding motif 405-WTAPE-409 in the C-lobe. Most notably, the coevolutionary distributions pointed out to a group of residues that may simultaneously display a strong coevolutionary signal, high centrality and present a binding hot spot for the ABL kinase inhibitors (Figs 14 and 15). In particular, high centrality Y253 and T315 residues showed high conservation and cMI score as well as a large number of coevolving interactions with other residues in different regions (Fig 14A-14C). Several other global mediating residues (H361, F382) exhibited similarly strong coevolutionary signatures. The integrating role of these residues in coevolutionary networks and residue interaction networks necessitates their high rigidity, which is determined by catalytic and structural stability requirements. Hence, binding preferences of Nilotinib may arise from energetic coupling with rigid networks of these highly connected residues that define stability of the specific kinase states (Fig 15B and 15C). Among energetic hot spots of Dasatinib binding, the gate-keeper T315 was the only functional residue .The third and fourth circles show the pMI (pMI) and cMI values as histograms, facing inwards and outwards respectively. Each node represents a position in the MSA and lines between nodes in the circle connect pairs of positions with MI > 6.5 as defined by MISTIC program [130]. Node color represents site conservation: red for highly conserved sites and blue for less conserved sites. that displayed a markedly high centrality in interaction networks and exceptionally high coevolutionary propensity (Fig 14). Accordingly, Dasatinib coupling with the networks of coevolving and interacting kinase residues may be almost exclusively controlled by the identity and structural stability of the gate-keeper residue. As a result, structural stability requirements for productive Dasatinib binding may be relatively relaxed and rely primarily on rigidity of the hinge region near the gate-keeper residue and integrity of the C-spine architecture. Mapping of coevolutionary residues onto the active ABL structure illustrated how Dasatinib can communicate with this network (Fig 15D and 15E). The inhibitor is coupled with the hinge region and a subnetwork that links the ATP binding site with the substrate region via the C-spine. We observed that important drug resistant sites could be identified based on high coevolutionary scores and a large number of coevolving residue interactions. Of special interest, the emergence of Y253 and T315 residues as highly coevolving and high centrality positions in the network distributions. This confirms the unique role of the P-loop and the gate-keeper residue in modulating kinase activity and binding. The vulnerability of the interaction networks in the ABL complexes to targeted perturbations of these residues may explain why T315I mutation and Ploop mutations could leverage a local perturbation to switch the global conformational equilibrium away from the specific inactive conformation and stabilize the active form [133].

Discussion
In this work, we have introduced a computational framework for explaining binding preferences of the ABL kinase inhibitors by unraveling a relationship between energetic hot spots of ligand binding and organization of the residue interaction networks in various kinase states. The approach is rooted in the energy landscape theory and invokes conformational selection mechanism, according to which the relative populations of preexisting kinase states can be modulated by inhibitor binding through selection and preferential stabilization of a particular functional form [134][135][136]. The main assertion of this study was that kinase inhibitors could intervene into organization of the residue interaction networks through binding hot spots that also function as global mediating sites in protein kinase structures. In this mechanism, by using local perturbations in the binding site, selective inhibitors may facilitate a global population shift of conformational ensembles towards a specific bound form. Structure-based network analysis has recapitulated various experimental observations by linking residue centrality and structural stability of the residue interaction networks with kinase-specific binding preferences to drug binding.
Modeling of the residue interaction networks and networks of coevolving residues in the ABL-drug complexes has determined a group of conserved functional sites that can simultaneously embody strong coevolutionary signal, high centrality, and the propensity to be energetic hot spots of inhibitor binding. These residues corresponded to functionally critical sites that regulate kinase activity, determine protein structure stability and control inhibitor binding. By probing ligand-induced changes in the residue interaction networks, we have found that inhibitor binding can change the local connectivity of the binding hot spot residues and rewire the global interaction network in the kinase complexes. In this mechanism, binding of selective inhibitors, such as Imatinib and Nilotinib, could enhance centrality of the hot spot residues that mediate structural environment favored by the specific inactive state. More promiscuous inhibitors Bosutinib and Dasatinib appeared to impose modest changes in the network organization, in which ligand binding is strongly coupled only with the identity and stability of the gatekeeper residue. These findings may have interesting implications, particularly in light of recent experiments suggesting that selectivity of kinase inhibitors is not associated with a specific type but may rather result from targeted exploitation of the kinase conformational landscape [55].
The results of this study may have a particular relevance in relating differences in the network organization of kinase complexes with drug resistance effects. A strong relationship between residue positions subjected to drug resistant mutations and their unique topological signatures emerged from the analysis of the residue interaction networks and networks of coevolving residues. We have found that residue centrality may be used as a robust metric for assessing severity of drug resistance mutations and differentiating between highly resistant and moderately resistant positions. Although many highly coevolving kinase residues can be associated with drug resistance, coevolutionary relationships could often hinder subtle differences in mutational sensitivities of binding site residues in response to drug binding. A consensus analysis that identifies energetic hot spots exhibiting both high coevolutionary and high centrality scores in the residue networks may present a simple and potentially more useful strategy to quantify drug resistant and disease-associated positions in protein kinases.
A combination of molecular simulations and network approaches could provide a useful perspective on structure and energetics of ligand binding via lenses of network efficiency and robustness to external perturbations and mutations. Systems properties such as modularity, bow-tie architectures, degeneracy and other topological features are often associated with robustness of protein systems [137]. Robustness is also supported through functional redundancy and response diversity. Our study has shown that selective and promiscuous inhibitor binding can be interrelated with the robustness of the residue networks in kinase structures. A resilient protein structure network should be able to cope with random and targeted attacks by providing efficient connectivity and interaction cooperativity required for proper function and signaling. According to our findings, binding selectivity in the ABL kinase complexes may be associated with the emergence of specific high betweenness nodes that can rapidly transmit allosteric signal in the network. However, this may result in the reduced robustness of the system and the increased dependency on a small number of critical mediating residues, making specific drug binding vulnerable to targeted mutations. In network terms, drug resistant mutations in ABL may be regarded as a protein response to regain its robustness against the drug through diversity of mutations that maintain the activity of tumor cells [137].
Robustness of complex systems may be drastically improved by either rewiring a fraction of the edges or adding a modest number of new edges at the suitable location of the network [138][139][140]. It has been discovered that a limited rewiring can create robust networks in which highly connected nodes are localized in the core surrounded by interconnected layers of less connected nodes [141]. This can improve network resilience by redistributing the network load among many high betweenness nodes without compromising the network topology and decreasing robustness against random failures [141,142]. Considering protein targets as robust interaction networks exhibiting points of fragility may provide insights for the development of new drugs using systems-based approaches [143,144]. The association of network properties with kinase regulation and binding suggests that residue interaction networks may be reorganized and specifically tailored through therapeutic intervention that is informed by knowledge of high centrality residue nodes. Structure-based network approach may thus provide a useful framework for rational design of therapies that target less essential nodes to increase synergetic performance and decrease side effects [145,146].
System-based approaches to drug discovery, such as network pharmacology, have exploited advances in chemical biology and network science to develop strategies for multitarget drugs and targeted drug combinations in the context of structural and biological networks [147][148][149]. These combinatorial therapies are considered as a promising alternative to combat drug resistance and side effects that arise from compensatory mechanisms and robustness of cellular systems to external perturbations. Rational design of drug combinations and multi-target drugs requires development of novel theoretical and experimental approaches that could bridge microscopic analysis of protein structure networks with macroscopic modeling of cellular networks and signaling pathways [137]. Dissecting the relationships between protein robustness, specificity, and drug binding within a unified framework of protein structure networks may prove to be a useful step in this direction.

MD Simulations and PCA Computations
Combinations of 200 ns MD simulations were carried out for the crystal structures of the kinase-inhibitor complexes studied in this work. For each structure, starting from the same initial conditions, we performed a single 200 ns simulation and 4 independent 50 ns simulations. The collected statistics used in the computational analysis was based on all MD trajectories. The following crystal structures were simulated: the inactive ABL structure bound with Nilotinib (pdb id 3CS9) and Ponatinib (pdb id 3OXZ); the active-like ABL structure bound with Bosutinib (3UE4), the active SRC structure bound with Bosutinib (pdb id 4NMXO); Dasatinib complexed with the active structures of EPHA4 (pdb id 2Y6O); LYN (pdb id 2ZVA), P38 (pdb id 3LFA), BTK (pdb id 3K54); Dasatinib binds with the inactive structures of BMX (pdb id 3SXR) and BTK kinases (pdb id 3OCT). An active DFG-in conformation is adopted in Dasatinib complexes with ABL, SRC, EPHA4, LYN, P38, and BTK kinases (pdb id 3K54). Dasatinib binds to a DFG-out conformation of the inactive BMX kinase (pdb id 3SXR). In the complex with an inactive BTK conformation (pdb id 3OCT) an intermediate DFG position is adopted which is between the fully DFG-in and DFG-out conformations. All crystal structures were obtained from the Protein Data Bank (RCSB PDB www.rcsb.org) [150]. All crystallographic water molecules were removed and missing hydrogen atoms of the protein were added. All ionizable residues were considered in the standard ionization state at neutral pH condition. The missing residues, unresolved structural segments and disordered loops were modeled and evaluated with the ModLoop server [151,152] and the ArchPRED server [153]. The unresolved portions of the P-loop and A-loop in the kinase complexes were assembled and energetically refined using the ArchPRED server.
MD simulations were carried out using NAMD 2.6 [154] with the CHARMM27 force field [155,156] and the explicit TIP3P water model. The employed MD protocol is consistent with the overall setup described in details in our earlier studies [157]. The initial structures were solvated in a water box with the buffering distance of 10 Å. MD simulations were run in the NPT ensemble at 1atm and 300K using extended system pressure algorithm [158] and Nosé-Hoover thermostat [159,160]. Long-range nonbonded van der Waals interactions were treated using an atom-based cutoff of 12Å with switching van der Waals potential beginning at 10Å. The smooth particle mesh Ewald (PME) method [161] was employed to treat the long-range electrostatics. Numerical integration was performed using the leap-frog Verlet algorithm with 2fs time step [162]. Covalent bond lengths involving hydrogen were constrained using the SHAKE algorithm [163].The following protocol preceded the production stage of MD simulations. All atoms of the complex were first restrained at their crystal structure positions with a force constant of 10 Kcal mol -1 Å -2 . The system was subjected to the following simulation annealing to ensure the proper equilibration. The temperature was increased from 0K to 500K at a rate of 1K per 1ps and was kept at 500K for 500ps. The temperature was then decreased from 500 K to 300K at a rate of 1K per 1ps and was kept at 300K for additional 500ps. An NPT production simulation was then run on each of the equilibrated structures keeping the temperature at 300 K and constant pressure of 1 atm.
To obtain charge parameters for the studied kinase inhibitors, the conformations of Nilotinib, Ponatinib, Bosutinib and Dasatinib were extracted from the respective crystal structures of the ABL complexes and then fully optimized at RHF/6-31++ level using the Gaussian 03 package [164]. The initial charge parameters for ligand atoms were obtained using the restrained electrostatic charge fitting procedure [165] from the electrostatic potential of RHF/6-31++ single point calculations on the optimized ligand structures. To ensure consistency of the ligand parameters with the CHARMM force field methodology, we also employed a VMD (Visual Molecular Dynamics) [166] plugin, ffTK (http://www.ks.uiuc.edu/Research/vmd/plugins/fftk/) that facilitates generation of CHARMM-compatible parameters for small molecules [167]. In this protocol and according to best CHARMM force field practices, partial atomic charges were optimized from water-interaction profiles to reproduce QM interactions with a TIP3P water molecule [168]. The adopted procedure follows an original approach for derivation of small molecule parameters in the CHARMM General Force Field (CGENFF) that covers a range of chemical groups in drug molecules, including heterocyclic scaffolds present in the studied ligands [169][170][171][172]. The bond lengths, bond angles, dihedral and van der Waals parameters along with the associated force constants were also consistent with a detailed force field parameterization of Imatinib [76].
PCA of the MD conformational ensembles was based on the set of backbone heavy atoms (N, Cα, Cβ, C, O) and on the Cα atoms only to determine the essential dynamics of the protein systems. The calculations were performed using the CARMA package [173].The frames are saved every 5 ps, and a total of 10,000 frames were used to compute the correlation matrices for each simulation. For comparison, we also employed the elastic network model (ENM) and computed ENM-based lowest normal modes using the Anisotropic Network Model web server [174].

Binding Free Energy Calculations
The binding free energy if inhibitor-kinase binding was calculated using MM-GBSA approach [87][88][89]. In this approach the binding free energy ΔG bind is written as the sum of the gas phase contribution ΔG MM , the solvation free energy ΔG solv , and an entropic contribution -TΔS The brackets <> denote an average of these contributions calculated over the MD trajectories. The gas-phase contribution < ΔG MM > to the binding free energy is the difference in the molecular mechanics energy of the complex and the isolated protein and ligand. This contribution is the sum of the differences in the internal energies ΔE intra , the van der Waals interaction energy ΔE vdw , and the electrostatic interaction energy ΔE elec : where E bond is the energy of the bonded terms (bonds, angles, dihedral angles, and improper angles) of a given molecule; E vdw is the van der Waals energy of the molecule; and E elec is the electrostatic energy of the molecule. These contributions are calculated according to the CHARMM22 molecular mechanics force field.
The solvation free energy ΔG solv is the difference between the solvation energy of the complex and solvation free energies of the isolated protein and ligand: The solvation free energy of a molecule is given as the sum of nonpolar and polar contributions. The nonpolar contribution is computed using the solvent accessible surface are (SASA) model and give as DG np solv ¼ s Ã SASA where the parameter σ = 0.0072 kcal/ (mol Ã Å 2 ). The electrostatic contribution to the solvation free energy DG elec solv was calculated using the analytical generalized Born (GB) model implemented in CHARMM. This model is known to accurately reproduce the solvation free energies calculated by solving the Poisson equations. All energy terms were calculated for 10,000 frames regularly separated by 20 ps along the 200ns trajectory performed for the complex.
The entropy contribution consists of translational ΔS trans , rotational ΔS rot and vibrational ΔS vib components: The vibrational entropy terms were computed using normal mode analysis that yields better convergence than the quasiharmomic analysis from MD trajectories. The VIBRAN module of the CHARMM program was used to calculate and diagonalize the force constant matrix for the normal mode vectors and frequencies determination. The normal modes were calculated on minimized average structures obtained from MD simulations. The minimization was performed using the Newton-Raphson minimization algorithm, using the same cutoff scheme and constraints as for the normal mode calculations. Since the normal mode analysis is computationally demanding, −TΔS was averaged over only 500 frames extracted from the 200 ns MD trajectories. All energy terms are calculated using single 200 ns trajectories of the kinase-inhibitor complexes. This is followed by separation of the complexes into isolated protein and ligand structures and subsequent minimization of the isolated molecules without conducting additional simulations of the individual protein. With this simplified approach, the difference in the internal energy of the isolated molecules upon complexation is neglected and generally yields better convergence of binding free energies.

Computational Alanine Scanning
Computational alanine scanning was performed by replacing the side chain of a given residue by an alanine and recalculating the absolute binding free energy for the mutated system [85,86]. The difference in the binding free energy of the wild type and alanine mutant ΔΔG bind may be evaluated as follows: The binding free energy of the alanine mutant is calculated using the MM-GBSA approach using the snapshots obtained from MD simulations of the WT complex. The side chain of the residue under investigation is truncated, replacing the Cγ atom by hydrogen. All energy terms were calculated for 10,000 snapshots along the 200ns trajectory performed for the WT complex. For each of these snapshots of the WT complex, the mutated side chain was minimized under the fixed position of the remaining system using 1000 steps of steepest decent and Newton-Raphson minimization before calculating the energy terms. A central assumption of computational alanine scanning approach is that mutations would introduce only local structural perturbations of the system that are sufficiently moderate that the effect on the binding free energy may be gleaned from MD simulations of the WT system. Due to computational cost, the contribution of the entropy term to the binding free energies upon alanine mutations was neglected. Although it is possible to run separate MD trajectories on both WT and mutant complexes, these simulations should be more warranted when the crystal structure of the alanine mutant is available [175]. In this case, the absolute ΔG values can be calculated for each complex by averaging the results over separate MD trajectories and their difference can be effectively considered as ΔΔG bind . These simulations performed for WT and mutant complexes should represent global structural changes and describe a more rigorous effect of a give residue on binding affinity. Given a significant number of studied complexes and an attempt to estimate local energetic changes and determine energetic hot spots, we opted out for the former approach of using only MD simulations of the WT complexes.

The Residue Interaction Networks and Topological Parameters
A graph-based representation of proteins was used in the protein structure network analysis, where amino acid residues were considered as nodes and edges correspond to the nonbonding residue-residue interactions. The pair of residues with the interaction strength I ij greater than a user-defined cut-off (I min ) are connected by edges and produce a protein structure network graph for a given interaction strength I min . The strength of interaction between two amino acid side chains is evaluated as follows: where n ij is number of distinct atom pairs between the side chains of amino acid residues i and j that lie within a distance of 4.5 Å. N i and N j are the normalization factors for residues i and j respectively [176,177]. The number of interaction pairs including main-chain and side-chain made by residue type i with all its surrounding residues in a protein k is evaluated. Similar to the arguments presented in our earlier studies [157], we considered any pair of residues to be connected if I min was greater than 3.0%. We treat protein-ligand complexes as interaction networks in which the nodes of the network are formed by both amino acid residues and ligand atoms. In the network model, the binding of a ligand introduces new edges in the protein network. A weighted network representation of the protein structure is adopted that includes non-covalent connectivity of side chains and residue cross-correlation fluctuation matrix [178]. In this model of a protein network, the weight w ij of an edge between nodes i and j is measured as w ij = − log(|C ij |) where C ij is the element of the covariance matrix measuring the cross-correlation between fluctuations of residues is i and j obtained from MD simulations. The shortest paths between two residues are determined using the Floyd-Warshall algorithm [179]. Network calculations were performed using the python module NetworkX (http:// networkx.github.io/). Using the constructed protein structure networks, we computed the residue-based betweenness parameter. The betweenness of residue i is defined to be the sum of the fraction of shortest paths between all pairs of residues that pass through residue i: where g jk denotes the number of shortest geodesics paths connecting j and k, and g jk (i) is the number of shortest paths between residues j and k passing through the node n i . Residues with high occurrence in the shortest paths connecting all residue pairs have a higher betweenness values. For each node n, the betweenness value is normalized by the number of node pairs excluding n given as (N -1)(N -2) / 2, where N is the total number of nodes in the connected component that node n belongs to. The normalized betweenness of residue i can be expressed as follows: g jk is the number of shortest paths between residues j and k; g jk (i) is the fraction of these shortest paths that pass through residue i.

Mutual Information Networks of Coevolving Protein Residues
Analysis of coevolving residues was carried out using mutual information (MI) between two positions in the multiple sequence alignment (MSA), which reflects the extent to which knowing the amino acid at one position can predict the amino acid identity at the other position. MI was calculated between pairs of columns in the MSA using MISTIC approach and web server [130]. MI is a nonlinear statistic that measures the information between two random and discrete variables. Given that any two positions in the MSA can be considered random variables x and y, the MI between is these two positions is given by the relationship: where P(a i , b j) is the frequency of amino acid a occurring at position i and amino acid b at position j of the same sequence. P(a i ) is the frequency of amino acid a occurring at position i and P (b j ) is the frequency of amino acid b at position j. The amino acid frequency pair P(a i , b j) is calculated as N(a i , b j ) / N where N(a i , b j ) is the number of times that an amino acid pair (a i , b j ) is observed at positions i and j respectively. N is the total number of sequences in the MSA. MI networks were defined as graphs where nodes were defined as positions in the MSA of the protein tyrosine kinase family and edges were defined between any pair of nodes of MI. MSA were obtained from Pfam [180][181][182], a database of protein families that includes their annotations and MSA information generated using hidden Markov models. The human tyrosineprotein kinase ABL1 sequence (UniProtKB accession id P00519) matches Pfam-A entry Pkina-se_Tyr (PF007714). The seed alignment of the Pkinase_Tyr family (145 sequences) was used as the curated alignment which contains a set of representative sequences of the tyrosine kinase family, from which a profile hidden Markov model (HMM) is generated using the HMMER3 program (http://hmmer.janelia.org/). Each profile is then searched against a primary sequence database based on UniProtKB [183,184] to create the full MSA profile (24806 sequences). The criteria parameters for generating a MSA are the E-value and the number of columns in the multiple-sequence alignment for which sufficient sequences can be found to infer evolutionary couplings. Based on these criteria, for sequences in the dataset, an E-value of 10 −2 and a column-inclusion threshold of 80% were used in the MSA generation. All sequences in the full MSA score above curated thresholds are included in the full alignments for that family. A statistically significant and diverse number of sequences along with a high quality of the MSA in the Pfam database is an important prerequisite for the calculation of MI. In the MISTIC approach, MSAs with more than 400 sequences and less than 62% identity typically yield good predictive performance values [105].
In the colevolutionary network analysis, we explored the following parameters: the Kullback-Leibler conservation score, the number of coevolving residue interactions per residue, a cumulative mutual information score (cMI), and the proximity mutual information score (pMI). The Kullback-Leibler (KL) divergence score KLConsScore is used for measuring sequence conservation in protein tyrosine kinase sequences. For each column of the MSA, the KL conservation is calculated according to the following formula: Here, P(i) is the frequency of amino acid i in that position and Q(i) is the background frequency of the amino acid in nature calculated using an amino acids background frequency distribution obtained from the UniProtKB database [183,184]. CMI is a mutual information score per residue position that characterizes the extent of mutual information in its physical neighborhood. This sequence-based parameter defines an extent to which a given amino acid residue contributes to a mutual information network. CMI is calculated as the sum of MI values above a certain threshold (MI > 6.5) for every amino acid pair in which a particular residue of interest appears. The MI threshold of 6.5 was shown to be an adequate and reliable lower boundary of mutual information interactions [120,121,130]. pMI score for each position is defined as the average of cMI scores of all the residues within a certain physical distance from a given residue when mapped on the protein kinase structure. The distance between each pair of residues in the structure was calculated as the shortest distance between any two atoms, other than hydrogen atoms, that belong to each of the two positions. The threshold distance of 5 Å defines structural proximity of each residue in defining pMI score [130]: A sequential circular representation of the MI network implemented in the MISTIC approach and is provided by the web server [130]. The labels in the first (outer) circle correspond to the alignment position and the amino acid code of the reference sequence. The colored square boxes of the second circle indicate the MSA position conservation (highly conserved positions are in red, while less conserved ones are in blue).The third and fourth circles show the pMI and cMI scores as histograms, facing inwards and outwards respectively. The connected lines in the center of the circle link pairs of positions with MI > 6.5. ). An active DFG-in conformation is adopted in Dasatinib complexes with EPHA4 (pdb id 2Y6O), LYN (pdb id 2ZVA), P38 (pdb id 3LFA), and BTK kinases (pdb id 3K54). Dasatinib binds with the inactive structures of BMX (pdb id 3SXR) and BTK kinases (pdb id 3OCT). In the complex with an inactive BTK conformation (pdb id 3OCT) an intermediate DFG position is adopted which is between the fully DFG-in and DFGout conformations. Conformational dynamics profiles were computed by projecting MD trajectories onto the essential space of the three lowest frequency modes. The color gradient from blue to red indicates the decreasing structural rigidity (or increasing conformational mobility) of the protein residues and refers to an average value over the backbone atoms in each residue. The R-spine residues are annotated in spheres and colored according to their degree of structural stability. Dasatinib binding modes are shown in sticks and atom-based color-coded. The Pymol program was used for visualization of the protein kinase structures (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, and LLC). In Dasatinib complexes with the DFG-out inactive conformations of BMX (D) and BTK (E), the R-spine is partially disassembled, the αC-helix moved away from the catalytically competent position, and the characteristic salt bridge K-E is broken. The binding mode of Dasatinib remained stable in MD simulations of the crystal structures in all Dasatinib-kinase complexes. The R-spine residues are annotated in spheres, the αC-helix is shown in ribbons, and the binding mode of Dasatinib is shown in sticks. The residues forming a catalytic K-E salt bridge in the active site are shown in sticks.