Evolutionary divergence in the conformational landscapes of tyrosine vs serine/threonine kinases

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

Inactive conformations of protein kinase catalytic domains where the DFG motif has a “DFG-out” orientation and the activation loop is folded present a druggable binding pocket that is targeted by FDA-approved ‘type-II inhibitors’ in the treatment of cancers. Tyrosine kinases (TKs) typically show strong binding affinity with a wide spectrum of type-II inhibitors while serine/threonine kinases (STKs) usually bind more weakly which we suggest here is due to differences in the folded to extended conformational equilibrium of the activation loop between TKs vs. STKs. To investigate this, we use sequence covariation analysis with a Potts Hamiltonian statistical energy model to guide absolute binding free-energy molecular dynamics simulations of 74 protein-ligand complexes. Using the calculated binding free energies together with experimental values, we estimated free-energy costs for the large-scale (~17–20 Å) conformational change of the activation loop by an indirect approach, circumventing the very challenging problem of simulating the conformational change directly. We also used the Potts statistical potential to thread large sequence ensembles over active and inactive kinase states. The structure-based and sequence-based analyses are consistent; together they suggest TKs evolved to have free-energy penalties for the classical ‘folded activation loop’ DFG-out conformation relative to the active conformation, that is, on average, 4–6 kcal/mol smaller than the corresponding values for STKs. Potts statistical energy analysis suggests a molecular basis for this observation, wherein the activation loops of TKs are more weakly ‘anchored’ against the catalytic loop motif in the active conformation and form more stable substrate-mimicking interactions in the inactive conformation. These results provide insights into the molecular basis for the divergent functional properties of TKs and STKs, and have pharmacological implications for the target selectivity of type-II inhibitors.

Editor's evaluation

This important paper provides a convincing mechanism for the relative binding specificity of Type II inhibitors to kinases. The combination of a sequence-derived Potts model with experimental dissociation constants and calculated free energies of binding to the DFG-out state is highly compelling and goes beyond the current state-of-the-art. Given the importance of kinases in pathophysiological processes, the results will be of interest to a broad audience and, in addition, the combination of computational methods can be applicable to a wide variety of other biophysical processes that involve conformational rearrangements.

https://doi.org/10.7554/eLife.83368.sa0

Introduction

The human genome contains approximately 500 eukaryotic protein kinases which coordinate signaling networks in cells by catalyzing the transfer of a phosphate group from ATP to serine, threonine, or tyrosine residues (Manning, 1995; Modi and Dunbrack, 2019a). The GO (gene ontology) database identifies 351 (~70%) of these enzymes as serine/threonine kinases (STKs) and 90 (~18%) as tyrosine kinases (TKs). STKs are an ancient class of protein kinases that predate the divergence of the three domains of life (bacteria, archaea, eukaryote) (Stancik et al., 2018), whereas TKs are a more recent evolutionary innovation, having diverged from STKs about 600 million years ago during early metazoan evolution (Miller, 2012; Sebé-Pedrós et al., 2016). Kinases are important therapeutic targets in a large number of human pathologies and cancers. Both TKs and STKs share a striking degree of structural similarity in their catalytic domains, owing to evolutionary selective pressures that preserve their catalytic function; in particular, the location and structure of the ATP binding site are highly conserved which raises significant challenges for the design of small-molecule ATP-competitive inhibitors that are both potent for their intended target(s) and have low off-target activity for unintended kinase targets. The latter is referred to as the ‘selectivity’ of an inhibitor, a property which is difficult to predict and control but is nonetheless very important for developing drugs with minimal harmful side effects.

A particular class of ATP-competitive kinase inhibitors which were proposed to have a high potential for selectivity are called ‘type-II inhibitors’ which only bind when the kinase adopts an inactive ‘DFG-out’ conformation. ‘DFG’ (Asp-Phe-Gly) refers to a conserved catalytic motif located at the N-terminus of an ~20 residue-long ‘activation loop’ that is highly flexible and controls the activation state of the kinase and the structure of the substrate binding surface. The precise arrangement of catalytic residues and the structural organization of large regulatory elements, such as the activation loop and nearby ‘αC-helix’, are strongly coupled to the conformation of the DFG motif and the DFG-1 residue preceding it, which is well described by regions on the Ramachandran map occupied by the Asp, Phe, and DFG-1 residues (beta-turn, right-handed alpha-helix, left-handed alpha-helix) and the χ₁ rotamer state of the DFG-Phe sidechain (trans, gauche-minus, gauche-plus). Recently, Dunbrack and co-workers identified eight major conformational states in the Protein Data Bank (PDB) based on these metrics (Modi and Dunbrack, 2019b). The most common state, which is evolutionarily conserved in all kinases, corresponds to the active ‘DFG-in’ conformation. In this conformation all structural requirements for catalysis are typically met, e.g., a complete hydrophobic spine, a salt bridge between the conserved β3-Lys and αC-Glu residues, and an extended activation loop which forms the substrate binding surface. Inactive kinases in the PDB are most frequently seen in an ‘Src-like inactive’ conformation where the DFG is ‘in’, but the αC-helix is swung outward, breaking the β3-Lys → αC-Glu salt bridge and disassembling the hydrophobic spine. Disassembly of the hydrophobic spine caused by αC-helix rotation increases the cavity volume around the DFG-Phe residue, allowing it to pass through the Src-like inactive conformation and completely ‘flip’ from DFG-in to DFG-out (Levinson et al., 2006; Shan et al., 2013). The classical DFG-out conformation, targeted by type-II inhibitors, displays a highly reorganized activation loop that is folded away from the αC-helix, projecting toward solvent or forming stable secondary structure and substrate-mimicking interactions. We refer to these states of the activation loop collectively as ‘folded’, to describe its ~17 Å reorganization relative to the active ‘extended’ conformation, wherein the substrate binding surface has been ‘folded up’ toward the kinase N-terminal lobe and away from the αC-helix.

In both TKs and STKs, the activation loop undergoes this large-scale conformational change when the DFG motif flips from the active ‘DFG-in’ conformation to the classical DFG-out conformation. The DFG flip swaps the positions of DFG-Phe and DFG-Asp, opening a hydrophobic ‘back pocket’ that is connected to the conserved ATP binding site through the ‘gatekeeper’ residue. Type-II inhibitors typically have a long chemical fragment that allows them to bind across the gatekeeper and form interactions with residues in the back pocket. In contrast, type-I inhibitors (the majority of kinase drugs) occupy the ATP pocket but not the back pocket and can bind to either DFG-in or DFG-out. For these reasons, it has been proposed that type-II inhibition holds greater potential for the design of highly selective drugs (Vijayan et al., 2015; Davis et al., 2011; Anastassiadis et al., 2011); it has been shown that different kinase sequences have different propensities to adopt DFG-out in the absence of inhibitor (Haldane et al., 2016; Hari et al., 2013), and the DFG-out back pocket has been suggested to have a lesser degree of sequence and structural homology between kinases (Liu and Gray, 2006). However, the notion that type-II inhibitors developed to-date are more selective than type-I inhibitors has been brought into question (Zhao et al., 2014; Klaeger et al., 2017), suggesting that further consideration of the energetic contributions described above is required.

In order to fully exploit the target-selective potential of type-II inhibitors it is necessary to understand the underlying sequence-dependent principles that control the conformational preferences of their kinase targets, and the extent to which this has been diversified by evolution. This can, in principle, be directly approached using free-energy simulations to estimate the reorganization free-energy required for different kinases to adopt DFG-out, although this is computationally very expensive and of uncertain reliability for conformational changes involving large-scale loop reorganizations, such as the ~17 Å ‘folding’ of the activation loop that accompanies the transition from active DFG-in to the inactive, classical DFG-out state. To accommodate this limitation, we employ modern sequence-based computational methods to characterize the conformational selection process over the entire kinome and combine the sequence-based results with structure-based free energy simulations with the goal of identifying evolutionarily divergent features of the energy landscape that control the preference of individual kinases for the active (DFG-in) vs inactive (DFG-out ‘folded activation loop’) states. To this end, we report evidence that TK catalytic domains have a molecular evolutionary bias that shifts their conformational equilibrium toward the inactive ‘folded activation loop’ DFG-out state in the absence of activation signals. In contrast, STKs as a class have a more stable active conformation which is favored over the DFG-out state due to sequence constraints in the absence of other signals.

As described below, our analysis of a previously published kinome-wide assay suggests that TKs have properties which privilege the binding of type-II inhibitors in comparison to STKs, which leads us to hypothesize an evolutionary divergence in their conformational energy landscapes. To investigate this, we used a Potts Hamiltonian statistical energy model derived from residue-residue covariation in a multiple sequence alignment (MSA) of protein kinase sequences to probe the active DFG-in ↔ classical DFG-out conformational equilibrium as previously described (Haldane et al., 2016). Using an approach that involves ‘threading’ a large number of kinase sequences onto ensembles of active DFG-in and classical DFG-out structures from the PDB and scoring them using the Potts Hamiltonian, we are able to view the evolutionary divergence in TK and STK conformational landscapes. This calculation only probes the free-energy difference between the active DFG-in and classical DFG-out conformations, and by construction does not consider alternative conformations (e.g. ‘Src-like inactive’) that might be important for analyzing the type-II binding pathway. As discussed below, the Potts calculations from this two-state model correlate well with the free-energy cost to adopt the classical DFG-out conformation.

To validate our results, we used the Potts statistical energy threading calculations to guide target selection for a set of more computationally intensive free-energy simulations. These simulations use type-II inhibitors as tools to probe kinase targets that have already reorganized to DFG-out, allowing us to estimate the free-energy of reorganization ( $Δ G_{r e o r g}$ ) as the excess between the absolute binding free-energy (ABFE) calculated from simulations and the standard binding free-energy measured experimentally in vitro, which already includes the cost to reorganize. Although our methods avoid sampling the conformational change directly, we show how important structural determinants of the conformational change can be identified by analyzing residue-pair contributions to the Potts threading calculations, enabling us to reason about the molecular evolutionary basis for the differences in conformational behavior observed for TKs and STKs.

Results

Insights into the sequence-dependent conformational free-energy landscape

The binding of type-II inhibitors is achieved once a protein kinase has reorganized to the DFG-out with activation loop folded conformation (classical DFG-out). We sought initial insight into the conformational equilibrium from type-II binding data available publicly in the form of literature-reported dissociation constants (K_d). From the binding assay reported by Davis et al., 2011, we report a ‘hit’ where an inhibitor binds to a kinase with K_d ≤10 μM. Using this criterion, a type-II inhibitor hit rate was calculated for each kinase. Analysis of the type-II hit rate distributions for STKs and TKs from the Davis assay (Figure 1A) indicates that STKs, on average, have an unfavorable contribution to the binding of type-II inhibitors relative to TKs.

Figure 1 with 3 supplements see all

Download asset Open asset

Viewing the conformational landscape of the human kinome.

(A) Hit rate distributions from kinome-wide experimental binding assays with type-II inhibitors for human serine/threonine kinases (STKs; blue, top) and tyrosine kinases (TKs; orange, bottom) with small gatekeepers (solid bars; sidechain volume <110 Å³) and large gatekeepers (hatched bars; sidechain volume >110 Å³). (B) PyMol (pymol, 2015) visualization of two conformational ensembles populated by Abl kinase from recent solution NMR (Nuclear Magnetic Resonance) experiments (Xie et al., 2020). The active DFG-in conformation where the activation loop is ‘extended’ (red, the Protein Data Bank [PDB]: 6XR6) and the inactive classical DFG-out conformation where the activation loop is ‘folded’ (blue, PDB: 6XRG) both exist in the absence of ligands, but there is a free-energy cost to transition between them (Xie et al., 2020). Type-II inhibitors preferentially bind to this folded DFG-out state. (C) Correlation between Potts DFG-out penalty (ΔE Potts) and hit rates for kinases with small gatekeepers only, to control for gatekeeper effects (Pearson correlation of –0.59, p<0.001). (D) Potts DFG-out penalties calculated for the human kinome and plotted using CORAL (Metz et al., 2018); the TK branch appears to have lower penalties relative to the rest of the kinome, which represent STKs. See Figure 1—source data 1 for values of the calculated type-II hit rates and Potts threaded-energy scores over the human kinome.

Figure 1—source data 1 Type-II inhibitor hit rates and Potts threaded-energy penalties for DFG-out calculated for tyrosine kinases and serine/threonine kinases from the human kinome.: https://cdn.elifesciences.org/articles/83368/elife-83368-fig1-data1-v2.csv
Download elife-83368-fig1-data1-v2.csv

