Classifying reliably active and inactive molecular conformations of wildtype (WT) and mutated oncogenic proteins is a key, ongoing challenge in molecular cancer studies. Here, we probe the GTP-bound K-Ras4B conformational dynamics using long-time atomistic molecular dynamics (MD) simulations. We extract and analyze the detailed underlying free energy landscape of WT K-Ras4B. We use two key reaction coordinates, labeled d1 and d2 (i.e., distances coordinating the Pβ atom of the GTP ligand with two key residues, T35 and G60), shown to correlate closely with activities of WT and mutated K-Ras4B. However, our new K-Ras4B conformational kinetics study reveals a more complex network of equilibrium Markovian states. We show that a new reaction coordinate is required to account for the orientation of acidic K-Ras4B sidechains such as D38 with respect to the interface with binding effector RAF1 and rationalize the activation/inactivation propensities and the corresponding molecular binding mechanisms. We use this understanding to unveil how a relatively conservative mutation (i.e., D33E, in the switch I region) can lead to significantly different activation propensities compared with WT K-Ras4B. Our study sheds new light on the ability of residues near the K-Ras4B—RAF1 interface to modulate the network of salt bridges at the binding interface with the RAF1 downstream effector and, thus, to influence the underlying GTP-dependent activation/inactivation mechanism. Altogether, our hybrid MD-docking modeling approach enables the development of new in silico methods for quantitative assessment of activation propensity changes (e.g., due to mutations or local binding environment). It also unveils the underlying molecular mechanisms and facilitates the rational design of new cancer drugs.

The gene products of Kristen rat sarcoma virus oncogene (KRAS) homologs are notorious key cancer-related proteins for almost four decades.1,2 K-Ras proteins are small GTPases that control cellular proliferation by playing an important role in many cancer signal transduction pathways.3 K-Ras acts as a molecular switch, flipping between an inactive guanine diphosphate (GDP)-bound state and an active form of guanine triphosphate (GTP)-bound state.3–6 There are many known K-Ras mutations that have been associated with many fatal cancers such as colorectal cancer, pancreatic ductal adenocarcinoma, and lung cancer [Fig. 1(b)].7–10 Many computational and experimental efforts have been made to understand the conformational dynamics of K-Ras, the effects of its mutations, and to find mutation-specific drugs.3,11–14 Understanding the detailed thermodynamic and equilibrium kinetic behavior of wildtype (WT) K-Ras4B proteins remains an important, outstanding aim due to its involvement in a variety of cancers.3,4,15–17

FIG. 1.

(a) K-Ras4B structural elements. The switch I and switch II regions are highlighted in yellow and magenta, respectively. The GTP ligand is shown in licorice and colored by atom type. The Cα (for residues T35 and G60) and Pβ (for GTP) atoms are shown as blue spheres. d1 is the distance (dashed black line) between the Cα atom of G60 and the Pβ atom of GTP. d2 is the distance between the Cα atom of T35 and the same Pβ atom of GTP (dashed black line). The D33 residue is circled in red. (b) Schematic representation of the GTP-dependent activation of K-Ras4B and the hypothesized relationship between its active/inactive states and the d1 and d2 distances.

FIG. 1.

(a) K-Ras4B structural elements. The switch I and switch II regions are highlighted in yellow and magenta, respectively. The GTP ligand is shown in licorice and colored by atom type. The Cα (for residues T35 and G60) and Pβ (for GTP) atoms are shown as blue spheres. d1 is the distance (dashed black line) between the Cα atom of G60 and the Pβ atom of GTP. d2 is the distance between the Cα atom of T35 and the same Pβ atom of GTP (dashed black line). The D33 residue is circled in red. (b) Schematic representation of the GTP-dependent activation of K-Ras4B and the hypothesized relationship between its active/inactive states and the d1 and d2 distances.

Close modal

K-Ras has two splice variants, namely K-Ras4A and K-Ras4B. Both variants are oncogenic but have a distinct mechanisms of membrane-related and subcellular trafficking. K-Ras4A has 189 amino acids, like H-Ras and N-Ras. However, K-Ras4B has 188 amino acids. The two variants of K-Ras have a high degree of sequence identity of their catalytic domains, but a very low degree of sequence identity of their hypervariable region. In this paper, K-Ras4B and K-Ras notations are used interchangeably, both referring to the K-Ras4B splice variant, which is the focus of our study. K-Ras when bound to GDP is inactive, as it makes an unstable complex with downstream effectors such as RAF1.18–20 However, GTP-bound K-Ras can bind to downstream effectors, leading to signal transduction via pathways such as the Ras-Raf-MEK-ERK pathway.21,22 A major difference between GDP-bound K-Ras and GTP-bound K-Ras is in the regions switch I (residues 30–38) and switch II (residues 59–76) [Fig. 1(a)]. In GDP-bound K-Ras, switches are more flexible, and their distance from the GDP ligand is large. Experimental observations suggest that switch regions in the GTP-bound K-Ras structures are less flexible and tighter packed to the ligand. Computational studies of K-Ras have shown the existence of GDP-bound K-Ras-like conformations in the GTP-bound structures.3 This suggests that even when K-Ras is bound to GTP, it can still adopt conformations that are not favorable for RAF1 binding, suggesting existence of both active and inactive states for GTP-bound K-Ras. Here, we probe the detailed atomistic conformational dynamics of GTP-bound K-Ras4B proteins solvated with explicit water molecules, and the GTP-dependent activation/inactivation mechanism of the wild-type K-Ras4B(WT), and of the K-Ras4B(D33E) mutation. We note that the rather conservative D33E mutation leads to non-trivial large effects, in agreement with previous studies.23 

