Conformational Selection and Induced Fit Mechanisms in the Binding of an Anticancer Drug to the c-Src Kinase

Understanding the conformational changes associated with the binding of small ligands to their biological targets is a fascinating and meaningful question in chemistry, biology and drug discovery. One of the most studied and important is the so-called “DFG-flip” of tyrosine kinases. The conserved three amino-acid DFG motif undergoes an “in to out” movement resulting in a particular inactive conformation to which “type II” kinase inhibitors, such as the anti-cancer drug Imatinib, bind. Despite many studies, the details of this prototypical conformational change are still debated. Here we combine various NMR experiments and surface plasmon resonance with enhanced sampling molecular dynamics simulations to shed light into the conformational dynamics associated with the binding of Imatinib to the proto-oncogene c-Src. We find that both conformational selection and induced fit play a role in the binding mechanism, reconciling opposing views held in the literature. Moreover, an external binding pose and local unfolding (cracking) of the aG helix are observed.

Values decrease from yellow to red (the information is further displayed by the thickness of the backbone which is inversely proportional to the degree of order). In both panels residues for which data is not available are displayed in gray.         S12. Carbonyl (blue) and C↵ (red) chemical shift deviations from random coil values measured for the catalytic domain of c-Src. The figure reports centered running average over 5 residues. Elements of secondary structure found in the crystal (pdb 2OIQ) are schematically displayed and an ideal deviation of 4 ppm is reported as a continuous shadowed area (positive for alpha helices, negative for beta strands) Figure S13. Residues appearing in the bound NMR spectra and belonging to the main (green) or secondary (red) binding site as represented on top of the canonical bound structure (left, purple, corresponding to structure A in Figure 4) and of the secondary bound structure (right, blue, structure C in Figure 4). The residues highlighted in red are far apart in the bound structure and come closer to outline the secondary binding cavity in the secondary bound pose. Figure S14. Surface Plasmon Resonance spectra of imatinib binding to Src WT at 20 C.

Molecular Dynamics Simulations
The kinases structures were retrieved from the Protein Data Bank (PDB entries 2SRC and 2OIQ). Missing residues in 2OIQ were added using the software Modeller [2], according to the respective Uniprot sequence and using 2SRC as a template. For the apo structure we used the Amber99SB*-ILDN [3,4] force field, including backbone corrections by Hummer and Best [5] with explicit TIP3P [6] water molecules. For the holo structure we used Amber99SB*-ILDN [3,4] force field for the protein part and GAFF [7] for the ligand. The ligand charges were calculated at the HF level using a 631-G* basis set with the Gaussian03 [8] package. Preliminary simulations indicated that the default GAFF parameters for the dihedrals ca ca n c and ca ca c3 n3 didn't reproduce the correct torsional profile. For this reason, QM-level torsional scans with a step of 10 degrees were carried out for these dihedrals. Torsional profiles with the GAFF force field in vacuum were then fitted to the QM ones in order to derive more accurate parameters. All the systems were minimized with 10000 steps of conjugated gradient and equilibrated in the isothermal-isobaric (NPT) ensemble for 10 ns, using a Berendsen barostat to keep the pressure at 1 atm. The temperature was kept at 305 K by a Langevin thermostat. A 1 µs production run was carried out for all the systems in the canonical (NVT) ensemble. All the simulations were carried out with the ACEMD program [9] running on a server with four Nvidia Fermi GPUs. The van der Waals interactions were cut o↵ at 1.0 nm, and long range electrostatic interactions were calculated by the particle mesh Ewald algorithm with a mesh spaced 0.1 nm. These values have been shown to perform optimally on GPUs. The parallel-tempering metadynamics (PT-MetaD) simulations were performed using Gromacs 4.5.5 [10] combined with the PLUMED 1.3 [11] plug-in. The particle mesh Ewald (PME)-Switch algorithm was used for electrostatic interactions with a cut-o↵ of 1 nm. A single cut-o↵ of 1.2 nm was used for Van der Waals interactions. Temperature coupling was done with the V-re-scale algorithm [12]. We performed PT-MetaD with 5 replicas using the Well Tempered Ensemble [13] with 3 Collective Variables (CVs), the 2 Path Collective Variables (PCVs) s and z [14], plus a CV counting the number of water molecules interacting with both the ligand and the cavity. The bias factor was set to 15.0 and the Gaussians height to 1.25 kJ/mol, with a deposition rate of 1/2000 steps. The Gaussian width was set to 0.1, 0.005 and 0.1 for the three CVs, respectively. As customary in the case of ligand binding [15,16], we performed a preliminary metadynamics run using a simple geometrical CV to obtain an initial pathway for the setup of the PCVs. We selected 27 frames from the lowest energy path obtained in the preliminary run and optimized this initial guess using the methodology described by Branduardi et al. [14]. As the preliminary metadynamics showed large rearrangements of the A-loop, we included C↵ atoms of the loop and of the ↵C-helix in the definition of the PCVs. Two extra CVs defined as the distance between Asp404 and Lys295 and the distance between Phe405 and Leu317 (as in Lovera et al. 2012) were used to describe the DFG-flip. The sampling convergence was checked by comparing the reconstructed free energy surfaces at di↵erent time intervals during the last 50 ns of the simulations. As the sampling on top of the rate-determining TS was less exhaustive than in the main minima, we also performed Multiple-Walkers metadynamics using 5 walkers starting from the bound state and the same setup of the PTMetaD run. The value of the multiple-walkers determined tranistion state energy is within the error bar of the PTmetaD one.
In order to assess the agreement between the simulations and the dynamics observed by NMR, we used the MD trajectories to calculate the order parameter S 2 from the rotational correlation function of the backbone amide vector, according to Lipari and Szabo. In our case, in the context of an extended model [17], we used the more general expression introduced by Bremi et al. [18], according to which the internal motion correlation function is fitted with a sum of five exponentials where A 0 is analogous to the order parameter S 2 , while the other A k>0 are indicative of motions on di↵erent (slower) time-scales. The internal NH motion was decoupled from the overall motion by rotational and translational fitting to a reference structure. As for the RMSF, the analysis was repeated in windows to assess the accuracy of the estimate.

