
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12192-5
10.1016/j.heliyon.2024.e36161
e36161
Research Article
Insight into structural dynamics involved in activation mechanism of full length KRAS wild type and P-loop mutants
Jani Vinod
Sonavane Uddhavesh uddhaveshs@cdac.in
⁎
Joshi Rajendra
Centre for Development of Advanced Computing (C-DAC), Panchavati, Pashan, Pune, India
⁎ Corresponding author. uddhaveshs@cdac.in
13 8 2024
30 8 2024
13 8 2024
10 16 e3616111 4 2023
6 8 2024
11 8 2024
© 2024 The Authors. Published by Elsevier Ltd.
2024

https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
KRAS protein is known to be frequently mutated in various cancers. The most common mutations being at position 12, 13 and 61. The positions 12 and 13 form part of the phosphate binding region (P-loop) of KRAS. Owing to mutation, the protein remains in continuous active state and affects the normal cellular process. Understanding the structural changes owing to mutations in GDP-bound (inactive state) and GTP-bound (active state) may help in the design of better therapeutics. To understand the structural flexibility due to the mutations specifically located at P-loop regions (G12D, G12V and G13D), extensive molecular dynamics simulations (24 μs) have been carried for both inactive (GDP-bound) and active (GTP-bound) structures for the wild type and these mutants. The study revealed that the local structural changes at the site of mutations allosterically guide changes in distant regions of the protein through hydrogen bond and hydrophobic signalling network. The dynamic cross correlation analysis and the comparison of the correlated motions among different systems manifested that changes in SW-I, SW-II, α3 and the loop preceding α3 affects the interactions of GDP/GTP with different regions of the protein thereby affecting its hydrolysis. Further, the Markov state modelling analysis confirmed that the mutations, especially G13D imparts rigidity to structure compared to wild type and thus limiting its conformational state in either intermediate state or active state. The study suggests that along with SW-I and SW-II regions, the loop region preceding the α3 helix and α3 helix are also involved in affecting the hydrolysis of nucleotides and may be considered while designing therapeutics against KRAS.

Keywords

KRAS
Molecular dynamics simulation
Markov state model
DCCM
Cancer
==== Body
pmc1 Introduction

The RAt Sarcoma (RAS), a GTPase is known to regulate various cellular processes like cell proliferation, cell development and differentiation, migration and apoptosis [1,2]. It has an intrinsic function as a molecular switch and transits between GTP-bound active state and GDP-bound inactive state [[3], [4], [5]]. The activation process involves transformation of GDP-bound form with GTP-bound form and this process is accelerated with the help of guanine-nucleotide exchange factor (GEF) [6,7]. The inactivation process is guided by the GTPase activating protein (GAP) which leads to transformation from GTP-bound form with GDP-bound form [6,7]. HRAS, KRAS, and NRAS are three RAS isoforms with overall high sequence similarity [[8], [9], [10]]. Within the three RAS isoforms, KRAS is the known to be most frequently mutated in different types of cancers [[9], [10], [11], [12], [13]]. The patients with the pancreatic cancer, colorectal tumours and lung cancer mostly harbours mutations in the KRAS protein as compared to HRAS and NRAS [12].

The KRAS protein comprises a G-domain and hyper variable C-terminus domain. The initial 1–166 residues form the G-domain and comprises six beta-strands and five alpha-helices [14]. The residue 167–188 forms the hyper variable region (HVR) which helps KRAS to anchor to the membrane [[15], [16], [17]]. The details of different regions of the KRAS protein is given in Fig. 1. The important functional regions of the KRAS includes switch-I (residue 30–38, SW-I), switch-II (residue 59–76, SW-II), phosphate binding loop (residues 10–17, P-loop), and nucleobase binding loops (residues 116–120 and 145–147). These regions play an important role in the activation mechanism. In most cancers the mutation maps to either the residues of P-loop or SW-II [18]. Across the GTPase superfamily the nucleotide-binding pocket is conserved and for RAS it comprises of P-loop, NKxD motif (residue 116 to 119), and ExSAK motif (residue 143 to 147) [19,20].Fig. 1 Structure of wild-type (WT) KRAS depicting different secondary structure elements and their sequence.

Fig. 1

Out of the 44 different point mutations that are known to occur in the KRAS protein, around 99 % of them are located at the codons 12, 13 and 61 [21]. These mutations affect intrinsic as well as GAP associated GTP hydrolysis leading to the continuous active state of KRAS leading to unregulated cell growth [[22], [23], [24], [25], [26], [27], [28], [29]]. It has been also observed by researchers that different mutations in the KRAS protein distinctly alter its GPTase activity or its affinity for downstream effector proteins [[30], [31], [32], [33]]. Several experimental and theoretical studies have been carried out to study the kinetics of GTP hydrolysis for the RAS protein. Most of the theoretical studies have been carried out on mutants of the RAS. As these mutations are known to impair the affinity of RAS for the GTP, researchers are trying to understand the structural changes owing to these mutations [[32], [33], [34], [35]]. Majority of studies have reported that SW-I and SW-II regions show major conformational changes during GTP hydrolysis [[32], [33], [34], [35]]. The study by Diaz et al. showed that the SW-II domain shows transition during GTP hydrolysis [36]. The study states that the removal of γ phosphate creates more space which allows the movement of residue Gly12 and Gly13, thereby allowing SW-II to adopt alternate conformations [36,37]. The studies on various mutants has helped in exploring various intermediates in the activation process [[38], [39], [40]].

In the NMR study, researchers have identified that SW-I tends to occur in two distinct conformations i.e in the open conformation (state 1) and in the closed conformation (state 2) [41]. The closed conformation is similar to that found in the RAS effector complex [[42], [43], [44]]. Initially the above finding was observed for HRAS, but later it was found in KRAS as well [45]. The open conformation (state 1) is characterized by the residues Tyr32 and Thr35 positioning away from the nucleotide and the Thr35 residue does not form contact with γ phosphate or Mg2+. The closed conformation (state II) is characterized by the Tyr32 residue going close towards the nucleotide and Thr35 sidechain forming contact with γ phosphate and Mg2+ [46,47]. Further, state II is known to exist in the two substates i.e R state and T state. The R state has a fully ordered active site and T state has disordered SW-II region [48]. Also, the water is known to play an important role in the hydrolysis of GTP and has been demonstrated by various researchers using classical molecular dynamics simulations [[49], [50], [51]].

Most of the studies in this regard to understand the effect of mutations have been done on HRAS isoform. Molecular dynamics study by Gorfe et al. has shown that the KRAS is more flexible as compared to HRAS and NRAS isoforms [26]. Also, Johnson et al. has shown that the three isoforms KRAS, HRAS and NRAS are biochemically distinct [52]. Since, the KRAS is mostly mutated in different cancer types, it becomes important to understand the structural changes in the KRAS owing to mutation in its GDP and GTP-bound states. The study by Lu et al. [53] on the full length KRAS oncogenic mutations showed that the mutants G12C, G12V and Q61H lead to weakened association between hypervariable region and G domain. The dynamics of only G domains of various oncogenic mutants viz. G12D, G12V and G13D for H-RAS and KRAS and heir comparative study revealed that the G12V and G13D mutants of H-Ras protein are significantly more flexible than their corresponding KRAS counterpart [54]. There have been studies on the conformational flexibility of the structure owing to various mutations [33,[54], [55], [56]], still there is scope to understand the local changes and how these changes allosterically regulate different regions of the protein. Vatansever et al. [56] carried out simulation studies to explore local events that affect the function of KRAS protein owing to G12D mutations. The study by Grudzien et al. showed that different mutants at residue position 12, namely G12D and G12V showed variation in conformational dynamics suggesting divergent activation mechanism by distinct mutants [57].

Various studies have suggested that the different mutants affect different downstream signaling pathways [[58], [59], [60], [61]]. The study by Ihle et al. [59] showed that in the non-small cell lung cancer (NSCLC) cell lines, owing to the G12C and G12D mutations the Ral signaling and growth factor-dependent Akt activation is affected. They also observed that the G12D mutation activates phosphatidylinositol 3-kinase (PI-3-K) and mitogen-activated protein/extracellular signal-regulated kinase (MEK) signalling. Hammond et al. [60] in their study showed that each of these mutations G12V, G12D and G13D have distinct effects on downstream signaling pathways. The review by Abdelkarim et al. suggest that the mutation at position 12 of KRAS may affect its interaction with calmodulin [61]. These studies clearly show that the mutation not only affects the region responsible for controlling the activation/deactivation of KRAS but also affects the regions responsible for interactions with downstream proteins.

Hence in the current work, molecular dynamics simulations have been performed for KRAS wild type, and three of its mutants i.e G12D, G12V and G13D in GDP and GTP-bound form. The study has been performed to explore the various conformational states of KRAS in GDP and GTP bound forms in order to get structural insight into the wild type and mutant protein due to mutations in the P-loop region.

2 Methodology

2.1 System preparation

The starting structures for the system preparation were downloaded from RCSB Protein Data Bank (PDB) with PDB id: 4DSU (GDP-bound) and PDB id: 4DSO (GSP bound) [62]. As most of the structures for KRAS have a truncated C-terminal end, hence above mentioned pdbs were considered as it had the coordinates for residues 1–180. The pdb structures had mutation at the position 12 i.e G12D, hence for obtaining wild type system (WT) residue Asp at position 12 was replaced by Gly using Chimera tool [63]. To obtain the G12V system, the residue position 12 was replaced by Val. Further the system G13D was obtained by replacing Asp at position 12 by Gly and Gly at position 13 with Asp. Also, to obtain GTP-bound systems, GSP was replaced by GTP.

2.2 Molecular dynamics simulations

The Gromacs v5x version was used to carry out for the simulations of all GDP and GTP bound KRAS systems [64]. The protein was parameterized using Amber99SB force field [65]. The parameters for the GTP and GDP were generated using prodrg2 server [66]. For each of the system of GDP-bound (WT, G12D, G12V and G13D) and of GTP-bound (WT, G12D, G12V and G13D), four sets of simulations (two 1 μs and two 500ns) were carried out. Details of the simulation lengths for each system has been given in Table 1. The systems were solvated in a cubic water box with TIP3P water model with distance of 10 Å between the closest protein fragment and edges of the box. The systems were neutralized by adding Na+ and Cl-ions. Each of the systems were energy minimized for around 5000 steps followed by two steps of equilibration. Initially using V-rescale thermostat [67], position restrained NVT equilibration was carried out for around 500 ps. This was followed by NPT equilibration for around 1 ns using Berendsen barostat [68]. The Particle-Mesh Ewald (PME) method was employed for calculation of long-range electrostatic calculations [69]. For the van der Waals's interactions and Coulombic interactions, a cut-off of 12 Å was set. The hydrogen bonds were constrained using the LINC algorithm [70]. The integration step of 2 fs was set for solving newton equations of motion. All the production runs were carried out using Parrinello-Rahman barostat [71].Table 1 Detail of molecular simulations of WT and Mutant systems.