To model wild-type K-Ras4B bound to GTP [notation K-Ras(WT).GTP], the structure of human Q61H mutant in complex with the GTP analog GNP [PDB ID 3GFT, K-Ras(Q61H).GNP] was used, with residue 61 mutated to Gln and GNP replaced with GTP. The system was solvated using 16 274 explicit water molecules (TIP3P water model24). Na+ and Cl ions were used to neutralize the system and achieve physiological conditions (0.15 mol/l). All molecular dynamics (MD) simulations used the Nanoscale Molecular Dynamics (NAMD) package (2.13 multicore CUDA version)25 with the CHARMM36 force field, as described below and in our recent work on a different system of similar size (i.e., Abl kinase).26–28 Here, the system was subjected to potential energy minimization, followed by runs in the isothermal-isobaric (NPT, 5 ns long) and canonical ensembles (NVT, 20 ns long) using periodic boundary conditions at a temperature of 300 K with a 1 fs integration time step. For data collection, 40 MD trajectories, each 40 ns long, were launched from the S1, S4, and S6 initial conformations (see Fig. 1). S1 corresponds to the crystal structure and has rather small d1 and d2 values. The S4 and S6 states are already observed in the first 25 ns of MD runs. Overall, a total of 4.8 µs of simulation data for K-Ras4B(WT).GTP were used for the final analysis. These data are available to download, as mentioned at the end of the article.

The K-Ras4B(D33E).GTP model was prepared by mutating the K-Ras(WT).GTP prepared above. The same minimization and equilibration protocol was followed as explained for K-Ras(WT).GTP. For data collection, 14 trajectories 100 ns long each were launched from the high energy area between S2 and S4 of the WT free energy map [see Fig. 2(a)]. We adopted this improved sampling approach, compared with previous MD studies of K-Ras, to extract a detailed and fully converged landscape in (d1, d2) coordinates.3,12 A total of 1.4 µs of MD simulation data were used for analysis for the K-Ras4B(D33E).GTP system.

FIG. 2.

(a) The detailed K-Ras4B wild-type (WT) free energy landscape (ΔG, kcal/mol) for its GTP-bound structure in (d1, d2) coordinates [see Fig. 1(a)]. The six main conformational states of K-Ras4B WT are labeled S1 to S6, (yellow). (b) The corresponding ΔG landscape was calculated for the GTP-bound K-Ras4B D33E mutant, with the new conformational basins labeled S1’ to S6’. (c) The normalized distributions of values, P(d), in the case of K-Ras4B WT for d1 and the corresponding distances between the neighboring residues Q61 and E62 and the Pβ atom of GTP. (d) The corresponding P(d) distributions was calculated for the GTP-bound K-Ras4B D33E mutant.

FIG. 2.

(a) The detailed K-Ras4B wild-type (WT) free energy landscape (ΔG, kcal/mol) for its GTP-bound structure in (d1, d2) coordinates [see Fig. 1(a)]. The six main conformational states of K-Ras4B WT are labeled S1 to S6, (yellow). (b) The corresponding ΔG landscape was calculated for the GTP-bound K-Ras4B D33E mutant, with the new conformational basins labeled S1’ to S6’. (c) The normalized distributions of values, P(d), in the case of K-Ras4B WT for d1 and the corresponding distances between the neighboring residues Q61 and E62 and the Pβ atom of GTP. (d) The corresponding P(d) distributions was calculated for the GTP-bound K-Ras4B D33E mutant.

Close modal

In a second step, the PatchDock program29 was used in the molecular docking trials to study the relative compatibility of the conformational states revealed from our detailed MD-based free energy map with RAF1 binding and the underlying activation/inactivation propensity. The experimental RAF1 conformation22 (i.e., the Ras binding domain, RBD, and the cysteine-rich domain (CRD); PDB ID 6XI7) was used to dock to MD-generated K-Ras4B.GTP structures. A list of potential binding sites residues on RAF1 and Ras was specified as a constraint (see Figs. S4 and S5).22 This list includes residues from both the RBD and CRD domains of RAF1, as CRD increases the binding affinity to K-Ras and the CRD–K-Ras interaction is necessary for RAF1 activation and, thus, for downstream signaling.22,30–34 Unlike RBD, which primarily binds with the switch regions of K-Ras, CRD interacts with inter-switch region of K-Ras and with the C-terminal α-helix via nonbonded interactions and hydrogen bonds. The root mean square deviation (RMSD) was calculated for the binding interface (i.e., defined using Cα atoms on RAF1, which are within 5 Å of K-Ras atoms) for the top five high-score structures, generated by the docking program. RMSD was calculated with respect to the 6XI7 crystal structure,22 after aligning K-Ras as a reference. The smallest RMSD value of these five values in each case is reported in Table I. Visual Molecular Dynamics (VMD) was used for generating and analyzing representative structures.35 