The size of the gatekeeper residue is important for type-II binding as it controls access to a hydrophobic pocket adjacent to the ATP binding site that is traversed by type-II inhibitors (Zuccotto et al., 2010; Bosc et al., 2015; Liu et al., 1998; Ghose et al., 2008; van Linden et al., 2014), and the size of the gatekeeper residue is thought to negatively affect type-II binding (Zuccotto et al., 2010; Azam et al., 2008; Lovera et al., 2015; Yun et al., 2008). Because TKs tend to have small gatekeepers in comparison to STKs (Zuccotto et al., 2010; Taylor and Kornev, 2011), we considered this as a possible explanation behind the bias for TKs to have larger type-II hit rates. By plotting the hit rate distributions for STKs and TKs where the gatekeeper is either small or large (Figure 1A) we confirm that gatekeeper size has an important influence on type-II binding for both STKs and TKs (Azam et al., 2008). However, the hit rate distribution for TKs appears more sensitive to gatekeeper size than STKs. Even with small gatekeepers, there is a significant fraction of STKs that have hit rates of zero compared with TKs, suggesting the difference in hit rates between TKs and STKs cannot be accounted for primarily by the size of the gatekeeper residue.

Recent solution NMR experiments with Abl kinase revealed two DFG-out conformational states (Xie et al., 2020) one where the DFG motif has flipped ‘DFG-in to DFG-out’ but the activation loop remains in a ‘minimally perturbed’ active-like conformation, and the other state is a classical ‘folded’ DFG-out conformation where the activation loop has moved ~17 Å away from the active conformation (Figure 1B), and the DFG motif is in a ‘classical DFG-out’ (Vijayan et al., 2015) or ‘BBAminus’ (Modi and Dunbrack, 2019b) state. Type-II inhibitors were shown to preferentially bind to this folded DFG-out state, confirming observations that Abl is almost always co-crystallized with type-II inhibitors in this conformation. This binding behavior is also exhibited by other kinases, suggested by the large number of activation loop folded DFG-out states seen in type-II bound co-crystal structures (Figure 1—figure supplement 1). Hence, the importance of large-scale activation loop conformational changes in type-II binding and the large number of residue-residue contact changes involved in this transition (Figure 1—figure supplement 2) suggests the sequence variation of the activation loop and the catalytic loop with which it interacts, might contour the conformational landscape differently for TKs compared with STKs. To investigate this, we used a Potts statistical energy model of sequence covariation to estimate the energetic cost of the active DFG-in (activation loop extended) → inactive DFG-out (activation loop folded) transition for human TKs and STKs (see Methods).

Patterns of coevolution of amino acids at different positions in an MSA are thought to largely reflect fitness constraints for fold stability and function between residues close in 3D space (Lapedes et al., 2012; Hopf et al., 2015; Hopf et al., 2017; Morcos et al., 2014), and these coevolutionary interactions can be successfully modeled by a Potts Hamiltonian (Weigt et al., 2009; Lunt et al., 2010) which we inferred using Mi3-GPU, an algorithm designed to solve ‘Inverse Ising’ problems for protein sequences with high accuracy (Haldane and Levy, 2021). The pairwise interactions from the Potts model can be used as a simple threaded energy function to estimate energetic differences between two conformations, based on changes in residue-residue contacts in the PDB (Haldane et al., 2016). We have calculated the threading penalty for all kinases in the human kinome. Our calculations show the Potts predicted DFG-out penalty ( $Δ E_{P o t t s}$ ), which is dominated by large-scale reorganization of the activation loop to the folded DFG-out state, is correlated with type-II hit rates (Figure 1C) when controlling for gatekeeper size. From this, we determine that sequence variation of the activation loop and the contacts broken/formed by its large-scale conformational change (Figure 1—figure supplement 2) makes an important contribution to the binding affinity of type-II inhibitors.

Notably, our calculations over the entire human kinome show that the large majority of kinases with large $Δ E_{P o t t s}$ (unfavorable conformational penalties) are STKs, and the large majority of low-penalty kinases are TKs (Figure 1D). To validate this finding, we next perform an independent and more computationally intensive prediction of the conformational reorganization energy of TKs and STKs for select kinase targets, chosen based on the kinome calculations of $Δ E_{P o t t s}$ and type-II hit rates shown in Figure 1, in which we use type-II inhibitors as probes in ABFE simulations as described in the following section. By comparing the conformational penalties predicted from these structure-based molecular dynamics (MD) free-energy simulations with the Potts conformational penalty scores, we also identify the scale of $Δ E_{P o t t s}$ in physical free-energy units. This allows us to predict physical conformational free energies based on Potts calculations which can be carried out at scale on large numbers of sequences, to evaluate the evolutionary divergence of the conformational penalty between STKs and TKs generally.

Structure-based free-energy simulations guided by the sequence-based Potts model

Relative binding free-energy simulations are now widely employed to screen potent inhibitors in large-scale drug discovery studies (Wang et al., 2015). These methods are used to determine the relative free energy of binding between ligands that differ by small substitutions, which permit one to simulate along an alchemical pathway that mutates one ligand to another. By leaving the common core scaffold unperturbed, the cost and difficulty of sampling the transition between unbound (apo) and bound (holo) states of the system are avoided (Wang et al., 2015; Hayes et al., 2022; Guest et al., 2022). Alternatively, alchemical methods to determine ABFEs, such as the ‘double decoupling’ method employed in this work, sample the apo → holo transition along a pathway that decouples the entire ligand from its environment. While more computationally expensive, the advantage of ABFE is that the computed $Δ G_{b i n d}$ can be directly compared with experimental binding affinities, and successful convergence does not rely on the structural similarity of compounds being simulated (Cournia et al., 2020; Li et al., 2020; Heinzelmann and Gilson, 2021; Lee et al., 2020; Gilson et al., 1997; Qian et al., 2019; Sun et al., 2022).

Our alchemical ABFE simulations of type-II inhibitors binding to TKs and STKs simulate the apo and holo states of the kinase domains in the classical DFG-out conformation with the activation loop folded, starting from the experimentally determined co-crystal structure of the holo state. The apo state remains DFG-out with the activation loop folded throughout the simulations, and therefore the calculated ABFE ( $Δ G_{b i n d}^{A B F E}$ ) excludes the cost to reorganize from DFG-in ( $Δ G_{r e o r g}$ ). On the other hand, standard binding free-energies ( $Δ G_{e x p}^{o}$ ) determined experimentally from inhibition or dissociation constants (Equation 1) implicitly include the free-energy cost to reorganize. Therefore given the experimentally determined total binding free energy, $Δ G_{e x p}^{o}$ , ABFE simulations can be used to separate the free energy of ligand-receptor association in the inactive state ( $Δ G_{b i n d}^{A B F E}$ ) from the cost to reorganize from the active to inactive state, ${Δ G}_{r e o r g}$ (Equation 3; Deng et al., 2011; Lin et al., 2014; Lin et al., 2013).

We calculated $Δ G_{e x p}^{o}$ (Equation 1) from literature reported IC₅₀ or $K_{d}$ values, where the standard concentration $C_{0}$ is set to 1 $M$

{Δ G}_{e x p}^{o} = k_{b} T l n (K_{d} / C_{0})

$Δ G_{e x p}^{o}$ can be expressed as the sum of the free-energy change to reorganize from the active to inactive state, $Δ G_{r e o r g}$ plus the free energy to bind to the inactive state $Δ G_{b i n d}^{A B F E}$ (Equation 2). $Δ G_{r e o r g}$ is therefore the excess free-energy difference between $Δ G_{e x p}^{o}$ and $Δ G_{b i n d}^{A B F E}$ (Equation 3).

{Δ G}_{e x p}^{o} = {Δ G}_{r e o r g} + {Δ G}_{b i n d}^{A B F E}_{}

{Δ G}_{r e o r g} = {Δ G}_{e x p}^{o} - {Δ G}_{b i n d}^{A B F E}

Type-II inhibitors generally bind when the activation loop is in a folded DFG-out conformation (Figure 1B), which presents major challenges for direct simulations to determine the free energy cost of the conformational change in contrast to the method employed here (Equation 3).

Because the type-II inhibitor imatinib is co-crystallized in a type-II binding mode with MAPK14 (p38α), an STK, and several other TKs (e.g. ABL1, DDR1, LCK, CSF1R, KIT, and PDGFRA), we chose this inhibitor as an initial probe of our hypothesis that TKs evolved to have lower $Δ G_{r e o r g}$ than STKs (Figure 2). In this example we note that TKs bind strongly to imatinib (‘STI’ in Figure 2) with an average $Δ G_{e x p}^{o}$ of –9.3 kcal/mol, in contrast to the STK MAPK14 which binds this drug very weakly ( $Δ G_{e x p}^{o} = - 6.1$ kcal/mol). At face value this appears consistent with our analysis from Figure 1D, where we calculated a large Potts DFG-out penalty for MAPK14 ( $Δ E_{P o t t s} = 5.2$ ) and low penalties for TKs, suggesting that the weak binding of imatinib to MAPK14 is due at least partially to large $Δ G_{r e o r g}$ . To confirm this, we used ABFE simulations with the imatinib: MAPK14 complex to evaluate Equation 3, confirming that MAPK14 incurs a large penalty to adopt the DFG-out conformation with the activation loop folded ( $Δ G_{r e o r g} = 5$ kcal/mol) (Figure 2).

Figure 2 with 1 supplement see all

Download asset Open asset

Overview of the conformational landscapes between serine/threonine kinases (STKs) and tyrosine kinases (TKs) from absolute binding free-energy simulations, where we compare ΔG_bind (hatched bars) from binding free-energy simulations with ΔG_exp (solid bars) for the type-II inhibitors imatinib (the Protein Data Bank [PDB] code: STI) and BIRB-976 (PDB code: B96) vs several TKs (orange) and STKs (blue).

Figure 2—source data 1 Absolute binding free-energy results for kinases bound to imatinib and BIRB-796.: https://cdn.elifesciences.org/articles/83368/elife-83368-fig2-data1-v2.csv
Download elife-83368-fig2-data1-v2.csv

Despite the large $Δ G_{r e o r g}$ predicted for MAPK14 by both the Potts model and the simulations with imatinib described above, highly potent type-II inhibitors have been successfully developed for this kinase. For example, BIRB-796 (Pargellis et al., 2002) binds to MAPK14 about 7 kcal/mol more strongly than imatinib. This stronger binding of BIRB-796 is captured by $Δ G_{b i n d}^{A B F E}$ from our simulations (Figure 2), and the calculated value of $Δ G_{r e o r g}$ for this complex ( $Δ G_{r e o r g} \approx 4$ kcal/mol) is very close to the corresponding estimate of $Δ G_{r e o r g}$ based on simulations with imatinib (Figure 2). Importantly, this result suggests that STKs can be potently inhibited by type-II inhibitors despite their large $Δ G_{r e o r g}$ . To support this, we performed additional ABFE simulations with BIRB-796 and calculated $Δ G_{r e o r g}$ for two additional STKs predicted to have large reorganization penalties (MAPK9 and BRAF, $Δ E_{P o t t s} \geq 4$ ). We calculated $Δ G_{r e o r g} > 8$ kcal/mol for MAPK9 and BRAF, which is consistent with predictions from the Potts model, and comparison of $Δ G_{e x p}^{o}$ and $Δ G_{b i n d}^{A B F E}$ in Figure 2 confirms that BIRB-796 is able to overcome the large $Δ G_{r e o r g}$ of certain kinases to attain high experimental potencies (e.g. MAPK14 and MAPK9). To further validate this result, we calculated $Δ G_{r e o r g}$ via ABFE simulations of BIRB-796 binding to a TK predicted by the Potts model to have a low penalty (PTK2B, $Δ E_{P o t t s} < 1$ ), which again shows consistency with our Potts prediction of the conformational landscape (Figure 2). The relatively weak value of $Δ G_{b i n d}^{A B F E}$ for this kinase compared with MAPK14 is also consistent with observations of the BIRB-796: PTK2B co-crystal structure (PDB: 3FZS), where the binding mode in PTK2B is more weakly associated with the ATP pocket in comparison with MAPK14 (Han et al., 2009).

The analysis above provides initial support for our hypothesis about the evolutionarily divergent STK and TK conformational landscapes. To further develop this approach, we identified five STKs and five TKs which are predicted by the Potts threading calculations to have large and small $Δ G_{r e o r g}$ , respectively, and for which there are sufficient experimental structural and inhibitory data (co-crystal structures and binding constants) to calculate an average $Δ G_{r e o r g}$ via Equation 3 for each target using multiple type-II inhibitor probes. For each of the TK and STK targets, these sets of calculations can be visualized as a linear regression of $Δ G_{e x p}^{o}$ vs $Δ G_{b i n d}^{A B F E}$ where the slope is constrained to one, consistent with Equation 2 (see Methods for details). We employed this workflow for the set of five TK and five STK targets by simulating 22 and 23 type-II inhibitor complexes, respectively.