Table 1S.No	System	Number of simulations	Time	
1	WT(GDP)	4 sets	Two sets of 500ns and two sets of 1μs	
2	G12D(GDP)	4 sets	Two sets of 500ns and two sets of 1μs	
3	G12V(GDP)	4 sets	Two sets of 500ns and two sets of 1μs	
4	G13D(GDP)	4 sets	Two sets of 500ns and two sets of 1μs	
5	WT(GTP)	4 sets	Two sets of 500ns and two sets of 1μs	
6	G12D(GTP)	4 sets	Two sets of 500ns and two sets of 1μs	
7	G12V(GTP)	4 sets	Two sets of 500ns and two sets of 1μs	
8	G13D(GTP)	4 sets	Two sets of 500ns and two sets of 1μs	

2.3 Analysis

The gromacs analysis utilities and amber cpptraj were used for the analysis of the all simulation trajectories. The interaction analysis was carried out using GetContacts tools [72]. The visualization and rendering of molecules were done using VMD [73] and Chimera.

2.4 Dynamic cross-correlation analysis

The dynamic cross correlation (DCC) analysis was carried out using the cpptraj module of AmberTools. It gives an idea about correlation of motions between different atoms of the system. equation (1) is used for calculation of the DCC between the ith and jth atoms of the protein.Eq. 1 Cij=⟨(ri−⟨ri⟩)(rj−⟨rj⟩)⟩(⟨ri2⟩−⟨ri⟩2)(⟨rj2⟩−⟨rj⟩2)

The Cij represents the cross correlation between the atom i and j where the ri and rj represent the coordinates at a given time point and the <ri> and <rj> represents average over the trajectory [74]. In-order to calculate the DCC, all the conformations are overlaid over the starting structure and atom-wise cross-correlation is calculated. The values for the Cij may vary between −1 and 1. The –ve values indicate the negative correlation i.e asynchronous motion between atoms and +ve values indicate the positive correlation i.e synchronous motion. Calculation of difference between correlated motions between the different systems and plotting them was done by using python.

2.5 Markov state models analysis (MSM)

Pyemma software package was used for carrying out MSM analysis [75]. MSM analysis helps in identifying important conformational states [76]. It also helps in determining stable intermediates states and the kinetics between these intermediates [[76], [77], [78]]. The first step in the MSM analysis is the choice of the collective variable (CV) i.e the features which can help to differentiate conformational states of the system [78,79]. In the current study, residue min_dist was used as featurizer. Following set of residues viz. 12, 29, 32, 34, 35, 36, 40, 48, 59, 60, 61, 62, 67, 105, 122, 126, 138 and GDP/GTP were considered for residue min_dist calculation. The K-means clustering algorithm with the number of initial clusters equal to 100 was chosen. From these clusters (i.e. microstates), macrostates are created using PCCA + algorithm. Finally, the transition matrix between these macrostates is computed in such a way that it follows the Markovian behaviour. The Chapman Kolmogorov test given by equation (2) is used for validation of the MSM [79].Eq. 2 T(nτ)=T(τ)n

where n is an integer number of steps, T(τ) is transition matrix at lag time τ.

3 Results and discussion

In current work, 24μs-long molecular dynamics simulations of KRAS GDP-bound wild type and three mutant variants and KRAS GTP-bound wild type and three mutant variants were carried out to study structural and dynamic features, which may be responsible for continuous activation state of mutants. Structural changes and dynamics of the systems are studied by dynamic cross correlation, interaction analysis and other structural analysis like MSM.

3.1 Root mean square deviation (RMSD)

The structural comparison with respect to start structure was carried out for each of the systems. Fig. 2(A–D) shows the RMSD distribution for GDP-bound systems. Fig. 2A shows RMSD distribution for the overall structure for GDP-bound systems, and it can be clearly seen that the mutation adds rigidity to the structure. The RMSD for WT ranges from 3 to 6 Å, while for the mutants it ranges from 3 to 4.5 Å. Further to observe the impact of mutations on important regions of KRAS viz. switch I (SW–I), switch II (SW-II) and P-loop regions, the RMSD distribution was plotted for all these regions (Fig. 2B–D). From the figure, it can be seen that the WT type shows major deviation in the SW-I region as compared to the mutants. The P-loop regions show deviation in the mutants system. The SW-II region showed fluctuation in all the systems with RMSD ranging from 2 to 5.5 Å.Fig. 2 Population distribution plot for root mean square deviation for GDP/GTP-bound systems for wild type (black color) and G12D (red color), G12V (green color) and G13D (blue color). Figure A to D are for GDP-bound system A) Complete protein. B) SW-I region. C) SW-II region. D) P-loop region. Figures E to H are for GTP-bound system E) Complete protein. F) SW-I region. G) SW-II region. H) P-loop region.

Fig. 2

Fig. 2E shows the RMSD distribution for the overall structure for GTP-bound systems. In GTP-bound systems, it was observed that mutant systems tend to show slightly higher fluctuations as compared to the WT. The average RMSD of mutant systems showed rise by 0.3, 0.18 and 0.4 Å as compared to WT for G12D, G12V and G13D systems respectively (Fig. 2E). Further it was observed that the G12D and G12V systems showed comparatively enhanced deviations in SW-I and P-loop regions (Fig. 2F & H). However, it was seen SW-II region showed comparatively higher deviation in WT (Fig. 2G). Thus introduction of mutations tends to make structure more rigid as compared to wild type in GDP-bound systems. However, in case of GTP-bound systems, the mutant systems G12D and G12V tend to become more flexible as compared to wild type. The switch regions are known to play an important role in the activation of the protein. The changes in the flexibility of the protein and these switch regions are known to have impact on binding of the effector proteins thereby affecting the activation process [80,81].

3.2 Root mean square fluctuation (RMSF)

The structural flexibility or rigidity owing to the mutations was assessed by computing the root mean square fluctuation (RMSF). Fig. 3A shows RMSF for GDP-bound systems. The SW-I region, loop between SW-I and SW-II (residue range 45–52), β4 (residue range 82–92) region, α4, β5 and loop region between α4 and β5 (residue range 125–150) majorly fluctuating in WT type system. The G12D and G12V systems showed higher fluctuations in the SW-I region, loop between SW-I and SW-II (residue range 45–52). In the case of the G13D system, major fluctuations were observed in SW-I region, loop between SW-I and SW-II (residue range 45–52) region, α4, β5 and loop region between α4 and β5 (residue range 125–150). In case of GTP-bound systems (Fig. 3B), the major fluctuation in WT systems were observed in the SW-I region, loop between SW-I and SW-II (residue range 45–52). Similar observations were seen for the G12D and G13D systems. The G12V system tends to show higher fluctuations as compared to other systems, indicating that the introduction of hydrophobic residue has a major impact. In the G12V systems, the SW-I region, loop between SW-I and SW-II (residue range 45–52), β4 (residue range 82–92), region α4, β5 and loop region between α4 and β5 (residue range 125–150) showed majorly fluctuation.Fig. 3 Root mean square fluctuations for all the four systems, WT (black), G12D (red), G12V (green) and G13D (blue). (A) GDP-bound systems. (B) GTP-bound systems.

Fig. 3

3.3 Distance based population landscape

In the KRAS, the active state of the protein is defined based on distance between the SW-I residue Thr35 and GDP/GTP and distance between the SW-II residue Gly60 and GDP/GTP. Hence the two dimensional population landscape is plotted between Thr35-GDP/GTP and Gly60-GDP/GTP. Fig. 4(A–D) shows that population landscape for GDP-bound systems. For the WT system, two major clusters were seen, cluster 1 has conformations with Thr35 and GDP distance around 11 Å and Gly60 and GDP distance around 7 Å. While cluster 2 has the population with Thr35 and GDP distance ranging from 5.5 to 11 Å and Gly60 and GDP distance around 14 Å. In the G12D system, one single cluster was observed with Thr35 and GDP distance around 6.5 Å and Gly60 and GDP distance around 5 Å. For the G12V system two clusters were observed, one with Thr35 and GDP distance around 6.5 Å and Gly60 and GDP distance around 4.5 Å. The other clusters have the Thr35 and GDP distance spread from 6 to 13 Å and Gly60 and GDP distance around 6 Å. In the case of the G13D system, two clusters were observed, one with Thr35 and GDP distance around 4 Å and Gly60 and GDP distance around 5 Å.Fig. 4 Population landscape based on distance between Thr35-GDP and Gly60-GDP. (A) Wild type system (B) G12D system. (C) G12V system. (B) G13D system. Population landscape based on distance between Thr35-GTP and Gly60-GTP. (E) Wild type system (F) G12D system. (G) G12V system. (H) G13D system.

Fig. 4

Similar to GDP for GTP-bound systems, the two dimensional population landscape plot is plotted between Thr35-GTP and Gly60-GTP. Fig. 4(E–H) shows the population landscape for GTP-bound systems. Unlike GDP-bound systems, in the GTP-bound systems a single major cluster was observed. For the WT the cluster has conformation with Thr35 and GTP distance around 7 Å and Gly60 and GTP distance around 7 Å. In the G12D system, the cluster was observed with Thr35 and GTP distance around 7 Å and Gly60 and GTP distance around 6.5 Å. For the G12V system, the cluster was observed with Thr35 and GTP distance around 6.5 Å and Gly60 and GTP distance showed spread from 6 to 8 Å. In the case of the G13D system, the cluster was observed as one with Thr35 and GTP distance spread from 5 to 8 Å and Gly60 and GTP distance around 6 Å.

It can be observed that the WT and mutant proteins had different distribution patterns, it may suggest that they explored different conformational states. On comparing the GDP-bound system with the GTP-bound system, it was observed that for WT the overall population distribution was similar. However in WT GDP-bound system the conformations in the major cluster were mostly open inactive state structures with higher Thr35 to GDP distance and in GTP-bound system the conformations in the major cluster were in intermediate state or closed active state with structure having smaller values of Thr35 to GTP distance. In the case of the G12D system, it was observed that population distribution was more dispersed along both axes in GTP-bound states as compared to GDP-bound. Unlike WT, in the G12D system major clusters in both GDP and GTP systems had conformations representing intermediate state or closed active state. For the G12V system, in GTP-bound state the population was more dispersed along both axes, however the conformation in major clusters in GTP-bound were representing mostly intermediate state or closed active state structure as compared to that of GDP-bound system.

It was observed that in all the GTP-bound systems, the population was more dispersed. The observations were similar to NMR studies which suggested that GTP-bound systems are more flexible and tend to explore different conformations [44,82]. In spite of the population being dispersed, the majority of the population was observed with a lower value of Thr35 and GTP distance as indicated by a single major cluster. This may suggest that in GTP-bound systems the structures try to converge towards an active closed state. However in the case of all the GDP-bound systems the population was not much dispersed but multiple clusters were observed. In the WT system the dominant clusters have conformations with larger values of Thr35 and GTP distance, while in mutant systems major clusters have mixed population, a fraction with lower values of Thr35 and GTP distance and a fraction with higher values of Thr35 and GTP distance. It may suggest that in GDP-bound systems the structure try to explore different intermediates in mutant systems and mostly remains in inactive open state structures in the WT system. For GTP-bound G12D and G12V systems, the population was more dispersed along SW-I. Also, it was observed that in the G13D system (GDP and GTP-bound), the population was concentrated. This shows that G13D mutations tend to be rigid system with respect to SW-I and SW-2.

