# Insight on mutation-induced resistance from molecular dynamics simulations of the native and mutated CSF-1R and KIT

Priscilla da Silva Figueiredo Celestino Gomes, Isaure Chauvot de Beauchêne, Nicolas Panel, Sophie Lopes, Paulo de Sepulveda, Pedro G Pascutti, E Solary, Luba Tchertanov

► To cite this version:

Priscilla da Silva Figueiredo Celestino Gomes, Isaure Chauvot de Beauchêne, Nicolas Panel, Sophie Lopes, Paulo de Sepulveda, et al.. Insight on mutation-induced resistance from molecular dynamics simulations of the native and mutated CSF-1R and KIT. PLoS ONE, 2016, 11 (7), e0160165. doi:10.1371/journal.pone.0160165

HAL Id: hal-01505869

https://hal.science/hal-01505869v1

Submitted on 4 Oct 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

CrossMark
click for updates

RESEARCH ARTICLE

# Insight on Mutation-Induced Resistance from Molecular Dynamics Simulations of the Native and Mutated CSF-1R and KIT

Priscila Da Silva Figueiredo Celestino Gomes¹,², Isaure Chauvot De Beauchêne¹*,
Nicolas Panel¹, Sophie Lopez³, Paulo De Sepulveda³, Pedro Geraldo Pascutti²,
Eric Solary⁴,⁵, Luba Tchertanov¹,⁶*

1 Laboratoire de Biologie et Pharmacologie Appliquée (LBPA), ENS Cachan, CNRS, Université, Paris-Saclay, 94235 Cachan, France,
2 Laboratório de Modelagem e Dinâmica Molecular, Instituto de Biofísica Carlos Chagas Filho Universidade Federal do Rio de Janeiro, 373 av. Carlos Chagas Filho, 21941-902, Rio de Janeiro, Brazil,
3 INSERM U1068, Centre de Recherche en Cancérologie de Marseille, 13009 Marseille, France,
4 Inserm UMR1170, Gustave Roussy, 114 rue Edouard Vaillant, 94805 Villejuif cedex, France,
5 Paris-Sud University, Faculty of Medicine, Le Kremlin-Bicêtre, France,
6 Centre de Mathématiques et de Leurs Applications, (CMLA), ENS Cachan, CNRS, Université, Paris-Saclay, 94235 Cachan, France

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Da Silva Figueiredo Celestino Gomes P, Chauvot De Beauchêne I, Panel N, Lopez S, De Sepulveda P, Geraldo Pascutti P, et al. (2016) Insight on Mutation-Induced Resistance from Molecular Dynamics Simulations of the Native and Mutated CSF-1R and KIT. PLoS ONE 11(7): e0160165. doi:10.1371/journal.pone.0160165

**Editor:** L. Michel Espinoza-Fonseca, University of Minnesota, UNITED STATES

**Received:** March 3, 2016

**Accepted:** July 14, 2016

**Published:** July 28, 2016

**Copyright:** © 2016 Da Silva Figueiredo Celestino Gomes et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** The authors have no support or funding to report.

**Competing Interests:** The authors have declared that no competing interests exist.

## Abstract

The receptors tyrosine kinases (RTKs) for the colony stimulating factor-1, CSF-1R, and for the stem cell factor, SCFR or KIT, are important mediators of signal transduction. The abnormal function of these receptors, promoted by gain-of-function mutations, leads to their constitutive activation, associated with cancer or other proliferative diseases. A secondary effect of the mutations is the alteration of receptors' sensitivity to tyrosine kinase inhibitors, compromising effectiveness of these molecules in clinical treatment. In particular, the mutation V560G in KIT increases its sensitivity to Imatinib, while the D816V in KIT, and D802V in CSF-1R, triggers resistance to the drug. We analyzed the Imatinib binding affinity to the native and mutated KIT (mutations V560G, S628N and D816V) and CSF-1R (mutation D802V) by using molecular dynamics simulations and energy calculations of Imatinib•target complexes. Further, we evaluated the sensitivity of the studied KIT receptors to Imatinib by measuring the inhibition of KIT phosphorylation. Our study showed that (i) the binding free energy of Imatinib to the targets is highly correlated with their experimentally measured sensitivity; (ii) the electrostatic interactions are a decisive factor affecting the binding energy; (iii) the most deleterious impact to the Imatinib sensitivity is promoted by D802V (CSF-1R) and D816V (KIT) mutations; (iv) the role of the juxtamembrane region, JMR, in the imatinib binding is accessory. These findings contribute to a better description of the mutation-induced effects alternating the targets sensitivity to Imatinib.

<footnote>ª Current address: Physik-Department (T38), Technische Universität München, James-Franck-Str. 1, 85748 Garching, Germany</footnote>

<footnote>* Luba.Tchertanov@ens-cachan.fr</footnote>
Introduction

Receptors tyrosine kinases (RTKs) act as primary mediators of the ligand-induced responses to
control cellular signaling. The type III RTKs, comprising the stem cell factor (SCF) receptor
KIT, the macrophage colony-stimulating factor-1 (CSF-1) receptor CSF-1R (or FMS), the
platelet-derived growth factor α and β (PDGFR-α and PDGFR-β) and the FMS-like tyrosine
kinase 3 (FLT3) receptors, are crucial for the development and physiology of different cells
under normal conditions, and are implicated in different diseases [1]. The ligand-induced
dimerization of the extracellular domain, leading to activation of the intracellular tyrosine
kinase domain (TKD), promotes a conformational switch of key regulatory elements—the acti-
vation (A-) loop, the Cα-helix and the juxtamembrane region (JMR)—from inactive to active
state (Fig 1A-1C) required for the phosphotransfer reaction [2]. Phosphorylation of specific
tyrosine residues in RTKs controls various inter- and intra-cellular signaling pathways. The
kinase activity and post-transduction processes are highly ordered and tightly regulated in nor-
mal cells [3]. Their constitutive activation promoted by mutations is associated with different
forms of cancer [3-6]. Conversely, constitutive loss of function CSF-1R mutations were
recently involved in severe neurodegenerative disorders [7]. RTKs are therefore crucial objects
for fundamental research in biology and important targets for drug development.