TABLE I.

The results of docking-based modeling. RMSD values (in Å) obtained for comparing the docking interface of K-Ras4B and RAF1 obtained experimentally (PDB code 6XI7) with the dimer structures corresponding to the two peaks (denoted here by θ1 and θ2, respectively) of the angle θ38 (i.e., the C-Cα-Cβ-Cγ dihedral angle for residue D38, in degrees) used for measuring the relative orientation (w.r.t. the local backbone) of the D38 side chain in the K-Ras4B WT structure. See text for discussion.

K-Ras4B WTK-Ras4B D33E
 θ1 θ2  θ1 θ2 
S1 5.0 2.8 S1’ 19.4 2.0 
S2 5.1 20.1 S2’ 21.3 3.0 
S3 6.8 4.9  ⋯  
 ⋯  S3’ 20.1 11.6 
S4 21.4 5.4 S4’ 9.6 19.8 
S5 7.5 7.6  ⋯  
S6 20.8 3.1 S5’ 20.7 8.9 
 ⋯  S6’ 4.2 20.3 
K-Ras4B WTK-Ras4B D33E
 θ1 θ2  θ1 θ2 
S1 5.0 2.8 S1’ 19.4 2.0 
S2 5.1 20.1 S2’ 21.3 3.0 
S3 6.8 4.9  ⋯  
 ⋯  S3’ 20.1 11.6 
S4 21.4 5.4 S4’ 9.6 19.8 
S5 7.5 7.6  ⋯  
S6 20.8 3.1 S5’ 20.7 8.9 
 ⋯  S6’ 4.2 20.3 

Figure 2(b) shows the corresponding free energy surface calculated for the mutated system K-Ras4B(D33E).GTP showing significant changes compared to the WT case [Fig. 2(b)], as discussed below.

First, to study the mechanistic role of switch 1, we check the suitability of using the d2 distance (see Figs. 1 and 2) as reaction coordinate (RC). We thus compute the corresponding distances between Cα atoms of several residues in the 32–40 sequence range and the Pβ atom of GTP from our long MD trajectories (Fig. S2). The distributions of these distances for residues 32–34 and 36–40 clearly show a single peak (located at ∼6.5, 8.5, 8, 10.2, 12, 14, 15, and 15 Å, respectively, see Fig. S2). However, the distribution for d2 (residue T35) captures clearly the most structural states (i.e., three peaks at ∼6.5, 8.2, and 12 Å, see Fig. S2). Clearly, the d2 distance can capture more conformational metastable states compared with other options.

Similarly, to study switch II and thus check the suitability of using d1 as an RC, residues G60, Q61, and E62 located on switch 2 are selected, and the corresponding distance [defined as the distance between the Cα atom of G60 and the Pβ atom of GTP, dashed black line, Fig. 1(a)] is computed. The corresponding distance distributions for residues Q61 and E62 have only a single peak (at ∼10.5 and ∼13 Å, respectively, see Fig. 2). However, the distance distribution for residue G60 shows at least two clearly defined peaks at ∼6.5 and ∼8.5 Å [see Fig. 2(c)]. Once again, the choice of using G60 and, thus, the d1 distance as an RC appears to be the better choice.

Thus, for the force field used here, d1 and d2 proposed are reasonable RCs (Figs. 2 and S2), allowing the construction of a detailed underlying free energy landscape by measuring the population density in the (d1, d2) space [Fig. 2(a) for WT]. We identified six representative conformations as the centers of free energy minima (highest population density), illustrated in Fig. S1. These six centers of energy minima represent states labeled, S1, S2, S3, S4, S5, and S6 [Fig. 2(a)]. Our detailed free energy map suggests that GTP bound K-Ras (WT) has a higher resolution than previously described.3,12

We have also measured the d1 and d2 distances for various GTP analog/GTP-bound K-Ras and GDP analog/GDP-bound K-Ras crystal structures (see Fig. S3). These corresponding distances for crystal structures show that most GTP bound K-Ras WT and mutant structures cluster closer to the S1 state on our computational free-energy map.

To improve on the (d1, d2) description, we note that the structural analysis of the RAF1 and K-Ras complex from the crystal structure PDB ID 6XI722 shows that the acidic residues of K-Ras (e.g., residue D38) interact with the basic residues of RAF1 creating strong salt bridges. To measure the orientation of D38, we computed the distribution of values of the dihedral angle θ38 (i.e., defined as the C-Cα-Cβ-Cγ dihedral angle for residue D38, in degrees, see Fig. 3). For the 6XI7 crystal structure,22 dihedral angle for residue D38 is 151.63°. The D38 dihedral angle distribution for MD-generated trajectories presents two clear peaks (Fig. 3). Interestingly, the analysis of values of the dihedral angle, θ38, in all conformational states S1, S2, S3, S4, S5, and S6 shows the presence of both theta peaks in all six (d1, d2) conformational states.