3.4 Dynamic cross correlation (DCC)

The relative motion between different regions of the protein may help in deducing the effect of change in one region of protein on a distant region of the protein. The dynamic cross correlation calculation helps to identify the synchronous and asynchronous motion in different regions of the protein. Fig. 5(A–D) shows the DCC motion for GDP-bound systems. For the WT system, (Fig. 5A), the c-terminal tail (residues 170–180) showed strong negative correlation with SW-2, β4 to β6 and α3 to α5 regions. The strong positive correlation was observed within the β4 to β6 and α3 to α5 region. This indicate synchronous motion in these regions. The SW-I region showed negative correlation with β4- α3- β5- α4 region. For the G12D system (Fig. 5B), very weak positive correlation was observed for different regions of the protein. The G12V system (Fig. 5C), showed similarly pattern as that of WT. For the G13D system (Fig. 5D), it was seen that the c-terminal tail (residues 170–180) showed strong negative correlation with α3 and α4 region.Fig. 5 Dynamic Cross Correlation Motion (DCCM) for GDP-bound systems. (A) Wild type system (B) G12D system. (C) G12V system. (B) G13D system. Dynamic Cross Correlation Motion (DCCM) for GTP-bound systems. (E) Wild type system (F) G12D system. (G) G12V system. (H) G13D system.

Fig. 5

Fig. 5(E–H) shows the DCC motion for the GTP-bound system. For the WT system (Fig. 8A), it can be seen that the c-terminal tail (residues 170–180) showed strong negative correlation with the entire protein except the SW-I region. The weak positive correlation was observed between β4 to β6 and α3 to α5 regions. In the G12D system (Fig. 5F), very strong negative correlation was observed between c-terminal tail and β1, β3, α3, α4, and α5 regions. The SW-2 region shows negative correlation with the SW-1 region and β6. The loop region between β5 and α4 shows strong negative correlation with β1, SW-2, β4 and α3. The α3 shows strong positive correlation with the SW-2 region. The G12V system (Fig. 5G), shows a similar pattern for c-terminal region. The strong positive correlation was seen within the α3- β5- α4- β6 region. In the case of the G13D system (Fig. 5H), the c-terminal region showed strong negative correlation with α2, α3, α4, and α5. The strong positive correlation was seen within the α3- β5- α4- β6 region. In all the systems, it was observed that the HVR region tend to strong correlation with rest of the protein. The full length KRAS studies have suggested the folding of the HVR on top of the switch regions [83,84]. The study by Lu et al. suggests that the mutants are known to affect the interactions between HVR and G domain [53].

3.4.1 Difference in correlated motion between GTP and GDP-bound system

Further to see how these motions differ between the different systems, the root mean square change in cross-correlation was calculated. Fig. 6 shows the root mean square (RMS) cross correlation difference between GTP and GDP-bound systems. For the WT system (Fig. 6A), it was observed that on an average the RMS difference was around three units. Maximum deviation of around 4.7 was observed for α3 and loop between α3 and β5. Besides these regions, other regions which showed higher deviation were SW-I region and α5 region. For the G12D system (Fig. 6B), the average RMS cross correlation difference was around four and half units. The α3 and loop between α3 and β5 showed maximum difference and was around six units. Both the SW-I and SW-2 regions showed higher RMS cross correlation differences of around five units. The G12V (Fig. 6C) showed least RMS cross correlation difference values with an average of around one and half units. The G13D system (Fig. 6D) showed an average difference of around three units. Maximum RMS cross correlation difference was observed in the SW-2 region and was around four units. The β6 also showed a higher difference of around four units. Thus it was observed that not only SW1 and SW2 showed major deviation on transition from GDP to GTP, but the α3 and loop between α3 and β5 showed deviation. This suggests that these regions may be playing an important role in the transition from GDP-bound state to GTP-bound state.Fig. 6 Delta Root mean Square Cross- Correlation difference between GDP and GTP-bound system. (A) Between GDP-bound WT and GTP-bound WT. (B) Between GDP-bound G12D and GTP-bound G12D. (C) Between GDP-bound G12V and GTP-bound G12V. (D) Between GDP-bound G13D and GTP-bound G13D.

Fig. 6

3.4.2 Difference in correlated motion between WT and mutants for GDP-bound system

The Supplementary Fig. S1A shows the root mean square (RMS) cross correlation difference in correlation between WT and G12D for GDP-bound systems. It can be seen that a maximum difference of around six units was observed in α3 and loop between α3 and β5. Besides this the P-loop, SW-I region and region spanning β5 to α4 showed RMS cross correlation difference of about around four units. The average RMS cross correlation difference between WT and G12V for the GDP-bound system was on the lower side with a maximum value of four units was observed for SW-2 region (Supplementary Fig. S1B). Besides SW-2, P-loop and SW-I showed an RMS difference of around three units. The RMS cross correlation difference between WT and G13D for GDP-bound systems is shown in Supplementary Fig. S1C. It can be seen that a difference of around three and half units was observed for a region spanning from residue range 105–150 consisting of α3- β5-α4-β6. The α5 showed the RMS cross correlation difference of around four units. The difference of three units was observed for P-loop, SW-I and SW-2. Thus the comparison between WT and mutant systems gives an idea about the region which shows more deviation owing to mutations and these regions are α3 and loop between α3 and β5 and β6.

3.4.3 Difference in correlated motion between WT and mutants for GTP-bound system

The root mean square (RMS) cross correlation difference in correlation between WT and G12D for GTP-bound systems is shown in Supplementary Fig. S2A. From the figure it can be seen that the maximum difference of around five units was observed for the α3 and α5 region. Besides these regions the difference of around four units was observed for the P-loop region, around three and half units for the SW-I and SW-2 region. The RMS cross correlation difference between WT and G12V for GTP-bound systems is shown in Supplementary Fig. S2B. The RMS cross correlation difference of around 3.8 units was observed for SW-1 region, three and half units for the SW-2 region and around four units for loop region between α3 and β5. The α5 region showed an RMS difference of around four units. The RMS difference between WT and G13D for GTP-bound systems is shown in Supplementary Fig. S2C. The maximum difference of around five units was seen for α5. The P-loop, SW-I and SW-2 region showed differences of around three units.

From the comparison of differences in dynamics correlation, it was observed that along with SW-I, SW-2 region the α3 region also plays an important role in the activation process and hence may be considered for drug designing study. The maximum deviation was observed in the G12D system and minimum for the G12V system. The study by Chenn et al. [85] has also shown that G12D mutants tend to show maximum fluctuation in molecular simulation study. Also, the dynamic correlation motion suggests not only the motions of individual domains plays an important role in the activation but the correlation in the motion between different regions viz. SW-I, SW-2 and α3 may play significant roles during the transition of KRAS from one conformational state to another. The importance of communication between SW-2 and α3 has also been shown by Khaled et al. [86]. Grant et al. has also suggested the importance of correlated motion between α2 (SW-2 region) and α3 in GTP-bound states [87].

3.5 Local conformational changes

As the backbone angle plays important role in maintaining the structure of the protein, any deviation in the backbone angle may give idea about flipping/motion of that residue/s with respect to rest of the protein. Hence, the backbone dihedral angle PHI/PSI have been calculated and the angles which are showing deviation among WT and mutant systems have been given in the Supplementary Fig. 3. From the figure it can be clearly seen that in GDP-bound state, the residues 28–32, 35–40, 57–68, 71, 72, 104, 105, 106, 108,149 show deviation in PHI/PSI angles (Supplementary Figs. S5A and S5B). Among the residues, the residue 29, 30, 32, 36, 37, 40, 59–64, 106 and 107 showed major difference among WT and mutant systems for PHI angle ((Supplementary Fig. S3A) and the residues 30, 31, 35–38, 60–62, 64, 67, 104, 108 and 149 showed major differences for PSI angle (Supplementary Fig. S3B). In the GTP-bound systems the number of residues with variation in torsion angle increased (Supplementary Figs. S3C and S3D). The residues 12–14, 28–40, 46, 47, 57–68, 101,103–109, 115, 122,123, 145,150, 154 showed deviation in PHI/PSI angles. Among these residues maximum deviation in PHI angle between WT and mutant system was observed 12, 14, 28, 29, 35, 38, 39, 40, 61–64, 105, 106,108,123, 154 (Supplementary Fig. S3C). In case of PSI angle maximum deviation was observed in 12, 13, 34, 37, 38, 60–64, 101,104,105,108, 121 and 154 (Supplementary Fig. S3D). It can be clearly seen that in both GDP and GTP-bound system the maximum deviation is observed in SW-I, SW-II and α3 regions, however the number of residues showing deviation in SW-I, SW-II and α3 regions increased in GTP-bound systems. Also it was observed that in the GTP-bound system the P-loop region and loop between β6 and α5 also showed changes. Thus, these changes in backbone dihedral of different residues may suggest that the GDP-bound mutant system tends to explore different intermediates conformations by fluctuation in the SW-I, SW-II and α3 regions. However the presence of GTP expands these conformation exploration in different regions thereby reaching the active state. These regions showing fluctuations in the mutant GTP-bound state may be of interest to target mutant KRAS systems.

3.6 Interactions of GTP/GDP with rest of protein

3.6.1 Hydrogen and van Der waals interactions

As various residues tend to form interactions with GDP/GTP thereby assisting in the hydrolysis or exchange, hence hydrogen bond and VDW interactions are calculated for GDP and GTP with the rest of the protein. All the interactions having occupancy more than 20 % have been shown. Fig. 7A and B shows the hydrogen bond and van der Waals interaction for GDP with the rest of the protein residues. The mutation in P-loop region affected the hydrogen bond interaction of P-loop region with GDP, the decrease in the hydrogen bond occupancy was observed for Gly15 and Lys16 residues for all the mutant systems and complete loss of hydrogen bond for Val14 residue was observed. However, the mutated Asp13 in G13D tends to form a hydrogen bond with GDP. Further it was observed that G13D mutation helps the GDP to form hydrogen bonds with Thr35 and thereby converge the SW-I region towards closed active state. Also the mutation G12D and G13D tend to affect the β5 (Asn116) and β6 (Ser145, Ala146 and Lys147) regions where the hydrogen bonding occupancy was reduced for both these mutants. The effect of mutation on van der Waals interaction was not prominent in the P-loop region, except the Val14 where the occupancy for the mutants decreased as compared to WT. The SW-I showed variations in the van der Waals interactions. The interaction with Asp30 was lost for G12V and G13D systems and Glu31 showed decrease in occupancy for the all mutant systems. However the Tyr32 showed increased in occupancy and for Pro34 and Thr35 no interactions were seen for WT type. The mutation G12D and G13D also tend to affect the β5 (Asn116, Lys117 and Leu120) and β6 (Ser145, Ala146 and Lys147) region where the van der Waals interaction occupancy was reduce for both this mutants.Fig. 7 (A) Hydrogen bonding and (B) hydrophobic interactions between GDP and residues of the KRAS for wild type and mutant systems. The X-axis represents the protein residues and Y-axis frequency of occurrence of the interaction during course of interactions.

Fig. 7

Fig. 8 (A) Hydrogen bonding and (B) hydrophobic interaction between GTP and the residues of KRAS for wild type and mutant systems. The X-axis represents the protein residues and Y-axis represent frequency of occurrence of the interaction during the course of interactions.

Fig. 8

For the GTP-bound systems, the hydrogen bond interactions are shown in Fig. 8. Mutations at 12 and 13 affects the hydrogen at position 13 (Fig. 8A). The occupancy of hydrogen bond showed decrease for residue 13. The Gly15 and Lys17 residues also showed decrease in the occupancy for the mutant systems. The decrease in occupancy was also seen for SW-I region residues Tyr32 and Thr35 for the mutant systems. The effect of mutation on van der Waals interaction (Fig. 8B) was seen in the G13D system in the P-loop region where small decrease in occupancy was seen for Asp13, Val14, Gly15, Lys16, Ser17 and Ala18. The increase in occupancy for Asp33 residue was seen for the G12D and G13D systems.

3.6.2 Water-based interactions

In the KRAS, it is known that water plays an important role in hydrolysis of GTP. Also water based bridge interactions may help the KRAS to adopt different conformational states. The water based interactions were calculated for both GDP and GTP. Fig. 9A represents the water based bridge interactions for GDP-bound systems. From the figure it can be seen that the G12D and G13D mutations tend to increase the occupancy of water based interaction in the P-loop region. In the G12D system Ala11, Asp12 and Val14 form water based interaction with GDP which are not present in any other GDP-bound systems. Similarly, in the G13D system the Ala18 interaction with GDP is not found in any system and also the occupancy of Ser17 interaction with GDP is increased. In the SW-I region the occupancy of Tyr32, Asp33 and Thr35 water mediated interaction with GDP is increased in mutant systems. Also in the G13D system Asp38 forms water mediated interaction with GDP. In the SW-2 region, Asp58 forms interaction with GDP only in the G13D system. In the G12D system the occupancy of Gly60 interaction is increased and Gln61 was found to form interaction. In the case of the GTP-bound system (Fig. 9B), the effect of mutations on water mediated interaction of protein residues with GTP is predominant in the P-loop region. The Asp12 and Val14 form interaction only in the G12D system and the occupancy of Asp13 interaction with GTP shows increase with respect to other residues at position 13 in other systems. In the SW-I region the occupancy of Asp33 interaction is increased for G12D and G13D systems.Fig. 9 Water mediated interaction between GDP/GTP and residues of KRAS for wild type and mutant systems. The X-axis represents the protein residues and Y-axis frequency of occurrence of the interaction during the course of interactions. A) Water mediated interaction between GDP and residues of KRAS. B) Water mediated interaction between GTP and residues of KRAS.

