THE EFFECTS OF ONCOGENIC G12D MUTATION ON K-RAS STRUCTURE, CONFORMATION AND DYNAMICS

K-Ras is the most frequently mutated protein in human tumors. Activating K-Ras mutations drive cancer initiation, progression and drug resistance, directly leading to nearly a million deaths per year. To understand the mechanisms by which mutations alter K-Ras function, we need to understand their effects on protein dynamics. However, despite decades of research, how oncogenic mutations in K-Ras alter its conformation and dynamics remain to be understood. Here, we present how the most recurrent K-Ras oncogenic mutation, G12D, leads to structural, conformational and dynamical changes that lead to constitutively active KRas. We have developed a new integrated MD simulation data analysis approach to quantify such changes in a protein and applied it to K-Ras. Our results show that G12D mutation induces strong negative correlations between the fluctuations of SII and those of the P-loop, Switch I (SI) and α3 regions in K-Ras G12D . Furthermore, characteristic decay times of SII fluctuations significantly increase after G12D mutation. We have further identified causal relationships between correlated residue pairs in K-Ras G12D and show that the correlated motions in K-Ras dynamics are driven by SII fluctuations, which have the strongest negative correlations with other protein parts and the longest characteristic decay times in mutant KRas. Ours is arguably the first study that shows the causal relationships between residue pairs in K-Ras G12D , relates them to the decay times and correlates their fluctuations.