FIG. 3.

(a) The positions of atoms defining the angle θ (i.e., the C-Cα-Cβ-Cγ dihedral angle for residue D38, in degrees) used for measuring the relative orientation (w.r.t. the local backbone) of the D38 side chain in the K-Ras4B WT structure. The GTP ligand is shown as licorice and the D38 atoms as balls-and-sticks. (b) The corresponding overall distribution of θ38 values for K-Ras4B WT. (c) Dihedral angle (θ38) distributions in each of the six states of K-Ras4B WT. (d) The corresponding dihedral angle (θ38) distributions in each of the six states of K-Ras4B D33E.

FIG. 3.

(a) The positions of atoms defining the angle θ (i.e., the C-Cα-Cβ-Cγ dihedral angle for residue D38, in degrees) used for measuring the relative orientation (w.r.t. the local backbone) of the D38 side chain in the K-Ras4B WT structure. The GTP ligand is shown as licorice and the D38 atoms as balls-and-sticks. (b) The corresponding overall distribution of θ38 values for K-Ras4B WT. (c) Dihedral angle (θ38) distributions in each of the six states of K-Ras4B WT. (d) The corresponding dihedral angle (θ38) distributions in each of the six states of K-Ras4B D33E.

Close modal

Based on these observations, which are in agreement with other studies that suggested that the orientations of interface residues are important descriptors of K-Ras-binding propensity to effector molecules,36,37 we hypothesize that for K-Ras4B GTP bound to be active, it should be able (i.e., have a high propensity) to bind to RAF1. To determine which possible combination of d1, d2, and θ38 triad would facilitate binding, we dock RAF1 conformations, extracted using the 6XI7 structure, to the peak 1 and peak 2 (i.e., of the θ38 distribution) representative structures for all six (S1 to S6) states. We can safely assume that if the high-scoring docked structure is close to the crystal structure orientation, then binding is facilitated, and thus, the K-Ras4B GTP conformation is active. The binding interface RMSD for the top five structures, generated by the docking program, was computed. RMSD values were calculated with respect to the 6XI7 crystal structure, after aligning K-Ras. A low RMSD value means closer to the crystal structure and a higher propensity of RAF1-Ras binding. The smallest of the five RMSD values, in each case, has been reported in Table I. Clearly, the S1 state conformation with the large dihedral angle (i.e., peak 2) binds to RAF1 better than S1 state conformation with the small dihedral angle (i.e., peak 1). Previous studies showed that if d1 and d2 distances are small, then K-Ras-GTP is active and binds to RAF1.3 Interestingly, our study finds that even for small values for d1 and d2, not all structures (including S1) are necessarily always active. We also note that, although key, the value of the new θ38 RC alone is not sufficient for establishing RAF1 binding propensity. For example, for states like S2, a large θ38 value does not result in a low RMSD for the binding interface compared with the experiment.

Importantly, from the docking calculations reported in Table I, we can conclude that structures with small d1 and d2, and large theta values would have the highest propensity to bind to RAF1. However, other RCs, besides d1, d2, and θ38 may be able to also correlate with the K-Ras(GTP)-RAF1 binding propensity. Interestingly, θ38 has the advantages over using the (d1, d2) combination alone as an RC that (i) it presents a largely bimodal distribution overall and (ii) one of its two modes (peak 2) is much more compatible with RAF1 binding, as it promotes interfacial salt bridges between the acidic D37 and D38 residues of K-Ras and the basic residues of RAF1such as R89.

Figure 4 illustrates important residues located at the binding interface between K-Ras4B and its RAF1 effector (from PDB ID 6XI7).22 Our analysis shows that a main stabilizing factor is the network of complementary acidic (red) and basic (blue) amino acids located at this interface [Fig. 4(a)]. Neutralizing mutations such as acting on D38 would have a destabilizing effect. However, our analysis of the θ38 angle shows that mutations, such as D33E, can in fact facilitate the K-Ras binding to RAF1 by allowing the D33 side chain (and thus, the switch 1 region) to interact more closely with RAF1.

FIG. 4.

The RAF1-RAS binding interface. (a) Crystal structure of K-Ras4B-GNP in complex with RAF1 (PDB ID 6XI7). The binding interface for RAF1 is shown in surface representation, colored by the residue type. Binding interface residues of K-Ras are shown as licorice. Acidic residues of K-Ras-like residue E37 and D38 (red, circled) form salt bridges with the basic residues (blue) on the RAF1 surface. (b) RAF1 (cyan) docked to the S2’ for the theta peak 1 structure and RAF1 (purple) docked to the S2’ peak 2. The two structures share only a small part of the binding interface involving D38 (red, circled in black). The RAF1 bound to peak 2 (cyan) interface is closer to the RAF1 crystal structure (see Table I).

FIG. 4.