Sample preparation
Expression and purification were performed following the procedure published elsewhere [19]. Protein concentrations were determined by absorbance spectroscopy using a NanoDrop ND -1000 (✏ 280 =52370 M 1 cm 1 ). Stock solutions of the proteins were stored at 4 C, for short-term use, or at -20 C.

NMR assignment
Protein samples were prepared in 20 mM sodium phosphate at pH 6.4, 0.5 M NaCl, 1mM MgCl2, 1mM TCEP, and 0.03% NaN 3 . The temperature was set to 293 K. The choice of the pH was dictated by a compromise between protein solubility (decreasing at lower pH) and the detectability of water-exchanging amide protons. NMR experiments were performed on Bruker Avance 600 MHz, 700 MHz and 1 GHz, the latter equipped with cryogenically-cooled triple resonance probe. Sample concentration was 300 µM. 1 H and 15 N backbone assignment of the catalytic domain of Src was performed using a series of three-dimensional experiments recorded at both 700 MHz and 1 GHz, namely: trHNCO , trHN(CA)CO, trHNCA, trHNCACB and 1 H, 15 N -NOESY-HSQC recorded on a doubly or triply labeled sample with deuteration > 70%. The assignment was aided by the comparison with the assignment of the Src-Imatinib complex [19], corrected for the di↵erent pH by a titration from pH 8 to 6.4. Chemical shift values were further checked by comparison with predictions from the crystallographic structure and random coil using SPARTA+ [20] and ncIDP programs [21]. Spectra were processed with XWINNMR 2.1 program and analysed with NMRViewJ (One Moon Scientific) [22]. Carbon Ca and Cb chemical shifts were corrected for the deuterium isotope e↵ect [23]. Chemical shift Index was calculated for C↵ , C and carbonyl chemical shifts as described by Wishart et al.

