Computational protocol: Towards a Molecular Understanding of the Link between Imatinib Resistance and Kinase Conformational Dynamics

Similar protocols

Protocol publication

[…] The kinase structures were retrieved from the Protein Data Bank (PDB entries 2G1T, 2SRC, 1PKG, 3KMM, 2WGJ and 2ITW). Missing residues were added using the software Modeller [], according to the respective Uniprot sequences. We used the Amber99SB*-ILDN [, ] force field, including backbone corrections by Hummer and Best [], with explicit TIP3P [] water molecules. The unbiased MD simulations were carried out with the ACEMD program [] running on GPUs. 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 at 1 atm. The temperature was kept at 305 K by a Langevin thermostat. A 400 ns production run was carried out for all the systems in the canonical (NVT) ensemble. The runs for Abl, Src and Kit were extended to a total length of 1 μs in both the DFG-in and DFG-out conformations. The simulations of the five Abl mutants (G250E, E279K, H396P, E450K and T315I) were carried out with the same setup in both the DFG-in and DFG-out conformation, after mutating in-silico the respective residue in the Abl structure. To capture the sub-μs dynamics, the RMSF of the C α were averaged on 30ns non-overlapping windows after discarding the first 100ns of each run. We also compared the diffusion in the space described by the first two principal component analysis vectors (see ).Since the conformational changes of TKs take place on time scales longer than those accessible by standard MD simulations, here we used a combination of enhanced sampling approaches. Parallel Tempering—Well Tempered Metadynamics [, ] in the well tempered ensemble (WTE) [, ] (PTmetaD) was chosen due to its proven ability to fully converge complex conformational free energy surfaces such as those relevant in kinases (including the DFG-flip). Indeed, the PTmetaD approach (both in the standard and WTE variants) has already successfully used to study the conformational dynamics and its associated free energy landscape in many kinases [, , , , ] and other flexible proteins [, ]. PTmetaD was performed using the software Gromacs 4.5 [] and the PLUMED plugin [], using an integration step of 2 fs. The particle mesh Ewald algorithm was used for electrostatic interactions. Temperature coupling was done with the V-rescale algorithm []. The WTE allowed the use of a reduced number of replicas compared to standard PTMetaD [, ]. An average exchange probability of 24% was obtained using 5 replicas in the temperature range 305–400 K. We used the same four collective variables (CVs) that were used to reconstruct the free energy surface (FES) associated with the DFG flip in Src and Abl [].They are shown in and defined in the following (SRC numbering): CV1 is the distance between the centre of mass of Asp404 (DFG motif) and Lys295. CV2 is the distance between the centre of mass of Phe405 and the Cβ of Ile293. CV3, is a function of 3 dihedral angles f ( ϕ 404, ψ 405, ψ 408 ) ranging from 3, when the three dihedral arguments correspond to the DFG-Asp-in position, to 0, when they are in the DFG-Asp-out conformation. CV4 is the distance between the centre of mass of residues Asn381, …, His384 and residues Ala408, …, Ile411 of the activation loop, known to form a β-sheet in the active conformation. The height of the Gaussians was set at 2.0 kJ/mol with a deposition rate of 1/2000 steps and a bias factor of 5. The Gaussian width used for the CVs was 0.1 for the dihedral similarity (CV3) and 0.3 Å for all the others.A minimum of 400 ns of sampling per replica in the NVT ensemble were needed to reach full convergence of the free energy. In the case of G250E and T315I, as the convergence was slower than in the other cases, we performed more than 600 ns and 1200 ns of sampling per replica, respectively. The total sampling time amounted to more than 14 μs across all mutants. The free energy surface reconstruction was obtained from the PT-metaD-WTE by reweighting the fixed potential energy bias and using two independent approaches: integrating the bias or using the time independent estimator of Ref. []. The convergence of the free energy reconstruction was monitored by integrating the cumulative added bias as a function of time (50 ns intervals) and comparing the reconstruction to that obtained by the time-independent estimator (as shown in , – Figs). Changes of all the CVs used were also monitored to guarantee that the system diffuses freely in the CV space and is able to visit all the basins several times (see ).The FES were also reprojected as a function of two other CVs describing the conformational change of the A-loop from open to closed. Two path collective variables [] were built by using the open and closed crystallographic structures of Src (CV1) and Abl (CV2). The reweight was performed by using our python implementation of the approach of Ref. [] (available on our homepage: https://www.ucl.ac.uk/chemistry/research/group_pages/prot_dynamics/. The entropic contribution to the DFG-in DFG-out free energy difference was computed by linear fitting of the free energy differences as a function of temperature. It has been shown that this approach is more accurate than competing ones. []. We also tested the robustness of the estimate with respect to the range of temperatures on both standard PT and WTE-PT (see ).The drug binding free energy surface was calculated using metadynamics with the same software described above. The General Amber Force Field (GAFF) was used for the ligand. The ligand charges were calculated at the HF level using a 631-G* basis set with the Gaussian03 [] package. QM-level torsional scans with a step of 10 degrees were carried out for the ca-ca-n-c and ca-ca-c3-n3 dihedrals which appeared to have wrong torsional profiles. The profiles with the GAFF force field in vacuum were then fitted to the QM ones to refine the dihedral parameters. As customary in the case of ligand binding [, ], we performed a preliminary metadynamics run using sub-optimal geometrical CV to obtain an initial pathway for the setup of the path collective variables (PCV) []. We selected 23 frames from the lowest free energy path obtained in the preliminary run and optimized this initial guess using the methodology described by Branduardi et al. []. To take into account possible rearrangements, we included the Cα atoms of the A-loop and of the αC-helix in the definition of the PCVs. The 2 PCVs s and z were used to run a 300 ns metadynamics. The height of the Gaussians was set at 2.0 kJ/mol with a deposition rate of 1/2000 steps and a bias factor of 10. The Gaussian width used for the CVs was 0.1 for s and 0.003 for z. The free energy corresponding to the last leg of the unbinding mechanism (from the external binding pose to a fully solvated state) was computed again by using Well-tempered metadynamics following the approach of Ref. []. […]

Pipeline specifications

Software tools MODELLER, GROMACS, PLUMED
Application Protein structure analysis
Organisms Homo sapiens
Diseases Leukemia, Myeloid, Neoplasms
Chemicals Tyrosine