The RAF1-RAS binding interface. (a) Crystal structure of K-Ras4B-GNP in complex with RAF1 (PDB ID 6XI7). The binding interface for RAF1 is shown in surface representation, colored by the residue type. Binding interface residues of K-Ras are shown as licorice. Acidic residues of K-Ras-like residue E37 and D38 (red, circled) form salt bridges with the basic residues (blue) on the RAF1 surface. (b) RAF1 (cyan) docked to the S2’ for the theta peak 1 structure and RAF1 (purple) docked to the S2’ peak 2. The two structures share only a small part of the binding interface involving D38 (red, circled in black). The RAF1 bound to peak 2 (cyan) interface is closer to the RAF1 crystal structure (see Table I).

Close modal

Finally, we have used all the MD trajectory data generated in this study to build the corresponding network of transitions between the coarse Markovian conformational basins corresponding to the states S1 to S6 illustrated in Fig. 2(a) for the GTP-bound KRas-4B WT system. Here, we used a similar approach as described in our recent studies on different proteins38,39 to calculate the equilibrium populations of these six coarse states and the corresponding transition rates between them as illustrated in Fig. S6. We note that a similar analysis based on Markov State Models has been published recently on the KRas-4B G12C mutant (for unveiling the effect of binding two different inhibitors)40 and also on the GDP-bound KRas-4B WT [for probing the role of its the C-terminal hypervariable region (HVR), in activation and membrane anchoring].41 Although both these studies focus on different systems and choose to describe the Markovian states using other reaction coordinates (i.e., RMSD of loops I and II,40 or principal components41), they unveil a similar level of complexity (i.e., in the overall number of Markovian states unveiled) of the underlying free energy landscape of K-Ras4B. We are currently pursuing a detailed, computationally intensive study using Markov models to rationalize systematically the effect of mutations on K-Ras4B’s activation/inactivation propensities.

To sum up, we have generated and used a new set of long-time MD trajectories, initialized from various regions of the (d1, d2) conformational space of K-Ras4B to extract high resolution, converged free energy maps of both the K-Ras4B WT protein and its K-Ras4B WT D33E mutant [Figs. 2(a) and 2(b)].

For the K-Ras4B WT molecule, the new map showcases six distinct (d1, d2) states that, however, cannot discriminate alone between active and inactive molecular conformations of K-Ras4B WT and mutated structures, demonstrating the need for additional reaction coordinates. The basins of the new ΔG map in (d1, d2) coordinates correlate well with experimental information available from the crystal structures that are either likely to be active (i.e., GTP or GTP-analog bound) or known as inactive (GDP or GDP-analog bound) as shown in Fig. S3. This observation is not unexpected as the d1, and d2 RCs were introduced to capture active-like conformations in regions such as S1 and S2. However, we find that some inactive-like conformations can also correspond to small values of (d1, d2) coordinates (e.g., labeled 8 in Fig. S3), whereas some active-like cases (e.g., labeled 7 and 9 in Fig. S3) correspond to rather large (d1, d2) values. To address this issue, we showed that the dihedral angle, θ38, of an acidic interface residue, D38, presents a bimodal distribution in all the six main conformational states evidenced in the (d1, d2) coordinates, with conformations corresponding to its two peaks discriminating well between active and inactive conformations (see Fig. S3b). This finding is also supported by our docking calculations showing that both high scores and small RMSD values for the docking interface are obtained by S1 and S2 structures that are docked to RAF1. Interestingly, including the θ38 conformations in the analysis showed that experimental-like (i.e., small RMSD values from the experimental interface structure in 6XI7, in the range of ∼2–3 Å) can be obtained for states S1 and S2 if the θ38 conformations are active-like conformations (i.e., peak 2, Fig. 3).

These observations lead us to suggest a new paradigm in analyzing the effect of K-Ras (WT) conformations and its mutant for inferring propensities of activation/inactivation, namely that a more complex collective variable such as obtained by monitoring all three RCs discussed (d1, d2, and θ38) simultaneously, is necessary. This ability would be particularly useful in quantifying the effects of mutations on K-Ras propensities of activation/inactivation. As shown here, using the (d1, d2, and θ38) triad, we can explain the otherwise puzzling observation that a rather conservative mutation, D33E, can have a rather dramatic effect by significantly increasing the K-Ras activation propensities for its S1 and S2 confirmational states, in agreement with previous studies.23 

Here, the (d1, d2, and θ38) triad provides a new and more detailed mechanistic RC and unveils an explanation of how Ras uses the same binding site to engage with multiple effectors forming diverse binding interfaces through modulation of intermolecular interactions (e.g., salt bridges). Indeed, binding affinities between Ras and effectors vary a lot,20 as do contributions of individual amino acid contacts (e.g., “hot-spots”) at the different Ras-effector interfaces.42,43 Our analysis also paves the way for a better mechanistic understanding of Ras oncogenic mutations that are suggested to differentially impact binding to effectors.20,44 Altogether, our hybrid MD-docking modeling approach enables the development of new in silico methods for quantitative assessment of activation propensity changes (e.g., due to mutations or local binding environment). It also unveils the underlying molecular mechanisms (e.g., in this case, the modulation of the network of salt bridges at the binding interface with the RAF1 downstream effector) and facilitates the rational design of new cancer drugs.