NMR Relaxation
Relaxation experiments for the mobility measurements of free c-Src were performed at 700 and 1000 MHz using pulse sequences for 15 N R 1 and R 2 rates, and heteronuclear NOE spectra as described in the literature [24][25][26][27][28][29]. Relaxation delays for T 1 were: 50, 200, 400, 600, 1000, 1500, 2000, 3000 ms while for T2 measurements we used the following delays: 0.0, 6.3, 12.6, 18.9, 25.2, 31.4 and 37.7 ms. Due to the short transversal relaxation times the pulse sequence was modified in order to minimize the CPMG delay to 6 ms. This allowed to sample more time points in T 2 decays. Measurements were repeated in the presence of Imatinib (1.5 equivalents) at 700 MHz. Addition of Imatinib in ratio protein:drug 1: 1.5 was achieved adding few micro-liters of a concentrated (typically 50 mM) DMSO solution. TENSOR2 [30], was used to analyze the data, yielding the isotropic rotational correlation time ⌧ c and the Lipari-Szabo sub-nanosecond order parameters (S 2 ) [31].
The isotropic rotational correlation time ⌧ c estimated from the R2/R1 ratio for a solution 300 µM of free c-Src at 700 MHz is slightly larger than expected (26 ns). This is probably due to the onset of aggregation, as in the case of Abl [32]. When we repeated the analysis on the data obtained on a 1GHz spectrometer on a less concentrated sample (< 100 µM) , it resulted in a value for ⌧ c of 21 ns, which is close to the value expected for a monomeric protein of this size [33] and in excellent agreement with the value estimated by Hydronmr [1] using the crystallographic c-Src PDB structure (2OIQ).
However, at the lower concentration the signal to noise ratio was too low to estimate the S 2 for many of the observed peaks, while for the well resolved and intense peaks, the estimated S 2 is in good agreement with that obtained at 700MHz on the more concentrated solution. Indeed, the chemical shifts of the protein remain unaltered up to concentrations of 600 µM, a further indication that the onset of aggregation phenomena at 300 µM should have limited e↵ects on the structure and local dynamics of the protein.
Thus, we used the relaxation data obtained with a concentration of 300 µM at 700 MHz to estimate the order parameter S 2 . Due to the large number of residues with exchange contributions to R2 (large R1R2 products) or fast mobility (low heteronuclear NOE values), a spherical model was used for the rotational di↵usion tensor, as in the case of Abl kinase [32].
Deuterium exchange measurements H/D exchange was measured following the disappearance of peaks in the 1 H, 15 N -HSQC spectrum. 2D HSQC spectra were recorded at increasing times after resuspension of the lyophilized protein in D 2 O-based bu↵er (20 mM phosphate at pD 6.4, NaCl 0.5 M, MgCl2 1mM, TCEP 1mM, NaN3 0.03%.), extending to a final exchange time of 52.5 days (1260 hours). A total of 26 2D spectra were recorded for each of the two di↵erent samples (in the presence and in the absence of Imatinib), with variable number of scans (increasing from 64 in the first time-point to 1024 for the last ones), according to which intensities were also corrected for quantification of exchange rates. Exchange rates were estimated fitting peak intensities as a function of time by regression analysis with a three parameters exponential curve of the type Ae Bt + C, where B is the exchange rate k exch . Uncertainties in the exchange rates were obtained from 250 MonteCarlo simulations at a confidence interval of 90% and using the standard deviation of the spectral noise to account for intensity uncertainties. Protection factors (pf ) were estimated as the ratio k intr /k exch [34,35] where k intr is the intrinsic exchange rate calculated from the primary sequence using the SPHERE program [34] and k exch is the experimental value. Protection factors were converted in free energies of local folding in the assumption of EX2 regime [36] (where both the intrinsic exchange rate and the local unfolding kinetic constants are much larger than the local kinetic constant of folding). Such regime is considered quite common in proteins especially for folded protein below pH 9 [37]. Under this assumption the free energy can be calculated as RT ln(pf ). Peaks already disappeared in the first spectrum after dissolution in D 2 O were given an exchange rate equal to k intr . Peaks not decreasing in intensity after 52 days were given an exchange rate of 10 4 h 1 .

Surface Plasmon Resonance
All measurements were performed with a Biacore X (GE Healthcare) system. The protein was covalently bound to the activated carboxyl groups of the dextran molecules of the Sensor Chip CM5 at 25 C with running bu↵er PBS-P (10 mM phosphate bu↵er, 2.7 mM KCl, 137 mM NaCl, 0.05% P20). The dextran matrix was activated with a solution of EDC/NHS. The protein was diluted to 50 µg/ml in 10 mM sodium acetate pH 5.0 and injected over the surface until the target level of 12300 RU was reached. The remaining activated dextran was inactivated by injection of 1 M ethanolamine. Flow rate was 10 µl/min throughout the immobilization. The cell placed upstream to the active one was left untreated and served as reference cell. The theoretical R max value for the binding of imatinib (MW 494) to Src (MW 35476) was calculated as 170 RU. Kinetic assays were performed at 20 C in running bu↵er supplemented with 5% DMSO. Imatinib was diluted into running bu↵er with 5% DMSO and a serial dilution series was created with a 2-fold dilution factor to give concentrations of 1.56, 3.25, 6.5, and 12.5 µM. The analyte was injected for 30 seconds at 100 µl/min and the dissociation was monitored for 70 seconds. Experimental data were fit with a heterogeneous ligand model, using one global R max parameter, one global k a parameter, one global k d parameter and local RI (bulk refractive index) parameters set constant to 0.