The result of this workflow for the set of five TKs and their type-II complexes revealed a low average ${Δ G}_{r e o r g}$ of <1 kcal/mol (Figure 3B), consistent with our initial predictions from Potts $Δ E s$ and type-II hit rates (Table 1). On the other hand, the binding free-energy simulations for the set of five STKs and their type-II complexes show an average of ~6 kcal/mol of ${Δ G}_{r e o r g}$ is required for these kinases to adopt DFG-out conformation, which is also consistent with our initial predictions from the Potts model (Table 1). To verify that the large $Δ G_{r e o r g}$ identified for STKs is a property of conformational selection for DFG-out rather than systematic overestimation of $Δ G_{b i n d}^{A B F E}$ for these kinases, we performed ABFE simulations of type-I inhibitors binding to the same set of STKs (an additional 23 complexes). For the binding of type-I inhibitors, we expect there to be no reorganization penalty due to the lack of DFG-out conformational selection in their binding mechanism. As anticipated, the calculated values of ${Δ G}_{b i n d}^{A B F E}$ for type-I inhibitors are very close to their experimental binding affinities ( $Δ G_{e x p}^{o}$ ) (Figure 3A).

Figure 3 with 1 supplement see all

Download asset Open asset

Using type-II inhibitors as tools to probe the conformational landscape of TKs and STKs.

(A) The average ΔG_reorg calculated via absolute binding free-energy (ABFE) simulations with 23 type-I (stars) and 23 type-II inhibitors (circles) complexes in the active DFG-in and inactive DFG-out state, respectively, computed from five serine/threonine kinase (STK) targets (Table 1) and (B) computed with 22 type-II inhibitors vs five tyrosine kinase (TK) targets in the DFG-out state (Table 1). (C) Kinome plot created with CORAL (Metz et al., 2018) illustrating the selection of five TKs and five STKs which are detailed in Table 1.

Figure 3—source data 1 Absolute binding free-energy results for five tyrosine kinases and five serine/threonine kinases bound to type-II and type-I inhibitors.: https://cdn.elifesciences.org/articles/83368/elife-83368-fig3-data1-v2.csv
Download elife-83368-fig3-data1-v2.csv

Table 1

Calculation of reorganization free-energy for five TKs and five STKs.

Type-II hit rates from Davis et al. and Potts threaded energy penalties from Figure 1 were used to guide the selection of five serine/threonine kinase (STK) and five tyrosine kinase (TK) targets for absolute binding free-energy simulations. For some kinases, the hit rate binary classifier captures a set of relatively weak hits with average binding which, in context with large Potts penalty (see Figure 1D), might be explained by a large ΔG_reorg incurred for the folded DFG-out state (Figure 1B). See Figure 3—source data 1 for detailed data and references for experimental binding affinities.

Kinase	Class	Hit rate^* (K_d <10 µM)	Potts penalty^† (ΔE_Potts)	Calculated ΔG_reorg^‡
MELK	STK	3	5.9	5.6±0.2
MAPK9	STK	5	4.7	6.9±0.3
CDK2	STK	2	5.3	7.7±0.2
IRAK4	STK	0	2.5	5.4±0.2
BRAF	STK	7	4.0	6.5±0.1
ABL1	TK	10	–1.0	1.3±0.3
LCK	TK	11	0.5	1.0±0.3
TIE2	TK	6	1.1	–0.3±0.2
NTRK2	TK	6	–1.0	1.7±0.1 §
DDR1	TK	11	0.4	0.3±0.3

*

Type-II inhibitors only, data from Davis et al., 2011.
†

Calculations from Figure 1d.
‡

ΔG_reorg was calculated from Equation 3. Reported standard deviations were calculated by propagating error from the simulations used in the calculation of average ΔG_reorg in units of kcal/mol (see Figure 3—source data 1 for statistics from individual simulations).
§

Average and standard deviation calculated from two simulated complexes.

We find that the set of type-II inhibitors complexed with STKs in this dataset tends to have more favorable binding free energies to the reorganized receptor ( ${Δ G}_{b i n d}^{A B F E}$ ) than type-II inhibitors complexed with TKs, as shown by their distribution along the horizontal axes in Figure 3. The reason for this can be understood as a consequence of selection bias. Our selection of STK complexes for this study usually involved lead compounds from the literature, which were designed for high on-target experimental potency and published for their pharmaceutical potential, similar to BIRB-796: MAPK14 which is a tightly bound complex with high experimental affinity despite the large $Δ G_{r e o r g}$ incurred by this kinase (Figure 2). This tight binding is reflected by the favorability of the $Δ G_{b i n d}^{A B F E}$ term which must be implicitly tuned by medicinal chemists to overcome the large $Δ G_{r e o r g}$ found in STKs. Meanwhile, the chemical space of type-II inhibitors studied against TKs appears to be privileged by their low $Δ G_{r e o r g}$ , judging by the comparably weak $Δ G_{b i n d}^{A B F E}$ for these complexes. This ultimately gives rise to similar experimental potencies for the binding of type II inhibitors to TKs and STKs plotted in Figure 3.

The results of the MD binding free-energy simulations when combined with experimental binding affinities, reveal significant differences in the conformational free-energy landscapes between STKs and TKs. The DFG-in (activation loop extended) to DFG-out (activation loop folded) reorganization penalties are strongly correlated with corresponding $Δ E s$ calculated from the Potts model ( $R^{2} = 0.75, {P \approx 10}^{- 3}$ ) emphasizing the connection between coevolutionary statistical energies in sequence space and physical free-energies in protein conformational space (Figure 4). From this relationship, we can approximate a scale for the Potts $Δ E$ scores in physical free-energy units which describe the conformational landscapes of folded proteins in a similar manner to that of an earlier study of protein folding landscapes (Morcos et al., 2014); we find that a Potts statistical energy difference ΔE of one unit corresponds approximately to 1.3 kcal/mol of $Δ G_{r e o r g}$ .

Figure 4

Download asset Open asset

Correlation between ΔE_Potts and averaged calculations of ΔG_reorg for five tyrosine kinases (TKs) and five serine/threonine kinases (STKs; Table 1).

Structural and evolutionary basis for the divergent TK and STK conformational landscapes

The consistency of the predictions between $Δ G_{r e o r g}$ and $Δ E_{P o t t s}$ identified from Potts-guided free-energy simulations (Figure 4) led us to investigate whether the observed difference in the free energy to reorganize from the active to inactive state is a more general feature that distinguishes TKs from STKs. To this end, we extracted ~200,000 STKs and ~10,000 TKs from the large MSA of Pfam sequences used in the construction of our Potts model, based on patterns of sequence conservation that clearly distinguish the two classes (see Methods). For each sequence, we calculated $Δ E_{P o t t s}$ threaded over the structural database (a total of 4268 active DFG-in and 510 classical DFG-out PDB structures) and plotted the distributions for TKs and STKs, revealing a bias for STKs toward larger Potts conformational penalties (Figure 5). The average difference between these distributions, $Δ Δ E = 3.2,$ is extremely unlikely to be observed by chance (p≤10⁻¹⁵, see Methods) and supports the hypothesis that TKs are evolutionarily biased toward a lower free-energy cost to adopt the classical ‘folded activation loop’ DFG-out conformation ( $Δ G_{r e o r g}$ ) compared to STKs. We estimate that $Δ Δ E = 3.2$ corresponds to ~4.3 kcal/mol based on the analysis summarized in Figure 4.

Figure 5

Download asset Open asset

The distributions of Potts conformational penalties for (orange) 10,345 tyrosine kinases (TKs) from 471 different species and (blue) 210,862 serine/threonine kinases (STKs) from 2713 different species, showing that TKs tend to have smaller energetic penalties on average.

The difference in averages between these distributions is shown ( $Δ Δ E = 3.2$ ), which we estimate to be ≈4.3 kcal/mol based on the analysis in Figure 4. The probability density was plotted after down-weighting each sequence by the number of times another sequence in the same class (e.g. TKs) is observed within 40% identity (see Methods for details). The effective number of TKs (N_eff =1,096) and STKs (N_eff =22,893) s the sum totals of their down-weights, which is an unbiased measure of the sequence diversity in each probability distribution.

To gain insight into the molecular basis for this effect which distinguishes the conformational landscape of TKs from STKs, we examined the residue-residue interactions that make the most significant contributions to the observed $Δ Δ E$ . The difference between average Potts threaded-energy penalties, $Δ Δ E = {⟨Δ E⟩}_{S T K s} - {⟨Δ E⟩}_{T K s}$ , can also be written as a sum over pairs of alignment positions $i$ and $j$ along length $L$ of the aligned kinase domains, $\sum_{| i - j | > 4}^{L} Δ Δ E_{i j}$ (see Methods for details). We find that ~75% of the total contribution to $Δ Δ E$ (approximately 3 kcal/mol) can be traced back to a small number (10) of residue-residue interactions involving the activation loop, suggesting that mutations within the activation loop are largely responsible for the evolutionary divergence between the conformational free-energy landscapes of TKs and STKs. These interactions occur between important structural motifs responsible for controlling the stability of the active ‘extended’ conformation of the activation loop (Figure 6A), especially the activation loop N-terminal and C-terminal ‘anchors’ (Nolen et al., 2004), and the regulatory ‘RD-pocket’ formed by the HRD motif of the catalytic loop which functions to stabilize or destabilize this conformation depending on the activation loop’s phosphorylation state (Nolen et al., 2004; Johnson et al., 1996). The remaining top $Δ Δ E_{i j}$ correspond to contacts that stabilize the DFG-out ‘folded’ conformation of the activation loop for TKs, wherein the kinase’s substrate binding site recognizes its own activation loop tyrosine (Figure 6C, right; Hubbard et al., 1994). We describe these interactions below, focusing on the strongest effects involving these structural motifs that lead to differences in the conformational free-energy landscapes of TKs and STKs. The residue nomenclature we use in our descriptions follows the format $M S A_{A B L 1}^{C D K 6}$ , which is the unique residue numbering in our MSA followed by Abl1 (TK) numbering in the subscript (active PDB: 2GQG, inactive PDB: 1IEP) and CDK6 (STK) numbering in the superscript (active PDB: 1XO2, inactive PDB: 1G3N), corresponding to the original PDB files used to generate Figure 6B and Figure 6C.

Figure 6 with 2 supplements see all

Download asset Open asset

Molecular basis for the evolutionary divergence between tyrosine kinases (TKs) and serine/threonine kinases (STKs).

Visualizing the top interaction pairs that contribute to the result in Figure 5, which is discussed in the main text. (A) Diagrams depicting general features of the active ‘extended’ conformation of the activation loop (left) and the primary structure of these motifs in our multiple sequence alignment (MSA; right) with the HRD (pink), DFG (red), and APE (teal) motifs color-coded for reference. (**B–C**) Structural examples of a representative STK (CDK6) and TK (Abl) in the active DFG-in conformation (left) and the folded DFG-out conformation (right). Residues are labeled according to their position in our MSA, and colored according to A. The inset (center) displays the ΔE_Potts of the reference kinase derived from Figure 5 as well as a cartoon depicting their location in the distributions. The diagrams of CDK6 in the active DFG-in conformation (PDB: 1XO2, chain B), CDK6 in the folded DFG-out conformation (PDB: 1G3N, chain A), Abl in the active DFG-in conformation (PDB: 2G2I, chain A), and Abl in the folded DFG-out conformation (PDB: 1IEP, chain A) were generated with PyMol (pymol, 2015). All ligands and some backbone atoms were hidden for clarity.