INTRODUCTION
K-Ras is the most frequently mutated oncoprotein in human cancers [1][2][3] . Patients with oncogenic K-Ras mutations have very poor response to standard therapies. Unfortunately, K-Ras mutations eventually emerge during the course of their treatment and drive acquired resistance [4][5][6] . K-Ras functions as a small GTPase and is critical in the regulation of intracellular signaling networks in cellular growth, proliferation and differentiation 7 . To perform its cellular roles, it switches between its inactive GDP-bound and active GTP-bound states 8,9 . Only active K-Ras (K-Ras-GTP) can bind and activate its downstream effector proteins 10 . Active K-Ras catalyzes GTP hydrolysis to become inactive (K-Ras-GDP). Intrinsic GTPase function of active K-Ras can be accelerated by the binding of GTPase-activating proteins (GAPs) 11,12 . However, certain mutations in K-Ras impair its intrinsic GTPase function and GAP binding and thereby GTP hydrolysis. Unable to switch to its GDP-bound inactive state, mutant K-Ras remains continuously active, causing prolonged activation of downstream pathways associated with oncogenic cell growth 10,[13][14][15] .
In cancer patients, oncogenic K-Ras mutations are recurrently observed at positions 12, 13 and 61. G12 is the most frequently mutated residue (89%), which most prevalently mutates to aspartate (G12D, 36%) followed by valine (G12V, 23%) and cysteine (G12C, 14%) 3,10 . This residue is located at the protein active site, which consists of a phosphate binding loop (P-loop, residues 10-17) and switch I (SI, residues [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] and II (SII, residues 60-74) regions. The active site residues are bound to the phosphate groups of GTP and are responsible for the GTPase function of K-Ras. In its side-chain, G12 has only a single hydrogen. However, the mutation to aspartate (G12D) leads to the projection of a bulkier side group into the active site, which causes a steric hindrance in GTP hydrolysis 16 , impairs the GTPase function and locks K-Ras in its active GTP-bound state 12 . There is strong evidence that blocking mutant K-Ras activity can be very effective in the treatment of cancer patients 17,18 . Yet, despite decades of research, there are still no drugs in the clinic today that can directly target mutant K-Ras 19,20 . Although the effects of G12D mutation on the structure, conformation and flexibility of K-Ras have been studied [21][22][23][24] , the relationship between its conformational and dynamical changes still remains to be understood. At the same time, there is increasing evidence that suggests that crystal structure studies alone may miss drugbinding pockets on mutant K-Ras surface 18,[25][26][27][28][29][30][31][32][33] . Studies that include dynamics information have recently achieved promising results, however, these are limited to K-Ras G12C mutant, and the mechanisms that regulate its dynamics remain unknown. Understanding the mechanisms of dynamic regulation of mutant K-Ras can present novel opportunities for identifying clinically actionable regulatory sites on its surface, making a profound impact in the treatment of millions of patients. mutation-based changes in protein structure and its consequences. Our approach is based on the conditional time-delayed correlations (CTC) method we have recently developed, which uses time delayed correlation functions from MD datasets to provide information on the correlation of two events, one taking place at time t and the other at a later time, t+. This approach has been remarkably accurate in predicting sites that control a protein's motions 45 , and has shown excellent agreement with experimental data when tested on wild-type K-Ras (K-Ras WT ).
We first performed MD-simulations of GTP-bound active forms of both wild type and G12D mutant K-Ras (hereafter K-Ras WT and K-Ras G12D ) and applied our CTC method to predict sites on mutant K-Ras G12D surface that control its activity. From the CTC results, we first studied the structural changes in both GTP and GDP-bound K-Ras upon G12D mutation, and discovered salt bridges that are either formed or destroyed upon mutation. Second, we evaluated the changes in the pair-wise distances between residues and quantified the local volume changes to identify changes in protein conformation. Third, we identified changes in protein dynamics through a multi-step process where we (i) quantified the residue fluctuations; (ii) evaluated correlation of residue fluctuations and identified lost or newly formed correlations upon mutation; (iii) calculated the characteristic decay times of residue fluctuations; and (iv) identified residue pairs that are causally related (i.e., residue pairs that show driver-responder behavior). Finally, we related the observed structural changes to the conformational and dynamical alterations, which enabled us to identify the important changes that affect protein function. Overall, our study identifies regulatory sites on K-Ras G12D , which enhance our understanding of its dynamics and can assist in the development of direct inhibitors.

MD SIMULATIONS
We performed all-atom MD simulations for both Mg +2 GTP-bound K-Ras WT and K-Ras G12D . We obtained the K-Ras-GTP WT structure from the final frame of the 300ns simulation of active state protein by Vatansever et al 45 . For constructing K-Ras G12D structure, we mutated glycine to aspartate at position 12 in K-Ras-GTP WT structure using Discovery Studio 4.5 software, (DS) 46 . To optimize the K-Ras G12D -GTP complex, we used Clean Geometry tool of DS. For MD simulations, we used NAMD 2.10 47 with AMBER ff99SB 48 and general amber force fields (GAFF) 49 . Briefly, we performed energy minimization of the initial model after we introduced the G12D mutation in K-Ras, and then ran MD simulations of each complex following the protocols from our previous study 45 , the details of which we provide in Supplementary Methods. During the simulations, we applied minimization for 10,000 steps and equilibration for 500,000 steps, after which we performed 1 microsecond MD simulations, and saved atomic coordinates ̂ of all atoms every 10ps. we used the last 900ns of the simulation trajectories in all computations described in this study. To eliminate all rotational and translational motions, we aligned the trajectories to the first frame using VMD software 1.9.2 50 . We visualized the trajectories with VMD. To identify salt bridges formed in the protein during the MD simulations, we used Salt Bridges Plugin, Version 1.1, of VMD.
To quantify the effect of the G12D mutation on the distances between K-Ras residue pairs, we developed a new computational algorithm detailed in Figure S1. Briefly, we first assumed K-Ras WT as the initial state and K-Ras G12D as the final state. Then, we calculated the distances between Cα atoms of two residues (i, j) as we previously described 45 . A widely used model to analyze protein dynamics is Gaussian network model (GNM). Studies that use GNM typically assume the maximum Cα-Cα distance for the separation between two contacting residues at ~7.2Å 51,52 , and label it the 'first coordination shell' 53,54 . We followed this same protocol and determined the first coordination shell around a selected residue by choosing its Cα as the center of a volume V with a radius of r1 ~7.2 Å 52 . However, because the contribution of non-bonded pairs to higher-order coordination shells may also be significant 54,55 , we also studied residue pairs that are within their 'second coordination shell' in K-Ras WT structure. We defined this second coordination shell at twice the volume of the first, with a radius of ~9.1 Å 54 .
For every residue pair (i, j) where j is in the second coordination shell of i, we first calculated its time-averaged distance in K-Ras WT ( ̅ ij WT ) and in K-Ras G12D ( ̅ ij G12D ). We then calculated the difference (∆ ̅ ) between ̅ ij WT and ̅ ij G12D , where ∆ ̅ = ̅ ij G12D -̅ ij WT . The magnitude of the difference is the degree of distortion resulting from the G12D mutation. We present ∆ ̅ values in the pairwise distances map ( Figure 1) where a positive value indicates that a residue pair moves apart upon G12D mutation, while a negative value indicates that the pair gets closer. Then, we identified the residue pairs (ij) significantly distorted by the G12D mutation. For this purpose, we selected the residue pairs that have the greatest (positive and negative) ∆ ̅ values. We assumed that the residue pairs whose ∆ ̅ values are greater than 2.75 or smaller than -1.35 showed the most significant distance changes. For those identified residue pairs, we drew the distribution graphs W(R ij ) of their distances (R ij ) during the simulations of K-Ras WT and K-Ras G12D (For details see Supplementary).
Then, to quantify the changes in local volumes upon G12D mutation, for each residue i we calculated the average of all ∆ ̅ ij values based on the formula 〈∆ ̅ 〉 = ∑ ∆ ̅ ⁄ , where N n is the number of residues j in the second coordination shell of residue i. In detail, for a residue i, at the center of a volume V with a radius of 9.1 Å (the second coordination shell) and we defined the residues j within this volume V as the neighbors of residue i. Then, we calculated the total change in the distance between residue i and its neighbors,∑ ∆ ̅ , and divided it by the number of neighbors ∑ ∆ ̅ ⁄ . The resulting 〈∆ ̅ 〉 value is a measure of the change in volume around residue i due to G12D mutation.

PAIRWISE CORRELATION, TIME DELAYED CORRELATION AND CHARACTERISTIC DECAY TIME CALCULATIONS
To show coupled motions in protein dynamics, we calculated the correlation coefficients between the fluctuations of residue pairs (C ij ). A correlation coefficient value of a residue pair ranges from -1 to 1, where for residue pair fluctuations that are not coupled Cij= 0; perfectly positively correlated Cij=1, and perfectly negatively correlated Cij=-1. Additionally, to identify the directionality (causality) in the correlated motions of residue pairs, we calculated the timedelayed correlations between them (C ij (τ)), where τ is the time-delay. C ij (τ) value represents the degree and manne in which the past fluctuations of residue i and the present fluctuations of residue j are coupled.
We further calculated autocorrelations C ij (τ) of each protein residue, which is the coupling of its past and present fluctuations. A slow decay of C ii (τ) function indicates that the present fluctuations of residue i tend to be strongly coupled with its own past fluctuations and thereby that residue i has a strong memory of its past fluctuations. To measure the memory length for a residue i, we calculated the time-delay value (τ) where the autocorrelation decay curve of the residue i reaches to 0.5 and we accepted it as the characteristic decay time of the residue i. We adopted the 0.5 criterion instead of the commonly used 1/e decay criterion because it was a better indicator of short time relaxation differences while the latter led to values too close to each other.
The the calculations of correlation coefficients, time-delayed correlations and autocorrelations were as described in our previous study 45 and are summarized in detail in Supplementary.

1.Structural Changes
New close-range electrostatic interactions (salt bridges) are formed in K-Ras G12D ., which are not present in K-Ras WT . Substitution of a non-polar, non-charged amino acid (glycine) with a negatively-charged amino acid (aspartate) triggers several conformational and dynamical changes in K-Ras G12D . With the plausible assumption that the negatively charged residue D12 may cause new electrostatic interactions within the protein and that those interactions can be the sources of conformational changes upon G12D mutation, we compared the close-range electrostatic interactions (i.e. salt bridges) in K-RAS WT vs. K-RAS G12D . In the mutated structure, D12 forms salt bridges with K16 (P-loop) and K88, and K16 forms a salt bridge with D57. However, none of these interactions are present in K-RAS WT . To identify the effects of the new electrostatic interactions caused by G12D mutation, we then investigated the conformational changes in K-Ras.

2.Conformational Changes
Pairwise Distance Calculations show conformational changes in K-Ras G12D vs. K-Ras WT . We compared the distances between residue pairs within their second coordination shell in K-Ras WT vs. K-Ras G12D , where the first coordination shell is the maximum Cα-Cα distance of the separation between two contacting residues ~7.2Å 51-54 and the second coordination shell is twice the volume of the first, at a radius of ~9.1Å 54 . In Figure 1A, we show the ∆ ̅ ij values for all residue pairs where K-Ras WT is the reference and K-Ras G12D is the final structure. As seen from the abundance of positive ∆ ̅ ij values, the dominant distortion of the protein upon mutation is expansion, where SII region (T58-T74) moves away from the phosphate binding loop (P-loop, G10-S17), SI (Q25-Y40), α3 (T87-K104) in K-Ras G12D . The SI region also moves away from the P-loop. In detail, all the P-loop residues move away from T58-R68 (SII) residues as a block. Second, the residues in the segment P34-Y40 of SI move away from N-terminal residues of SII (D57-Q61). Third, the α2 helix of SII moves away from α3, D92-R102. Fourth, SI residues move away from the P-loop residues. On the other hand, SI residues (F28-Y32, Y32-I36) assume a closer conformation where an H-bond between D33 and I36 is established. This causes the T35-E37 part of SI to move away from the Ploop residues, G13-K16. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. We observed the most significant changes in the distances between residue pairs Q61 (SII)-A11 (P-loop), Q61 (SII)-G12D (P-loop), Q61 (SII)-G13 (P-loop), E37(SI)-D57 (SII) and S65 (SII)-H95(α3) ( Figure 1B). Residues R68, M72 and Y96 were also further apart in K-Ras G12D . On the other hand, we observed that some residue pairs got significantly closer in K-Ras G12D , including E63-R68, M72-G75 and R73-K104, demonstrated by their negative ∆ ̅ ij values in Figure 1A. K-Ras G12D has broader distance distributions than K-Ras WT . To uncover the conformational differences between wild type and mutant K-Ras and to thereby better understand the effects of the G12D mutation, we calculated the probability distributions of the distances between pairs of residues that exhibited the largest changes in the distance calculations. Distance distributions of the residue pairs which underwent the largest changes due to G12D mutation are shown in Figure 2. These are between the alpha carbons of residue pairs A11-Q61, G12D-Q61 and G13-Q61, respectively, as shown in Fig. 2A-C. Their distribution patterns exhibit marked differences between WT and mutant K-Ras. In WT K-Ras, they exhibit a narrow distance distribution with smaller peaks compared to the mutant form. Residues A11, G12 and G13 are located at the P-loop. The P-loop has an omega shape and forms a turn at the C-terminal neighborhood of G12. As shown in Figure 2A-B, A11 and G12 are located within H-bond distance from Q61, while the distance between G13 and Q61 is increased. Furthermore, we observed that in K-Ras WT , G12 (side chain H atom) forms an H-bond with G60 (backbone O atom), the neighbor of G60. However, in K-Ras G12D , this H-bond disappears because of the bulkier side chain of D12. In the absence of this H-bond, D12 and its neighbors move away from Q61. The broadened distance distribution values of the mutant protein are indicators of lost hydrogen bonds.
In K-Ras WT simulations, we observed an H-bond between D38 (SI) and D57 (SII). However, this H-bond dissappers in K-Ras G12D simulations, leading D57 to move away from D38 and more remarkably from E37. The absence of this H-bond in K-Ras G12D is shown in Figure 2D where the peak distance value is greater in the K-Ras G12D .
Among all residue-residue pairs, S65-H95 undergoes the largest conformational change upon mutation. S65 is on α2 and moves away from H95, which is on α3 in K-Ras G12D ( Figure  2E). The simulations show that this significant change due to the brakeage of the salt bridge between α2 and α3 helices. In K-Ras WT , these two helices interact through the salt bridge between R68 and D92. However, the G12D mutation breaks the R68-D92 salt bridge and causes the α2 and α3 helices to move away from each other. This conformation of α2 and α3 in K-Ras G12D can be deduced from Figure 2E, which shows that the distance distribution curve of S65 (α2)-H95 (α3) has a higher peak. Finally, the distances between the residues R68, M72 and Y96 are also more stable in K-Ras WT as shown in Figure 2F-H.
Overall, we observed that distance distribution curves for K-Ras WT are characterized by Gaussian-shaped, narrow dispersion curves. We use the term 'stable' in the sense that the distribution has a sharp peak, with a small dispersion around it. However, mutant K-Ras showed significant deviations from the Gaussian, except for the residue pair E37(SI)-D57 (SII).
We show the distance distributions of the residue pairs which became significantly closer upon mutation in Figure S2. During the MD simulations, E63-R68 pair is in a closer conformation in K-Ras G12D ( Figure S2A). This conformation may have resulted from the residues adjacent to E63, Y64 and S65, making H-bonds with R68 and D69, respectively. We observed these H-bonds between the backbone H atoms of Y64 & S65 and the backbone O atoms of R68 & D69 only in K-Ras G12D . Additionally, during the K-Ras G12D simulations R73-K104 pair also assumed a closer conformation which may have resulted from the H-bonds between R73-D105 and G75-K104 ( Figure S2B).
On the other hand, in Figure S2C, we show that the residue pair M72-G75 switches between two conformations in K-Ras WT , where one has a closer conformation around 5.5 Å and the other one has a distant conformation around 8.5 Å. However, their distance distribution has a single peak at 5.5 Å in K-Ras G12D that is similar to the closer conformation in K-Ras WT . We have observed that the residue M72 forms an H-bond with G75 only in K-Ras G12D . This Hbond may have caused their single-peaked distance distribution. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Figure 2. Distribution of distances between residue pairs in K-Ras WT (black) and K-Ras G12D (red). Distance distribution of residue pairs
The copyright holder for this preprint (which was this version posted August 19, 2017. ; https://doi.org/10.1101/178483 doi: bioRxiv preprint G12D mutation causes changes in local volume. According to the GNM model, a residue typically fluctuates within its first or second coordination shell 51,52 . Within these fluctuation volumes, there are several other residues, which are either near-neighbors along the chain or are spatially distant. As has been shown by the GNM model, each residue has a different number of neighbors 52 . A residue with a smaller number of neighbors will show larger fluctuations than a residue with a larger number of neighbors. Therefore, the neighborhood of a given residue significantly affects its fluctuations, and as we will show here there is also a significant effect on its dynamics. In order to understand which parts of K-Ras move away from its neighbors and which parts move closer upon mutation, we have calculated the average of all ∆ ̅ ij values for each residue, 〈∆ ̅ 〉. Figure 1C shows that most of the protein parts, especially the P-loop and SII, move away from their neigbors after G12D mutation. In summary, we identified the residue pairs that underwent the largest distance change due to G12D mutation. In addition, for each residue in the identified pairs, we showed the extent of its deviation from its neighbors. For this purpose, we have compared the individual ∑ ∆ ̅ values of the residues in the identified pairs (Table S1). We discovered that the residues which move further away from each other in K-Ras G12D also move away from their neighbors. On the other hand, the residues that move closer to each other in K-Ras G12D have different conformations relative to their neighbors. G12D . To understand how the flexibility of K-Ras changes upon G12D mutation, we calculated the rootmean-square fluctuations (RMSF) of each residue in both wild-type and mutant protein, where RMSF is a measure of the average atomic fluctuations of a residue. In Figure 3A, we display our results, which show that the fluctuations of central residues of SII are increased in K-Ras G12D .

The correlated motions of residues are markedly increased in K-Ras G12D in comparison to K-Ras WT . Regulation of protein dynamics is strictly coordinated by the correlations of residue fluctuations. Figures 3C-D present pairwise correlations of residue fluctuations ( ),
where we show the results for K-Ras WT and K-Ras G12D in the left and right panels, respectively. Specifically, β3-SII residues become negatively correlated with the residues of P-loop, SI, and β4-α3.

The characteristic decay times of residues in P-loop, SI and SII are longer in K-Ras G12D .
Fluctuations of each residue has a characteristic decay time that corresponds to the memory of its past 56 . In other words, characteristic decay time is the time-delay for which the present fluctuations of residue i becomes decoupled from its past fluctuations. Figure 3B shows that the characteristic decay times of residue fluctuations. Moreover, SII residues show the longest correlation decay times within the K-Ras G12D residues.  G12D . We analyzed mutant K-Ras dynamics in depth by applying time-delayed correlation analysis to MD simulation data of K-Ras G12D . The time-delayed correlation of residues i and j ( ( )) is the correlation of the fluctuations of residue i with the later fluctuations of residue j. Through this analysis, we identified driver-responder residue pairs in K-Ras G12D motions. Then, we compared the driver-responder pairs for K-Ras G12D that we have calculated here with those for K-Ras WT that we have previously published 45 .

SII motions are the main drivers of correlated motions in K-Ras
Our causality analysis calculations show that in K-Ras G12D , the motions of SII are the main drivers, driving the motions of both P-loop and β3. We present the time delayed correlation plots of Q70 (SII) with V9 (P-loop) ( Figure 4A) and C80 (β4) ( Figure 4B) and D69 (SII) with not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted August 19, 2017. ; https://doi.org/10.1101/178483 doi: bioRxiv preprint D54 (β3) ( Figure 4C) for K-Ras G12D . These are newly formed driver-responder interactions which do not exist in K-Ras WT .  G12D . Our results show that upon K-Ras G12D mutation, salt bridges are formed between the sidechain atoms of the residues D12 (P-loop)-K16 (P-loop) and K16-D57 (SII). Furthermore, the SII region moves away from the P-loop, and they both move away from their neighbors. MD simulations of K-Ras G12D show that the sidechains of the K16-D57 residue pair approach each other and form salt bridges, while their backbone Cα atoms move away from each other ( Figure S1). Therefore, the new salt bridges formed between the P-loop and SII residues may cause the conformations of these two regions to move away from each other and from their neighbors. As a consequence of these structural and conformational changes in SII region, SII also moves away from SI and α3 regions. G12D . Since K-Ras G12D simulations show the formation of a salt bridge between of D12 (P-loop) and K16 (P-loop), we investigated whether D12 and K16 moving closer due to the salt bridge causes changes in distal regions of K-Ras G12D . For this purpose, we first defined a connectivity vector, ΔR 12→16 , between D12 and K16 based on the definition in our previous paper 45 , as summarized in Supplementary. We then defined connectivity vectors between the correlated residue pairs in K-Ras G12D . Finally, we calculated the correlations of ΔR 12→16 with the other connectivity vectors. We discovered that that ΔR 12→16 is significantly correlated with ΔR 60→70 , ΔR 61→75 and ΔR 60→82 . Consequently, D12-K16 pair moving closer as a result of salt bridge formation affects the dynamics of distant residue pairs such as G60-Q70, Q61-G75 and G60-F82.

Relationships Between Conformational and Dynamic Changes in K-Ras G12D
Fluctuations of dilated regions become negatively correlated in K-Ras G12D . We observed that not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted August 19, 2017. ; https://doi.org/10.1101/178483 doi: bioRxiv preprint after G12D mutation, negative correlations occur between the regions that move away from each other and from their neighbors. Combination of distance and correlation calculations gives us the relation between conformational and dynamic changes in the protein as a result of the mutation in its structure ( Figure 5).  G12D . Our results show that upon K-Ras G12D mutation, the fluctuations of the SII region become negatively correlated with fluctuations of the P-loop, SI and α3 regions. Moreover, fluctuations of these regions have the longest correlation decay times in the protein ( Figure  6). To investigate whether increased negative correlations between the residue fluctuations slow down the autocorrelation decay of the residue fluctuations, we calculated the average of negative correlation values for each residue of K-Ras G12D . In Figure 7, we show that the residues whose fluctuations are more negatively correlated with other parts of the protein also have longer characteristic decay times.

Relationships Between RMSF Values and Characteristic Decay Times in K-RAS G12D .
Autocorrelations of residue fluctuations tend to decay slowly as fluctuation magnitudes increase. Since our results showed that the residues in the SII region have the longest characteristic decay times and the greatest RMSF values, to investigate whether there is a correlation between these RMSF and autocorrelation decay time values, we plot RMSF vs. decay times per residue in Figure S3. Our results show that autocorrelation decay times become longer as fluctuations increase. G12D . Negatively correlated residue pairs show causality. All driverresponder residue pairs we identified in K-Ras G12D correspond to regions that are negatively correlated during MD simulations.

Relationships Between Negative Correlations, Characteristic Decay Times and Causality in K-RAS
Driver residues exhibit slow autocorrelation decays. A slow autocorrelation decay indicates that the residue has a strong memory of its past fluctuations. To compare the memory lengths of the residues in K-Ras G12D driver-responder pairs, we plot their autocorrelation decay curves (Figure 8). These plots show slower decay in drivers than responders, suggesting that they have longer memory lengths.

DISCUSSION
K-Ras is an important GTPase in cellular signaling, and is only active in its GTP-bound state 8,9 . Structurally, in active K-Ras, the P-loop, SI and SII regions are bound to the phosphate groups of GTP and are responsible for its GTPase function. However, when there is a G12D mutation in the P-loop, GTP hydrolysis is impaired and K-Ras freezes in its active state 9 , causing uncontrollable cellular growth and evasion of apoptotic signals 24,57,58 . In order to quantify the conformational and dynamic changes caused by the G12D mutation, we followed a new integrated MD data analysis approach. Our method enabled us to relate the structural and conformational changes in the residues of a protein to changes in the regulation of the protein's dynamics. Using our method, we first discovered that the G12D mutation causes the formation of new salt bridges between the residues of K-Ras, which alter the conformations of those at the active site and α3, potentially affecting GTP hydrolysis, and freezing K-Ras in its active state. We then correlated the conformational changes in the residues at the active site and α3 to the dynamic changes in these regions. Our results show that the motions of SII become coupled to the motions of other residues at the active site and α3.
Next, to understand the time-dependent regulation of K-Ras G12D motions, we measured the 'memory length' of each residue, and discovered that residues whose motions show strong coupling to the rest of the protein have a longer memory of their past. Specifically, the fluctuations of SII residues show the highest negative correlations with other parts of the protein and have the strongest memory. Most importantly, our results reveal the regulatory mechanism of K-Ras G12D dynamics, where not only the residues in SII exhibit correlated fluctuations with other residues, they also drive them. These memory lengths of the residues not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted August 19, 2017. ; https://doi.org/10.1101/178483 doi: bioRxiv preprint determine their driver-responder roles, and our results show that among all residues, those on SII have the strongest memory and are drivers in driver-responder pairs. In this study, our first important finding is that the broad distribution of Q61-P-loop distances affects GTP hydrolysis function of K-Ras G12D . Because in the G12D mutation, glycine residue is substituted by aspartate, which has a bulkier side group, it causes a structural change in the P-loop which affects the conformations of the other active site regions SI and SII, as observed in previous studies [22][23][24]59 . Our distance calculations show that after G12D mutation, the P-loop residues 11, 12, 13 move away from SII residue Q61. Moreover, the distances between Q61-A11, Q61-G12D and Q61-G13 display broad distributions (Figure 2). Since Q61 is a known catalytic residue that plays a critical role in both intrinsic and GAP-mediated GTP hydrolysis 10 , the highly variable nature of the P-loop-Q61 distance may affect the GTP hydrolysis in K-Ras G12D . Furthermore, our distance calculations also show that the residues in SII move away from some of the residues in SI, and α3, and distances between those residue pairs also display broad distributions in K-Ras G12D (Figure 2). Considering the increased fluctuations of SII, these broad distance distributions with larger peak values between SII and the other parts of the protein may be arising from the increased flexibility of SII in K-Ras G12D .
Our second important finding is that the residue E37 moves away from residue D57 in K-Ras G12D and alters Raf effector protein binding. The residues E37 and D57 are known hotspots for the interaction of K-Ras with Raf 60 . This interaction has been previously shown to be impaired upon G12D mutation 61,62 . Our distance analyses describe the possible reason for this impartment is that these two hot-spot residues move away from each other and their neighbors in K-Ras G12D .
Third, residues R73 and K104 move closer upon mutation and alter PI3Kγ effector binding. Pairwise distance calculations indicate that the distances of residue pairs 63-68, 72-75 and 73-104 significantly decrease upon G12D mutation. The 73-104 pair is particularly important since R73 is critical for interaction with the effector protein PI3Kγ 63 . Experiments have shown that PI3Kγ is preferably activated by K-Ras G12D with higher binding affinity 61,62 . Additionally, the backbone carbonyls of R73 interact with the amino sidechain of K104 64 . Therefore, the decreased distance between R73 and K104 in K-Ras G12D ( Figure 2D) may affect the interaction of these residues that allow PI3Kγ binding with higher affinity.
Fourth, active site residues that move apart in K-Ras G12D may alter its GTPase function. K-Ras undergoes conformational changes when it binds to GTP. The P-loop, SI and SII constitute the active site of the protein that binds to phosphate groups of GTP and participates in GTP hydrolysis 10 . SI and SII are also responsible for controlling the binding to effector molecules. Conformational changes in the active site affect K-Ras interactions with the GAPs which amplify the GTPase activity of K-Ras 65 . Our distance calculations showed that the active site residues of K-Ras G12D move away from their neighbors ( Figure 1C) due to the structural change in the P-loop. SII fluctuations also increase ( Figure 3A). These observations are consistent with previous studies that showed that the larger distances between active site residues and increased SII fluctuations in K-Ras G12D24 . Since conformations of the active site play important roles in GTPase activity of K-Ras and its binding ability to GAPs, we postulate that the deviation of active site residues may impair the GTP hydrolysis and also GAP binding ability which leads to the constitutively active K-Ras G12D .
And fifth, G12D mutation augments the correlations between SII and the regions by increasing the flexibility of SII. As seen in the pairwise correlation maps, K-Ras G12D differs markedly from K-Ras WT and displays increased correlations between the fluctuations of SII residues and other parts of the protein. This is consistent the increased amplitude of the SII fluctuations in K-Ras G12D , as shown in Figure 3. These results are consistent with a previous study that showed that SII fluctuations display increased level of fluctuations and negative correlations 24 .
Overall, our results provide a mechanistic understanding of the dynamic regulatory mechanisms of K-Ras G12D as well as predictions of sites that regulate its motions (such as SII). Targeting these sites with small molecules can be an effective strategy for the allosteric inhibition of oncogenic K-Ras G12D for the treatment of several cancers.