The supplementary material consists of a PDF file including six supplementary figures (Figs. S1 to S6).

We wish to thank the DJEI/DES/SFI/HEA Irish Center for High-End Computing (ICHEC) and the Sonic HPC cluster (UCD Research IT) for the provision of computational facilities and support. We acknowledge financial support by Science Foundation Ireland under Project No. 16/FRL/3886 (to C.K.). B.N. and N.-V.B. acknowledge the financial support received from the Thomas Preston Ph.D. Scholarship Fund (University College Dublin) and from the European Union’s Horizon 2020 research and innovation program “NanoinformaTIX” (H2020-NMBP-14-2018, Grant No. 814426).

We also wish to thank Walter Kolch for helpful discussions, and Shaoyong Lu, Hyunbum Jang, Ruth Nussinov, and Jian Zhang for kindly providing PDB structures that represent active and inactive KRas-4B conformations found in their previous simulations,45 which helped us with the initial conceptualization of our study.

The authors have no conflicts to disclose.

Brajesh Narayan: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Christina Kiel: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Nicolae-Viorel Buchete: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Methodology (lead); Project administration (equal); Resources (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead).

All molecular dynamics and molecular docking data presented, analyzed, and discussed in this paper is available upon reasonable requests submitted to the corresponding author. The following data is also available online at https://zenodo.org/record/7672195: (i) file KRas4B_pdbs.zip, a compressed archive including coordinate (.pdb) and structure (.psf) files for KRas-4B WT and D33E proteins, (ii) KRas4B_WT_traj.zip, a compressed archive including trajectory files (.trr) for the WT KRas-4B runs, including 120 trr trajectories each corresponding to 40 ns of MD simulation time, and (iii) a sample Python script to generate a free energy plot as shown in Fig. 2.

1.
N.
Tsuchida
,
T.
Ryder
, and
E.
Ohtsubo
, “
Nucleotide sequence of the oncogene encoding the p21 transforming protein of Kirsten murine sarcoma virus
,”
Science
217
,
937
939
(
1982
).
2.
A.
Fernández-Medarde
,
J.
De Las Rivas
, and
E.
Santos
, “
40 Years of RAS-A historic overview
,”
Genes
12
,
681
(
2021
).
3.
S.
Lu
,
H.
Jang
,
S.
Muratcioglu
,
A.
Gursoy
,
O.
Keskin
,
R.
Nussinov
, and
J.
Zhang
, “
Ras conformational ensembles, allostery, and signaling
,”
Chem. Rev.
116
,
6607
6665
(
2016
).
4.
A.
Fernandez-Medarde
and
E.
Santos
, “
Ras in cancer and developmental diseases
,”
Genes Cancer
2
,
344
358
(
2011
).
5.
W.
Kolch
,
D.
Berta
, and
E.
Rosta
, “
Dynamic regulation of RAS and RAS signaling
,”
Biochem. J.
480
,
1
23
(
2023
).
6.
A.
Kapoor
and
A.
Travesset
, “
Differential dynamics of RAS isoforms in GDP- and GTP-bound states
,”
Proteins: Struct., Funct., Bioinf.
83
,
1091
1106
(
2015
).
7.
M.
Cefalì
,
S.
Epistolio
,
M. C.
Palmarocchi
,
M.
Frattini
, and
S.
De Dosso
, “
Research progress on KRAS mutations in colorectal cancer
,”
J. Cancer Metastasis Treat.
7
,
26
(
2021
).
8.
N.
Umetani
,
S.
Sasaki
,
T.
Masaki
,
T.
Watanabe
,
K.
Matsuda
, and
T.
Muto
, “
Involvement of APC and K-ras mutation in non-polypoid colorectal tumorigenesis
,”
Br. J. Cancer
82
,
9
15
(
2000
).
9.
M.
Ando
,
M.
Maruyama
,
M.
Oto
,
K.
Takemura
,
M.
Endo
, and
Y.
Yuasa
, “
Higher frequency of point mutations in the c-K-ras2 gene in human colorectal adenomas with severe atypia than in carcinomas
,”
Jpn. J. Cancer Res.
82
,
245
249
(
1991
).
10.
M.
Gerber
,
S.
Goel
, and
R.
Maitra
, “
In silico comparative analysis of KRAS mutations at codons 12 and 13: Structural modifications of P-Loop, switch I&II regions preventing GTP hydrolysis
,”
Comput. Biol. Med.
141
,
105110
(
2022
).
11.
E.
Kempf
,
B.
Rousseau
,
B.
Besse
, and
L.
Paz-Ares
, “
KRAS oncogene in lung cancer: Focus on molecularly driven clinical trials
,”
Eur. Respir. Rev.
25
,
71
76
(
2016
).
12.
X.
Wang
, “
Conformational fluctuations in GTP-bound K-Ras: A metadynamics perspective with harmonic linear discriminant analysis
,”
J. Chem. Inf. Model.
61
,
5212
5222
(
2021
).
13.
M.
Eren
,
N.
Tuncbag
,
H.
Jang
,
R.
Nussinov
,
A.
Gursoy
, and
O.
Keskin
, “
Normal mode analysis of KRas4B reveals partner specific dynamics
,”
J. Phys. Chem. B
125
,
5210
5221
(
2021
).
14.
S. M.
Park
,
C. Y.
Hwang
,
S. H.
Cho
,
D.
Lee
,
J. R.
Gong
,
S.
Lee
,
S.
Nam
, and
K. H.
Cho
, “
Systems analysis identifies potential target genes to overcome cetuximab resistance in colorectal cancer cells
,”
FEBS J.
286
,
1305
1318
(
2019
).
15.
S.
Lu
,
H.
Jang
,
J.
Zhang
, and
R.
Nussinov
, “
Inhibitors of Ras-SOS interactions
,”
ChemMedChem
11
,
814
821
(
2016
).
16.
S.
Catozzi
,
C.
Ternet
,
A.
Gourrege
,
K.
Wynne
,
G.
Oliviero
, and
C.
Kiel
, “
Reconstruction and analysis of a large-scale binary Ras-effector signaling network
,”
Cell Commun. Signaling
20
,
24
(
2022
).
17.
D. K.
Simanshu
,
D. V.
Nissley
, and
F.
McCormick
, “
RAS proteins and their regulators in human disease
,”
Cell
170
,
17
33
(
2017
).
18.
C.
Kiel
,
D.
Filchtinski
,
M.
Spoerner
,
G.
Schreiber
,
H. R.
Kalbitzer
, and
C.
Herrmann
, “
Improved binding of raf to Ras.GDP is correlated with biological activity
,”
J. Biol. Chem.
284
,
31893
31902
(
2009
).
19.
M.
Ilter
,
R.
Kaşmer
,
F.
Jalalypour
,
C.
Atilgan
,
O.
Topcu
,
N.
Karakaş
 et al “