The ‘RD-pocket’ (Nolen et al., 2004) is a conserved basic pocket formed by the Arg and Asp residues of the HRD motif (Johnson et al., 1996) ( $R 123_{361}^{144}$ and $D 124_{362}^{145}$ ) and a positively charged Lys or Arg that is often present in the N-terminal anchor of the activation loop ( $147_{386}^{168}$ ). $R 123_{361}^{144}$ in the HRD motif and Lys/Arg $147_{386}^{168}$ in the N-terminal anchor form an unfavorable like-charge interaction when the activation loop is in the active, extended conformation (Nolen et al., 2004). Kinase activation is typically a complex process involving many layers of regulation from other protein domains, cofactors, and phosphorylation events (Endicott et al., 2012) however, a general activation mechanism that applies to the majority of kinases involves quenching the net-charge of the RD-pocket by addition of a negatively charged phosphate group to a nearby residue in the activation loop, stabilizing the active conformation. The conservation of this regulatory mechanism in most protein kinases, particularly those bearing the HRD-Arg residue (termed ‘RD-kinases’ Johnson et al., 1996), explains why Lys or Arg is frequently observed at position $147_{386}^{168}$ of the N-terminal anchor of the activation loop. However, RD-TKs prefer Arg at this position (78%) which the Potts model suggests has a greater destabilizing effect on the active conformation than Lys (9%) due to interactions with HRD-Arg. The activation loop Arg also forms part of an electrostatic interaction network that stabilizes the ‘Src-like inactive’ conformation in TKs (Ozkirimli and Post, 2006; Wu et al., 2020), a conformation with a ‘partially’ folded activation loop (Figure 1—figure supplement 1) that is suggested to be an intermediate state along the transition to DFG-out (Levinson et al., 2006; Shan et al., 2013). On the other hand, RD-STKs display $K 147$ more frequently (26%), which the Potts model suggests interacts more favorably with the HRD-Arg independently of activation loop phosphorylation, thus contributing to greater stabilization of the active DFG-in conformation for RD-STKs in comparison with TKs. Additionally, the Potts statistical energy analysis suggests that packing interactions between HRD-Arg and $V 160_{400}^{180}$ located near the activation loop C-terminal also contributes to phosphorylation-independent stabilization of the active conformation in RD-STKs (Figure 6B left), appearing in 28% of RD-STKs and only 2% of RD-TKs. In RD-TKs, the residue at this position is usually Arg or Lys which appears to be repelled from the RD-pocket (Figure 6C left) and is suggested by the Potts model to result in a less stable active conformation.

The largest $Δ Δ E_{i j}$ term, which contributes 16% of the total difference in average Potts conformational penalties between STKs and TKs, comes from an interaction pair in the C-terminal of the activation loop ( $162_{402}^{182}$ ) and the C-terminal of the catalytic loop ( $126_{364}^{147}$ ). These residues form part of the ‘C-terminal anchor’ (Nolen et al., 2004) which is important for creating a suitable binding site for the substrate peptide. The C-terminal anchor residue $162_{402}^{182}$ is Pro in TKs and typically Ser or Thr in STKs (Taylor et al., 1995). In STKs, the sidechain hydroxyl of this residue forms a hydrogen bond with ${K 126}_{364}^{147}$ in the catalytic loop, creating a stable binding site for substrate phosphoacceptor residues and stabilizing the C-terminal anchor (Taylor et al., 1995). $K 126$ is also directly involved in catalysis by interacting with and stabilizing the gamma phosphate of ATP (Zheng et al., 1993); hence, it is often referred to as the ‘catalytic lysine’. The hydrogen bond between $K 126$ and $S$ or $T 162$ is almost always formed in the active DFG-in conformation, and we observe breakage of this hydrogen bond in many STKs crystallized in the DFG-out/activation loop folded conformation (e.g. CDK6, Figure 6B), suggesting that deformation of the C-terminal anchor contributes an energetic penalty for the active → inactive conformational change. In TKs, however, the catalytic lysine is almost always replaced with Ala, with the exception of a few TKs (e.g. c-Src) which have instead adopted Arg at this position. The C-terminal anchor of TKs containing $A 126$ and $P 162$ is less stable in comparison to STKs containing ( $K 126, T 162$ ) or ( $K 126, S 162$ ) for which the Potts coupling is very favorable, consistent with the structural observation that ( $A 126, P 162$ ) forms weak interactions (Figure 6C). Our analysis suggests the interaction pair ( $A 126, P 162$ ) weakens the C-terminal anchor, leading to a less stable active conformation in TKs as compared with STKs. Another significant contribution to the stabilization of the C-terminal anchor in the active DFG-in conformation for STKs comes from interactions between the residue pair ( $161_{401}^{181}$ , $166_{406}^{186}$ ) which are both located within the activation loop. We observe ( $G 161, M 166$ ) at this position pair in 33% of STKs, but never in TKs (Figure 6—figure supplement 1D). The Potts coupling between these residues is highly favorable. In contrast, we observe ( $L 161, M 166$ ) in 40% of TKs but never in STKs (Figure 6—figure supplement 1H), which have weaker coupling. The bulky sidechains of ( $L 161, M 166$ ) observed in TKs cause the activation loop to ‘bulge’ in this C-terminal region which has been previously identified as a feature of TKs that helps shape the substrate binding site to accommodate Tyr residues (Nolen et al., 2004). In addition to this paradigm, our analysis suggests that the C-terminal bulge results in weaker structural constraints on the active conformation relative to STKs.

In summary, TKs are suggested by the Potts statistical energy model which is based on sequence covariation, to have on average, weaker N-terminal anchor, RD-pocket, and C-terminal anchor interactions than STKs. This mechanism of shifting the TK conformational equilibrium away from the active DFG-in/extended activation loop conformation can explain 7 of the top 10 $Δ Δ E_{i j}$ , accounting for ~80% of these contributions to the divergence between the STK and TK conformational landscapes. The remaining ~20% of the top contributions can be attributed to residue-residue interactions that occur within the folded DFG-out conformation wherein the activation loop of TKs binds to the kinase’s own active site as though it was engaging a peptide substrate in trans (Figure 6C right) (Nolen et al., 2004). STKs, however, rarely adopt a folded DFG-out conformation with this property, and instead the activation loop is typically found to be unresolved and/or projecting outward toward solvent (Figure 6B right). The Potts model suggests that this substrate mimicry of the folded DFG-out activation loop observed in TKs is highly dependent on the presence of a Tyr phosphorylation site at position in the activation loop (Figure 6C). In the active conformation, the anionic $p Y 154$ stabilizes the active conformation by binding to the basic RD-pocket (Figure 6C left, phosphate not shown). However, in the (unphosphorylated) folded DFG-out conformation, this $Y 154$ mimics a substrate by stacking against the TK-conserved ${P 162}_{402}^{182}$ residue (Figure 6C right) (Nolen et al., 2004). The substrate mimicking nature of this binding mode is demonstrated by the autophosphorylation dimer structure of TK FGFR3 (PDB: 6PNX) solved recently (Figure 6—figure supplement 1K; Chen et al., 2020).

The striking connection between the ability of TKs to phosphorylate tyrosine substrates and their enhanced access to the DFG-out conformation via substrate-competitive contacts from their own activation loop described above suggests an evolutionary model for the TK conformational behavior characterized in this work. In this model, the coevolution of residues that form substrate-competitive contacts in folded DFG-out appears to be a byproduct of the evolutionary pressure for TKs to phosphorylate other TKs on their activation loop tyrosine residues. STKs, on the other hand, have optimized the binding of Ser and Thr substrates via a different binding mode (Endicott et al., 2012) which does not have the same energetic feedback with the stability of the folded DFG-out conformation. Additionally, the catalytic domains of TKs appear to have a less energetically stable ‘extended’ activation loop conformation than STKs, which may have encouraged the evolution of more complex mechanisms of allosteric regulation and autophosphorylation which are highly important regulatory mechanisms in TKs (Chen et al., 2020; Beenstock et al., 2016; Lemmon and Schlessinger, 2010). The combined effect of these two TK phenotypes, the former favoring the stabilization of folded DFG-out and the latter favoring destabilization of active DFG-in, may explain their low free-energy cost for the DFG-in → DFG-out conformational change compared with STKs.

Discussion

In this work, we have combined sequence and structure-based approaches to analyze the conformational free energy difference between active DFG-in and inactive DFG-out kinase states. Using a Potts statistical energy model derived from residue-residue covariation in a kinase family multiple sequence alignment, we first threaded all human STKs and TKs onto large ensembles of active DFG-in and classical DFG-out structures from the PDB. We found distinctly different distributions of threading scores for STKs compared with TKs, with STKs having a significant conformational reorganization penalty compared with TKs. The molecular basis for the evolutionary divergence in the conformational landscapes was analyzed; a substantial contribution to the difference is associated with sequence position pairs that couple the N and C terminal anchor residues of the activation loop to N and C terminal residues in the catalytic loop, according to the Potts statistical energy analysis. We then used the Potts statistical energy model to guide the selection for structure-based MD binding free energy simulations of 74 protein-ligand complexes; using the calculated binding free-energy estimates together with experimental values, we were able to estimate free-energy costs for the large-scale (~17–20 Å) conformational change of the activation loop by an indirect approach. The structure-based estimates of the reorganization free-energy penalties are consistent with the sequence-based estimates. Additionally, the strong correlation between $Δ G_{r e o r g}$ and $Δ E_{P o t t s}$ identified in this study reveals that the conformational landscape has a strong sequence dependence with STKs having an ~4 kcal/mol conformational free energy bias favoring the active state over the inactive state relative to TKs (Figure 4). We note that the most potent type-II inhibitors from the literature which target STKs bind with nanomolar K_ds, similar to that for TKs, despite the substantial additional reorganization penalty that STKs must overcome. This suggests that medicinal chemists have implicitly been able to exploit particularly favorable characteristics of the type-II binding pocket to design inhibitors with extremely strong affinities to the DFG-out (activation loop folded) receptor conformations of STKs, and that further analysis of the molecular basis for this tight binding could provide a basis for designing more selective inhibitors.

Materials and methods

Multiple sequence alignment (MSA) and classification of serine/threonine vs tyrosine kinases

Request a detailed protocol

An MSA of 236,572 protein kinase catalytic domains with 259 columns was constructed as previously described (McGee et al., 2021). STKs and TKs were classified based on patterns of sequence conservation previously identified by Taylor and co-workers (Taylor et al., 1995); characteristic sequence features of TKs and STKs which form their respective phosphoacceptor binding pockets are found at the HRD +2 (Ala or Arg in TKs, Lys in STKs) and HRD +4 (Arg in TKs, variable in STKs) in the catalytic loop, as well as the APE-2 residue (Trp in TKs, variable in STKs) in the activation loop. These residues correspond to positions 126, 128, and 165 in our MSA, respectively. Kinases which satisfy the conditions for TKs at all three positions were classified as TKs (10,345 raw sequences, 1069 effective sequences), and those that satisfy the condition for STKs at position 126 and are non-overlapping with the TK condition at position 128 were classified as STKs (210,862 raw sequences, 22,893 effective sequences). The effective number of sequences in each class was calculated by summing over sequence weights, where each sequence was assigned a weight defined as the fraction of the number of sequences in the same class that are within 40% identity. In this way, we correct for the effects of phylogenetics in the calculation of sample size as well as other quantities (see below).

Human kinome dataset

Request a detailed protocol

497 human kinase catalytic domain sequences were acquired from Modi and Dunbrack, 2019b (excluding atypical kinases). These sequences were aligned using a Hidden Markov Model (HMM) that contains 259 columns (L=259), which was constructed from the same MSA used to derive our Potts model. 447 human kinases remained after filtering-out sequences with 32 or more gaps.

Classification of gatekeeper size

Request a detailed protocol

The designation of gatekeeper residues as ‘large’ or ‘small’ was based on sidechain van der Waals volumes (Miller et al., 1987; Figure 1—figure supplement 3), where small gatekeepers have a volume of <110 Å³ (Gly, Ala, Ser, Pro, Thr, Cys, Val), and large gatekeepers have a volume of >110 Å³ (Asn, His, Ile, Leu, Met, Lys, Phe, Glu, Tyr, Gln, Trp, Arg).

PDB dataset and conformational states

X-ray crystal structures of tyrosine, serine/threonine, and dual-specificity eukaryotic protein kinases in the PDB were collected from http://rcsb.org on July 30, 2020. The protein sequences of 7919 chains were extracted from 6805 PDB files by parsing the SEQRES record and aligned to the MSA used to construct our Potts model, using an HMM.

Contact frequency differences between ensembles of active/DFG-in and inactive/DFG-out (classical DFG-out) PDB structures (Figure 1—figure supplement 2) are incorporated into the calculation of Potts threaded energies, which are central to this work. Our classification of the active DFG-in and classical DFG-out conformational states is based on ref 6, which we describe in further detail here:

Active DFG-in (BLAminus)

Request a detailed protocol