The high sequence similarity and structural conservation of the kinase ATP-binding pocket
—the site targeted by most kinase inhibitors—is the main challenge in the development of
selective inhibitors of tyrosine kinase activity. Nevertheless, many of such ATP competitive
inhibitors have been conceived and are currently used in the clinic or are under clinical trials.
Debates on development of either selective or multi-targets inhibitors are still ongoing. On the
one side, the use of a selective inhibitor of a particular kinase for a given tumor diminishes the
possible non-desirable side effects. On the other side, a multi-targeted inhibitor inhibits divers
signal transduction pathways that coexist in an apt tumor, thus providing a higher therapeutic
efficacy. The archetype of such multi-kinase inhibitors is Imatinib (marketed by Novartis as
Gleevec or Glivec), a rationally designed signal transduction inhibitor (antineoplastic agent),
that selectively binds to a limited number of targets (KIT, BCR/ABL, RET proto-oncogene,
CSF-1R, PDGFR-α and PDGFR-β) [8] and acts either as inhibitor (BCR/ABL, RET) or antago-
nist (KIT, CSF-1R, PDGFR-α and PDGFR-β) (DrugBank, http://www.drugbank.ca/).

Imatinib is successfully used for the treatment of several malignancies, such as Philadelphia
chromosome-positive chronic myelogenous leukaemia (CML) [9], acute lymphoblastic leukaem-
ia (ALL) [10] and KIT-positive gastrointestinal stromal tumors (GISTs) [11]. However, after
several years of clinical use, it becomes clear that Imatinib-based treatment constantly results
in acquired resistance coming from appearances of point mutations in the target kinase. Those
mutations usually occur within amino-acid sequences that encode crucial structural (and func-
tional) elements of the kinase, such as the ‘gatekeeper’ residues [12–14], the P-loop [15–17],
the JMR and the A-loop amino acids [18,19]. In fact, therapeutic inactivation of protein creates
selective pressures for tumor cells, analogous to those in natural selection, to evolve towards
resistance phenomenon, mainly through the emergence of minor preexisting resistant clones
[20,21]. The molecular mechanisms of the target resistance to Imatinib are not yet fully under-
stood. The description of such mechanisms at the atomic scale will lead to the straightforward
development of new generation of tyrosine kinase inhibitors, effective against the mutated
targets.

Despite of numerous attempts, the structural-based explaining of Imatinib selectivity by the
targets is only partial [22–25]. Extensive computational studies shed light on different aspects
of Imatinib recognition by the native targets [26–31]. In particular, molecular dynamics (MD)
simulations and free energy calculations on kinase-Imatinib complexes indicated comparable
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Structure of cytoplasmic region in RTKs. (A-C) Crystallographic structures (referenced with pdb code) of KIT in the inactive autoinhibited (1T45, left), inactive non-autoinhibited (1T46, middle) and active (1PKG, right) states are presented as cartoon. Zoomed views (bottom panel) show the active site of receptor with DFG motif (D810-F811-G812) represented in sticks. (D) Point mutations D816, V560 and S628 are showed in KIT (inactive state) as balls. In (A-D) the key structural fragments are highlighted in color—the αC-helix is in green, the activation loop (A-loop) is in red, the juxtamembrane region (JMR) is in orange. (E) Superposition of the Imatinib-binding site in structures of KIT (green) and CSF-1R (blue) in their auto-inhibited state non-bound to imatinib (PDB codes: 1T45 and 2OGV, respectively) and the structures of KIT (orange) and CSF-1R (yellow) bound to imatinib (1T46 and 4R7I, respectively). While a part of the JMR is absent in the crystallographic data, the residue W550/557 is not shown in the structures bound to imatinib.</figcaption></figure>

doi:10.1371/journal.pone.0160165.g001
binding free energies of Imatinib to different kinases—Abl, KIT, Lck and Src [26]. Computing the binding free energy of Imatinib to Abl with the hybrid method QM-MM/PBSA differentiated the forces contributing to the inhibitor binding and showed that the van der Waals energy was the main component favoring the binding of a non-protonated imatinib [32]. Another detailed computational study applying FEP/MD simulations and umbrella sampling [33] showed that the selectivity of Imatinib to Abl over Src is a consequence of distinct stabilities of the targeted conformation of the kinase domain in the two receptors, which apparently affects the conformational selection by the drug [27]. A similar conclusion was delivered for the Imatinib binding by the wild-type targets KIT and Lck [27].

The mutation effects on the energy of Imatinib binding to the kinases were computationally studied for Bcr-Abl [29,34–36], FGFR [37,38] and KIT [39]. Exploring the X-ray structure of Imatinib—KIT complex by diverse online servers lead to the prediction of the impact of KIT mutations on the stability and dynamics of the non-bound target and to the evaluation of the role of particular residues to Imatinib binding [39].

We have previously proposed mechanisms of the mutation-induced constitutive activation of two RTKs, KIT and CSF-1R, having the equivalent D816V/D802V mutation [40,41]. These two mutants are highly resistant to Imatinib [18,42]. We have shown that the structural and dynamical effects induced by the KIT oncogenic mutations V506D/G, D816V/Y/N/H and S628N (Fig 1D), correlate with the auto-activation rates of KIT mutants measured in vivo and in vitro and their Imatinib sensitivity [43,44].

To examine rigorously the mutation-induced effects on Imatinib affinity towards these RTKs, we studied the structural, dynamical and thermodynamic properties of molecular complexes formed by their clinically relevant mutants and Imatinib. Two RTKs, KIT and CSF-1R, each in the native and mutated states, with the point mutations either alternating sensitivity to Imatinib or not, were considered as targets. As in our previous studies of non-bound KIT and CSF-1R, we used all-atom MD simulations and binding free energy (MM/PBSA) calculation. We aimed to distinguish between the binding affinities of Imatinib in the different complexes and to identify the main factors driving the binding of Imatinib to the native and mutated targets. Our study demonstrated that the Imatinib binding energy to the native and mutated KIT and CSF-1R targets is quite different. The electrostatic interactions between the protonated Imatinib and the targets residues were identified as a main factor contributing to such difference, and their alternation leads either to Imatinib resistance or to increasing of the target sensitivity to this inhibitor. Further, we experimentally evaluated the relative sensitivity of the various KIT receptors to Imatinib by the inhibition of KIT phosphorylation. The two types of data, obtained *in silico* and measured *in vitro*, were found to be highly correlated, strongly suggesting that the chosen computational methods are satisfying to produce biologically interpretable results.

## Materials and Methods

### Targets selection

Structural models of the inactive state of the cytoplasmic domain of CSF-1R and KIT in their native or wild-type (WT) (CSF-1R<sup>WT</sup> and KIT<sup>WT</sup>) and mutated forms (CSF-1R<sup>D802V</sup>, KIT<sup>V560G</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup>) were selected from our previous MD simulations [40–44]. Representative conformations over the MD trajectories were extracted by a fast clustering, based on convergence analysis [45]. Briefly, a first *reference* structure is initially picked up randomly and associated with a bin of MD conformations distant by less than an arbitrary RMSD cut-off (here we chose 2.0 Å); then other *reference* structures are picked up randomly among the remaining MD conformations, in an iterative way until no MD conformations are
remaining. Then each MD conformation in the total pool is associated to its closest (RMSD) *reference* structure, to form as many *reference* clusters. Considering only the Cα atoms, we applied this technique to a concatenated trajectory of KIT, containing two MD replicas of 50 ns for each form of the protein (KIT<sup>WT</sup>, KIT<sup>V560G</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup>), and a concatenated trajectory of CSF-1R, containing two replicas of 50 ns for each form of the protein (CSF-1R<sup>WT</sup> and CSF-1R<sup>D802V</sup>). Each MD replica contained 5 000 frames, but 500 frames generated during the first 5 ns were discarded. The total number of conformations/frames are 36 000 for KIT and 18 000 for CSF-1R. Each concatenated trajectory was spliced into the ‘sections’ describing each form of the receptor. Generally, a good convergence is reached when each *reference structure* is represented by conformations from all the parts of MD trajectories. In our case, instead of taking this well-represented *reference structures*, we have chosen the ones that were more represented for a particular target, instead of been visited by all forms, WT or mutants, of CSF-1R and KIT. Only one structure for each form of KIT or CSF-1R was selected for further calculations. Since the N-terminal portion of the JMR in these selected structures is buried into the ATP-binding site provoking a clash with the inhibitor, the JMR residues 543–581 in CSF-1R and 547–588 in KIT were manually removed from each selected structure, except for KIT<sup>V560G</sup> in which the JMR was truncated at position 558 to conserve the mutated residue. Further, we have prepared the other KIT targets (KIT<sup>WT</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup>) with the JMR, truncated as in KIT<sup>V560G</sup>.

## Molecular docking

The preparation of the receptors and imatinib structures, as well as the docking runs, were performed using the Schrödinger suite Maestro (Schrödinger Release 2014). The module *protein preparation wizard* was used to re-assign hydrogens, charges and to minimize the structures of WT CSF-1R/KIT and mutants, using the default parameters. Imatinib coordinates were retrieved from the crystallographic structure of KIT in complex with this inhibitor–1T46 [25]. Imatinib was prepared using *LigPrep* (Schrödinger Release 2014) in the environment at pH 7.0. From the possible protonation states of Imatinib, generated by *LigPrep*, we have chosen the protonated state (+1), which seems to be the correct state for this inhibitor in complex with kinases, according to its predicted theoretical value [26]. The receptor structures were superposed to the crystal structure 1T46 to center the docking workspace on the ligand. The docking calculations were performed by using the *Induced Fit Docking* (IFD) [46] with the *Extended Sampling* protocol in which all targets residues positioned at a distance inferior to 5 Å from Imatinib were considered flexible. The binding affinity of Imatinib to the targets was characterized by docking scores (GLIDE and IFD). Each docking model was superposed with the crystal structure of the Imatinib-bound KIT (1T46) using *Maestro* from the Schrödinger suite, and the RMSD values were calculated in respect to 1T46. The best generated models of Imatinib•target complexes (target is CSF-1R<sup>WT</sup>, CSF-1R<sup>D802V</sup>, KIT<sup>WT</sup>, KIT<sup>V560G</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup> with the entirely truncated JMR), reconstituting the Imatinib position in 1T46, were used for MD simulations. Imatinib was manually placed into KIT<sup>WT</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup> structures bearing the partially truncated JMR. This was easily done by superposition with the successfully docked complexes for their truncated forms. Possible steric clashes were eliminated by energy minimization of the complexes.

## Molecular dynamics simulations

**Imatinib preparation.** To generate the Imatinib topology parameters compatible with CHARMM all-atoms force field [47], an Imatinib structure was retrieved from the Zinc Database (http://zinc.docking.org/) and used as input for the web server Swissparam [48].
**Set up of the systems.** The set-up of MD simulations was performed with the CHARMM27 all-atom force field integrated in GROMACS package 4.6.5 [49]. For simulations, each model was placed in a cubic box and centred with a 1.2 nm distance to the faces under periodic boundary conditions and solvated with explicit TIP3P water molecules [50]. Cl⁻ counter ions were added when necessary to neutralize the overall charge. The minimization procedure consisted of two steps: (i) steepest descent energy minimization (EM) with the solute atoms restrained and (ii) EM with all atoms free. The equilibration procedure was performed on the solvent, keeping the solute atoms (except H-atoms) restrained for 500 ps at 310 K and a constant volume (canonical NVT ensemble).

**Production of MD trajectories.** Two runs of 50-ns MD simulations were carried out for each complex ‘Imatinib•target’ (target is CSF-1R<sup>WT</sup>, CSF-1R<sup>D802V</sup>, KIT<sup>WT</sup>, KIT<sup>V560G</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup> with the entirely and partially cleaved JMR) and for each cleaved wild-type target in absence of Imatinib, further called ‘non-bound’ or ‘apo’ form. The temperature of solute (Imatinib•target) and solvent (water and ions) was separately coupled to the velocity rescale thermostat [51] at 310 K with relaxation time of 0.1 ps. The pressure was maintained at 1 atm by isotropic coordinate scaling (NPT ensemble) with relaxation time of 1 ps using Berendsen thermostat [52]. A time step of 2 fs was used to integrate the equations of motion based on the Leap-Frog algorithm [53]. The Lennard-Jones interactions were shifted to a cut-off 1.4 nm, and the Particle Mesh Ewald (PME) method [54] was used to treat long-range electrostatic interactions. The list for the electrostatic interactions was updated every 5 steps. All bonds were controlled using the Linear Constraint Solver (P-LINCS) algorithm [55]. The SETTLE algorithm [56] was used to constrain the water molecules geometry. Coordinates files were recorded every 10 ps.

**Analysis of the trajectories.** The MD trajectories were analysed (RMSDs and RMSFs computation and Principle Component Analysis) with tools included in the GROMACS package. The first 5 ns of each trajectory (equilibration time), were removed prior to analysis. A convergence analysis was performed using an ensemble-based approach described in details in [40–41,44]. Each analysis was performed on 9 000 conformations generated over two MD runs of 45 ns per simulation. Visual inspection of the MD data was made with PyMOL [57] and VMD [58]. Graphs were generated using Grace (http://plasma-gate.weizmann.ac.il/Grace/).

## Free energy calculation

Standard MM-PBSA (Molecular Mechanics—Poisson Boltzmann Surface Area) method [59] has been used for calculation of thermodynamic parameters and free energy of binding. In MM-PBSA, the free energy difference of binding (ΔG<sub>bind</sub>) between a ligand (L is Imatinib) and a receptors (R) to form a complex RL is defined as:

$$ \Delta G_{bind} = \Delta H - T\Delta S \approx \Delta E_{MM} + \Delta G_{sol} - T\Delta S \quad (1) $$

$$ \Delta E_{MM} = \Delta E_{int} + \Delta E_{elect} + \Delta E_{vdw} \quad (2) $$

where ΔE<sub>MM</sub>, ΔG<sub>sol</sub> and –TΔS are the changes of the gas phase MM energy, the solvation free energy and the conformational entropy upon binding, respectively.

ΔE<sub>MM</sub> includes ΔE<sub>int</sub> (covalent bonding including the bond, angle and dihedral energies), electrostatic ΔE<sub>elect</sub> and van der Waals ΔE<sub>vdw</sub> energies. ΔG<sub>sol</sub> is the sum of the electrostatic solvation energy (polar contribution) and non-electrostatic solvation component (non-polar contribution). The polar contribution was calculated using PB model, while the non-polar free energy estimated by solvent accessible surface area (SASA).

ΔG<sub>bind</sub> was computed for each generated Imatinib•target complex over the merged trajectories containing 9 000 frames (and also for the individual MD simulation replicas over a
trajectory containing 4500 frames). We used the g_mmpbsa module of GROMACS, which combines subroutines from GROMACS and APBS [60] to calculate the enthalpy contribution (ΔH) to the Gibbs free energy and the per-residue energy decomposition. The energy components ΔE<sub>MM</sub> and the polar and non-polar components of ΔG<sub>sol</sub> were calculated in the bound and unbound form, and subsequently their contribution to the binding energy ΔR<sub>x</sub><sup>BE</sup> of residue x [61] was calculated as follows:

$$ \Delta R_x^{BE} = \sum_{i=0}^{n} (A_i^{bound} - A_i^{free}) \quad (3) $$

Where A<sub>i</sub><sup>bound</sup> and A<sub>i</sub><sup>free</sup> are the energy of ith atom from x residue in bound and free forms, respectively, and n is the total number of atoms in the residue.

The g_mmpbsa module does not include the calculation of entropic terms and therefore it gives the relative binding energy. Since we have used the same compound (Imatinib) interacting with proteins of similar structures and showing the same binding mode, the entropy contribution was neglected [62].

## Data simulation analysis

**Electrostatic potential surface.** The electrostatic potential surface was computed with the APBS software running through the PDB2PQR webserver (http://nbcr-222.ucsd.edu/pdb2pqr/). As the input data, the equilibrated conformations of the Imatinib-target complexes prior to MD simulations, were used.

**Salt bridges.** The salt bridges occurrences were calculated for the merged trajectories containing 9 000 frames using the Timeline plugin available with VMD. The residues-pair selection was made within a cutoff distance of 10 Å from Imatinib. Plots were generated using R software [63], discarding interactions occurring less than 10% of the simulation time.

## Inhibition of KIT phosphorylation

**Cell transfection.** COS-7 cells were transfected with 1 µg of KIT<sup>WT</sup>, KIT<sup>D816V</sup>, KIT<sup>V560G</sup> or KIT<sup>S628N</sup> expression vectors with FuGENE 6 (Promega). Twenty four hours following transfection, cells were split in four 60-mm plates and serum-starved overnight in Dulbecco's modified Eagle's medium with 0.5% FBS. Cells were then treated for 90 minutes with Imatinib at concentrations ranging from concentration 0 to 5µM. Cells transfected with hKIT WT were stimulated for 5 min with 250 ng/ml murine SCF at 37°C.

**Immunoblotting.** Cells were washed in ice-cold phosphate-buffered saline (PBS), pelleted, and lysed in HNTG buffer (50mM HEPES [pH 7], 50mM NaF, 1mM EGTA, 150mM NaCl, 1% Triton X-100, 10% glycerol, 1.5 mM MgCl2) containing the protease inhibitor mix Complete EDTA-free (Roche Applied Science) and 100 µM Na<sub>3</sub>VO<sub>4</sub>. Proteins were separated on SDS-polyacrylamide gels and transferred to Immobilon-P membrane (Millipore). Membranes were saturated with 5% bovine serum albumin (Euromedex) and probed with phospho-c-Kit (Tyr719) (Cell Signaling Technology) and c-Kit (D13A2) antibodies (Cell Signaling Technology). Signals were revealed using horseradish peroxidase-conjugated secondary antibodies and SuperSignal West Pico Chemiluminescent Substrate (ThermoScientific).

## Results and Discussions

### Imatinib docking into the targets

For the study of stability of the Imatinib•target complexes and the comparison of their free binding energy, we used the representative conformations of the inactive form of native (or wild-type, WT) and mutated kinase domain of two proteins, CSF-1R (CSF-1R<sup>WT</sup> and
CSF-1R<sup>D802V</sup>) and KIT (KIT<sup>WT</sup>, KIT<sup>D816V</sup>, KIT<sup>V560G</sup> and KIT<sup>S628N</sup>), selected from our previous MD simulations [40–41, 43–44]. Our earlier attempts to dock Imatinib into the inactive form of the kinase domain of KIT and CSF-1R have failed [64], probably due to the buried configuration of the JMR. To clarify the reason of this failed assay, we compared the crystallographic structures of KIT and CSF-1R bound and non-bound to Imatinib-1T46 [25], 4R7I [65], 2OGV [66] and 1T45 [25], respectively. We observed that the JMR buried configuration and the side-chain orientation of residues W550/557 and F798/811 (CSF-1R/KIT) in both receptors, impairs the binding of Imatinib into the ATP-binding site (Fig 1E). Repositioning of the side chain of F811/797 in KIT/CSF-1R bound to Imatinib, evidenced by crystallography (1T46/4R7I), confirms its inhibitor-induced (or inhibitor-selected) relocation. The Imatinib binding effect on the JMR residue W550/557 (CSF-1R/KIT) inserted in the binding pocket remains obscure, because the N-extremity of the JMR is not resolved in the crystallographic structures 1T46 and 4R7I, suggesting that this portion might be very flexible and dislocated from the kinase domain. To avoid steric hindrance, we removed the JMR from the selected MD conformations, producing suitable targets for the docking procedure. In the mutant KIT<sup>V560G</sup>, to conserve the point mutation, the JMR was truncated at position 558, yielding a target, longer by 30 residues (V558-R588) than the other targets. Expecting a tiny impact of this JMR partial fragment on Imatinib binding, we have accepted such exception to compare the Imatinib binding to all studied targets.

The optimised structure of the protonated Imatinib [26], generated with Ligprep, was successfully docked into each prepared target (see Methods). The best docking poses of Imatinib show a striking similarity in all studied proteins (Fig 2D, insert) and all are very close to that observed in the crystallographic structure of Imatinib•KIT (1T46) and Imatinib•CSF-1R (4R7I) complexes, in which the position of Imatinib and its conformation is virtually identical (Fig A in S1 File). Taking 1T46 as a reference, the docking poses, showing Root Mean Square Deviation (RMSD) values less than 1 Å and reproducing five H-bonds described in the crystallographic structure, were preselected. The best poses, ranked by RMSDs and GLIDE scores, were retained for further analysis. The GLIDE scores of the best poses, ranging from -9.2 to -6.3 kcal/mol over the generated complexes, showed a common tendency to be diminished (higher affinity) in the native proteins and KIT<sup>V560G</sup>, and worse in the resistant mutants (Table 1). Although GLIDE docking scores show only qualitative correlation with the experimentally determined Imatinib affinity, their values help the choice of ‘Imatinib•target’ complexes further studied by MD simulations. The IFD scores, distributed in a narrower range, from -12.78 to -10.96 kcal/mol, do not significantly distinguish the binding affinity of Imatinib to the different targets.

## Molecular dynamic simulations of CSF-1R and KIT complexes

To get further insight into the molecular mechanisms of the resistance to Imatinib, it is of primary importance first to investigate the structural and dynamical properties of the Imatinib-bound targets and second, to define the energy-related parameters driving the molecular complex formation and stability. For these purposes, we performed MD simulations on Imatinib-bound receptors, KIT and CSF-1R, in the native and mutated forms.

**Stability of the target•Imatinib complexes.** The RMSDs calculated on the backbone atoms from the initial structures of CSF-1R and KIT in their free and bound states are relatively stable over MD simulations, except the very flexible residues located at the C-terminal loop of both receptors, which are significantly contributing to the RMSD values increase. These residues (7/6 residues from the C-terminal of CSF-1R/KIT) were excluded from the RMSD analysis. Overall, the RMSD values in almost all simulated complexes vary around 0.15–0.20 nm (Fig 2 and Fig B in S1 File), while one of the two MD simulations replica of each receptor
<figure><img src="image_5.png" /><figcaption>**Fig 2.** MD simulations of the RTK-Imatinib complexes. (A-B) The Root Mean Square Deviation (RMSD) values of the native non-bonded receptors (CSF-1R<sup>WT-APO</sup> and KIT<sup>WT-APO</sup>) and the complexes formed by Imatinib and the native receptors (CSF-1R<sup>WT</sup> and KIT<sup>WT</sup>) and their mutants (CSF-1R<sup>D802V</sup>, KIT<sup>D816V</sup>, KIT<sup>V560G</sup> and KIT<sup>S628N</sup>) were calculated for the backbone atoms from trajectories 1 and 2 (orange and blue respectively) of MD simulations of all studied systems. (C) Simulated complexes were obtained by the docking of Imatinib (in green) to the active site of each receptor. The JMR was removed prior the docking procedure. The superimposed receptors CSF-1R/KIT are shown as grey/cyan cartoons. The key structural fragments of receptors, the αC-helix and the activation loop (A-loop) are highlighted in cyan/brown (CSF-1R/KIT) and in red/rose (CSF-1R/KIT) color respectively. **Insert in C:** The best docking poses of Imatinib bound to CSF-1R<sup>WT</sup>(blue), CSF-1R<sup>D802V</sup>(orange), KIT<sup>WT</sup>(cyan), KIT<sup>V560G</sup>(green), KIT<sup>S628N</sup>(pink) and KIT<sup>D816V</sup>(magenta) superposed with the crystallographic structure 1T46 (grey). The protein backbone and Imatinib are represented in cartoon and sticks, respectively.</figcaption></figure>

doi:10.1371/journal.pone.0160165.g002

possessing the equivalent mutation, D802V (CSF-1R) or D816V (KIT), shows slightly increased RMSD values (0.25 and 0.35nm, respectively).

To identify the residues responsible for this increased structural instability, we computed the Root Mean Square Fluctuations (RMSFs) of all amino acids in each protein. The MD
<table><caption>**Table 1.** Imatinib docking into the native targets (CSF-1R<sup>WT</sup> and KIT<sup>WT</sup>) and their mutants (CSF-1R<sup>D802V</sup>, KIT<sup>V560G</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup>). The two energy values, the GLIDE and Induced Fit Docking (IFD) scores, were produced with MAESTRO (Schrödinger release, 2014). The RMSD values were calculated for the predicted docking poses taking crystallographic structure 1T46 as a reference. The IC<sub>50</sub> values (µM) are taken from the literature.</caption><thead><tr><th>Target/Parameter</th><th>Score Glide (kcal/mol)</th><th>Score IFD (kcal/mol)</th><th>RMSD (Å)</th><th>IC<sub>50</sub> Reference (µM)</th></tr></thead><tbody><tr><td>CSF-1R<sup>WT</sup></td><td>-7.9</td><td>-10.9</td><td>0.5</td><td>0.3 [42]</td></tr><tr><td>CSF-1R<sup>D802V</sup></td><td>-6.3</td><td>-12.4</td><td>0.5</td><td>>4 [42]</td></tr><tr><td>KIT<sup>WT</sup></td><td>-7.8</td><td>-12.2</td><td>0.5</td><td>0.1-1 [18,19]</td></tr><tr><td>KIT<sup>V560G</sup></td><td>-9.2</td><td>-12.1</td><td>0.6</td><td>0.025 [18]</td></tr><tr><td>KIT<sup>S628N</sup></td><td>-8.6</td><td>-12.3</td><td>1.0</td><td>not determined</td></tr><tr><td>KIT<sup>D816V</sup></td><td>-7.0</td><td>-12.1</td><td>0.4</td><td>>10 [19]</td></tr></tbody></table>

doi:10.1371/journal.pone.0160165.t001

RMSFs, usually coupled with the 3D interpretation of the temperature factors (beta-values) from crystallographic data, show that the highest atomic fluctuations are principally observed, as was expected, in some particular loop regions, as we established previously in the non-bonded targets [41,43,44]. Comparing the native non-bound targets, CSF-1R and KIT, we note some differences in the atomic fluctuations of residues from the molecular fragments involved in the activation process—the catalytic (C-) loop, the activation (A-) loop and the loop preceding αC-helix, annotated on Fig 3 as 6, 7 and 3, respectively. When analysing the Imatinib-bound complexes, we remarked that binding of Imatinib influences differently the atomic fluctuations of the native and mutated targets.

In both Imatinib•CSF-1R complexes, formed by the native and mutated receptors, we observed an alteration of atomic fluctuations respectively to those in the non-bound target—higher for residues of the loop preceding the αC-helix (3 in Fig 3), and diminished for residues of the 3<sup>10</sup>-helix of A-loop (7 in Fig 3). The RMSFs of the N-lobe loops residues (2 and 4 in Fig 3), showing similar values in the non-bound target (CSF-1R<sup>WT/APO</sup>) and in Imatinib•CSF-1R<sup>WT</sup> complex, are significantly diminished in the complex formed by the mutated target CSF-1R<sup>D802V</sup>. The P-loop residues (annotated as 1 in Fig 3) show diminished RMSFs in complexes compared to non-bound receptor.

In KIT complexes, the RMSFs pattern generally corresponds to that in CSF-1R's (Fig 3). Nevertheless, some features, characterizing the atomic fluctuations in KIT, mutated or not, bound or not to Imatinib, were different from those in CSF-1R's. In particular, in the non-bound native KIT, the RMSFs of residues from the loop (3) preceding the C-helix are increased respective to those in CSF-1R, whereas they are significantly diminished in the A-loop 3<sup>10</sup>-helix (7) and in the catalytic (C-) loop (6). When KIT is bonded to Imatinib, the RMSF values in the N-lobe's loop preceding the C-helix are diminished in the native target—the effect opposite to those observed in CSF-1R-, while in complexes formed by the mutants theirs values are conserved. For the Imatinib-bound KIT<sup>V560G</sup> and KIT<sup>S628N</sup> mutants, specifically, the increasing of RMSFs is observed for residues from the A-loop (annotated as 7' in Fig 3) and the C-lobe loop preceding G-helix (9 in Fig 3). One of the differences between the two receptors, KIT and CSF-1R, is the atomic fluctuations of the catalytic (C-) loop residues. These values are diminished in both CSF-1R complexes compared to the non-bound target—significantly in Imatinib•CSF-1R<sup>WT</sup> and moderately in Imatinib•CSF-1R<sup>D802V</sup>. Apparently, the increase of the C-loop atomic fluctuations in Imatinib•CSF-1R<sup>D802V</sup> complex is directly related with the alteration in the H-bonds pattern involving residues of the C- and A-loops (Table 2). In all KIT targets, wild-type and mutated, bound or not, the RMSFs of the catalytic (C-) loop are similarly very small.

**The H-bonds pattern stabilizing the Imatinib binding in the complexes.** The five H-bonds observed in the crystallographic structures (1T46 and 4R7I for KIT and CSF-1R,
respectively) and described as the principal non-covalent contacts stabilising the binding of Imatinib to the native targets (KIT<sup>WT</sup> and CSF-1R<sup>WT</sup>), are well preserved along MD trajectories (Fig 4A and 4C, scheme on the left). These H-bonds are formed by functionally significant amino acids in KIT<sup>WT</sup>/CSF-1R<sup>WT</sup>—the gatekeeper residue T670/663, the highly conserved

<figure><img src="image_6.png" /><figcaption>**Fig 3.** MD simulations of the RTK-Imatinib complexes. (A) The Root Mean Square Fluctuations (RMSFs) computed on the backbone atoms averaged over the total production simulation time of CSF-1R<sup>WT</sup> and CSF-1R<sup>D802V</sup> (blue and red, respectively) (top panel) and of KIT<sup>WT</sup>, KIT<sup>D816V</sup>, KIT<sup>V560G</sup> and KIT<sup>S628N</sup> (blue, red, yellow and green, respectively) complexes (bottom panel) were compared to those in the non-bound receptors (dotted black lines). (B) The average conformations of CSF-1R (top panel) and KIT (bottom panel) are presented as tubes. The tube size is proportional to the residue atomic fluctuations computed on the backbone atoms. The highly fluctuating regions are colored in red and numerated for CSF-1R in (A) and (B).</figcaption></figure>

doi:10.1371/journal.pone.0160165.g003
<table><caption>**Table 2.** H-bonds between the A- and C-loops residues in CSF-1R and KIT. Occurrences (in %) of the H-bonds are averaged over the two MD replicas.</caption><thead><tr><th>Contact/Target</th><th>CSF-1R<sup>WT</sup></th><th>Im•CSF-1R<sup>WT</sup></th><th>Im•CSF-1R<sup>D802V</sup></th></tr></thead><tbody><tr><td>D778...N783</td><td>92</td><td>100</td><td>70</td></tr><tr><td>D778...H776</td><td>69</td><td>60</td><td>67</td></tr><tr><td>D778...R801</td><td>93</td><td>78</td><td>1</td></tr><tr><td>D778...Y809</td><td>18</td><td>39</td><td>78</td></tr><tr><td>D778...R782</td><td>2</td><td>7</td><td>94</td></tr></tbody><thead><tr><th>Contact/Target</th><th>KIT<sup>WT</sup></th><th>Im•KIT<sup>WT</sup></th><th>Im•KIT<sup>D816V</sup></th></tr></thead><tbody><tr><td>D792...N797</td><td>50</td><td>51</td><td>89</td></tr><tr><td>D792...H790</td><td>23</td><td>34</td><td>43</td></tr><tr><td>D792...R815</td><td>100</td><td>100</td><td>74</td></tr><tr><td>D792...Y823</td><td>47</td><td>45</td><td>45</td></tr></tbody></table>

doi:10.1371/journal.pone.0160165.t002

residue E640/633, the DFG residue D810/796, the hinge residue C673/666 and the C-loop residue I789/775. Over MD simulations of Imatinib•KIT<sup>WT</sup>/CSF-1R<sup>WT</sup> complexes, these H-bonds are characterized by high, moderate or low occurrence values—98% (C673/666), 47/38% (T670/663), 42/37% (D810/796), 12/34% (E640/633) and 96/98% (I789/775), respectively (Table 3, Scheme A in S1 File).

This analysis shows that in KIT and CSF-1R complexes formed by the native targets, the H-bonds of Imatinib with the hinge residue C673/666 and with the hydrophobic residue I789/775 are constantly and similarly maintained during the MD simulations with an occurrence close to 100%. The other H-bonds, characterised by relatively lower occurrences, show some difference between the two targets. In particular, the H-bond of Imatinib with D810 in KIT shows higher occurrence (42%) than compared to E640 (12%), while the corresponding H-bonds in CSF-1R have similar occurrence. The H-bonds involving the residue T670/663 have close occurrences in both native targets.

As it was evidenced by the analysis of different X-ray structures, all these H-bonds stabilizing Imatinib in the binding site are observed in complexes formed by KIT, CSF-1R and other tyrosine kinases, while the author's interpretation of the observed contacts may differ. For instance, in the crystal structure of KIT complexed with Imatinib (1T46) [25], the Imatinib interactions with I789 through the N-methylpiperazine fragment was described as non-specific. In the recently solved structure of Imatinib in complex with CSF-1R (4R7I) [3], this contact is considered as alternative to that with H776. In other structures of kinases complexed with Imatinib, such as Abl (2HYY) [67], c-Src (2OIQ) [22] and DDR1 (4BKJ) [68], the N-methylpiperazine fragment of Imatinib makes direct contacts with a hydrophobic residue (isoleucine or valine) and with histidine. Such interpretation of the N-methylpiperazine moiety's stabilising contacts is probably related to considering the Imatinib as a neutral molecule over the structure refinement in the X-ray analysis. Therefore, our interpretation of Imatinib contact with I789 as a strong and stable H-bond interaction should be considered for the improving or the development of kinase inhibitors.

In CSF-1R<sup>WT</sup> and CSF-1R mutant, the H-bond interactions stabilizing Imatinib in its binding site are overall conserved; however, their occurrence values vary significantly among the different mutants (Table 3). In particular, an essential decrease of the H-bond occurrence is detected for the interaction of Imatinib with I775 in CSF-1R<sup>D802V</sup> (18% against 96% for CSF-1R<sup>WT</sup>. In KIT<sup>D816V</sup> the occurrences for the same H-bonds were surprisingly maintained. I775/789 (CSF-1R/KIT) is placed in the C-loop, and the binding of this residue with Imatinib correlates directly with the C-loop RMSFs observed in the related proteins.
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Imatinib binding to CSF-1R and KIT. (A) 3D structure of Imatinib•target complex illustrated by using KIT<sup>WT</sup> bound with Imatinib. Imatinib and interacting residues of receptor are shown by sticks. (B) Graphical representation of the free binding energy (ΔG) of Imatinib•target complexes. The computed total ΔG energy value (this work) and the experimentally measured affinity (the literature data and the present data) are shown for each complex. (C) The H-bonds (dashed lines) pattern stabilizing Imatinib in complexes formed by the native KIT and its mutants KIT<sup>D816V</sup> and KIT<sup>V560G</sup> (left), and in mutant KIT<sup>S628N</sup> (right).</figcaption></figure>

doi:10.1371/journal.pone.0160165.g004
<table><caption>**Table 3. Imatinib binding to the targets.** Occurrences of the H-bonds between Imatinib and the targets residues were averaged over each of the two MD replicas. Imatinib atoms participating in the interactions with the targets are represented in Scheme A in S1 File. The atom pairs for donor-acceptor interactions are depicted in the Table.</caption><thead><tr><th>Target/Atoms</th><th>N3-(N)C666</th><th>O1-(N)D796</th><th>N7-(Oδ)D796</th><th>N7-(O)I775</th><th>N5-(Oε)E633</th><th>N4-(Oγ)T663</th></tr></thead><tbody><tr><td>CSF-1R<sup>WT</sup></td><td>98</td><td>37</td><td>0</td><td>96</td><td>34</td><td>38</td></tr><tr><td>CSF-1R<sup>D802V</sup></td><td>99</td><td>26</td><td>0</td><td>18</td><td>34</td><td>51</td></tr><tr><td></td><td>N3-(N)C673</td><td>O1-(N)D810</td><td>N7-(Oδ)D810</td><td>N7-(O)I789</td><td>N5-(Oε)E640</td><td>N4-(Oγ)T670</td></tr><tr><td>KIT<sup>WT</sup></td><td>98</td><td>42</td><td>0</td><td>98</td><td>12</td><td>47</td></tr><tr><td>KIT<sup>V560G</sup></td><td>98</td><td>73</td><td>0</td><td>95</td><td>46</td><td>34</td></tr><tr><td>KIT<sup>S628N</sup></td><td>93</td><td>68</td><td>56</td><td>33</td><td>34</td><td>27</td></tr><tr><td>KIT<sup>D816V</sup></td><td>98</td><td>60</td><td>0</td><td>99</td><td>27</td><td>51</td></tr></tbody></table>

doi:10.1371/journal.pone.0160165.t003

As expected, given its higher sensitivity to Imatinib compared to KIT<sup>WT</sup>, KIT<sup>V560G</sup> binds Imatinib through H-bonds characterised by higher occurrences than in KIT<sup>WT</sup>, particularly for interactions with D810 (73% vs 42%) and D640 (46% vs 12%). Some change in occurrences of the H-bond interactions was also evidenced in the Imatinib•KIT<sup>S628N</sup> complex, formed by a mutant recently found in a metastatic melanoma and described as still sensitive to Imatinib [43]. In this complex, mutation affects mainly the inhibitor interactions with residues I789 and D810, which are reduced for I789 (33% vs 98%) and increased for D810 (68% vs 42%). Furthermore, Imatinib interacts with residue D810 by H-bond with the δ-oxygen atom in KIT<sup>S628N</sup>, a contact that was not observed in complexes formed by the other targets (Table 3, Scheme A in S1 File). Visual inspection of the Imatinib•KIT<sup>S628N</sup> MD conformations showed that this H-bond occurs due to a flip of the N-methylpiperazin fragment of Imatinib and a change on the side chain orientation of D810, favouring the interaction between the inhibitor and the D810 δ-oxygen atom (Fig 4, scheme on the right). Such conformational adaptation of both interacting partners, target and inhibitor, induces a significant diminishing or entire loss of the inhibitor interaction with I789 (33% KIT<sup>S628N</sup> vs 98% in KIT<sup>WT</sup>) but stabilizes the bifurcate H-bond with D810 side chain. For each complex, the H-bonds patterns stabilizing the Imatinib-target binding and the H-bonds occurrences were similar when either the individual or concatenated MD trajectories were analyzed (Table 3, Table A in S1 File).

## The free binding energy in Imatinib•target complexes and its relation with targets sensitivity to Imatinib

The conformational change in the ATP-binding site is a footprint of the induced-fit mechanism according of which the protein adapts its active site conformation to accommodate the ATP or inhibitor. Previously, based on normal modes analysis, it was concluded that KIT mutants (including D816V) had a decrease in flexibility of the majority of residues that participated in the binding of Imatinib [39]. In our MD simulations, the H-bond pattern found for KIT<sup>D816V</sup> indicates an increased rigidity of the ATP-binding site. Due to the reduced flexibility, residues forming the ATP-binding site in KIT<sup>D816V</sup>, would be less able to adapt their conformation over the recognition process, possibly reducing the affinity for the inhibitor. In contrast, the conformational adaptation of mutated target (in Imatinib•KIT<sup>V560G</sup>) or together with the inhibitor (in Imatinib•KIT<sup>S628N</sup>) may be sufficient to reach a good affinity between two partners.

**Experimental evaluation of KIT sensitivity to Imatinib.** To clarify the role of the S628N substitution in KIT, we determined experimentally the relative sensitivity of the various KIT (KIT<sup>WT</sup>, KIT<sup>D816V</sup> KIT<sup>V560G</sup> and KIT<sup>S628N</sup>) receptors. COS7 cells were transfected with the corresponding expression vectors, and were then treated with various concentrations of
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Inhibition of KIT phosphorylation by Imatinib. Cells expressing KIT<sup>WT</sup>, KIT<sup>D816V</sup>, KIT<sup>V560G</sup> or KIT<sup>S628N</sup> were treated with the indicated concentrations of Imatinib for 90 minutes. Protein lysates were analyzed by western blotting to reveal KIT protein expression (KIT) and KIT phosphorylation (Phospho-KIT) as readout of Imatinib inhibition. Two independent experiments were performed.</figcaption></figure>

doi:10.1371/journal.pone.0160165.g005

Imatinib ranging from 0.01 to 5 µM. The inhibition of KIT phosphorylation was then evaluated by western blotting (Fig 5). As previously described, KIT<sup>D816V</sup> was not affected even at 5 µM, while KIT<sup>WT</sup> was partially inhibited at 0.1 and fully inhibited at 1 µM. As expected also, KIT<sup>V560G</sup> was inhibited at 0.1µM. Finally, KIT<sup>S628N</sup> had an inhibition comparable to KIT<sup>WT</sup>. Consequently, the experimental data indicates the following order of KIT sensitivity to Imatinib: KIT<sup>V560G</sup> > KIT<sup>WT</sup> = KIT<sup>S628N</sup> > KIT<sup>D816V</sup>.

**The binding free energy of Imatinib•target complexes.** Since the H-bonds pattern does not explain exhaustively a difference of the Imatinib binding to the sensitive and resistant targets, we further analysed the binding free energy, (ΔΔG), [69] of Imatinib to each target using the MM/PBSA approach [59]. The binding free energy values are considerably different in the studied complexes (Fig 4B). The lowest (-94 kcal/mol) and the biggest (-42/-43 kcal/mol) values are found in complexes formed by the most sensitive (KIT<sup>V560G</sup>) and the most resistant (CSF-1R<sup>D802V</sup> and KIT<sup>D816V</sup>) targets to Imatinib, respectively. Based on these binding free energy values, the *in silico* estimated sensitivity of the studied targets to Imatinib may be described by the following order—KIT<sup>V560G</sup> > CSF-1R<sup>WT</sup>/KIT<sup>WT</sup> ≥ KIT<sup>S628N</sup> ≥ CSF-1R<sup>D802V</sup>/KIT<sup>D816V</sup> and CSF-1R<sup>WT</sup>/KIT<sup>WT</sup> > CSF-1R<sup>D802V</sup>/KIT<sup>D816V</sup>, which is qualitatively very consistent with our *in vitro* data.

In general, the experimentally measured sensitivity of the different KIT receptors to Imatinib shows the same tendency as predicted by the free energy of binding calculations. Only the KIT<sup>S628N</sup> energy profile could not be significantly distinguished *in silico* from those of KIT<sup>WT</sup> and KIT<sup>D816V</sup>. Such divergence between the theoretical prediction and the experimentally measured values may be attributed to the limitation of the MM/PBSA method, in comparison to more accurate approaches such as alchemical free energy perturbation molecular dynamics (FEP/MD) [70].

However, high absolute quantitative accuracy was not in the scope of the present study that focuses on the discrimination of Imatinib binding by the native and mutated targets (RTKs). To our knowledge, this demonstration of a difference in the energy of binding of Imatinib to the native and resistant forms of these two kinase targets, evidenced computationally, is new. The binding energy alteration in Imatinib-targets complexes together with the earlier described mutation-induced perturbation of the Imatinib-specific inactive conformation as an important step for transition toward the constitutively active form [41,44], fully explain the resistant profile of the studied CSF-1R<sup>D802V</sup>/KIT<sup>D816V</sup> mutants.
<table><caption>**Table 4. The free energy of binding in the Imatinib-targets complexes.** The free binding energy (ΔG<sub>bind</sub>) and contributions of electrostatic (ΔG<sub>elect</sub>), van der Waals (ΔG<sub>vdw</sub>) and solvation (ΔG solv) energies were calculated using MM-PBSA. All energies are in kcal/mol.</caption><thead><tr><th>RTK</th><th>ΔG<sub>vdw</sub></th><th>ΔG<sub>elect</sub></th><th>ΔG<sub>sol</sub></th><th>ΔG<sub>bind</sub></th></tr></thead><tbody><tr><td>CSF-1R<sup>WT</sup></td><td>-66</td><td>-76</td><td>73</td><td>-69</td></tr><tr><td>CSF-1R<sup>D802V</sup></td><td>-66</td><td>-54</td><td>78</td><td>-42</td></tr><tr><td>KIT<sup>WT</sup></td><td>-67</td><td>-75</td><td>81</td><td>-61</td></tr><tr><td>KIT<sup>D816V</sup></td><td>-66</td><td>-56</td><td>80</td><td>-43</td></tr><tr><td>KIT<sup>V560G</sup></td><td>-64</td><td>-106</td><td>76</td><td>-94</td></tr><tr><td>KIT<sup>S628N</sup></td><td>-67</td><td>-67</td><td>83</td><td>-51</td></tr></tbody></table>

doi:10.1371/journal.pone.0160165.t004

Since Imatinib is held in the kinase pocket via an extended network of hydrogen bonding and charge-charge interactions, the decomposition of the binding free energy on the different terms contributing to the binding affinity of Imatinib helps to identify the primary forces driving such binding and to analyze their variations among different complexes. The electrostatic energy (ΔG<sub>elect</sub>) shows more variation (from -106 kcal/mol in KIT<sup>V560G</sup> to -54 kcal/mol in CSF-1R<sup>D802V</sup>) than the others terms and appears to be the principal force contributing to the difference in the final binding energy values (ΔG<sub>bind</sub>) in the studied complexes (Table 4, Table B in S1 File). The correlation coefficient between ΔG<sub>elect</sub> and ΔG<sub>bind</sub> is 0.89. The solvation energy (ΔG<sub>sol</sub>) and the van der Waals dispersion energy (ΔG<sub>vdw</sub>) show only minor variations among the studied complexes (ΔΔG<sub>sol</sub> and ΔΔG<sub>vdw</sub> is less than 10 and 3 kcal/mol, respectively). The significant role of the electrostatic component in the Imatinib binding affinity to KIT was reported earlier [27]. Recently, the analysis of the reported literature binding energies of Imatinib to tyrosine kinases, obtained by molecular mechanics (MM) and those found by quantum mechanics (QM), suggested that solvation energies are a major component of the overall binding energy [71].

In our quite simple analysis, we have shown that KIT<sup>V560G</sup> has the most favorable protein-Imatinib electrostatic interactions, greater than the native targets, KIT and CSF-1R, suggesting that mutation V560G promotes the largest electrostatic stabilization of the Imatinib binding. On the contrary, mutation D802V in CSF-1R and D816V in KIT, diminishes the electrostatic interactions with Imatinib.

**The per residue energy contributions to the binding energy.** To identify the key molecular determinants responsible for Imatinib binding to the targets, we calculated per residue energy contribution of all protein residues to the overall binding energy. As was expected, the free energy binding value is mostly influenced by the charged residues,–either favorably (impact of the negatively charged amino acids, asparagine and glutamic acid)–or unfavorably (impact of the positively charged residues, arginine and lysine) demonstrating nearly similar contributions from the corresponding residues in CSF-1R and KIT (Tables C and D in S1 File). The substitution of an aspartate for a valine (D→V) in CSF-1R<sup>D802V</sup> and KIT<sup>D816V</sup> reduces significantly the residue energy contribution—from -8 kcal/mol (CSF-1R<sup>WT</sup>/KIT<sup>WT</sup>) to 0 kcal/mol (CSF-1R<sup>D802V</sup>/KIT<sup>D816V</sup>). In KIT<sup>V560G</sup>, the residue V560G, distant by ~8 Å from the Imatinib, does not contribute considerably (-0.3 kcal/mol) to the binding energy. Similarly, in KIT<sup>S628N</sup>, the native or mutated residue affects lowly the binding energy. This estimation also evidenced that the equivalent Imatinib-binding site residues (in CSF-1R and KIT), making direct contact with the inhibitor, contribute almost similarly to the binding energy and show only small differences between the native and mutated targets.

**The JMR role in the Imatinib binding.** Since the ΔG<sub>bind</sub> energy associated with the binding of Imatinib to KIT<sup>V560G</sup> is low compared to the native protein and to the other studied targets, we evaluated the role of the partial JMR fragment, considered for the calculations in this
mutant. The newly generated MD trajectories of the Imatinib•KIT complexes (where
KIT = KITWT, KITS628N and KITD816V, having the same portion of the JMR (partially cleaved)
as the mutant KITV560G), were used for the ΔGbind energy calculation.

Comparing the data (MM-PBSA) obtained for Imatinib•KIT complexes formed by (I) the
targets having the partial JMR (Fig C in S1 File) and (II) the targets with the entirely cleaved
JMR (Fig 4B), we observed that (i) the tendency of the computed free energy of binding in the
two types studied complexes, formed by either the targets with totally cleaved JMR or having
partially cleaved JMR, is maintained; (ii) the free binding energy values are systematically
lower in all complexes of type (I) in respect to those of type (II); (iii) the difference in free
energy binding between the imatinib-bound complexes formed by KITD816V and KITWT con-
sists of 20–30% in the two types of complexes (I) and (II); (iv) the values of ΔGbind in com-
plexes formed by KITS628N and KITWT targets, having the partial JMR, are close. All these
observations indicate that the length of receptor, of type (I) or (II), influences the free energy
values but has only a tiny impact on the qualitative results distinguishing the energy binding to
the different targets. The last observation, (iv), corresponds well to the experimentally mea-
sured data that indicate the similar sensitivity of these proteins to Imatinib.

**Contributions of residues positioned out of the ATP-binding site to the binding energy.** Considering the inhibitor-target recognition as a cooperative process, we showed that the protein residues positioned at the vicinity of the ATP-binding site might affect the global binding energy. We have subtracted the energy contributions of the residues from mutant with those of the native protein, and the obtained increments (ΔΔG), ‘positive’ or ‘negative’, are illustrated in Fig 6.

In CSF-1R, except the most important impact of D802, several residues participate to either
increase the unfavorable contribution (K616) or reduce the favorable contribution (E626, I775
and D837) to the binding energy. With exception of I775, all these residues are not bound
directly to Imatinib (Fig 4). As a ‘positive’ effect on the Imatinib-CSF-1R recognition, we noted
only a diminishing of the unfavorable contribution (R777) and a small increasing of the favor-
able contribution (E628) in the mutant respectively to the native protein. The difference in the
residual contributions (ΔΔG) in KITD816V, in general, is close to those in CSF-1RD802V. In par-
ticular, the residues contributing negatively to the Imatinib binding, besides of the important
impact of the D816V substitution, are D810, making H-bond with Imatinib, and K623, not
interacting directly with Imatinib.

Such repartition of the favorable and non-favorable contributions resulted in a global
diminishing of Imatinib affinity to CSF-1RD802V and KITD816V mutants, explaining its resis-
tance to the drug. Curiously, all H-bonds stabilizing Imatinib binding by KIT are preserved in
KITD816V mutant. Moreover, D810 forms more prevalent H-bonds (60%) in Imatinib•
KITD816V complex compared to Imatinib• KITWT. It is not excluded that in KIT, together with
electrostatic interactions, the dispersion contribution is also important in the Imatinib binding.
According to the literature, in KITWT, the disperse interaction energy was estimated to be -7.1
kcal/mol, while the electrostatic contribution is 4.7 kcal/mol [27]. Moreover, the water mole-
cules localized in the binding pockets may contribute significantly to the Imatinib affinity, as
was reported in this study. Despite the H-bond pattern conservation, the overall energy in the
Imatinib complex formed by the KIT mutant is considerably changed as we demonstrated
using MM/PBSA method. A more detailed computation based on alchemical FEP/MD simula-
tions could be useful to investigate further the molecular determinants of the KITD816V
resistance.

Comparing KITWT and KITV560G, we noted that the ‘positive’ increment in the binding
energy of the mutant is furnished by residues E633, E640, D792, and R815 that contribute sig-
nificantly to the Imatinib affinity either by impact of the negatively charged residues E633,
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Binding energy of the Imatinib•target complexes. Difference in the per residues contribution to the binding energy (ΔΔG) between the native and mutated targets—(top) CSF-1R<sup>WT</sup> and CSF-1R<sup>D802V</sup>(blue), and (bottom) KIT<sup>WT</sup>, KIT<sup>V560G</sup> (orange), KIT<sup>S628N</sup> (green) and KIT<sup>D816V</sup> (red). The labeled residues are characterized by ΔΔG values more than 1 or less than -1 kcal/mol.</figcaption></figure>

doi:10.1371/journal.pone.0160165.g006

E640 and D792 or by diminishing of the unfavorable interaction by R815. Such increase of the favorable interactions could be a main factor contributing to the higher affinity of Imatinib to KIT<sup>V560G</sup> compare to the native target (Fig 4B). In this protein, a partial JMR may affect the global free binding energy, contributing favorably to the ΔΔG value in the mutant (Table C and Table D in S1 File). The energy contribution of the non-ATP binding site residues increases the global energy binding of all KIT targets, decreasing the ΔΔG value, in particular for the sensitive mutants KIT<sup>WT</sup> and KIT<sup>S628N</sup>.

In Imatinib•KIT<sup>S628N</sup> complex, two residues, I789 and R791, located in the C-loop, supply the most significant positive decrement to the binding energy, while residue D810 contributes to the ‘negative’ increment (Fig 6). The impact of these two opposite effects is compensated completely. Moreover, we observed that during the MD simulations of this complex, the methyl-piperazinyl moiety of Imatinib abolished its interaction with I789 and established a bifurcated H-bonding with D810 (Fig 4C). As already discussed, such alternation of the H-bond pattern either is induced by a flipping of the methyl-piperazinyl fragment inside the
ATP-binding site or promotes such flipping, suggesting the induced fit of imatinib by the ATP-binding site adjustment.

**Formation of the salt-bridges in Imatinib•target complexes.** Since we observed that the charged residues have the highest impact on the binding energy contributions in all studied complexes, we suggested a possible alteration of the salt-bridges in the mutants, especially in the vicinity of the Imatinib-binding site. The salt-bridges, a peculiar type of non-covalent interaction, combines two main components—hydrogen bond and electrostatic interaction—between two ionized sites.

Comparing the salt-bridges formed in the binding site of CSF-1R<sup>WT</sup> and CSF-1R<sup>D802V</sup>, we observe that the salt-bridges D778•••R801 and D796•••K616, established in the native protein, were completely disrupted in the mutant (Fig 7). As was identified through the MM-PBSA computing, the residues establishing salt-bridges in CSF-1R are important contributors to the total free binding energy with Imatinib, especially K616 (Table A in S1 File). All these residues are located at proximity of the Imatinib binding pockets, in particular, D778 is in the C-loop, D796 and R801 in the A-loop, and K616 in the β3 turn. The perturbed orientation of the side chains of R801 and D796 in CSF-1R<sup>D802V</sup> mutant, correlated to the changed atomic fluctuations (RMSFs), promotes their participation in an alternative salt-bridges pattern, forming in particular, the bifurcate salt-bridges D778•••H776 and D796••• H776 (Fig 7).

In the native KIT, the salt-bridge D792•••R815 is equivalent to that found in CSF-1R<sup>WT</sup>. In addition, D810 is bound to R791 (R791•••D810), completing a pattern of strong and paired salt-bridge interactions stabilizing the KIT binding pocket, which is smaller than that in CSF-1R (Fig 7). In KIT mutants, the salt-bridge pattern is considerably changed. In KIT<sup>D816V</sup> mutant, the most resistant to Imatinib, the salt-bridge D792•••R815 is absent. In KIT<sup>V560G</sup>, sensitive to Imatinib, residue D792 interacts either with R791 or with R815. In addition, in this mutant we observed the salt-bridge interaction D810•••K623 found in CSF-1R<sup>WT</sup>. The Imatinib complex formed by mutant KIT<sup>S628N</sup> shows salt-bridges pattern similar to KIT<sup>WT</sup>, nevertheless the interactions stability is different, reflected in a lower occurrence when compared to the native target.

Such change of the salt-bridges pattern in the mutants compared to the native proteins is prompted by an alternation of the atomic fluctuations (measured by the RMSFs) observed principally in molecular fragments containing the residues involved in the salt-bridges formation. Our results suggest that in Imatinib•target complexes formed by the native receptors (CSF-1R<sup>WT</sup>/KIT<sup>WT</sup>) and by their mutants, sensitive (KIT<sup>S628N</sup>) and very sensitive (KIT<sup>V560G</sup>) to Imatinib, the charges of residues localized at proximity of the binding pocket are equilibrated through alternative salt-bridges which apparently contributes to the stabilization of the binding pocket, providing the high affinity of Imatinib. In complexes formed by the resistant mutants, CSF-1R<sup>D802V</sup> and KIT<sup>D816V</sup>, the limited and locally positioned salt-bridges favor an excess of positive charges on residues non-involved in salt-bridges at proximity of the binding pocket (Fig D in S1 File) that would interfere with the Imatinib binding.

## Conclusions

The gain-of-function mutations in RTKs alter not only their tight factor-depending regulation but also influence their sensitivity to inhibitors, inducing drug resistance. In particular, it was demonstrated that the V560G substitution in KIT promotes higher sensitivity to Imatinib, while the substitution D816V in KIT and D802V in CSF-1R induces resistance to this drug either systematically (in KIT), or occasionally (in CSF-1R).

Our *in silico* calculation performed on molecular dynamics trajectories demonstrate that the native proteins and their mutants show different binding energy values with Imatinib, correlating with the experimental data. The main factor that drives the targets responsiveness to
<figure><img src="image_10.png" /><figcaption>**Fig 7.** Salt-bridges at proximity of the binding pocket calculated over the MD simulations of the Imatinib•target complexes. Representative conformations of the proteins are shown as cartoon, Imatinib (in blue) and residues forming salt-bridges (in orange) are drawn in stick. Non-covalent interactions are shown by dashed lines.</figcaption></figure>

doi:10.1371/journal.pone.0160165.g007
drug, sensitivity or resistance, is the electrostatic interactions between the protonated inhibitor and the negatively charged residues in the ATP-binding site or in the proximity of this site. The per residue energy decomposition indicates that the Asp to Val substitution in CSF-1R<sup>D802V</sup> and KIT<sup>D816V</sup> contributes the most to falling of the binding energy compared to the other KIT mutants. The energy calculations of complexes formed with the KIT targets having either the partially or entirely cleaved JMR showed that the JMR role in the Imatinib binding is accessory. The salt-bridges pattern suggests that in the resistant mutants, CSF-1R<sup>D802V</sup> and KIT<sup>D816V</sup>, a charge redistribution within residues at the vicinity of the ATP-binding site favors the repulsion of the positively charged Imatinib.

**Supporting Information**

**S1 File.** **Scheme A. Structural formula of Imatinib.** Atoms participating in the H-bonding with the targets are shown. **Fig A. Imatinib docking (best poses) into the targets, CSF-1R and KIT.** Structural formula of Imatinib is shown in the first box. Imatinib and protein residues that interact directly with Imatinib are represented in sticks and the protein backbone is shown in grey as cartoon. H-bonds between the protein and the ligand are shown as dotted lines. **Fig B.** **MD simulations of the KIT-Imatinib complexes.** (A) The Root Mean Square Deviation (RMSD) values of the complexes formed by Imatinib and KIT<sup>WT</sup> (top) and its mutant KIT<sup>D816V</sup> (bottom) were calculated for the backbone atoms from replicas 1 and 2 separately (black and red) of MD simulations of the studied systems. (B) Root Mean Square Fluctuations (RMSFs) computed on the backbone atoms from replicas 1 and 2 separately (black and red). **Fig C.** **Binding energy of imatinib•KIT targets.** Graphical representation of the free binding energy (ΔG) of Imatinib•target complexes. The MM-PBSA calculations were performed for KIT<sup>WT</sup>, KIT<sup>S628N</sup> and KIT<sup>D816V</sup> containing the partial JMR, identic to that in the mutant KIT<sup>V560G</sup>. The total ΔG energy value is shown for each complex. The ΔG value for KIT<sup>V560G</sup> is reproduced for comparison. **Fig D. Electrostatic potential (EP) surface of the Imatinib•target complexes.** EP surface for each complex was obtained on a representative equilibrated conformation before the MD simulations. EP calculations on the Connolly solvent-accessible surfaces of the receptors were performed with the APBS software. The color scale ranges from red (electronegative potential) through white (neutral) to blue (electropositive potential). **Table A. Imatinib binding to the targets.** Occurrences of the H-bonds between Imatinib and the targets residues were calculated for each of the two MD replicas (1 and 2). Imatinib atoms participating in the interactions with the targets are represented in scheme A in **S1 File**. The atom pairs for donor-acceptor interactions are depicted in the Table. **Table B. The free energy of binding in the Imatinib-RTK complexes (RTK = KIT<sup>WT</sup> and KIT<sup>D816V</sup>)** The free binding energy (ΔG<sub>bind</sub>) and contributions of electrostatic (ΔG<sub>elect</sub>), van der Waals (ΔG<sub>vdw</sub>) and solvation (ΔG<sub>solv</sub>) energies were calculated using MM-PBSA for each of the two MD replicas (1 and 2). All energies are in kcal/mol. **Table C. Per residue contribution to the final binding energy in Imatinib•CSF-1R complexes.** Contribution of all residues from the ATP binding pockets and the point mutation (distinguished by yellow and orange background, respectively), and all residues contributing with ΔΔG > 4 kcal/mol were considered. The favorable and unfavorable contributions of residues are highlighted in green and red colour, respectively. Std is standard deviation. **Table D. Per residue contribution to the final binding energy in Imatinib•KIT complexes.** Contribution of all residues from the ATP binding pockets and the point mutation (distinguished by yellow and by orange background, respectively) and all charged residues contributing with ΔΔG > 4 kcal/mol were considered. The favorable and unfavorable contributions of residues are highlighted in green and red colour, respectively. Std is standard deviation.
(DOCX)
Acknowledgments

The authors thanks Dr. Sjoerd de Vries J. for useful discussions. We also are grateful to the LMDM together with LabFisbio, at UFRJ, for attribution of the computer time on their in-house cluster.

Author Contributions

Conceived and designed the experiments: LT ES. Performed the experiments: PG NP SL. Analyzed the data: PG ICB NP PS PP LT. Contributed reagents/materials/analysis tools: SL PS. Wrote the paper: PG ICB LT.

References

1. Schlessinger J (2000) Cell signaling by receptor tyrosine kinases. Cell 103: 211-225. S0092-8674(00) 00114-8 [pii]. PMID: 11057895

2. Nolen B, Taylor S, Ghosh G (2004) Regulation of protein kinases; controlling activity through activation segment conformation. Mol Cell 15: 661-675. doi: 10.1016/j.molcel.2004.08.024; S1097276504004800 [pii]. PMID: 15350212

3. Lemmon MA, Schlessinger J (2010) Cell Signaling by Receptor Tyrosine Kinases. Cell 141: 1117-1134. doi: 10.1016/j.cell.2010.06.011 PMID: 20602996

4. Blume-Jensen P, Hunter T (2001) Oncogenic kinase signalling. Nature 411: 355-365. doi: 10.1038/35077225; 35077225 [pii]. PMID: 11357143

5. Dibb NJ, Dilworth SM, Mol CD (2004) Switching on kinases: oncogenic activation of BRAF and the PDGFR family. Nat Rev Cancer 4: 718-727. doi: 10.1038/nrc1434; nrc1434 [pii]. PMID: 15343278

6. Roussel MF, Sherr CJ (2003) Oncogenic potential of the c-FMS proto-oncogene (CSF-1 receptor). Cell Cycle 2: 5-6. 295 [pii]. PMID: 12695672

7. Rademakers R, Baker M, Nicholson AM, Rutherford NJ, Finch N, Soto-Ortolaza A, et al. (2012) Mutations in the colony stimulating factor 1 receptor (CSF1R) gene cause hereditary diffuse leukoencephalopathy with spheroids. Nat Genet 44: 200-205. ng.1027 [pii];10.1038/ng.1027.

8. Kitagawa D, Yokota K, Gouda M, Narumi Y, Ohmoto H, Nishiwaki E, et al. (2013) Activity-based kinase profiling of approved tyrosine kinase inhibitors. Genes Cells 18: 110-122. doi: 10.1111/gtc.12022 PMID: 23279183

9. Cohen MH, Williams G, Johnson JR, Duan J, Gobburu J, Rahman A, et al. (2002) Approval summary for imatinib mesylate capsules in the treatment of chronic myelogenous leukemia. Clin Cancer Res 8: 935-942.

10. Piccaluga PP, Rondoni M, Paolini S, Rosti G, Martinelli G, Baccarani M (2007) Imatinib mesylate in the treatment of hematologic malignancies. Expert Opin Biol Ther 7: 1597-1611. doi: 10.1517/14712598.7.10.1597 PMID: 17916051

11. Dagher R, Cohen M, Williams G, Rothmann M, Gobburu J, Robbie G, et al. (2002) Approval summary: imatinib mesylate in the treatment of metastatic and/or unresectable malignant gastrointestinal stromal tumors. Clin Cancer Res 8: 3034-3038. PMID: 12374669

12. Tamborini E, Bonadiman L, Greco A, Albertini V, Negri T, Gronchi A, et al. (2004) A new mutation in the KIT ATP pocket causes acquired resistance to imatinib in a gastrointestinal stromal tumor patient. Gastroenterology 127: 294-299. S0016508504002203 [pii]. PMID: 15236194

13. Gorre ME, Mohammed M, Ellwood K, Hsu N, Paquette R, Rao PN, et al. (2001) Clinical resistance to STI-571 cancer therapy caused by BCR-ABL gene mutation or amplification. Science 293: 876-880. doi: 10.1126/science.1062538; 1062538 [pii]. PMID: 11423618

14. Cools J, DeAngelo DJ, Gotlib J, Stover EH, Legare RD, Cortes J, et al. (2003) A tyrosine kinase created by fusion of the PDGFRA and FIP1L1 genes as a therapeutic target of imatinib in idiopathic hypereosinophilic syndrome. N Engl J Med 348: 1201-1214. doi: 10.1056/NEJMoa025217; 348/13/1201 [pii]. PMID: 12660384

15. Jabbour E, Kantarjian H, Jones D, Talpaz M, Bekele N, O'Brien S, et al. (2006) Frequency and clinical significance of BCR-ABL mutations in patients with chronic myeloid leukemia treated with imatinib mesylate. Leukemia 20: 1767-1773. 2404318 [pii]; doi: 10.1038/sj.leu.2404318 PMID: 16855631

16. Griswold IJ, MacPartlin M, Bumm T, Goss VL, O'Hare T, Lee KA, et al. (2006) Kinase domain mutants of Bcr-Abl exhibit altered transformation potency, kinase activity, and substrate utilization, irrespective of sensitivity to imatinib. Mol Cell Biol 26: 6082-6093. 26/16/6082 [pii]; doi: 10.1128/MCB.02202-05 PMID: 16880519
17. Cang S, Liu D (2008) P-loop mutations and novel therapeutic approaches for imatinib failures in chronic myeloid leukemia. J Hematol Oncol 1: 15. 1756-8722-1-15 [pii]; doi: 10.1186/1756-8722-1-15 PMID: 18828913

18. Frost MJ, Ferrao PT, Hughes TP, Ashman LK (2002) Juxtamembrane mutant V560GKit is more sensitive to Imatinib (STI571) compared with wild-type c-kit whereas the kinase domain mutant D816VKit is resistant. Mol Cancer Ther 1: 1115-1124. PMID: 12481435

19. Ma Y, Zeng S, Metcalfe DD, Akin C, Dimitrijevic S, Butterfield JH, et al. (2002) The c-KIT mutation causing human mastocytosis is resistant to STI571 and other KIT kinase inhibitors; kinases with enzymatic site mutations show different inhibitor sensitivity profiles than wild-type kinases and those with regulatory-type mutations. Blood 99: 1741-1744. PMID: 11861291

20. Diaz-Cano SJ (2012) Tumor heterogeneity: mechanisms and bases for a reliable application of molecular marker design. Int J Mol Sci 13: 1951-2011. doi: 10.3390/ijms13021951; ijms-13-01951 [pii]. PMID: 22408433

21. Diaz-Cano SJ (2015) Pathological bases for a robust application of cancer molecular classification. Int J Mol Sci 16: 8655-8675. ijms16048655 [pii]; doi: 10.3390/ijms16048655 PMID: 25898411

22. Seeliger MA, Nagar B, Frank F, Cao X, Henderson MN, Kuriyan J (2007) c-Src binds to the cancer drug imatinib with an inactive Abl/c-Kit conformation and a distributed thermodynamic penalty. Structure 15: 299-311. S0969-2126(07)00070-6 [pii]; doi: 10.1016/j.str.2007.01.015 PMID: 17355866

23. Seeliger MA, Ranjitkar P, Kasap C, Shan Y, Shaw DE, Shah NP, et al. (2009) Equally potent inhibition of c-Src and Abl by compounds that recognize inactive kinase conformations. Cancer Res 69: 2384-2392. 0008-5472.CAN-08-3953 [pii]; doi: 10.1158/0008-5472.CAN-08-3953 PMID: 19276351

24. Jacobs MD, Caron PR, Hare BJ (2008) Classifying protein kinase structures guides use of ligand-selectivity profiles to predict inactive conformations: structure of Ick/imatinib complex. Proteins 70: 1451-1460. doi: 10.1002/prot.21633 PMID: 17910071

25. Mol CD, Dougan DR, Schneider TR, Skene RJ, Kraus ML, Scheibe DN, et al. (2004) Structural basis for the autoinhibition and STI-571 inhibition of c-Kit tyrosine kinase. J Biol Chem 279: 31655-31663. doi: 10.1074/jbc.M403319200; M403319200 [pii]. PMID: 15123710

26. Aleksandrov A, Simonson T (2010) Molecular Dynamics Simulations Show That Conformational Selection Governs the Binding Preferences of Imatinib for Several Tyrosine Kinases. Journal of Biological Chemistry 285: 13807-13815. doi: 10.1074/jbc.M110.109660 PMID: 20200154

27. Lin YL, Roux B (2013) Computational Analysis of the Binding Specificity of Gleevec to Abl, c-Kit, Lck, and c-Src Tyrosine Kinases. J Am Chem Soc 135: 14741-14753. doi: 10.1021/ja405939x PMID: 24001034

28. Vinteler E, Stan NF, Luchian R, Cainap C, Ramalho JP, Chis V (2015) Conformational landscape and low lying excited states of imatinib. J Mol Model 21: 84. doi: 10.1007/s00894-015-2639-8 PMID: 25764326

29. Pereira EG, Moreira MA, Caffarena ER (2012) Molecular interactions of c-ABL mutants in complex with imatinib/nilotinib: a computational study using linear interaction energy (LIE) calculations. J Mol Model 18: 4333-4341. doi: 10.1007/s00894-012-1436-x PMID: 22570081

30. Lovera S, Sutto L, Boubeva R, Scapozza L, Dolker N, Gervasio FL (2012) The different flexibility of c-Src and c-Abl kinases regulates the accessibility of a druggable inactive conformation. J Am Chem Soc 134: 2496-2499. doi: 10.1021/ja210751t PMID: 22280319

31. Cowan-Jacob SW, Mobitz H, Fabbro D (2009) Structural biology contributions to tyrosine kinase drug discovery. Current Opinion in Cell Biology 21: 280-287. doi: 10.1016/j.ceb.2009.01.012 PMID: 19208462

32. Dubey KD, Ojha RP (2011) Binding free energy calculation with QM/MM hybrid methods for Abl-Kinase inhibitor. J Biol Phys 37: 69-78. doi: 10.1007/s10867-010-9199-z; 9199 [pii]. PMID: 22210962

33. Deng Y, Roux B (2009) Computations of standard binding free energies with molecular dynamics simulations. J Phys Chem B 113: 2234-2246. doi: 10.1021/jp807701h; 10.1021/jp807701h [pii]. PMID: 19146384

34. Lovera S, Morando M, Pucheta-Martinez E, Martinez-Torrecuadrada JL, Saladino G, Gervasio FL (2015) Towards a Molecular Understanding of the Link between Imatinib Resistance and Kinase Conformational Dynamics. PLoS Comput Biol 11: e1004578. doi: 10.1371/journal.pcbi.1004578; PCOMP-BIOL-D-15-00577 [pii]. PMID: 26606374

35. Lee TS, Potts SJ, Albitar M (2009) Basis for resistance to imatinib in 16 BCR-ABL mutants as determined using molecular dynamics. Recent Pat Anticancer Drug Discov 4: 164-173. PMID: 19519539

36. Lee TS, Ma W, Zhang X, Giles F, Cortes J, Kantarjian H, et al. (2008) BCR-ABL alternative splicing as a common mechanism for imatinib resistance: evidence from molecular dynamics simulations. Mol
Cancer Ther 7: 3834–3841. 1535-7163.MCT-08-0482 [pii]; doi: 10.1158/1535-7163.MCT-08-0482
PMID: 19056677

37. Bunney TD, Wan S, Thiyagarajan N, Sutto L, Williams SV, Ashford P, et al. (2015) The Effect of Mutations on Drug Sensitivity and Kinase Activity of Fibroblast Growth Factor Receptors: A Combined Experimental and Theoretical Study. EBioMedicine 2: 194–204. doi: 10.1016/j.ebiom.2015.02.009
PMID: 26097890

38. Rajasekaran R, Sethumadhavan R (2010) In silico identification of significant detrimental missense mutations of EGFR and their effect with 4-anilinoquinazoline-based drugs. Appl Biochem Biotechnol 160: 1723–1733. doi: 10.1007/s12010-009-8662-4 PMID: 19455431

39. Rajasekaran R, Sethumadhavan R (2010) Exploring the cause of drug resistance by the detrimental missense mutations in KIT receptor: computational approach. Amino Acids 39: 651–660. doi: 10.1007/s00726-010-0486-6 PMID: 20140688

40. Da Silva Figueiredo Celestino Gomes, Panel N, Laine E, Pascutti PG, Solary E, Tchertanov L (2014) Differential effects of CSF-1R D802V and KIT D816V homologous mutations on receptor tertiary structure and allosteric communication. PLoS One 9: e97519. doi: 10.1371/journal.pone.0097519; PONE-D-14-07239 [pii]. PMID: 24828813

41. Laine E, Chauvot de Beauchêne I, Perahia D, Auclair C, Tchertanov L (2011) Mutation D816V Alters the Internal Structure and Dynamics of c-KIT Receptor Cytoplasmic Region: Implications for Dimerization and Activation Mechanisms. PLoS Comput Biol 7: e1002068. doi: 10.1371/journal.pcbi.1002068; PCOMPBIOL-D-10-00316 [pii]. PMID: 21698178

42. Taylor JR, Brownlow N, Domin J, Dibb NJ (2006) FMS receptor for M-CSF (CSF-1) is sensitive to the kinase inhibitor imatinib and mutation of Asp-802 to Val confers resistance. Oncogene 25: 147–151. 1209007 [pii]; doi: 10.1038/sj.onc.1209007 PMID: 16170366

43. Vita M, Tisserand JC, Chauvot de Beauchêne I, Panel N, Tchertanov L, Agopian J, et al. (2014) Characterization of S628N: a novel KIT mutation found in a metastatic melanoma. JAMA Dermatol 150: 1345–1349. 1913663 [pii]; doi: 10.1001/jamadermatol.2014.1437 PMID: 25317746

44. Chauvot de Beauchêne I, Allain A, Panel N, Laine E, Trouve A, Dubreuil P, et al. (2014) Hotspot mutations in KIT receptor differentially modulate its allosterically coupled conformational dynamics: impact on activation and drug sensitivity. PLoS Comput Biol 10: e1003749. doi: 10.1371/journal.pcbi.1003749; PCOMPBIOL-D-13-02035 [pii]. PMID: 25079768

45. Lyman E, Zuckerman DM (2006) Ensemble-based convergence analysis of biomolecular trajectories. Biophysical Journal 91: 164–172. doi: 10.1529/biophysj.106.082941 PMID: 16617086

46. Sherman W, Day T, Jacobson MP, Friesner RA, Farid R (2006) Novel procedure for modeling ligand/receptor induced fit effects. Journal of Medicinal Chemistry 49: 534–553. doi: 10.1021/jm050540c
PMID: 16420040

47. Mackerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, et al. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 102: 3586–3616. doi: 10.1021/jp973084f PMID: 24889800

48. Zoete V, Cuendet MA, Grosdidier A, Michielin O (2011) SwissParam: a fast force field generation tool for small organic molecules. J Comput Chem 32: 2359–2368. doi: 10.1002/jcc.21816 PMID: 21541964

49. van der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJ (2005) GROMACS: fast, flexible, and free. J Comput Chem 26: 1701–1718. doi: 10.1002/jcc.20291 PMID: 16211538

50. Jorgensen WL, Jenson C (1998) Temperature dependence of TIP3P, SPC, and TIP4P water from NPT Monte Carlo simulations: Seeking temperatures of maximum density. Journal of Computational Chemistry 19: 1179–1186.

51. Bussi G, Donadio D, Parrinello M (2007) Canonical sampling through velocity rescaling. J Chem Phys 126: 014101. doi: 10.1063/1.2408420 PMID: 17212484

52. Berendsen HJC, Postma JPM, Vangunsteren WF, Dinola A, Haak JR (1984) Molecular-Dynamics with Coupling to an External Bath. Journal of Chemical Physics 81: 3684–3690.

53. Van Gunsteren WF, Berendsen HJC (1988) A Leap-frog Algorithm for Stochastic Dynamics. Molecular Simulation 1: 173–185.

54. Darden T, York D, Pedersen L (1993) Particle Mesh Ewald—An N.Log(N) Method for Ewald Sums in Large Systems. Journal of Chemical Physics 98: 10089–10092.

55. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS: A linear constraint solver for molecular simulations. Journal of Computational Chemistry 18: 1463–1472.

56. Miyamoto S, Kollman PA (1993) Absolute and relative binding free energy calculations of the interaction of biotin and its analogs with streptavidin using molecular dynamics/free energy perturbation approaches. Proteins 16: 226–245. doi: 10.1002/prot.340160303 PMID: 8346190
57. DeLano WL, Lam JW. PyMOL: A Communications Tool for Computational Models. Abstr. Pap. Am. Chem. Soc. 2005, 230, U1371–U1372.

58. Humphrey W, Dalke A, Schulten K (1996) VMD: Visual molecular dynamics. Journal of Molecular Graphics 14: 33–38. PMID: 8744570

59. Lee MR, Duan Y, Kollman PA (2000) Use of MM-PB/SA in estimating the free energies of proteins: application to native, intermediates, and unfolded villin headpiece. Proteins 39: 309–316. PMID: 10813813

60. Konecny R, Baker NA, McCammon JA (2012) iAPBS: a programming interface to Adaptive Poisson-Boltzmann Solver (APBS). Comput Sci Discov 5. doi: 10.1088/1749-4699/5/1/015005

61. Kumari R, Kumar R, Lynn A (2014) g_mmpbsa, a GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf Model 54: 1951–1962. doi: 10.1021/ci500020m PMID: 24850022

62. Brown S. P. & Muchmore S. W. (2009). Large-scale application of high-throughput molecular mechanics with Poisson-Boltzmann surface area for routine physics-based scoring of protein-ligand complexes. *J. Med. Chem.* 52, 3159–3165. doi: 10.1021/jm801444x PMID: 19385614

63. Ihaka R, Gentleman R (1996) A Language for Data Analysis and Graphics. J Comput Graph Stat 5: 299–314. doi: 10.1080/10618600.1996.10474713

64. Chauvot de Beauchêne I (2013) Etude des mécanismes d'activation et de résistance du récepteur tyrosine kinase KIT sauvage et mutant. Ph.D thesis, Ecole Normale Supérieure de Cachan, France. Available: https://drive.google.com/open?id=0BwzanqmJYoa0SER6VWJaX29zMTQ.

65. Tap WD, Wainberg ZA, Anthony SP, Ibrahim PN, Zhang C, Healey JH, et al. (2015) Structure-Guided Blockade of CSF1R Kinase in Tenosynovial Giant-Cell Tumor. N Engl J Med 373: 428–437. doi: 10.1056/NEJMoa1411366 PMID: 26222558

66. Walter M, Lucet IS, Patel O, Broughton SE, Bamert R, Williams NK, et al. (2007) The 2.7 A crystal structure of the autoinhibited human c-Fms kinase domain. J Mol Biol 367: 839–847. S0022-2836(07)00076-9 [pii]; doi: 10.1016/j.jmb.2007.01.036 PMID: 17292918

67. Cowan-Jacob SW, Fendrich G, Floersheimer A, Furet P, Liebetanz J, Rummel G, et al. (2007) Structural biology contributions to the discovery of drugs to treat chronic myelogenous leukaemia. Acta Crystallogr D Biol Crystallogr 63: 80–93. S0907444906047287 [pii]; doi: 10.1107/S0907444906047287 PMID: 17164530

68. Canning P, Tan L, Chu K, Lee SW, Gray NS, Bullock AN (2014) Structural mechanisms determining inhibition of the collagen receptor DDR1 by selective and multi-targeted type II kinase inhibitors. J Mol Biol 426: 2457–2470. S0022-2836(14)00198-3 [pii]; doi: 10.1016/j.jmb.2014.04.014 PMID: 24768818

69. Kollman PA, Massova I, Reyes C, Kuhn B, Huo SH, Chong L, et al. (2000) Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Accounts of Chemical Research 33: 889–897. doi: 10.1021/ar000033j PMID: 11123888

70. Jiang W, Roux B (2010) Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations. J Chem Theory Comput 6: 2559–2565. doi: 10.1021/ct1001768 PMID: 21857813

71. Fong CW (2015) Binding energies of tyrosine kinase inhibitors: Error assessment of computational methods for imatinib and nilotinib binding. Comput Biol Chem 58: 40–54. S1476-9271(15)30007-4 [pii]; doi: 10.1016/j.compbiolchem.2015.05.002 PMID: 26025598