Inhibition of mutant RAS-RAF interaction by mimicking structuraland dynamic properties of phosphorylated RAS
,”
eLife
11
,
e79747
(
2022
).
20.
C.
Kiel
,
D.
Matallanas
, and
W.
Kolch
, “
The ins and outs of RAS effector complexes
,”
Biomolecules
11
,
236
(
2021
).
21.
L.
Huang
,
Z.
Guo
,
F.
Wang
, and
L.
Fu
, “
KRAS mutation: From undruggable to druggable in cancer
,”
Signal Transduction Targeted Ther.
6
,
386
(
2021
).
22.
T. H.
Tran
,
A. H.
Chan
,
L. C.
Young
,
L.
Bindu
,
C.
Neale
,
S.
Messing
,
S.
Dharmaiah
,
T.
Taylor
,
J.-P.
Denson
,
D.
Esposito
 et al, “
KRAS interaction with RAF1 RAS-binding domain and cysteine-rich domain provides insights into RAS-mediated RAF activation
,”
Nat. Commun.
12
,
1176
(
2021
).
23.
J.
Chen
,
S.
Zhang
,
W.
Wang
,
L.
Pang
,
Q.
Zhang
, and
X.
Liu
, “
Mutation-induced impacts on the switch transformations of the GDP- and GTP-bound K-Ras: Insights from multiple replica Gaussian accelerated molecular dynamics and free energy analysis
,”
J. Chem. Inf. Model.
61
,
1954
1969
(
2021
).
24.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
,
926
935
(
1983
).
25.
J. C.
Phillips
,
R.
Braun
,
W.
Wang
,
J.
Gumbart
,
E.
Tajkhorshid
,
E.
Villa
,
C.
Chipot
,
R. D.
Skeel
,
L.
Kalé
, and
K.
Schulten
, “
Scalable molecular dynamics with NAMD
,”
J. Comput. Chem.
26
,
1781
1802
(
2005
).
26.
B.
Narayan
,
A.
Fathizadeh
,
C.
Templeton
,
P.
He
,
S.
Arasteh
,
R.
Elber
,
N.-V.
Buchete
, and
R. M.
Levy
, “
The transition between active and inactive conformations of Abl kinase studied by rock climbing and Milestoning
,”
Biochim. Biophys. Acta, Gen. Subj.
1864
,
129508
(
2020
).
27.
B.
Narayan
,
Y.
Yuan
,
A.
Fathizadeh
,
R.
Elber
, and
N.-V.
Buchete
, “
Chapter five—Long-time methods for molecular dynamics simulations: Markov State Models and Milestoning
,” in
Progress in Molecular Biology and Translational Science
, edited by
B.
Strodel
and
B.
Barz
(
Academic Press
,
2020
), pp.
215
237
.
28.
B.
Narayan
,
N.-V.
Buchete
, and
R.
Elber
, “
Computer simulations of the dissociation mechanism of gleevec from Abl kinase with milestoning
,”
J. Phys. Chem. B
125
,
5706
5715
(
2021
).
29.
D.
Schneidman-Duhovny
,
Y.
Inbar
,
R.
Nussinov
, and
H. J.
Wolfson
, “
PatchDock and SymmDock: Servers for rigid and symmetric docking
,”
Nucleic Acids Res.
33
,
W363
W367
(
2005
).
30.
R.
Thapar
,
J. G.
Williams
, and
S. L.
Campbell
, “
NMR characterization of full-length farnesylated and non-farnesylated H-Ras and its implications for Raf activation
,”
J. Mol. Biol.
343
,
1391
1408
(
2004
).
31.
J. G.
Williams
,
J. K.
Drugan
,
G.-S.
Yi
,
G. J.
Clark
,
C. J.
Der
, and
S. L.
Campbell
, “
Elucidation of binding determinants and functional consequences of Ras/Raf-cysteine-rich domain interactions
,”
J. Biol. Chem.
275
,
22172
22179
(
2000
).
32.
T.
Bondeva
,
A.
Balla
,
P.
Várnai
, and
T.
Balla
, “
Structural determinants of Ras-Raf interaction analyzed in live cells
,”
Mol. Biol. Cell
13
,
2323
2333
(
2002
).
33.
Z.-L.
Li
,
P.
Prakash
, and
M.
Buck
, “
A ‘Tug of War’ maintains a dynamic protein-membrane complex: Molecular dynamics simulations of C-Raf RBD-CRD bound to K-Ras4B at an anionic membrane
,”
ACS Cent. Sci.
4
,
298
305
(
2018
).
34.
V. P.
Mysore
,
Z.-W.
Zhou
,
C.
Ambrogio
,
L.
Li
,
J. N.
Kapp
,
C.
Lu
,
Q.
Wang
,
M. R.
Tucker
,
J. J.
Okoro
,
G.
Nagy-Davidescu
 et al, “