Fig. 9

The pocket of KRAS where GDP/GTP binds is overall positively charged and it helps to counter balance the negative charge of GDP/GTP phosphate groups. Along with Tyr32, Thr35, Gly60 and Gln61, the P-loop residues play an important role in stabilizing the position of GDP/GTP in the cavity thereby stabilizing the GDP/GTP interaction with various residues. The mutation of Gly to Asp at 12 and 13th position leads to charge reversal on the binding cavity and may be responsible for the breaking/formation of the important hydrogen bonds and water mediated interaction important for GTP hydrolysis and thereby promoting the KRAS in continuous activation state. Tyr32 is known to play an important role in GTP hydrolysis, as it known to show major reorientation during nucleotide exchange [88,89]. The increase in hydrogen bond occupancy for 32, 35 and 60 residues for GTP-bound systems can be seen and this may help in stabilizing the active state of KRAS. However, residues from P-loop have formed less stable water mediated H-bonds in GTP bound systems as compared to GDP bound systems. The SW-I and SW-II residues were preferred for water mediated H-bonds in all the WT and mutants in GTP bound systems.

3.7 Markov state model analysis

Markov state models help to explore the dynamics of the protein at a larger timescale which is not accessible by traditional molecular dynamics simulations. It helps to identify the highly stable or intermediate protein conformational states and kinetics between them by integrating the data from multiple simulations. The most important aspect while building the Markov state model is the choice of the reaction coordinates or feature which can help in discretizing the different conformational states of the protein with reduced dimensionality. In the current study, Markov state model is build using the residue_min_dist as the featurizer considering following residues 12, 29, 32, 34, 35, 36, 40, 48, 59, 60, 61, 62, 67, 105, 122, 126, 138 and GDP for GDP-bound system and GTP for GTP-bound systems. The following residues were chosen as these residues play an important role in the GTP hydrolysis and also were found to show maximum deviation in the current study. The two MSMs were constructed, one for GDP-bound systems and another for GTP-bound systems. All the GDP-bound systems viz. WT, G12D, G12V and G13D system trajectories were considered for building the first MSM. Similarly all the GTP-bound systems were considered for the second MSM. In case of both the MSMs, four metastable states were observed, however the distribution of the conformational states differs.

For the GDP-bound MSM construct, the four metastable states along with the representative structures generated in MSM are shown in Fig. 10A. The cluster 1 (C1) consists of the conformation belonging to WT system trajectories. The cluster 2 (C2) consists of the conformation mainly from G12V system trajectories. A small population of the C2 also consists of the population from G13D. The cluster 3 (C3) consists of the conformation from all the four systems, however it is mainly dominated by the population of the G12D and G13D systems. The cluster 4 (C4) consists of the conformation belonging to WT system trajectories. From the cluster distribution it can be seen that the WT system showed a more dispersed population spread across different clusters. Further to check the properties of these clusters, in terms of which state they belong i.e closed active or open inactive, distribution of the important distance and angle was seen for these clusters and has been represented in Fig. 10B. The cluster 1 mainly consists of the conformation belonging to inactive state and some of the intermediate I population, where the distance between Thr35 and GDP fluctuate in the range of the 5 to 14 Å, with maximum conformation having distance around 10 Å. In cluster 1 the Gly60 and GDP distance was quite stable with maximum conformation having distance around 9 Å. The PHI 61 angle showed variation of around 120° ranging from −50 to −170 with maximum population around −100°, while PSI 61 angle had values concentrated around 150°. The cluster 2 which consists of population from G12V and G13D showed conformation belonging to the intermediate state, where the distance between Thr35 and GDP fluctuated in the range of 4–12 Å, with maximum conformation having distance around 7 Å. The Gly60 and GDP distance fluctuated in the range of 4–12 Å, with maximum conformation having a distance around 9 Å. The PHI 61 angle showed maximum variation among all the clusters with range from 0 to −180 with maximum population around −50°, while PSI 61 angle had values concentrated around 130°. The cluster 3 consists of conformation belonging to intermediate as well as closed active state. The distance between Thr35 and GDP was concentrated around 5 Å, Gly60 and GDP distance fluctuated in the range of 4–12 Å, with maximum conformation having a distance around 10 Å. The majority of conformation have PHI 61 value around −60° and PSI 61 value around 140°. The cluster 4 consists of the population belonging to inactive open conformations. The distance between Thr35 and GDP was concentrated around 15 Å, Gly60 and GDP distance fluctuated in the range of 3–16 Å, with maximum conformation having a distance around 9 Å. The majority of conformation have PHI 61 value around −150° and PSI 61 value around 30°. Further to these distances and angles, RMSF was plotted for all the four clusters, the most fluctuating region was c-terminal hyper variable region with RMSF value of around 10 Å. The SW-I region showed fluctuations of around 4 Å, 3 Å, 2 Å and 1.5 Å in cluster C1, C2, C3 and C4 respectively. The SW-II region showed fluctuation of around 4.5 Å in cluster C1 to C3 and around 1.5 Å in C4. The region connecting SW-1 and SW-2 showed fluctuation of around 3 Å in C1 and C3 and around 1.7 Å in cluster C2 and C4.The α3 helix also showed higher fluctuation of around 3 Å. In case of GDP most of the WT conformation were in inactive open state, while all the mutants were in either intermediate state where the Thr35 tends to move close to GDP and Gly60 was away from GDP.Fig. 10 A Free energy distribution map along time-lagged independent components (ICs) 1 and 2 for GDP-bound systems. Metastable states (clusters C1 to C4) of KRAS are shown in the figures. The metastable states are illustrated with ten representative structures.

Fig. 10 A

Fig. 10 B Violin plots showing the distribution of various properties of metastable state obtained from Markov state model for GDP-bound systems.

a) Distance distribution between Thr35 and GDP.

b) Distance distribution between Gly60 and GDP.

c) PHI 61 angle distribution.

d) PSI 61 angle distribution.

e) RMSF plot showing fluctuations in various regions of KRAS.

Fig. 10 B

For the GTP-bound MSM construct, the four metastable states along the representative structure generated in MSM is shown in Fig. 11A. The cluster 1 (C1) consists of the conformation belonging to G12V system trajectories. The cluster 2 (C2) is one of the dominant clusters consisting of the conformation mainly from the G13D system, beside this it also consists of considerable conformation of other three systems. The cluster 3 (C3) consists of the conformation from only WT system trajectories. The cluster 4 (C4) consists of the conformation belonging from all four systems. From the cluster distribution it can be seen that owing to G13D mutation the system becomes more rigid. Further to check the properties of these clusters, in terms of which state they belong i.e closed active or open inactive, distribution of the important distance and angle was seen for these clusters and has been represented in Fig. 11B. The conformation in cluster 1 belongs mainly to the intermediate state or active closed state. In cluster 1, the conformations had a distance between Thr35 and GTP around 6 Å, while the distance between Gly60 and GTP was around 8 Å. For the majority of conformation the PHI 61 angle value was around −125° and PSI 61 angle value was around 150°. The cluster 2 which consist of population from all the systems consist of conformation mainly from intermediate state or active closed state. The distance between Thr35 and GTP was around 6 Å, and distance between Gly60 and GTP was around 8 Å. PHI 61 values were concentrated around −100° and PSI 61 values around 150. The difference between cluster 1 and cluster 2 was population with more distributed PHI and PSI angle values. The cluster 3 consists of conformation from open inactive and intermediate state populations. The distance between Thr35 and GTP ranges from 4 to 14 Å, with maximum around 8 Å and the distance between Gly60 and GTP have population mainly around 8 Å. The PHI 61 angle ranges from −175 to −50 with most of the population concentrated around −125° and PSI 61 value mainly concentrated around 140°. The cluster 4 mainly consists of intermediate state populations with Thr35 and GTP distance around 6 Å, Gly60 and GTP distance around 8 Å. The PHI 61 value mainly concentrated around −150° and PSI value around 130°. The RMSF was plotted for all the four clusters, similar to the GDP-bound system; the most fluctuating region was c-terminal hyper variable region with RMSF value more than 10 Å. The SW-I region showed fluctuation of around 3 Å, 2 Å, 5 Å and 2 Å in cluster C1, C2, C3 and C4 respectively. The SW-II region showed fluctuation of around 4.5 to 5 Å in all four clusters. The region connecting SW-1 and SW-2 showed fluctuation of around 3.5 Å in all four clusters. The α3 helix also showed higher fluctuation of around 3 Å.Fig. 11A Free energy distribution map along time-lagged independent components (ICs) 1 and 2 for GTP-bound systems. Metastable states (clusters C1 to C4) of KRAS are shown in the figures. The metastable states are illustrated with ten representative structures.

Fig. 11A