The BLAminus state of the DFG motif (active DFG-in) is the only active conformation of the activation loop, in contrast to several other inactive DFG-in states (ABAminus, BLAplus, BLBplus, BLBminus, BLBtrans). In the active conformation, all structural requirements for catalytic activity are typically met, e.g., a complete hydrophobic spine, a salt bridge between β3-Lys → αC-Glu, and an extended activation loop that ensures unobstructed substrate-binding, all of which have high correspondence with the BLAminus state (Modi and Dunbrack, 2019b). Our analysis using the software provided in Modi and Dunbrack, 2019b identifies 3643 structures in this conformation belonging to STKs and 625 structures belonging to TKs (4268 structures in total).

Classical DFG-out (BBAminus)

Request a detailed protocol

The DFG-out state is characterized by a ‘flip’ of the conserved Phe to occupy the ATP binding pocket which is otherwise occupied by the conserved Asp in the active conformation. The classical DFG-out conformation is associated with the binding of type-II inhibitors which occupy the back pocket region opened up by the DFG-flip (Vijayan et al., 2015). This conformation corresponds to the BBAminus rotamer state of these residues, and it is the dominant DFG-out conformation associated with the binding mode of type-II inhibitors. This classical DFG-out or BBAminus conformation is correlated with a larger-scale conformational change of the activation loop that involves an ~17 Å ‘folding’ transition with respect to the active conformation. This conformational change is usually accompanied by the formation of secondary structure that obstructs the typical substrate binding surface (e.g. PDB ID: 2HIW, Figure 1—figure supplement 1B). Our conformational analysis via Modi and Dunbrack, 2019b identifies 224 structures in this conformational state belonging to STKs and 286 structures belonging to TKs (510 structures in total).

Contact frequency differences

Request a detailed protocol

Each of the clustered PDB structures were converted to an adjacency matrix of binary contacts (1 for ‘in-contact’, 0 otherwise). A contact between residues $i$ and $j$ in structure $n$ was assigned when their nearest sidechain atoms (excluding hydrogen) were detected within a distance $r_{i j} (n) < 6 \overset{˚}{A}$ . The contact frequency $c_{i j}^{A}$ in cluster $A$ (e.g. active DFG-in) was calculated for each residue pair ( $i, j$ ) by taking a weighted average over all instances of a contact in that cluster –

c_{i j}^{A} = \frac{1}{\sum_{n}^{N_{A}} w_{n}^{A B}} \sum_{n}^{N_{A}} w_{n}^{A B} δ_{i j} (n)