A structural model of a Ras-Raf signalosome
,”
Nat. Struct. Mol. Biol.
28
,
847
857
(
2021
).
35.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD: Visual molecular dynamics
,”
J. Mol. Graphics
14
,
33
38
(
1996
).
36.
Y.
Li
,
Y.
Zhang
,
F.
Großerüschkamp
,
S.
Stephan
,
Q.
Cui
,
C.
Kötting
,
F.
Xia
, and
K.
Gerwert
, “
Specific substates of Ras to interact with GAPs and effectors: Revealed by theoretical simulations and FTIR experiments
,”
J. Phys. Chem. Lett.
9
,
1312
1317
(
2018
).
37.
S.
Vatansever
,
B.
Erman
, and
Z. H.
Gümüş
, “
Oncogenic G12D mutation alters local conformations and dynamics of K-Ras
,”
Sci. Rep.
9
,
11730
(
2019
).
38.
C. T.
Leahy
,
R. D.
Murphy
,
G.
Hummer
,
E.
Rosta
, and
N.-V.
Buchete
, “
Coarse master equations for binding kinetics of amyloid peptide dimers
,”
J. Phys. Chem. Lett.
7
,
2676
2682
(
2016
).
39.
B.
Narayan
,
C.
Herbert
,
B. J.
Rodriguez
,
B. R.
Brooks
, and
N.-V.
Buchete
, “
Replica exchange molecular dynamics of diphenylalanine amyloid peptides in electric fields
,”
J. Phys. Chem. B
125
,
5233
5242
(
2021
).
40.
G.
Tu
,
Q.
Liu
,
Y.
Qiu
,
E. L.
Leung
, and
X.
Yao
, “
In silico study of the acquired resistance caused by the secondary mutations of KRAS G12C protein using long time molecular dynamics simulation and Markov state model analysis
,”
Int. J. Mol. Sci.
23
,
13845
(
2022
).
41.
H.
Zhang
,
D.
Ni
,
J.
Fan
,
M.
Li
,
J.
Zhang
,
C.
Hua
,
R.
Nussinov
, and
S.
Lu
, “
Markov state models and molecular dynamics simulations reveal the conformational transition of the intrinsically disordered hypervariable region of K-Ras4B to the ordered conformation
,”
J. Chem. Inf. Model.
62
,
4222
4231
(
2022
).
42.
D.
Berta
,
S.
Gehrke
,
K.
Ny̚ri
,
B. G.
Vértessy
, and
E.
Rosta
,“
Redesigned GAP to activate oncogenic GAP
,”
ChemRxiv
(
2023
).
43.
C.
Kiel
,
L.
Serrano
, and
C.
Herrmann
, “
A detailed thermodynamic analysis of ras/effector complex interfaces
,”
J. Mol. Biol.
340
,
1039
1058
(
2004
).
44.
T. L.
Yuan
,
A.
Amzallag
,
R.
Bagni
,
M.
Yi
,
S.
Afghani
,
W.
Burgan
,
N.
Fer
,
L. A.
Strathern
,
K.
Powell
,
B.
Smith
 et al, “
Differential effector engagement by oncogenic KRAS
,”
Cell Rep.
22
,
1889
1902
(
2018
).
45.
S.
Lu
,
H.
Jang
,
R.
Nussinov
, and
J.
Zhang
, “
The structural basis of oncogenic mutations G12, G13 and Q61 in small GTPase K-Ras4B
,”
Sci. Rep.
6
,
21949
(
2016
).

Supplementary Material