Fig. 11B Violin plots showing the distribution of various properties of metastable state obtained from Markov state model for GTP-bound systems.

a) Distance distribution between Thr35 and GTP.

b) Distance distribution between Gly60 and GTP.

c) PHI 61 angle distribution.

d) PSI 61 angle distribution.

e) RMSF plot showing fluctuations in various regions of KRAS.

Fig. 11B

It can be clearly seen that in the GDP-bound population, the distance distribution for Thr35 and Gly60 with GDP was more dispersed and in the case of GTP-bound systems the distance between Thr35 and Gly60 with GTP was more concentrated towards lower values. Further it was observed that in GDP-bound systems, the SW-I region showed more flexibility and in GTP-bound systems the SW-II along with the region connecting the SW-I to SW-II showed higher fluctuations. Also α3 helix showed higher fluctuation in both the GDP and GTP-bound system. Such observations where multiple conformational states are explored by KRAS has been reported in the earlier experimental studies. Our observation of multiple conformational states for active KRAS is consistent with previous experimental studies, where two states viz. catalytic competent state (R state) and catalytic incompetent state (T state) has been reported[45, 48 89].

3.8 Conformational dynamics in KRAS mutants and its inhibitors

The cryptic nature of KRAS makes it difficult to understand the different conformations that it may adopt in the presence of varied cancerous mutations, thereby posing a challenge for inhibitor designing against them. It is already known through literature that different KRAS mutants show distinct activation mechanisms [45,48,57,90]. The study by Grudzien et al. [57] also has indicated divergent activation mechanism by distinct mutants. Thus, understanding of varied structural changes in distinct KRAS mutants may provide some insights in the designing of inhibitors. The KRAS was once considered as un-targetable, but the study by Shokat and his team in year 2013 helped to identify new allosteric pocket in G12C-mutant of KRAS [90].They identify new allosteric pocket beneath SW-II region and used covalent binding approach for targeting G12C-mutant KRAS [90]. Consequently, this study guided several studies and which led to development of different G12C-mutant KRAS inhibitors and some of them have are in the different phases of clinical study for example ARS-1620 [91,92]. Some of the drugs which are studied for selectively targeting the KRAS mutant are Adagrasib (MRTX849), Lumakras (Sotorasib; AMG 510), ARS-3248: ARS-3248, MRTX1133 and many more [93]. All the above mentioned drugs targets G12C, except the MRTX1133 which targets selectively to KRAS G12D mutant [94]. The MRTX1133 is the non-covalent inhibitor and provides an alternative strategy for targeting KRAS, as most of drug targeting G12C mutants are covalent inhibitors. Few of recent reviews give details about various inhibitors in different stages of drug development and the role of molecular simulations for their development [[95], [96], [97]]. The molecular dynamics simulations have played important role in identifying novel cryptic pockets [98,99]. Recent advances in molecular simulation techniques have enabled the exploration of the structural and dynamic properties of KRAS, providing new insights into its function and mechanisms of action. Thus, the knowledge of the regions playing important role in activation mechanism of KRAS in different mutants may guide in designing better and alternative strategy for targeting KRAS.

4 Conclusion

Detailed structural analysis of GDP and GTP bound systems helped to understand dynamics of multiple KRAS regions like SW-I, SW-II, α3, loops etc. Most clusters in WT GDP-bound system had conformations belonging to inactive state. However, the introduction of mutation shifted the population towards the intermediate state (showing smaller distance for Thr35 or Gly60 with GDP). In GTP-bound mutant systems, the major cluster had a population belonging to the active state. However, in the GTP-bound WT system along with active state conformations, clusters with intermediate state were also observed. The reason for such a shift in conformation in mutant systems may be attributed to interaction network formed by P-loop, SW-I, SW-II regions along with various water based interactions. The other regions like α3 and loops have indirect role in these activation/hydrolysis processes via allosteric effects. The P-loop residue 12, SW-I residues 29, 32, 34, 35, 36, 40, 48, SW-II residues 59, 60, 61, 62, 67, α3 residue 105, and loop region residues 122, 126 and 138 may be playing important role in transmitting the allosteric signals. The HVR region also tend to influence the activation process by having strong negative correlation with these regions. The analysis of clusters obtained by Markov state model analysis suggested that in case of GDP bound systems, the difference in motion of the SW-I region in WT and mutant systems observed. The presence of bulky Valine and Aspartic acid affect the motion of SW-I region leading to locking the conformations in mostly intermediate state with SW-I in closed conformation. In the case of GTP bound systems, the presence of bulky mutant residues and bulky GTP together affect the motion of SW-I and SW-II regions. Thus the introduction of bulky residue (Val) and the charged residue (Asp) may be playing a role in activation/hydrolysis by changing the local charge distribution and altering the dynamic interaction network of mentioned residues bringing SW-I and SW-II in closed active conformation. In the current study, it was also observed that along with SW-I and SW-II region, the region connecting SW-I and SW-II, the α3 region and the loop preceding α3 affects the interactions of GDP/GTP with different regions of the protein plays important in activation mechanism. These structural insights of WT and mutants in the GDP or GTP form may help for designing better therapeutics.

Supporting Information

The supporting information files contains Delta Root mean Square Correlation, PHI and PSI angle distribution for various residues.

Data availability statement