δ_{i j} (n) = {\begin{cases} 1, r_{i j} (n) < 6 \overset{˚}{A} \\ 0, r_{i j} (n) \geq 6 \overset{˚}{A} \end{cases}

Where $N_{A}$ is the number of PDB chains in cluster $A$ , and weights were calculated with $w_{n}^{A B} = \frac{1}{u_{n}^{A B}}$ , where $u_{n}^{A B}$ is the number of times the UniProt ID of structure $n$ is found within either cluster $A$ (i.e. active) or cluster $B$ (i.e. DFG-out). In this way, we have down-weighted contributions to the contact differences ${Δ c_{A B}^{i j} = c}_{A}^{i j} - c_{B}^{i j}$ that are due to overrepresentation of specific kinases in the PDB clusters, with the goal of using contact differences to represent conserved features of the conformational transition across many different kinases. Alignment gaps and unresolved residues were accounted for by excluding these counts in the summations. Only $| i - j | > 4$ were included in the calculation. The PDB clusters used to calculate these contact differences are described above. The contact frequency differences for both STKs and TKs were plotted on a contact map for visualization (Figure 1—figure supplement 2).

Potts model and threaded-energy calculation

Request a detailed protocol

Our Potts Hamiltonian was constructed from an MSA of protein kinase catalytic domains as previously described (McGee et al., 2021). The Potts Hamiltonian $H (S)$ takes the form –

H (S) = \sum_{i < j}^{L} J_{S_{i} S_{j}}^{i j} + \sum_{i}^{L} h_{S_{i}}^{i}

where $L$ is the number of columns in the MSA ( $L = 259$ ), $h$ is a matrix of self-interactions or ‘fields’, and $J$ is the coupling matrix which has the interpretation of co-evolutionary interactions between residues.

The Potts threaded-energy penalty $Δ E (S)$ for sequence $S$ to undergo the conformational transition $A \to B$ is calculated using contact frequency differences between the two conformational ensembles (Haldane et al., 2016) –

Δ E (S \in X) = - \sum_{i < j}^{L} J_{S_{i} S_{j}}^{i j} Δ c_{A B}^{i j} (X),

where $X$ represents a class or family of sequences for which sequence $S$ has membership, and $Δ c_{A B}^{i j} (X)$ represents the contact frequency difference between conformations $A$ and $B$ observed only for other sequences belonging to class $X$ (e.g. $X \equiv T K s$ or $X \equiv S T K s$ ; upper and lower triangle of Figure 1—figure supplement 2, respectively). As described previously (Haldane et al., 2016) the couplings ( $J_{S_{i} S_{j}}^{i j}$ ) and fields ( $h_{S_{i}}^{i}$ ) were transformed to the ‘zero-gauge’ prior to calculating $Δ E (S)$ .

Contributions to average shift in $Δ E_{P o t t s}$ between STKs and TKs ( $Δ Δ E$ ). Where $Δ E (S)$ is the Potts conformational penalty for sequence S to undergo the conformational change $A \to B$ , we define $Δ Δ E$ as the difference in average $Δ E$ between two groups of sequences $X$ and $Y$

Δ Δ E = {⟨Δ E⟩}_{X} - {⟨Δ E⟩}_{Y} .

To help interpret $Δ Δ E$ in a structural and coevolutionary context, we can write $Δ Δ E$ as a sum over position pairs $(i, j)$

Δ Δ E = \sum_{i < j}^{L} Δ Δ E^{i j} .

To evaluate this, we note that average $Δ E$ can be expressed as a sum over position pairs

{⟨ Δ E ⟩}_{X} = \sum_{i < j}^{L} {⟨ Δ E^{i j} ⟩}_{X}

where $⟨ Δ E^{i j} ⟩$ is calculated as follows

{⟨Δ E^{i j}⟩}_{X} = - \sum_{α} \sum_{β} f_{α β}^{i j} (X) J_{α β}^{i j} Δ c_{A B}^{i j} (X) .

$f_{α β}^{i j} (X)$ is the frequency (bivariate marginal) of residues $α$ and $β$ at positions $i$ and $j$ for sequences in the MSA which belong to group $X$ , which we calculate after applying the MSA-derived phylogenetic weights described above. Finally, by substituting Equation 10 back into Equation 8, we show how $Δ Δ E$ can be decomposed into contributions from individual residue pairs

Δ Δ E^{i j} = - \sum_{α} \sum_{β} J_{α β}^{i j} (f_{α β}^{i j} (X) Δ c_{A B}^{i j} (X) - f_{α β}^{i j} (Y) Δ c_{A B}^{i j} (Y)) .

By viewing the largest (most positive) $Δ Δ E^{i j}$ terms, where $X \equiv S T K s$ and $Y \equiv T K s$ in Equation 12, we are identifying position pairs that cause STKs to have higher penalties than TK in our Potts threading calculations for the active DFG-in to DFG-out conformational change (Figure 6—figure supplement 2).

Calculation of p-value for $Δ Δ E$

Request a detailed protocol

The quantity $Δ Δ E$ is a difference between two averages, ${⟨Δ E⟩}_{S T K} - {⟨Δ E⟩}_{T K}$ . Hypothesis testing to determine the statistical significance of this quantity was performed with respect to a null model where the populations of $Δ E s$ for STKs and TKs, from which our samples were drawn, are indistinguishable. To this end, a p-value was calculated for a t-statistic derived from Welch’s t-test (Welch, 1947), where $s$ is the standard error of average $Δ E$ –

\tilde{t} = \frac{{⟨Δ E⟩}_{S T K} - {⟨Δ E⟩}_{T K}}{\sqrt{s_{S T K}^{2} + s_{T K}^{2}}}

where the averages and standard errors are calculated after down-weighting each sequence as described above. This was done to lessen the effects of phylogenetic sampling bias from our MSA and ensure that $Δ Δ E$ captures general differences between TKs and STKs, rather than specific TK or STK families.

The one-tailed p-value was calculated using the cumulative t-distribution function generated in python using the SciPy package (Virtanen et al., 2020),

p = 1 - t_{c d f} (\tilde{t}, ν)

where the degrees of freedom for the t-distribution describing the combined population, $ν$ , was estimated via the Welch-Satterthwaite equation (Welch, 1947) from the degrees of freedom of the two samples $ν_{S T K}$ and $ν_{T K}$

ν = \frac{{(s_{S T K}^{2} N_{S T K}^{- 1} + s_{T K}^{2} N_{T K}^{- 1})}^{2}}{s_{S T K}^{4} N_{S T K}^{- 2} ν_{S T K}^{- 1} + s_{T K}^{4} N_{T K}^{- 2} ν_{T K}} .

where $N$ represents the effective number of STKs or TKs, which is an unbiased count of sequences in each dataset that can be obtained by summing the sample weights ( $N_{S T K} = 22,893$ , $N_{T K} = 1069$ ). From the calculation of $Δ Δ E = {⟨Δ E⟩}_{S T K} - {⟨Δ E⟩}_{T K} = 3.2$ , we determine the corresponding p-value to be less than $10^{- 15}$ , meaning it is highly unlikely for this large of a difference to be observed if the $Δ E s$ for TKs and STKs were randomly drawn from the same distribution rather than distinct distributions.

Enumeration of absolute binding free-energy simulations

Request a detailed protocol

In this work, we have performed all-atom MD simulations in explicit solvent for a total of 94 different kinase-inhibitor complexes to calculate ABFEs via the alchemical DDM. 74 of these free-energy calculations were guided by insights from the Potts model, specifically the patterns of Potts conformational penalties plotted in Figure 1D. 68 of these complexes correspond to the ABFEs plotted in Figure 3, of which 23 type-II inhibitors and 23 type-I inhibitors ABFEs are plotted for STKs in Figure 3A, and 22 type-II inhibitors ABFEs for TKs are plotted in Figure 3B. Six additional Potts-guided ABFEs corresponding to CSF1R, KIT, PDGFRA, MAPK14, and PTK2B are included in Figure 2. An additional 20 type-I and type-I ½ ABFEs were calculated as part of our benchmarking procedure described in the Methods section.

Double decoupling method setup

Request a detailed protocol

The double decoupling method (DDM), also known as an ‘alchemical’ method, was applied to compute ABFE ( ${Δ G}_{b i n d}^{o}$ ), as shown in Equation 16 (Deng et al., 2018; Sakae et al., 2020). This method computes the free energies of decoupling the inhibitor from the bulk solvent in the presence and absence of a receptor via a nonphysical thermodynamic cycle where the two end states are connected via the alchemical pathway. The starting holo-structures for ABFE calculations were taken from the available crystal structure. The absence of crystal structure prompted us to model the structure of the ligand into the active site of the kinase by superimposing over the binding pose of the available holo crystal structure.

{Δ G}_{b i n d}^{o} = {- Δ G}_{r e s t r a i n}^{b o u n d} {- Δ G}_{d e c o u p l e}^{b o u n d} {+ Δ G}_{r e s t r a i n}^{g a s} {+ Δ G}_{d e c o u p l e}^{b u l k}

Decoupling of the ligand was achieved by first turning off the coulombic intermolecular interactions followed by Lennard-Jones intermolecular interactions from both the legs. This allows DDM to estimate the free energy, i.e., in the presence of protein ( ${Δ G}_{d e c o u p l e}^{b o u n d}$ ) and absence of protein, i.e., in the bulk solvent, ( ${Δ G}_{d e c o u p l e}^{b u l k}$ ) as shown in Equations 17; 18.

{Δ G}_{d e c o u p l e}^{b o u n d} = {Δ G}_{d e c o u p l e - C o u l o m b}^{b o u n d} {+ Δ G}_{d e c o u p l e - L J}^{b o u n d}

{Δ G}_{d e c o u p l e}^{b u l k} = {Δ G}_{d e c o u p l e - C o u l o m b}^{b u l k} {+ Δ G}_{d e c o u p l e - L J}^{b u l k}

Substituting Equations 17; 18 into Equation 19 yields the estimated ( ${Δ G}_{b i n d}^{o}$ ) from DDM

{Δ G}_{b i n d}^{o} = {- Δ G}_{r e s t r a i n}^{b o u n d} + {{Δ Δ G}_{C o u l o m b} + {Δ Δ G}_{L J} + Δ G}_{r e s t r a i n}^{g a s}

where ${Δ Δ G}_{C o u l o m b}$ is the electrostatic energy contribution toward the total ABFE, and ${Δ Δ G}_{L J}$ is the non-polar energy contribution.

In this study, depending on the system’s convergence, either 20 or 31 total λs were used for decoupling the ligand from bulk solvent. For instance, either 5 λs with Δλ=0.5 or 11 λs with Δλ=0.1 were used for coulombic decoupling and 15 λs with Δλ=0.1 or 20 λs with Δλ=0.05 were used for decoupling Lennard-Jones interactions in the bulk solvent.

Similarly, depending on the convergence, either 30 or 42 total λs were used for decoupling ligand bound to protein. For instance, 11 or 12 non-uniformly distributed λs were used to restrain the ligand. Decoupling the coulombic interactions between ligand and protein was achieved by either using 4 λs with Δλ=0.25 or 10 λs with Δλ=0.1, whereas a large number of λs were used for decoupling Lennard-Jones interactions, i.e., 15 λs with Δλ=0.1 or 20 λs with Δλ=0.05 were used. The correction term developed by Rocklin and coworkers for treating charged ligands during DDM simulations was adopted (Rocklin et al., 2013). In this regard, it is well documented that the use of a finite-sized periodic solvent box during DDM simulations can lead to non-negligible electrostatic energy contribution toward the calculated total ABFE. Thus, calculated ( ${Δ G}_{b i n d}^{o}$ ) for charged ligand after addition of electrostatic correction term can be expressed as:

{Δ G}_{b i n d}^{o} = {- Δ G}_{r e s t r a i n}^{b o u n d} + {{Δ Δ G}_{C o u l o m b} + {Δ Δ G}_{L J} + {Δ G}_{e l e c t r o s t a t i c_c o r r e c t i o n}^{f i n i t e_s i z e} + Δ G}_{r e s t r a i n}^{g a s}

For a proper convergence during DDM simulations, the application of restrains is crucial. Herein, we have used six relative orthogonal restrains with harmonic potentials that include one distance, two angles, and three dihedral angles restrain between the ligand and the protein with a force constant of 10 kcal mol⁻¹ Å⁻² [deg⁻²]. At each λ, 10–30 ns of decoupling simulation via replica-exchange (Affentranger et al., 2006) were obtained to compute the ${Δ G}_{b i n d}^{o}$ over a well-converged trajectory.

Molecular dynamics setup

Request a detailed protocol

In this study, MD simulations were applied to compute the binding free-energy simulations via DDM. GROMACS-2018.8 (Abraham et al., 2015) was used as an MD engine for all simulations. The tleap module of AMBER16 was used to add the missing hydrogen atoms to the kinase enzymes. The system was solvated explicitly using TIP3P water boxes (Jorgensen et al., 1983) that extended at least 10 Å from the center of the system in each direction. The topology file for the kinase enzyme was created using the amber forcefield ff14SB (Maier et al., 2015). The AM1-BCC charge model (Jakalian et al., 2002) and general amber force field 2 (GAFF2) (Wang et al., 2004) were employed to parametrize different inhibitors used in this study. The overall charge of the system was maintained by adding a suitable number of counterions in each system. During the simulations, electrostatic interactions were computed using the particle mesh Ewald method (Essmann et al., 1995) with a cutoff and grid spacing of 10.0 and 1.0 Å, respectively. The NPT (constant Number of particles, Pressure, and Temperature) ensembles with a time step of 2 fs was used in the simulations.

Benchmarking calculations for absolute binding free-energy simulations

Request a detailed protocol

Accurate prediction of the ABFE difference between the Apo and Holo state of a protein is extremely important to achieve from force field-based MD simulations (Lin et al., 2013; Lovera et al., 2012). In this study, we used type-II inhibitors as probes to estimate ΔG_reorg via ABFE simulations, which is the excess free-energy between experimentally determined binding affinity and ΔG_bind calculated from ABFE. Target kinases, i.e., MAPK14, CDK2, and JNK1 bound with type-I and I 1/2 inhibitors in the active conformation states (where ΔG_reorg is expected to be close to zero) have been regularly used by the computational community as benchmark systems for absolute or relative binding free energy calculations (Wang et al., 2015; Lee et al., 2020; Goel et al., 2021; Kuhn et al., 2020; Gapsys et al., 2019; Khalak et al., 2021). Benchmarking calculations over multiple protein-ligand complexes show close agreement between calculated (ΔG_bind) and experimental (ΔG_exp) terms (Table 2, Figure 2—figure supplement 1). In the later part of benchmarking studies, we have included ABL1 bound to type-II inhibitors in the DFG-out/folded activation loop state. NMR studies have shown experimentally that ABL1 has a ΔG_reorg of 1.2 kcal/mol (Xie et al., 2020), which is consistent with our range of estimates of ΔG_reorg for this kinase (Table 2).

Table 2

Summary of benchmarking results.

All reported values are in units of kcal/mol. The root mean square error (RMSE) and mean unsigned error (MUE) were calculated with respect to the linear model $Δ G_{b i n d}^{A B F E} + Δ G_{r e o r g}$ , where $Δ G_{r e o r g}$ is calculated as $⟨ Δ G_{e x p} - Δ G_{b i n d}^{A B F E} ⟩$ . The average difference between $Δ G_{e x p}$ and $Δ G_{b i n d}^{A B F E}$ is shown in the last column.

Kinase	# Compounds	MUE	RMSE	$\| ⟨ Δ G_{e x p} - Δ G_{b i n d}^{A B F E} ⟩ \|$
CDK2	6 (type-I)	0.71	1.03	1.89
JNK1	5 (type-I)	0.37	0.44	1.29
MAPK14	9 (type-1¹/₂)	0.47	0.64	0.2
ABL1	6 (type-II)	1.47	1.57	1.26

Potts-guided target selection for absolute binding free-energy simulations

Request a detailed protocol

Our Potts threaded-energy calculations were used alongside experimental type-II binding data from the large-scale assay by Davis et al., 2011 to identify kinase targets that are likely to have very large or very small ΔG_reorg. As described in the main text, all-atom MD simulations to calculate ABFEs of type-II inhibitors can be used alongside experimental binding affinities to calculate the free-energy cost for a kinase to reorganize to the DFG-out/folded activation loop conformation. This relation, ΔG_reorg = ΔG_exp − ΔG_bind(ABFE), gives the free-energy cost to reorganize in physical energy units (kcal/mol) and can be used to approximate a scale for the Potts statistical energy differences provided one sample a sufficient range of ΔG_reorg and ΔE_Potts. However, ABFE simulations are much more computationally demanding than the Potts threading calculation, which we sought to mitigate by choosing kinase simulation targets which are likely to provide a strong signal, guided by the Potts model. We direct the reader to Table 1, which contains the Potts penalties and type-II hit rates for the targets of interest. For comparison, Figure 1A and Figure 5 provide the overall distributions of hit rates and Potts penalties for TKs and STKs.

A significant challenge for our target selection was the limited availability of type-II inhibitors co-crystallized against STKs which have (a) Potts penalties and type-II hit rates that predict very high ΔG_reorg, (b) experimental binding affinities available in the literature in the form of IC50, Ki, or Kd, (c) availability of protein-ligand co-crystallized structure(s), and (d) type-II inhibitor complex systems where the activation loop appears to have undergone a large-scale ‘folding’ conformational change relative to the active ‘extended’ conformation. STK complexes that satisfy all four criteria appear to be sparse, which is consistent with the notion that kinases with large reorganization penalties are more difficult to crystallize in the classical DFG-out conformation (Haldane et al., 2016). However, for some STKs with very high Potts threaded-energy penalties (e.g., MELK) there has been significant medicinal chemistry efforts to design potent type-II inhibitors and structurally characterize their complexes using x-ray crystallography. Using type-II co-crystal structures that cover five different STKs with high Potts penalties (Table 1) and five different TKs with low Potts penalties, we were able to sample a wide range of ΔG_reorg from a total of 45 ABFE simulations covering 45 type-II inhibitor complexes and 10 different kinase targets. These simulations and subsequent calculations of ΔG_reorg for each kinase resulted in a strong correlation with the Potts threaded energy scores (Figure 5), allowing us to establish a scale for the Potts energies in kcal/mol. We have provided detailed results from the ABFE simulations of these 10 kinase targets in the form of supplementary figures (Figure 3—figure supplement 1) where the average ΔG_reorg for each kinase is visualized as the y-intercept of a linear regression with the slope constrained to one.

Data availability

Request a detailed protocol

Values of $Δ G_{b i n d}^{A B F E}$ from all ABFE simulations described in this work, including benchmarking calculations (94 simulations in total), are provided in the form of supplementary tables (Figure 2—source data 1 and Figure 3—source data 1). Potts $Δ E s$ , type-II hit rates computed from Shan et al., 2013, the identity of gatekeeper residues and corresponding van der Waals volumes in Å (Stancik et al., 2018), and the classification of human kinases as TKs or STKs were provided in a separate supplementary table (Figure 1—source data 1).

Data availability

Our computational study makes use of experimental data from the literature, which we extracted and curated manually rather than relying on any specific database. Any experimental data used can be found in our supporting information in the form of tables alongside appropriate citations. A large set of experimental "hit rates" were derived from binding affinities available from Davis et al., 2011. The data used to generate various plots in the main text can be found in tables throughout the supporting information, as well as a distinct "supplementary table" which we provide. The Mi3-GPU (Haldane and Levy, 2021) source code required to reproduce the Potts model employed in this manuscript can be found at the following link: https://github.com/ahaldane/Mi3-GPU (Haldane and avikbiswas, 2021, v1.1, copy archived at swh:1:rev:b8fd4aa67bb2531fdc60e3a00fed6f80c8aceb49).

References

1. Abraham MJ
2. Murtola T
3. Schulz R
4. Páll S
5. Smith JC
6. Hess B
7. Lindahl E
(2015) GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers
SoftwareX 1–2:19–25.

https://doi.org/10.1016/j.softx.2015.06.001
- Google Scholar
(2006) A novel hamiltonian replica exchange MD protocol to enhance protein conformational space sampling
Journal of Chemical Theory and Computation 2:217–228.

https://doi.org/10.1021/ct050250b
- PubMed
- Google Scholar
(2011) Comprehensive assay of kinase catalytic activity reveals features of kinase inhibitor selectivity
Nature Biotechnology 29:1039–1045.

https://doi.org/10.1038/nbt.2017
- PubMed
- Google Scholar
1. Apperley JF
(2007) Part I: mechanisms of resistance to imatinib in chronic myeloid leukaemia
The Lancet. Oncology 8:1018–1029.

https://doi.org/10.1016/S1470-2045(07)70342-X
- PubMed
- Google Scholar
1. Azam M
2. Seeliger MA
3. Gray NS
4. Kuriyan J
5. Daley GQ
(2008) Activation of tyrosine kinases by mutation of the gatekeeper threonine
Nature Structural & Molecular Biology 15:1109–1118.

https://doi.org/10.1038/nsmb.1486
- PubMed
- Google Scholar
(2008) Dasatinib treatment can overcome imatinib and nilotinib resistance in CML patient carrying F359I mutation of bcr-abl oncogene
Journal of Applied Genetics 49:201–203.

https://doi.org/10.1007/BF03195613
- PubMed
- Google Scholar
(2016) How do protein kinases take a selfie (autophosphorylate)?
Trends in Biochemical Sciences 41:938–953.

https://doi.org/10.1016/j.tibs.2016.08.006
- PubMed
- Google Scholar
1. Bosc N
2. Wroblowski B
3. Aci-Sèche S
4. Meyer C
5. Bonnet P
(2015) A proteometric analysis of human kinome: insight into discriminant conformation-dependent residues
ACS Chemical Biology 10:2827–2840.

https://doi.org/10.1021/acschembio.5b00555
- PubMed
- Google Scholar
1. Chen L
2. Marsiglia WM
3. Chen H
4. Katigbak J
5. Erdjument-Bromage H
6. Kemble DJ
7. Fu L
8. Ma J
9. Sun G
10. Zhang Y
11. Liang G
12. Neubert TA
13. Li X
14. Traaseth NJ
15. Mohammadi M
(2020) Molecular basis for receptor tyrosine kinase A-loop tyrosine transphosphorylation
Nature Chemical Biology 16:267–277.

https://doi.org/10.1038/s41589-019-0455-7
- PubMed
- Google Scholar
1. Cournia Z
2. Allen BK
3. Beuming T
4. Pearlman DA
5. Radak BK
6. Sherman W
(2020) Rigorous free energy simulations in virtual screening
Journal of Chemical Information and Modeling 60:4153–4169.

https://doi.org/10.1021/acs.jcim.0c00116
- PubMed
- Google Scholar
1. Davis MI
2. Hunt JP
3. Herrgard S
4. Ciceri P
5. Wodicka LM
6. Pallares G
7. Hocker M
8. Treiber DK
9. Zarrinkar PP
(2011) Comprehensive analysis of kinase inhibitor selectivity
Nature Biotechnology 29:1046–1051.

https://doi.org/10.1038/nbt.1990
- PubMed
- Google Scholar
1. Deng N
2. Zhang P
3. Cieplak P
4. Lai L
(2011) Elucidating the energetics of entropically driven protein-ligand association: calculations of absolute binding free energy and entropy
The Journal of Physical Chemistry. B 115:11902–11910.

https://doi.org/10.1021/jp204047b
- PubMed
- Google Scholar
1. Deng N
2. Cui D
3. Zhang BW
4. Xia J
5. Cruz J
6. Levy R
(2018) Comparing alchemical and physical pathway methods for computing the absolute binding free energy of charged ligands
Physical Chemistry Chemical Physics 20:17081–17092.

https://doi.org/10.1039/c8cp01524d
- PubMed
- Google Scholar
(2012) The structural basis for control of eukaryotic protein kinases
Annual Review of Biochemistry 81:587–613.

https://doi.org/10.1146/annurev-biochem-052410-090317
- PubMed
- Google Scholar
1. Essmann U
2. Perera L
3. Berkowitz ML
4. Darden T
5. Lee H
6. Pedersen LG
(1995) A smooth particle mesh ewald method
J Chem Phys 103:8577–8593.

https://doi.org/10.1063/1.470117
- Google Scholar
(2019) Large scale relative protein ligand binding affinities using non-equilibrium alchemy
Chemical Science 11:1140–1152.

https://doi.org/10.1039/c9sc03754c
- PubMed
- Google Scholar
(2008) Knowledge based prediction of ligand binding modes and rational inhibitor design for kinase drug discovery
Journal of Medicinal Chemistry 51:5149–5171.

https://doi.org/10.1021/jm800475y
- PubMed
- Google Scholar
(1997) The statistical-thermodynamic basis for computation of binding affinities: a critical review
Biophysical Journal 72:1047–1069.

https://doi.org/10.1016/S0006-3495(97)78756-3
- PubMed
- Google Scholar
1. Goel H
2. Hazel A
3. Ustach VD
4. Jo S
5. Yu W
6. MacKerell AD
(2021) Rapid and accurate estimation of protein-ligand relative binding affinities using site-identification by ligand competitive saturation
Chemical Science 12:8844–8858.

https://doi.org/10.1039/d1sc01781k
- PubMed
- Google Scholar
(2022) Alchemical free energy methods applied to complexes of the first bromodomain of BRD4
Journal of Chemical Information and Modeling 62:1458–1470.

https://doi.org/10.1021/acs.jcim.1c01229
- PubMed
- Google Scholar
1. Haldane A
2. Flynn WF
3. He P
4. Vijayan RSK
5. Levy RM
(2016) Structural propensities of kinase family proteins from a potts model of residue co-variation
Protein Science 25:1378–1384.

https://doi.org/10.1002/pro.2954
- PubMed
- Google Scholar
Software
1. Haldane A
2. avikbiswas
(2021) Mi3-GPU, version swh:1:rev:b8fd4aa67bb2531fdc60e3a00fed6f80c8aceb49
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:5904df828b1f44cc243930567e984689b6c05e4e;origin=https://github.com/ahaldane/Mi3-GPU;visit=swh:1:snp:c9df01efd95662be6d37f7be55b23cce9eadd64c;anchor=swh:1:rev:b8fd4aa67bb2531fdc60e3a00fed6f80c8aceb49
1. Haldane A
2. Levy RM
(2021) Mi3-GPU: MCMC-based inverse Ising inference on GPUs for protein covariation analysis
Computer Physics Communications 260:107312.

https://doi.org/10.1016/j.cpc.2020.107312
- PubMed
- Google Scholar
1. Han S
2. Mistry A
3. Chang JS
4. Cunningham D
5. Griffor M
6. Bonnette PC
7. Wang H
8. Chrunyk BA
9. Aspnes GE
10. Walker DP
11. Brosius AD
12. Buckbinder L
(2009) Structural characterization of proline-rich tyrosine kinase 2 (PYK2) reveals a unique (DFG-out) conformation and enables inhibitor design
Journal of Biological Chemistry 284:13193–13201.

https://doi.org/10.1074/jbc.M809038200
- PubMed
- Google Scholar
(2013) Sequence determinants of a specific inactive protein kinase conformation
Chemistry & Biology 20:806–815.

https://doi.org/10.1016/j.chembiol.2013.05.005
- PubMed
- Google Scholar
1. Hauser K
2. Negron C
3. Albanese SK
4. Ray S
5. Steinbrecher T
6. Abel R
7. Chodera JD
8. Wang L
(2018) Predicting resistance of clinical abl mutations to targeted kinase inhibitors using alchemical free-energy calculations
Communications Biology 1:70.

https://doi.org/10.1038/s42003-018-0075-x
- PubMed
- Google Scholar
(2022) Addressing intersite coupling unlocks large combinatorial chemical spaces for alchemical free energy methods
Journal of Chemical Theory and Computation 18:2114–2123.

https://doi.org/10.1021/acs.jctc.1c00948
- PubMed
- Google Scholar
1. Heinzelmann G
2. Gilson MK
(2021) Automation of absolute protein-ligand binding free energy calculations for docking refinement and compound evaluation
Scientific Reports 11:1116.

https://doi.org/10.1038/s41598-020-80769-1
- PubMed
- Google Scholar
Preprint
1. Hopf TA
2. Ingraham JB
3. Poelwijk FJ
4. Springer M
5. Sander C
6. Marks DS
(2015) Quantification of the Effect of Mutations Using a Global Probability Model of Natural Sequence Variation
arXiv.

https://arxiv.org/abs/1510.04612
- Google Scholar
1. Hopf TA
2. Ingraham JB
3. Poelwijk FJ
4. Schärfe CPI
5. Springer M
6. Sander C
7. Marks DS
(2017) Mutation effects predicted from sequence co-variation
Nature Biotechnology 35:128–135.

https://doi.org/10.1038/nbt.3769
- PubMed
- Google Scholar
(1994) Crystal structure of the tyrosine kinase domain of the human insulin receptor
Nature 372:746–754.

https://doi.org/10.1038/372746a0
- Google Scholar
(2002) Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. parameterization and validation
Journal of Computational Chemistry 23:1623–1641.

https://doi.org/10.1002/jcc.10128
- PubMed
- Google Scholar
(1996) Active and inactive protein kinases: structural basis for regulation
Cell 85:149–158.

https://doi.org/10.1016/s0092-8674(00)81092-2
- PubMed
- Google Scholar
(1983) Comparison of simple potential functions for simulating liquid water
J Chem Phys 79:926–935.

https://doi.org/10.1063/1.445869
- Google Scholar
1. Khalak Y
2. Tresadern G
3. Aldeghi M
4. Baumann HM
5. Mobley DL
6. de Groot BL
7. Gapsys V
(2021) Alchemical absolute protein-ligand binding free energies for drug design
Chemical Science 12:13958–13971.

https://doi.org/10.1039/d1sc03472c
- PubMed
- Google Scholar
1. Klaeger S
2. Heinzlmeir S
3. Wilhelm M
4. Polzer H
5. Vick B
6. Koenig PA
7. Reinecke M
8. Ruprecht B
9. Petzoldt S
10. Meng C
11. Zecha J
12. Reiter K
13. Qiao H
14. Helm D
15. Koch H
16. Schoof M
17. Canevari G
18. Casale E
19. Depaolini SR
20. Feuchtinger A
21. Wu Z
22. Schmidt T
23. Rueckert L
24. Becker W
25. Huenges J
26. Garz AK
27. Gohlke BO
28. Zolg DP
29. Kayser G
30. Vooder T
31. Preissner R
32. Hahne H
33. Tõnisson N
34. Kramer K
35. Götze K
36. Bassermann F
37. Schlegl J
38. Ehrlich HC
39. Aiche S
40. Walch A
41. Greif PA
42. Schneider S
43. Felder ER
44. Ruland J
45. Médard G
46. Jeremias I
47. Spiekermann K
48. Kuster B
(2017) The target landscape of clinical kinase drugs
Science 358:eaan4368.

https://doi.org/10.1126/science.aan4368
- PubMed
- Google Scholar
1. Kuhn M
2. Firth-Clark S
3. Tosco P
4. Mey A
5. Mackey M
6. Michel J
(2020) Assessment of binding affinity via alchemical free-energy calculations
Journal of Chemical Information and Modeling 60:3120–3130.

https://doi.org/10.1021/acs.jcim.0c00165
- PubMed
- Google Scholar
Preprint
(2012) Using Sequence Alignments to Predict Protein Structure and Stability With High Accuracy
arXiv.

https://arxiv.org/abs/1207.2484
- Google Scholar
1. Lee TS
2. Allen BK
3. Giese TJ
4. Guo Z
5. Li P
6. Lin C
7. McGee TD
8. Pearlman DA
9. Radak BK
10. Tao Y
11. Tsai HC
12. Xu H
13. Sherman W
14. York DM
(2020) Alchemical binding free energy calculations in AMBER20: advances and best practices for drug discovery
Journal of Chemical Information and Modeling 60:5595–5623.

https://doi.org/10.1021/acs.jcim.0c00613
- PubMed
- Google Scholar
1. Lemmon MA
2. Schlessinger J
(2010) Cell signaling by receptor tyrosine kinases
Cell 141:1117–1134.

https://doi.org/10.1016/j.cell.2010.06.011
- PubMed
- Google Scholar
1. Levinson NM
2. Kuchment O
3. Shen K
4. Young MA
5. Koldobskiy M
6. Karplus M
7. Cole PA
8. Kuriyan J
(2006) A src-like inactive conformation in the abl tyrosine kinase domain
PLOS Biology 4:e144.

https://doi.org/10.1371/journal.pbio.0040144
- PubMed
- Google Scholar
1. Li Z
2. Li X
3. Huang YY
4. Wu Y
5. Liu R
6. Zhou L
7. Lin Y
8. Wu D
9. Zhang L
10. Liu H
11. Xu X
12. Yu K
13. Zhang Y
14. Cui J
15. Zhan CG
16. Wang X
17. Luo Hb
(2020) Identify potent SARS-cov-2 main protease inhibitors via accelerated free energy perturbation-based virtual screening of existing drugs
PNAS 117:27381–27387.

https://doi.org/10.1073/pnas.2010470117
- PubMed
- Google Scholar
1. Lin YL
2. Meng Y
3. Jiang W
4. Roux B
(2013) Explaining why gleevec is a specific and potent inhibitor of abl kinase
PNAS 110:1664–1669.

https://doi.org/10.1073/pnas.1214330110
- PubMed
- Google Scholar
1. Lin YL
2. Meng Y
3. Huang L
4. Roux B
(2014) Computational study of gleevec and G6G reveals molecular determinants of kinase inhibitor selectivity
Journal of the American Chemical Society 136:14753–14762.

https://doi.org/10.1021/ja504146x
- PubMed
- Google Scholar
1. Liu Y
2. Shah K
3. Yang F
4. Witucki L
5. Shokat KM
(1998) A molecular gate which controls unnatural ATP analogue recognition by the tyrosine kinase v-src
Bioorganic & Medicinal Chemistry 6:1219–1226.

https://doi.org/10.1016/s0968-0896(98)00099-6
- PubMed
- Google Scholar
1. Liu Y
2. Gray NS
(2006) Rational design of inhibitors that bind to inactive kinase conformations
Nature Chemical Biology 2:358–364.

https://doi.org/10.1038/nchembio799
- PubMed
- Google Scholar
1. Lovera S
2. Sutto L
3. Boubeva R
4. Scapozza L
5. Dölker N
6. Gervasio FL
(2012) The different flexibility of c-src and c-abl kinases regulates the accessibility of a druggable inactive conformation
Journal of the American Chemical Society 134:2496–2499.

https://doi.org/10.1021/ja210751t
- PubMed
- Google Scholar
(2015) Towards a molecular understanding of the link between imatinib resistance and kinase conformational dynamics
PLOS Computational Biology 11:e1004578.

https://doi.org/10.1371/journal.pcbi.1004578
- PubMed
- Google Scholar
Book
1. Lunt B
2. Szurmant H
3. Procaccini A
4. Hoch JA
5. Hwa T
6. Weigt M
(2010)
Inference of Direct Residue Contacts in Two-Component Signaling (1st ed)

Elsevier Inc.
- Google Scholar
(2015) Ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99sb
Journal of Chemical Theory and Computation 11:3696–3713.

https://doi.org/10.1021/acs.jctc.5b00255
- PubMed
- Google Scholar
1. Manning G
(1995)
The protein kinase complement of the human genome - supplemental information

Bulletin of the World Health Organization 73:7–14.
- Google Scholar
1. McGee F
2. Hauri S
3. Novinger Q
4. Vucetic S
5. Levy RM
6. Carnevale V
7. Haldane A
(2021) The generative capacity of probabilistic protein sequence models
Nature Communications 12:6302.

https://doi.org/10.1038/s41467-021-26529-9
- PubMed
- Google Scholar
(2018) Coral: clear and customizable visualization of human kinome data
Cell Systems 7:347–350.

https://doi.org/10.1016/j.cels.2018.07.001
- PubMed
- Google Scholar
1. Miller S
2. Janin J
3. Lesk AM
4. Chothia C
(1987) Interior and surface of monomeric proteins
Journal of Molecular Biology 196:641–656.

https://doi.org/10.1016/0022-2836(87)90038-6
- PubMed
- Google Scholar
1. Miller WT
(2012) Tyrosine kinase signaling and the emergence of multicellularity
Biochimica et Biophysica Acta 1823:1053–1057.

https://doi.org/10.1016/j.bbamcr.2012.03.009
- PubMed
- Google Scholar
1. Modi V
2. Dunbrack RL
(2019a) A structurally-validated multiple sequence alignment of 497 human protein kinase domains
Scientific Reports 9:19790.

https://doi.org/10.1038/s41598-019-56499-4
- PubMed
- Google Scholar
1. Modi V
2. Dunbrack RL
(2019b) Defining a new nomenclature for the structures of active and inactive kinases
PNAS 116:6818–6827.

https://doi.org/10.1073/pnas.1814279116
- PubMed
- Google Scholar
(2014) Coevolutionary information, protein folding landscapes, and the thermodynamics of natural selection
PNAS 111:12408–12413.

https://doi.org/10.1073/pnas.1413575111
- PubMed
- Google Scholar
1. Nolen B
2. Taylor S
3. Ghosh G
(2004) Regulation of protein kinases; controlling activity through activation segment conformation
Molecular Cell 15:661–675.

https://doi.org/10.1016/j.molcel.2004.08.024
- PubMed
- Google Scholar
1. Ozkirimli E
2. Post CB
(2006) Src kinase activation: a switched electrostatic network
Protein Science 15:1051–1062.

https://doi.org/10.1110/ps.051999206
- PubMed
- Google Scholar
1. Pargellis C
2. Tong L
3. Churchill L
4. Cirillo PF
5. Gilmore T
6. Graham AG
7. Grob PM
8. Hickey ER
9. Moss N
10. Pav S
11. Regan J
(2002) Inhibition of P38 MAP kinase by utilizing a novel allosteric binding site
Nature Structural Biology 9:268–272.

https://doi.org/10.1038/nsb770
- PubMed
- Google Scholar
Book
1. pymol
(2015)
The Pymol Molecular Graphics System, Version 2.0

Schrödinger, LLC.
- Google Scholar
(2019) Absolute free energy of binding calculations for macrophage migration inhibitory factor in complex with a druglike inhibitor
The Journal of Physical Chemistry. B 123:8675–8685.

https://doi.org/10.1021/acs.jpcb.9b07588
- PubMed
- Google Scholar
(2013) Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects
J Chem Phys 139:184103.

https://doi.org/10.1063/1.4826261
- PubMed
- Google Scholar
1. Sakae Y
2. Zhang BW
3. Levy RM
4. Deng N
(2020) Absolute protein binding free energy simulations for ligands with multiple poses, a thermodynamic path that avoids exhaustive enumeration of the poses
Journal of Computational Chemistry 41:56–68.

https://doi.org/10.1002/jcc.26078
- PubMed
- Google Scholar
(2016) High-throughput proteomics reveals the unicellular roots of animal phosphosignaling and cell differentiation
Developmental Cell 39:186–197.

https://doi.org/10.1016/j.devcel.2016.09.019
- PubMed
- Google Scholar
1. Shan Y
2. Arkhipov A
3. Kim ET
4. Pan AC
5. Shaw DE
(2013) Transitions to catalytically inactive conformations in EGFR kinase
PNAS 110:7270–7275.

https://doi.org/10.1073/pnas.1220843110
- PubMed
- Google Scholar
(2018) Serine/threonine protein kinases from bacteria, archaea and eukarya share a common evolutionary origin deeply rooted in the tree of life
Journal of Molecular Biology 430:27–32.

https://doi.org/10.1016/j.jmb.2017.11.004
- PubMed
- Google Scholar
1. Sun Q
2. Biswas A
3. Vijayan RSK
4. Craveur P
5. Forli S
6. Olson AJ
7. Castaner AE
8. Kirby KA
9. Sarafianos SG
10. Deng N
11. Levy R
(2022) Structure-based virtual screening workflow to identify antivirals targeting HIV-1 capsid
Journal of Computer-Aided Molecular Design 36:193–203.

https://doi.org/10.1007/s10822-022-00446-5
- PubMed
- Google Scholar
(1995) How do protein kinases discriminate between serine/threonine and tyrosine? structural insights from the insulin receptor protein‐tyrosine kinase
The FASEB Journal 9:1255–1266.

https://doi.org/10.1096/fasebj.9.13.7557015
- PubMed
- Google Scholar
1. Taylor SS
2. Kornev AP
(2011) Protein kinases: evolution of dynamic regulatory proteins
Trends in Biochemical Sciences 36:65–77.

https://doi.org/10.1016/j.tibs.2010.09.006
- PubMed
- Google Scholar
1. Tokarski JS
2. Newitt JA
3. Chang CYJ
4. Cheng JD
5. Wittekind M
6. Kiefer SE
7. Kish K
8. Lee FYF
9. Borzillerri R
10. Lombardo LJ
11. Xie D
12. Zhang Y
13. Klei HE
(2006) The structure of dasatinib (BMS-354825) bound to activated abl kinase domain elucidates its inhibitory activity against imatinib-resistant abl mutants
Cancer Research 66:5790–5797.

https://doi.org/10.1158/0008-5472.CAN-05-4187
- PubMed
- Google Scholar
(2014) KLIFS: a knowledge-based structural database to navigate kinase-ligand interaction space
Journal of Medicinal Chemistry 57:249–277.

https://doi.org/10.1021/jm400378w
- PubMed
- Google Scholar
1. Vijayan RSK
2. He P
3. Modi V
4. Duong-Ly KC
5. Ma H
6. Peterson JR
7. Dunbrack RL
8. Levy RM
(2015) Conformational analysis of the DFG-out kinase motif and biochemical profiling of structurally validated type II inhibitors
Journal of Medicinal Chemistry 58:466–479.

https://doi.org/10.1021/jm501603h
- PubMed
- Google Scholar
1. Virtanen P
2. Gommers R
3. Oliphant TE
4. Haberland M
5. Reddy T
6. Cournapeau D
7. Burovski E
8. Peterson P
9. Weckesser W
10. Bright J
11. van der Walt SJ
12. Brett M
13. Wilson J
14. Millman KJ
15. Mayorov N
16. Nelson ARJ
17. Jones E
18. Kern R
19. Larson E
20. Carey CJ
21. Polat İ
22. Feng Y
23. Moore EW
24. VanderPlas J
25. Laxalde D
26. Perktold J
27. Cimrman R
28. Henriksen I
29. Quintero EA
30. Harris CR
31. Archibald AM
32. Ribeiro AH
33. Pedregosa F
34. van Mulbregt P
35. SciPy 1.0 Contributors
(2020) SciPy 1.0: fundamental algorithms for scientific computing in python
Nature Methods 17:261–272.

https://doi.org/10.1038/s41592-019-0686-2
- PubMed
- Google Scholar
1. Wang J
2. Wolf RM
3. Caldwell JW
4. Kollman PA
5. Case DA
(2004) Development and testing of a general amber force field
Journal of Computational Chemistry 25:1157–1174.

https://doi.org/10.1002/jcc.20035
- PubMed
- Google Scholar
1. Wang L
2. Wu Y
3. Deng Y
4. Kim B
5. Pierce L
6. Krilov G
7. Lupyan D
8. Robinson S
9. Dahlgren MK
10. Greenwood J
11. Romero DL
12. Masse C
13. Knight JL
14. Steinbrecher T
15. Beuming T
16. Damm W
17. Harder E
18. Sherman W
19. Brewer M
20. Wester R
21. Murcko M
22. Frye L
23. Farid R
24. Lin T
25. Mobley DL
26. Jorgensen WL
27. Berne BJ
28. Friesner RA
29. Abel R
(2015) Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field
Journal of the American Chemical Society 137:2695–2703.

https://doi.org/10.1021/ja512751q
- PubMed
- Google Scholar
1. Weigt M
2. White RA
3. Szurmant H
4. Hoch JA
5. Hwa T
(2009) Identification of direct residue contacts in protein-protein interaction by message passing
PNAS 106:67–72.

https://doi.org/10.1073/pnas.0805923106
- PubMed
- Google Scholar
1. Welch ABL
(1947) The generalization of ` student ’ s ’ problem when several different population variances are involved published by: biometrika trust stable
Biometrika 34:28–35.

https://doi.org/10.1093/biomet/34.1-2.28
- PubMed
- Google Scholar
1. Wu H
2. Huang H
3. Post CB
(2020) All-atom adaptively biased path optimization of src kinase conformational inactivation: switched electrostatic network in the concerted motion of αC helix and the activation loop
J Chem Phys 153:175101.

https://doi.org/10.1063/5.0021603
- PubMed
- Google Scholar
1. Xie T
2. Saleh T
3. Rossi P
4. Kalodimos CG
(2020) Conformational states dynamically populated by a kinase determine its function
Science 370:29–37.

https://doi.org/10.1126/science.abc2754
- PubMed
- Google Scholar
1. Yun CH
2. Mengwasser KE
3. Toms AV
4. Woo MS
5. Greulich H
6. Wong KK
7. Meyerson M
8. Eck MJ
(2008) The T790M mutation in EGFR kinase causes drug resistance by increasing the affinity for ATP
PNAS 105:2070–2075.

https://doi.org/10.1073/pnas.0709662105
- PubMed
- Google Scholar
1. Zhao Z
2. Wu H
3. Wang L
4. Liu Y
5. Knapp S
6. Liu Q
7. Gray NS
(2014) Exploration of type II binding mode: a privileged approach for kinase inhibitor focused drug discovery?
ACS Chemical Biology 9:1230–1241.

https://doi.org/10.1021/cb500129t
- PubMed
- Google Scholar
1. Zheng J
2. Knighton DR
3. ten Eyck LF
4. Karlsson R
5. Xuong N
6. Taylor SS
7. Sowadski JM
(1993) Crystal structure of the catalytic subunit of camp-dependent protein kinase complexed with mgatp and peptide inhibitor
Biochemistry 32:2154–2161.

https://doi.org/10.1021/bi00060a005
- PubMed
- Google Scholar
(2010) Through the “gatekeeper door”: exploiting the active kinase conformation
Journal of Medicinal Chemistry 53:2681–2694.

https://doi.org/10.1021/jm901443h
- PubMed
- Google Scholar

Article and author information

Author details

Joan Gizzio
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Chemistry, Temple University, Philadelphia, United States
Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Writing - original draft, Writing - review and editing

Contributed equally with
Abhishek Thakur

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4127-2294
Abhishek Thakur
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Chemistry, Temple University, Philadelphia, United States
Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Writing - original draft, Writing - review and editing

Contributed equally with
Joan Gizzio

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4827-7602
Allan Haldane
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Physics, Temple University, Philadelphia, United States
Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Investigation, Methodology, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8343-1994
Ronald M Levy
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Chemistry, Temple University, Philadelphia, United States
Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing - review and editing

For correspondence
ronlevy@temple.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8696-5177

Funding

National Institutes of Health (R35-GM132090)

Joan Gizzio
Abhishek Thakur
Allan Haldane
Ronald M Levy

National Institutes of Health (OD020095)

Joan Gizzio
Abhishek Thakur
Allan Haldane
Ronald M Levy

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This research was supported by National Institutes of Health grant number R35-GM132090, and by NIH Computer Equipment Grant (OD020095). Gratitude is also expressed to the OWLSNEST high performance cluster at Temple University for its computing support in this project. We thank Shima Arasteh for helpful discussions related to kinase conformational states and alchemical free-energy simulations.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.