Molecular dynamics systems preparation.
In the simulation system of GDP-bound KRASG12D with MRTX1133, the initial structures of GDP-bound KRASG12D were obtained from Protein Data Bank (PDB, entry codes: 7RPZ, 7RT1, 7RT2, 7RT3, 7RT4, 7RT5) [28]. All these structures of MRTX1133 or its analogs bound to KRASG12D are used for the simulations to sample a comprehensive conformational space for the subsequent cluster analysis and Markov state model (MSM).
In the simulation system of GTP-bound KRASG12D with MRTX1133, the GTP-bound KRASG12D structures were prepared by converting the GDP molecule in GDP-bound KRASG12D with a GTP molecule.
In the simulation system of GTP-bound KRASG12D without MRTX1133, the crystal structure of GppNHp-bound KRAS4BQ61H were selected (PDB ID: 3GFT). Residue H61 was mutated back to Q61, residue G12 was mutated to D12. GppNHp was reverted into a GTP molecule.
These structure models were prepared by the Protein Preparation Wizard module in Schrödinger 2015 suite software package [33, 34] including filling in the missing residues, adding the hydrogen atoms and optimizing proteins at physiological pH. All protonated amino acids are visually inspected in detail to confirm the formation of favorable interaction with the surrounding environment. All protein structures were solvated in pre-equilibrated TIP3P water cubic boxes. And the sodium and chloride ions were added to maintain each system electrically neutral at physiological concentrations.
Molecular dynamics simulation.
We carried out MD simulations for each system. All the simulations were conducted by Desmond [34] with OPLS3 force field [35]. After the energy minimizations, 200 ns MD simulations in NPT ensembles were performed for each system. The simulations are carried out for three times in parallel with random initial velocities. The temperature was set at 300 K by using the Nosé-Hoover method and the pressure was maintained at 1.01325 bar by using Martyna-Tobias Klein method.
Cluster analysis.
To obtain a sufficiently overall structural dynamics of GDP-bound and GTP-bound KRASG12D, all snapshots of KRASG12D proteins from trajectories were clustered based on the root-mean-square deviations (RMSD) values of Cα backbone using the k-means algorithm [36]. The RMSD of Cα backbone of KRASG12D proteins were calculated by the MDAnalysis [37–39]. A total of 719,928 snapshots were divided into 12 clusters (C1-C12) according to the silhouette coefficient and the sum of the squared errors (Fig. S1). The most representative structures of each cluster were extracted from the MD trajectories, separately.
Markov state model construction.
The Markov state model, built from the clustering results, were used to compute transition probability, stationary probability [40]. The transition probability is defined as:
where any value n = 0, 1, 2, 3, ..., the times t indexed up to this value of n: t0, t1, t2, ... all states i: i0, i1, i2, i3,… and pij is the solution of the forward equation.