The protein structure of KRAS (PDB ID: 4DSU and 4DSO) were retrieved from RCSB protein data bank (https://www.rcsb.org/). The residues were mutated using Chimera software a free software (https://www.cgl.ucsf.edu/chimera). Molecular dynamic simulations were carried out on GROMACS, a free software (https://www.gromacs.org/). The analysis was performed using AmberTools and the gromacs analysis utilities. The additional raw data can be made available upon request to corresponding author.

CRediT authorship contribution statement

Vinod Jani: Writing – original draft, Visualization, Validation, Methodology, Data curation, Conceptualization. Uddhavesh Sonavane: Writing – review & editing, Conceptualization. Rajendra Joshi: Writing – review & editing, Funding acquisition, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is the supplementary data to this article:Multimedia component 1

Multimedia component 1

Acknowledgments

The authors gratefully acknowledge National Supercomputing Mission (NSM) project under the 10.13039/501100008628 Ministry of Electronics and Information Technology (MeitY) , and 10.13039/501100001409 Department of Science and Technology, Government of India, New Delhi , for providing financial support. The autor would like to acknowledge “Bioinformatics Resources and Applications Facility (BRAF)” and PARAM supercomputing facility for providing computing resources.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e36161.
==== Refs
References

1 Neves S.R. Ram P.T. Iyengar R. G protein pathways Science 296 5573 2002 May 31 1636 1639 10.1126/science.1071550 PMID: 12040175 https://pubmed.ncbi.nlm.nih.gov/12040175/ 12040175
2 White M.A. Nicolette C. Minden A. Polverino A. Van Aelst L. Karin M. Wigler M.H. Multiple Ras functions can contribute to mammalian cell transformation Cell 80 4 1995 Feb 24 533 541 10.1016/0092-8674(95)90507-3 PMID: 7867061 https://pubmed.ncbi.nlm.nih.gov/7867061/ 7867061
3 Khrenova M.G. Mironov V.A. Grigorenko B.L. Nemukhin A.V. Modeling the role of G12V and G13V Ras mutations in the Ras-GAP-catalyzed hydrolysis reaction of guanosine triphosphate Biochemistry 53 45 2014 Nov 18 7093 7099 10.1021/bi5011333 Epub 2014 Nov 4. PMID: 25339142 https://pubmed.ncbi.nlm.nih.gov/25339142/ 25339142
4 Spoerner M. Hozsa C. Poetzl J.A. Reiss K. Ganser P. Geyer M. Kalbitzer H.R. Conformational states of human rat sarcoma (Ras) protein complexed with its natural ligand GTP and their role for effector interaction and GTP hydrolysis J. Biol. Chem. 285 51 2010 Dec 17 39768 39778 10.1074/jbc.M110.145235 Epub 2010 Oct 11. PMID: 20937837; PMCID: PMC3000958 https://pubmed.ncbi.nlm.nih.gov/20937837/ 20937837
5 Ma J. Karplus M. Molecular switch in signal transduction: reaction paths of the conformational changes in ras p21 Proc Natl Acad Sci U S A 94 22 1997 Oct 28 11905 11910 10.1073/pnas.94.22.11905 PMID: 9342335; PMCID: PMC23651 https://pubmed.ncbi.nlm.nih.gov/9342335/ 9342335
6 Bos J.L. Rehmann H. Wittinghofer A. GEFs and GAPs: critical elements in the control of small G proteins Cell 129 5 2007 Jun 1 865 877 10.1016/j.cell.2007.05.018 Erratum in: Cell. 2007 Jul 27;130(2):385 https://pubmed.ncbi.nlm.nih.gov/17540168/ 17540168
7 Cherfils J. Chardin P. GEFs: structural basis for their activation of small GTP-binding proteins Trends Biochem. Sci. 24 8 1999 Aug 306 311 10.1016/s0968-0004(99)01429-2 PMID: 10431174 https://pubmed.ncbi.nlm.nih.gov/10431174/ 10431174
8 Goitre L. Trapani E. Trabalzini L. Retta S.F. Ras Signaling Methods and Protocols 2014 Humana Press Totowa, NJ 10.1007/978-1-62703-791-4
9 Castellano E. Santos E. Functional specificity of ras isoforms: so similar but so different Genes Cancer 2 3 2011 Mar 216 231 10.1177/1947601911408081 PMID: 21779495; PMCID: PMC3128637 https://pubmed.ncbi.nlm.nih.gov/21779495/ 21779495
10 Lu S. Jang H. Muratcioglu S. Gursoy A. Keskin O. Nussinov R. Zhang J. Ras conformational ensembles, allostery, and signaling Chem. Rev. 116 11 2016 Jun 8 6607 6665 10.1021/acs.chemrev.5b00542 Epub 2016 Jan 27. PMID: 26815308 https://pubmed.ncbi.nlm.nih.gov/26815308/ 26815308
11 van 't Veer L.J. Burgering B.M. Versteeg R. Boot A.J. Ruiter D.J. Osanto S. Schrier P.I. Bos J.L. N-ras mutations in human cutaneous melanoma from sun-exposed body sites Mol. Cell Biol. 9 7 1989 Jul 3114 3116 10.1128/mcb.9.7.3114-3116.1989 PMID: 2674680; PMCID: PMC362784 https://pubmed.ncbi.nlm.nih.gov/2674680/ 2674680
12 Bos J.L. Ras oncogenes in human cancer: a review Cancer Res. 49 17 1989 Sep 1 4682 4689 Erratum in: Cancer Res 1990 Feb 15;50(4):1352. PMID: 2547513 https://pubmed.ncbi.nlm.nih.gov/2547513/ 2547513
13 Riely G.J. Marks J. Pao W. KRAS mutations in non-small cell lung cancer Proc. Am. Thorac. Soc. 6 2 2009 Apr 15 201 205 10.1513/pats.200809-107LC PMID: 19349489 https://pubmed.ncbi.nlm.nih.gov/19349489/ 19349489
14 Wittinghofer A. Vetter I.R. Structure-function relationships of the G domain, a canonical switch motif Annu. Rev. Biochem. 80 2011 943 971 10.1146/annurev-biochem-062708-134043 PMID: 21675921 https://pubmed.ncbi.nlm.nih.gov/21675921/ 21675921
15 Hancock J.F. Paterson H. Marshall C.J. A polybasic domain or palmitoylation is required in addition to the CAAX motif to localize p21ras to the plasma membrane Cell 63 1 1990 Oct 5 133 139 10.1016/0092-8674(90)90294-o PMID: 2208277 https://pubmed.ncbi.nlm.nih.gov/2208277/ 2208277
16 Abraham S.J. Muhamed I. Nolet R. Yeung F. Gaponenko V. Expression, purification, and characterization of soluble K-Ras4B for structural analysis Protein Expr. Purif. 73 2 2010 Oct 125 131 10.1016/j.pep.2010.05.015 Epub 2010 Jun 8. PMID: 20566322 https://pubmed.ncbi.nlm.nih.gov/20566322/ 20566322
17 Welman A. Burger M.M. Hagmann J. Structure and function of the C-terminal hypervariable region of K-Ras4B in plasma membrane targetting and transformation Oncogene 19 40 2000 Sep 21 4582 4591 10.1038/sj.onc.1203818 PMID: 11030147 https://pubmed.ncbi.nlm.nih.gov/11030147/ 11030147
18 Prior I.A. Lewis P.D. Mattos C. A comprehensive survey of Ras mutations in cancer Cancer Res. 72 10 2012 May 15 2457 2467 10.1158/0008-5472.CAN-11-2612 PMID: 22589270; PMCID: PMC3354961 https://pubmed.ncbi.nlm.nih.gov/22589270/ 22589270
19 Bourne H.R. Sanders D.A. McCormick F. The GTPase superfamily: conserved structure and molecular mechanism Nature 349 6305 1991 Jan 10 117 127 10.1038/349117a0 PMID: 1898771 https://pubmed.ncbi.nlm.nih.gov/1898771/ 1898771
20 Valencia A. Chardin P. Wittinghofer A. Sander C. The ras protein family: evolutionary tree and role of conserved amino acids Biochemistry 30 19 1991 May 14 4637 4648 10.1021/bi00233a001 PMID: 2029511 https://pubmed.ncbi.nlm.nih.gov/2029511/ 2029511
21 Muñoz-Maldonado C. Zimmer Y. Medová M. A comparative analysis of individual RAS mutations in cancer biology Front. Oncol. 9 2019 Oct 18 1088 10.3389/fonc.2019.01088 PMID: 31681616; PMCID: PMC6813200 https://pubmed.ncbi.nlm.nih.gov/31681616/ 31681616
22 Bos J.L. Rehmann H. Wittinghofer A. GEFs and GAPs: critical elements in the control of small G proteins Cell 129 5 2007 Jun 1 865 877 10.1016/j.cell.2007.05.018 Erratum in: Cell. 2007 Jul 27;130(2):385. PMID: 17540168 https://pubmed.ncbi.nlm.nih.gov/17540168/ 17540168
23 Glennon T.M. Villà J. Warshel A. How does GAP catalyze the GTPase reaction of Ras? A computer simulation study Biochemistry 39 32 2000 Aug 15 9641 9651 10.1021/bi000640e PMID: 10933780 https://pubmed.ncbi.nlm.nih.gov/10933780/ 10933780
24 Cherfils J. Chardin P. GEFs: structural basis for their activation of small GTP-binding proteins Trends Biochem. Sci. 24 8 1999 Aug 306 311 10.1016/s0968-0004(99)01429-2 PMID: 10431174 https://pubmed.ncbi.nlm.nih.gov/10431174/ 10431174
25 Kapoor A. Travesset A. Differential dynamics of RAS isoforms in GDP- and GTP-bound states Proteins 83 6 2015 Jun 1091 1106 10.1002/prot.24805 Epub 2015 Apr 22. PMID: 25846136 https://pubmed.ncbi.nlm.nih.gov/25846136/ 25846136
26 Gorfe A.A. Grant B.J. McCammon J.A. Mapping the nucleotide and isoform-dependent structural and dynamical features of Ras proteins Structure 16 6 2008 Jun 885 896 10.1016/j.str.2008.03.009 PMID: 18547521; PMCID: PMC2519881 https://pubmed.ncbi.nlm.nih.gov/18547521/ 18547521
27 Raimondi F. Portella G. Orozco M. Fanelli F. Nucleotide binding switches the information flow in ras GTPases PLoS Comput. Biol. 7 3 2011 Mar e1001098 10.1371/journal.pcbi.1001098 Epub 2011 Mar 3. PMID: 21390270; PMCID: PMC3048383 https://pubmed.ncbi.nlm.nih.gov/21390270/
28 Hall B.E. Bar-Sagi D. Nassar N. The structural basis for the transition from Ras-GTP to Ras-GDP Proc Natl Acad Sci U S A 99 19 2002 Sep 17 12138 12142 10.1073/pnas.192453199 Epub 2002 Sep 4. PMID: 12213964; PMCID: PMC129411 https://pubmed.ncbi.nlm.nih.gov/12213964/ 12213964
29 Lukman S. Grant B.J. Gorfe A.A. Grant G.H. McCammon J.A. The distinct conformational dynamics of K-Ras and H-Ras A59G PLoS Comput. Biol. 6 9 2010 Sep 9 e1000922 10.1371/journal.pcbi.1000922 PMID: 20838576; PMCID: PMC2936511 https://pubmed.ncbi.nlm.nih.gov/20838576/
30 Smith M.J. Neel B.G. Ikura M. NMR-based functional profiling of RASopathies and oncogenic RAS mutations Proc Natl Acad Sci U S A 110 12 2013 Mar 19 4574 4579 10.1073/pnas.1218173110 Epub 2013 Mar 4. PMID: 23487764; PMCID: PMC3607025 https://pubmed.ncbi.nlm.nih.gov/23487764/ 23487764
31 Stolze B. Reinhart S. Bulllinger L. Fröhling S. Scholl C. Comparative analysis of KRAS codon 12, 13, 18, 61, and 117 mutations using human MCF10A isogenic cell lines Sci. Rep. 5 2015 Feb 23 8535 10.1038/srep08535 PMID: 25705018; PMCID: PMC4336936 https://pubmed.ncbi.nlm.nih.gov/25705018/ 25705018
32 Hunter J.C. Manandhar A. Carrasco M.A. Gurbani D. Gondi S. Westover K.D. Biochemical and structural analysis of common cancer-associated KRAS mutations Mol. Cancer Res. 13 9 2015 Sep 1325 1335 10.1158/1541-7786.MCR-15-0203 Epub 2015 Jun 2. PMID: 26037647 https://pubmed.ncbi.nlm.nih.gov/26037647/ 26037647
33 Lu S. Jang H. Nussinov R. Zhang J. The structural basis of oncogenic mutations G12, G13 and Q61 in small GTPase K-Ras4B Sci. Rep. 6 2016 Feb 23 21949 10.1038/srep21949 PMID: 26902995; PMCID: PMC4763299 https://pubmed.ncbi.nlm.nih.gov/26902995/
34 Pantsar T. The current understanding of KRAS protein structure and dynamics Comput. Struct. Biotechnol. J. 18 2019 Dec 26 189 198 10.1016/j.csbj.2019.12.004 PMID: 31988705; PMCID: PMC6965201 https://pubmed.ncbi.nlm.nih.gov/31988705/ 31988705
35 Sharma N. Sonavane U. Joshi R. Comparative MD simulations and advanced analytics based studies on wild-type and hot-spot mutant A59G HRas PLoS One 15 10 2020 Oct 16 e0234836 10.1371/journal.pone.0234836 PMID: 33064725; PMCID: PMC7567374 https://pubmed.ncbi.nlm.nih.gov/33064725/
36 Diaz J.F. Wroblowski B. Schlitter J. Engelborghs Y. Calculation of pathways for the conformational transition between the GTP- and GDP-bound states of the Ha-ras-p21 protein: calculations with explicit solvent simulations and comparison with calculations in vacuum Proteins 28 3 1997 Jul 434 451 PMID: 9223188 https://pubmed.ncbi.nlm.nih.gov/9223188/ 9223188
37 Díaz J.F. Escalona M.M. Kuppens S. Engelborghs Y. Role of the switch II region in the conformational transition of activation of Ha-ras-p21 Protein Sci. 9 2 2000 Feb 361 368 10.1110/ps.9.2.361 PMID: 10716188; PMCID: PMC2144537 https://pubmed.ncbi.nlm.nih.gov/10716188/ 10716188
38 Spoerner M. Hozsa C. Poetzl J.A. Reiss K. Ganser P. Geyer M. Kalbitzer H.R. Conformational states of human rat sarcoma (Ras) protein complexed with its natural ligand GTP and their role for effector interaction and GTP hydrolysis J. Biol. Chem. 285 51 2010 Dec 17 39768 39778 10.1074/jbc.M110.145235 Epub 2010 Oct 11. PMID: 20937837; PMCID: PMC3000958 https://pubmed.ncbi.nlm.nih.gov/20937837/ 20937837
39 Lu S. Jang H. Muratcioglu S. Gursoy A. Keskin O. Nussinov R. Zhang J. Ras conformational ensembles, allostery, and signaling Chem. Rev. 116 11 2016 Jun 8 6607 6665 10.1021/acs.chemrev.5b00542 Epub 2016 Jan 27. PMID: 26815308 https://pubmed.ncbi.nlm.nih.gov/26815308/ 26815308
40 Marcus K. Mattos C. Direct attack on RAS: intramolecular communication and mutation-specific effects Clin. Cancer Res. 21 8 2015 Apr 15 1810 1818 10.1158/1078-0432.CCR-14-2148 PMID: 25878362 https://pubmed.ncbi.nlm.nih.gov/25878362/ 25878362
41 Spoerner M. Hozsa C. Poetzl J.A. Reiss K. Ganser P. Geyer M. Kalbitzer H.R. Conformational states of human rat sarcoma (Ras) protein complexed with its natural ligand GTP and their role for effector interaction and GTP hydrolysis J. Biol. Chem. 285 51 2010 Dec 17 39768 39778 10.1074/jbc.M110.145235 Epub 2010 Oct 11. PMID: 20937837; PMCID: PMC3000958 https://pubmed.ncbi.nlm.nih.gov/20937837/ 20937837
42 Fetics S.K. Guterres H. Kearney B.M. Buhrman G. Ma B. Nussinov R. Mattos C. Allosteric effects of the oncogenic RasQ61L mutant on Raf-RBD Structure 23 3 2015 Mar 3 505 516 10.1016/j.str.2014.12.017 Epub 2015 Feb 12. PMID: 25684575; PMCID: PMC7755167 https://pubmed.ncbi.nlm.nih.gov/25684575/ 25684575
43 Johnson C.W. Mattos C. The allosteric switch and conformational states in Ras GTPase affected by small molecules Enzymes 33 Pt A 2013 41 67 10.1016/B978-0-12-416749-0.00003-8 Epub 2013 Aug 8. PMID: 25033800 https://pubmed.ncbi.nlm.nih.gov/25033800/ 25033800
44 Kearney B.M. Johnson C.W. Roberts D.M. Swartz P. Mattos C. DRoP: a water analysis program identifies Ras-GTP-specific pathway of communication between membrane-interacting regions and the active site J. Mol. Biol. 426 3 2014 Feb 6 611 629 10.1016/j.jmb.2013.10.036 Epub 2013 Nov 2. PMID: 24189050 https://pubmed.ncbi.nlm.nih.gov/24189050/ 24189050
45 Parker J.A. Volmar A.Y. Pavlopoulos S. Mattos C. K-ras populates conformational states differently from its isoform H-ras and oncogenic mutant K-RasG12D Structure 26 6 2018 Jun 5 810 820.e4 10.1016/j.str.2018.03.018 Epub 2018 Apr 26. PMID: 29706533 https://pubmed.ncbi.nlm.nih.gov/29706533/ 29706533
46 Johnson C.W. Mattos C. The allosteric switch and conformational states in Ras GTPase affected by small molecules Enzymes 33 Pt A 2013 41 67 10.1016/B978-0-12-416749-0.00003-8 Epub 2013 Aug 8. PMID: 25033800 https://pubmed.ncbi.nlm.nih.gov/25033800/ 25033800
47 Fetics S.K. Guterres H. Kearney B.M. Buhrman G. Ma B. Nussinov R. Mattos C. Allosteric effects of the oncogenic RasQ61L mutant on Raf-RBD Structure 23 3 2015 Mar 3 505 516 10.1016/j.str.2014.12.017 Epub 2015 Feb 12. PMID: 25684575; PMCID: PMC7755167 https://pubmed.ncbi.nlm.nih.gov/25684575/ 25684575
48 Buhrman G. Holzapfel G. Fetics S. Mattos C. Allosteric modulation of Ras positions Q61 for a direct role in catalysis Proc Natl Acad Sci U S A 107 11 2010 Mar 16 4931 4936 10.1073/pnas.0912226107 Epub 2010 Mar 1. PMID: 20194776; PMCID: PMC2841912 https://pubmed.ncbi.nlm.nih.gov/20194776/ 20194776
49 Sharma N. Sonavane U. Joshi R. Probing the wild-type HRas activation mechanism using steered molecular dynamics, understanding the energy barrier and role of water in the activation Eur. Biophys. J. 43 2–3 2014 Mar 81 95 10.1007/s00249-014-0942-4 Epub 2014 Jan 20. PMID: 24442446 https://pubmed.ncbi.nlm.nih.gov/24442446/ 24442446
50 Prakash P. Sayyed-Ahmad A. Gorfe A.A. The role of conserved waters in conformational transitions of Q61H K-ras PLoS Comput. Biol. 8 2 2012 e1002394 10.1371/journal.pcbi.1002394 Epub 2012 Feb 16. PMID: 22359497; PMCID: PMC3280954 https://pubmed.ncbi.nlm.nih.gov/22359497/
51 Pai E.F. Krengel U. Petsko G.A. Goody R.S. Kabsch W. Wittinghofer A. Refined crystal structure of the triphosphate conformation of H-ras p21 at 1.35 A resolution: implications for the mechanism of GTP hydrolysis EMBO J. 9 8 1990 Aug 2351 2359 10.1002/j.1460-2075.1990.tb07409.x PMID: 2196171; PMCID: PMC552258 https://pubmed.ncbi.nlm.nih.gov/2196171/ 2196171
52 Johnson C.W. Reid D. Parker J.A. Salter S. Knihtila R. Kuzmic P. Mattos C. The small GTPases K-Ras, N-Ras, and H-Ras have distinct biochemical properties determined by allosteric effects J. Biol. Chem. 292 31 2017 Aug 4 12981 12993 10.1074/jbc.M117.778886 Epub 2017 Jun 19. PMID: 28630043; PMCID: PMC5546037 https://pubmed.ncbi.nlm.nih.gov/28630043/ 28630043
53 Lu S. Banerjee A. Jang H. Zhang J. Gaponenko V. Nussinov R. GTP binding and oncogenic mutations may attenuate hypervariable region (HVR)-Catalytic domain interactions in small GTPase K-Ras4B, exposing the effector binding site J. Biol. Chem. 290 48 2015 Nov 27 28887 28900 10.1074/jbc.M115.664755 Epub 2015 Oct 9. PMID: 26453300; PMCID: PMC4661403 https://pubmed.ncbi.nlm.nih.gov/26453300/ 26453300
54 Prakash P. Hancock J.F. Gorfe A.A. Binding hotspots on K-ras: consensus ligand binding sites and other reactive regions from probe-based molecular dynamics analysis Proteins 83 5 2015 May 898 909 10.1002/prot.24786 Epub 2015 Mar 25. PMID: 25740554; PMCID: PMC4400267 https://pubmed.ncbi.nlm.nih.gov/25740554/ 25740554
55 Chen C.C. Er T.K. Liu Y.Y. Hwang J.K. Barrio M.J. Rodrigo M. Garcia-Toro E. Herreros-Villanueva M. Computational analysis of KRAS mutations: implications for different effects on the KRAS p.G12D and p.G13D mutations PLoS One 8 2 2013 e55793 10.1371/journal.pone.0055793 Epub 2013 Feb 20. PMID: 23437064; PMCID: PMC3577811 https://pubmed.ncbi.nlm.nih.gov/23437064/
56 Vatansever S, Erman B, Gümüş ZH. Oncogenic G12D mutation alters local conformations and dynamics of K-Ras. Sci. Rep. doi: 10.1038/s41598-019-48029-z. PMID: 31409810; PMCID: PMC6692342.(https://pubmed.ncbi.nlm.nih.gov/31409810/).
57 Grudzien P. Jang H. Leschinsky N. Nussinov R. Gaponenko V. Conformational dynamics allows sampling of an "Active-like" state by oncogenic K-Ras-GDP J. Mol. Biol. 434 17 2022 Sep 15 167695 10.1016/j.jmb.2022.167695 Epub 2022 Jun 23. PMID: 35752212
58 Hobbs G.A. Der C.J. Rossman K.L. RAS isoforms and mutations in cancer at a glance J. Cell Sci. 129 7 2016 Apr 1 1287 1292 10.1242/jcs.182873 Epub 2016 Mar 16. PMID: 26985062; PMCID: PMC4869631 https://pubmed.ncbi.nlm.nih.gov/26985062/ 26985062
59 Ihle N.T. Byers L.A. Kim E.S. Saintigny P. Lee J.J. Blumenschein G.R. Tsao A. Liu S. Larsen J.E. Wang J. Diao L. Coombes K.R. Chen L. Zhang S. Abdelmelek M.F. Tang X. Papadimitrakopoulou V. Minna J.D. Lippman S.M. Hong W.K. Herbst R.S. Wistuba I.I. Heymach J.V. Powis G. Effect of KRAS oncogene substitutions on protein behavior: implications for signaling and clinical outcome J. Natl. Cancer Inst. 104 3 2012 Feb 8 228 239 10.1093/jnci/djr523 Epub 2012 Jan 13. PMID: 22247021; PMCID: PMC3274509 https://pubmed.ncbi.nlm.nih.gov/22247021/ 22247021
60 Hammond D.E. Mageean C.J. Rusilowicz E.V. Wickenden J.A. Clague M.J. Prior I.A. Differential reprogramming of isogenic colorectal cancer cells by distinct activating KRAS mutations J. Proteome Res. 14 3 2015 Mar 6 1535 1546 10.1021/pr501191a Epub 2015 Feb 4. PMID: 25599653; PMCID: PMC4356034 https://pubmed.ncbi.nlm.nih.gov/25599653/ 25599653
61 Abdelkarim H. Leschinsky N. Jang H. Banerjee A. Nussinov R. Gaponenko V. The dynamic nature of the K-Ras/calmodulin complex can be altered by oncogenic mutations Curr. Opin. Struct. Biol. 71 2021 Dec 164 170 10.1016/j.sbi.2021.06.008 Epub 2021 Jul 24. PMID: 34311289 34311289
62 Maurer T. Garrenton L.S. Oh A. Pitts K. Anderson D.J. Skelton N.J. Fauber B.P. Pan B. Malek S. Stokoe D. Ludlam M.J. Bowman K.K. Wu J. Giannetti A.M. Starovasnik M.A. Mellman I. Jackson P.K. Rudolph J. Wang W. Fang G. Small-molecule ligands bind to a distinct pocket in Ras and inhibit SOS-mediated nucleotide exchange activity Proc Natl Acad Sci U S A 109 14 2012 Apr 3 5299 5304 10.1073/pnas.1116510109 Epub 2012 Mar 19. PMID: 22431598; PMCID: PMC3325706 https://pubmed.ncbi.nlm.nih.gov/22431598/ 22431598
63 Pettersen E.F. Goddard T.D. Huang C.C. Couch G.S. Greenblatt D.M. Meng E.C. Ferrin T.E. UCSF Chimera--a visualization system for exploratory research and analysis J. Comput. Chem. 25 13 2004 Oct 1605 1612 10.1002/jcc.20084 PMID: 15264254 https://pubmed.ncbi.nlm.nih.gov/15264254/ 15264254
64 Abraham M.J. Murtola T. Schulz R. Páll S. Smith J.C. Hess B. Lindah E. Gromacs: high performance molecular simulations through multi-level parallelism from laptops to supercomputers SoftwareX 1–2 2015 19 25 10.1016/j.softx.2015.06.001
65 Hornak V. Abel R. Okur A. Strockbine B. Roitberg A. Simmerling C. Comparison of multiple Amber force fields and development of improved protein backbone parameters Proteins 65 3 2006 Nov 15 712 725 10.1002/prot.21123 PMID: 16981200; PMCID: PMC4805110 https://pubmed.ncbi.nlm.nih.gov/16981200/ 16981200
66 Schüttelkopf A.W. van Aalten D.M. PRODRG: a tool for high-throughput crystallography of protein-ligand complexes Acta Crystallogr D Biol Crystallogr 60 Pt 8 2004 Aug 1355 1363 10.1107/S0907444904011679 Epub 2004 Jul 21. PMID: 15272157 https://pubmed.ncbi.nlm.nih.gov/15272157/ 15272157
67 Bussi G. Donadio D. Parrinello M. Canonical sampling through velocity rescaling J. Chem. Phys. 126 1 2007 Jan 7 014101 10.1063/1.2408420 PMID: 17212484 https://pubmed.ncbi.nlm.nih.gov/17212484/
68 Berendsen H.J. Postma J.V. Van Gunsteren W.F. DiNola A.R. Haak J.R. Molecular dynamics with coupling to an external bath J. Chem. Phys. 81 8 1984 Oct 15 3684 3690 10.1063/1.448118
69 Essmann U. Perera L. Berkowitz M.L. Darden T. Lee H. Pedersen L.G. A smooth particle mesh Ewald method J. Chem. Phys. 103 19 1995 Nov 15 8577 8593 10.1063/1.470117
70 Hess B. P-LINCS: a parallel linear constraint solver for molecular simulation J. Chem. Theor. Comput. 4 1 2008 Jan 116 122 10.1021/ct700200b PMID: 26619985 https://pubmed.ncbi.nlm.nih.gov/26619985/
71 Parrinello M. Rahman A. Polymorphic transitions in single crystals: a new molecular dynamics method J. Appl. Phys. 52 12 1981 Dec 7182 7190 10.1063/1.328693
72 https://getcontacts.github.io/.
73 Humphrey W. Dalke A. Schulten K. VMD: visual molecular dynamics J. Mol. Graph. 14 1 1996 Feb 33 38 10.1016/0263-7855(96)00018-5 PMID: 8744570 https://pubmed.ncbi.nlm.nih.gov/8744570/ 8744570
74 Hünenberger P.H. Mark A.E. van Gunsteren W.F. Fluctuation and cross-correlation analysis of protein motions observed in nanosecond molecular dynamics simulations J. Mol. Biol. 252 4 1995 Sep 29 492 503 10.1006/jmbi.1995.0514 PMID: 7563068 https://pubmed.ncbi.nlm.nih.gov/7563068/ 7563068
75 Scherer M.K. Trendelkamp-Schroer B. Paul F. Pérez-Hernández G. Hoffmann M. Plattner N. Wehmeyer C. Prinz J.H. Noé F. PyEMMA 2: a software package for estimation, validation, and analysis of Markov models J. Chem. Theor. Comput. 11 11 2015 Nov 10 5525 5542 10.1021/acs.jctc.5b00743 Epub 2015 Oct 14. PMID: 26574340 https://pubmed.ncbi.nlm.nih.gov/26574340/
76 Noé F. Schütte C. Vanden-Eijnden E. Reich L. Weikl T.R. Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations Proc Natl Acad Sci U S A 106 45 2009 Nov 10 19011 19016 10.1073/pnas.0905466106 Epub 2009 Nov 3. PMID: 19887634; PMCID: PMC2772816 https://pubmed.ncbi.nlm.nih.gov/19887634/ 19887634
77 Chodera J.D. Noé F. Markov state models of biomolecular conformational dynamics Curr. Opin. Struct. Biol. 25 2014 Apr 135 144 10.1016/j.sbi.2014.04.002 Epub 2014 May 16. PMID: 24836551; PMCID: PMC4124001 https://pubmed.ncbi.nlm.nih.gov/24836551/ 24836551
78 Chodera J.D. Singhal N. Pande V.S. Dill K.A. Swope W.C. Automatic discovery of metastable states for the construction of Markov models of macromolecular conformational dynamics J. Chem. Phys. 126 15 2007 Apr 21 155101 10.1063/1.2714538 PMID: 17461665 https://pubmed.ncbi.nlm.nih.gov/17461665/
79 Prinz J.H. Wu H. Sarich M. Keller B. Senne M. Held M. Chodera J.D. Schütte C. Noé F. Markov models of molecular kinetics: generation and validation J. Chem. Phys. 134 17 2011 May 7 174105 10.1063/1.3565032 PMID: 21548671 https://pubmed.ncbi.nlm.nih.gov/21548671/
80 Riely G.J. Marks J. Pao W. KRAS mutations in non-small cell lung cancer Proc. Am. Thorac. Soc. 6 2 2009 Apr 15 201 205 10.1513/pats.200809-107LC PMID: 19349489 https://pubmed.ncbi.nlm.nih.gov/19349489/ 19349489
81 Ostrem J.M. Shokat K.M. Direct small-molecule inhibitors of KRAS: from structural insights to mechanism-based design Nat. Rev. Drug Discov. 15 11 2016 Nov 771 785 10.1038/nrd.2016.139 Epub 2016 Jul 29. PMID: 27469033 https://pubmed.ncbi.nlm.nih.gov/27469033/ 27469033
82 Ito Y. Yamasaki K. Iwahara J. Terada T. Kamiya A. Shirouzu M. Muto Y. Kawai G. Yokoyama S. Laue E.D. Wälchli M. Shibata T. Nishimura S. Miyazawa T. Regional polysterism in the GTP-bound form of the human c-Ha-Ras protein Biochemistry 36 30 1997 Jul 29 9109 9119 10.1021/bi970296u PMID: 9230043 https://pubmed.ncbi.nlm.nih.gov/9230043/ 9230043
83 Chavan T.S. Jang H. Khavrutskii L. Abraham S.J. Banerjee A. Freed B.C. Johannessen L. Tarasov S.G. Gaponenko V. Nussinov R. Tarasova N.I. High-affinity interaction of the K-Ras4B hypervariable region with the ras active site Biophys. J. 109 12 2015 Dec 15 2602 2613 10.1016/j.bpj.2015.09.034 PMID: 26682817; PMCID: PMC4699860 https://pubmed.ncbi.nlm.nih.gov/26682817/ 26682817
84 Jang H. Banerjee A. Chavan T.S. Lu S. Zhang J. Gaponenko V. Nussinov R. The higher level of complexity of K-Ras4B activation at the membrane Faseb. J. 30 4 2016 Apr 1643 1655 10.1096/fj.15-279091 Epub 2015 Dec 30. Erratum in: FASEB J. 2018 Sep;32(9):5209. PMID: 26718888; PMCID: PMC4799498 https://pubmed.ncbi.nlm.nih.gov/26718888/ 26718888
85 Chen C.C. Er T.K. Liu Y.Y. Hwang J.K. Barrio M.J. Rodrigo M. Garcia-Toro E. Herreros-Villanueva M. Computational analysis of KRAS mutations: implications for different effects on the KRAS p.G12D and p.G13D mutations PLoS One 8 2 2013 e55793 10.1371/journal.pone.0055793 Epub 2013 Feb 20. PMID: 23437064; PMCID: PMC3577811 https://pubmed.ncbi.nlm.nih.gov/23437064/
86 Khaled M. Gorfe A. Sayyed-Ahmad A. Conformational and dynamical effects of Tyr32 phosphorylation in K-ras: molecular dynamics simulation and Markov state models analysis J. Phys. Chem. B 123 36 2019 Sep 12 7667 7675 10.1021/acs.jpcb.9b05768 Epub 2019 Aug 30. PMID: 31419909; PMCID: PMC7020251 https://pubmed.ncbi.nlm.nih.gov/31419909/ 31419909
87 Grant B.J. Gorfe A.A. McCammon J.A. Ras conformational switching: simulating nucleotide-dependent conformational transitions with accelerated molecular dynamics PLoS Comput. Biol. 5 3 2009 Mar e1000325 10.1371/journal.pcbi.1000325 Epub 2009 Mar 20. PMID: 19300489; PMCID: PMC2651530 https://pubmed.ncbi.nlm.nih.gov/19300489/
88 Sayyed-Ahmad A. Prakash P. Gorfe A.A. Distinct dynamics and interaction patterns in H- and K-Ras oncogenic P-loop mutants Proteins 85 9 2017 Sep 1618 1632 10.1002/prot.25317 Epub 2017 May 31. PMID: 28498561; PMCID: PMC5568977 https://pubmed.ncbi.nlm.nih.gov/28498561/ 28498561
89 Ma J. Karplus M. Molecular switch in signal transduction: reaction paths of the conformational changes in ras p21 Proc Natl Acad Sci U S A 94 22 1997 Oct 28 11905 11910 10.1073/pnas.94.22.11905 PMID: 9342335; PMCID: PMC23651 https://pubmed.ncbi.nlm.nih.gov/9342335/ 9342335
90 Ostrem J.M. Peters U. Sos M.L. Wells J.A. Shokat K.M. K-Ras(G12C) inhibitors allosterically control GTP affinity and effector interactions Nature 503 7477 2013 Nov 28 548 551 10.1038/nature12796 Epub 2013 Nov 20. PMID: 24256730; PMCID: PMC4274051 24256730
91 Zhao H. Li L. Liu J. Mai R. Chen J. Chen J. Discovery of ARS-1620 analogs as KRas G12C inhibitors with high in vivo antitumor activity Bioorg. Chem. 121 2022 Apr 105652 10.1016/j.bioorg.2022.105652 Epub 2022 Feb 10. PMID: 35158284
92 Li L. Liu J. Yang Z. Zhao H. Deng B. Ren Y. Mai R. Huang J. Chen J. Discovery of Thieno[2,3-d]pyrimidine-based KRAS G12D inhibitors as potential anticancer agents via combinatorial virtual screening Eur. J. Med. Chem. 233 2022 Apr 5 114243 10.1016/j.ejmech.2022.114243 Epub 2022 Mar 3. PMID: 35276423
93 Kwan A.K. Piazza G.A. Keeton A.B. Leite C.A. The path to the clinic: a comprehensive review on direct KRASG12C inhibitors J. Exp. Clin. Cancer Res. 41 1 2022 Jan 19 27 10.1186/s13046-021-02225-w PMID: 35045886; PMCID: PMC8767686 35045886
94 Wang X. Allen S. Blake J.F. Bowcut V. Briere D.M. Calinisan A. Dahlke J.R. Fell J.B. Fischer J.P. Gunn R.J. Hallin J. Laguer J. Lawson J.D. Medwid J. Newhouse B. Nguyen P. O'Leary J.M. Olson P. Pajk S. Rahbaek L. Rodriguez M. Smith C.R. Tang T.P. Thomas N.C. Vanderpool D. Vigers G.P. Christensen J.G. Marx M.A. Identification of MRTX1133, a noncovalent, potent, and selective KRASG12D inhibitor J. Med. Chem. 65 4 2022 Feb 24 3123 3133 10.1021/acs.jmedchem.1c01688 Epub 2021 Dec 10. PMID: 34889605 34889605
95 Wang H. Chi L. Yu F. Dai H. Gao C. Si X. Wang Z. Liu L. Zheng J. Shan L. Liu H. Zhang Q. Annual review of KRAS inhibitors in 2022 Eur. J. Med. Chem. 249 2023 Mar 5 115124 10.1016/j.ejmech.2023.115124 Epub 2023 Jan 16. PMID: 36680986
96 Parikh K. Banna G. Liu S.V. Friedlaender A. Desai A. Subbiah V. Addeo A. Drugging KRAS: current perspectives and state-of-art review J. Hematol. Oncol. 15 1 2022 Oct 25 152 10.1186/s13045-022-01375-4 PMID: 36284306; PMCID: PMC9597994 36284306
97 Punekar S.R. Velcheti V. Neel B.G. Wong K.K. The current state of the art and future trends in RAS-targeted cancer therapies Nat. Rev. Clin. Oncol. 19 10 2022 Oct 637 655 10.1038/s41571-022-00671-9 Epub 2022 Aug 26. PMID: 36028717; PMCID: PMC9412785 36028717
98 Bowman G.R. Geissler P.L. Equilibrium fluctuations of a single folded protein reveal a multitude of potential cryptic allosteric sites Proc. Natl. Acad. Sci. U.S.A. 109 29 2012 Jul 17 11681 11686 10.1073/pnas.1209309109 Epub 2012 Jul 2. PMID: 22753506; PMCID: PMC3406870 22753506
99 Oleinikovas V. Saladino G. Cossins B.P. Gervasio F.L. Understanding cryptic pocket formation in protein targets by enhanced sampling simulations J. Am. Chem. Soc. 138 43 2016 Nov 2 14257 14263 10.1021/jacs.6b05425 Epub 2016 Oct 20. PMID: 27726386 27726386
