Ras Conformational Switching: Simulating Nucleotide-Dependent Conformational Transitions with Accelerated Molecular Dynamics

1Ras proteins are guanine nucleotide-dependent conformational switches that couple cell-surface receptors to signaling pathways that mediate cell proliferation, growth and development. Signal propagation through Ras is mediated by a regulated GTPase cycle that induces distinct conformations with different affinities for downstream effectors. The binding of GTP switches Ras to an “active” effector interacting form. Subsequent GTP hydrolysis returns Ras to the “inactive” GDP-bound form. Two types of regulatory proteins enhance the intrinsically low rates of these processes. GTPase activating proteins (GAPs) promote GTP hydrolysis, whilst guanine nucleotide exchange factors (GEFs) promote GDP release and regeneration of the active GTP-bound state. Mutations that lead to deregulated Ras activity are found in over 25% of human tumors. These autonomously active variants are insensitive to the action of GAPs resulting in uncontrolled cell growth. The current work aims to better understand the basis of these oncogenic transformations by deciphering how Ras changes its structure as it executes its enzymatic cycle and how the fidelity of this process is affected by mutations.

2Conformational changes and oncogenic mutations are largely concentrated in the vicinity of the nucleotide binding site, including the so-called switch regions SI (residues 25–40) and SII (residues 57–75). Of particular note are the conserved SI threonine (residue 35) and SII glycine (residue 60) which converge to form hydrogen bonds with the γ-phosphate of GTP (effectively ‘closing’ the nucleotide binding pocket). In the absence of the γ-phosphate (or a suitable analogue such as aluminium fluoride (AlF3)) the switch regions display fewer structural contacts to the nucleotide and reside in a more ‘open’ conformation. This observation has been likened to a loaded spring, where release of the γ-phosphate after GTP hydrolysis allows the switch regions to relax into their ‘open’ GDP-bound conformations. Crystallographic and spectroscopic studies have indicated that the switch regions exhibit substantial mobility both within and between different nucleotide states –. However, a detailed sequence of events and the mechanism by which certain mutations affect key dynamic rearrangements remains elusive.

3In the present study we employ simulation approaches to perform a detailed characterization of the dynamics of nucleotide-dependent conformational transitions. Previous unbiased molecular dynamics (MD) simulations were restricted to characterizing fluctuations within individual nucleotide states,. As a result, external biasing forces (e.g., via targeted molecular dynamics), were required to characterize the conformational exchange between nucleotide states. However, sampling by these methods is biased along narrow channels whose transition states may be of unrealistically high energy (100 and 120 kcal/mol compared to estimates of ∼22 kcal/mol ). These findings highlight the need for new simulation approaches to probe bias-free transitions. Recently we reported the observation of spontaneous nucleotide-dependent transition during unbiased MD simulation of the oncogenically active G12V variant. This study suggested the existence of a lower thermally accessible energetic barrier between inactive and active states of this variant that renders it prone to adopt an active conformational state. Here we extend this work with a multi-scale simulation approach employing classical and accelerated molecular dynamics (cMD and aMD), along with normal mode analysis (NMA) to study both wild-type and mutant transitions. Multiple-barrier-crossing aMD trajectories are used to characterize the reaction paths of the transitions. Simulated conformations are evaluated by comparison to the distribution of Ras's crystallographic conformers. The application of aMD allows the observation of nucleotide-dependent conformational transitions for wild-type Ras that are practically impossible to see with cMD. Furthermore, NMA indicates that the dominant collective motion that occurs during nucleotide-dependent conformational exchange, and captured in aMD simulations, is a low-frequency motion intrinsic to the structure. A number of highly populated intermediate conformations are characterized and their relation to available experimental structures discussed. Finally, the pattern of correlated motions in the current simulations reveals nucleotide-dependent differences of possible functional significance for membrane association.

4Simulations were conducted with starting structures corresponding to wild-type GDP, wild-type GTP and mutant GDP states of Ras. Each of these systems was simulated with a bound GDP and GTP. Both classical and accelerated MD (cMD and aMD) simulations were performed with explicit solvent for 60 nanoseconds. In addition to conventional structural analysis, which assessed the stability of Ras during the various simulations (see Text S1, Figure S1, and Table S1), principal component analysis (PCA) was used to relate the conformational sampling in the various simulations to available crystallographic structures. To this end we employ our previously reported PCA based mapping of available crystallographic structures. This analysis clearly distinguished distinct conformational clusters representing wild-type GDP, wild-type GTP and a number of minor clusters populated by switch and P-loop mutants (Figure S2). Together, these crystallographic structures provide invaluable landmarks against which simulation results may be compared.

5Figure 1 shows the conformational space sampled by simulations projected onto the lowest PCs determined from the distribution of X-ray structures. These projections illustrate the correspondence of the crystal structure distribution and those adopted during simulations. Also shown, in terms of collective coordinates, is the time evolution of the distances between GDP and GTP conformation cluster centroids for each trajectory. All aMD simulations were found to have significantly larger fluctuations in the plane defined by the lowest PCs, indicating that the main collective displacements in the crystal structure cluster are accessible during aMD simulations (Figure 1C, 1D, 1G, 1H, 1K, and 1L and Table S2). In contrast, cMD simulations display a more limited sampling (Figure 1A, 1B, 1E, 1F, and 1J). A notable exception is the cMD simulation of the oncogenic G12V variant (Figure 1I), where the introduction of GTP into an initially GDP-bound structure induced evolution towards the cluster populated by the experimental GTP-like conformers. Additional cMD simulations of this varient with different starting conditions consistently sampled two distinct regions: one near the cluster containing the inactive form of the G12V variant, and another close to the cluster of the active structures, suggesting a spontaneous inactive to active transition.

6Interestingly, aMD simulations carried out to further probe the apparent low activation barrier of the G12V variant (Figure 1K and 1L) display a similar range of sampling for both GDP and GTP-bound systems. Note that cMD simulations of this system with a bound GDP (Figure 1J) display a similar pattern of sampling to all other cMD simulations (with simulated conformers being largely restricted to the vicinity of the starting structure (Figure 1A, 1B, 1E, and 1F)). Implementing a range of boost-valued aMD simulations (see Figure S3) indicated that employing half the boost required for transition of wild-type systems (discussed bellow) was sufficient to achieve similar G12V transitions. This is consistent with the existence of a comparatively low activation barrier, even in the presence of GDP, and again highlights the intrinsic susceptibility for activation of this oncogenic mutant. Furthermore, the overall agreement between the current and the previous cMD simulations, carried out with the AMBER and CHARMM force fields respectively, demonstrates that these results are robust and independent of simulation details.

7In contrast to the absence of evident transitions during cMD for wild-type Ras, notable transitions are sampled during aMD for systems with a swapped nucleotide (i.e. GTP inserted into a starting GDP structure and vice versa (Figure 1D, 1G, and 1K)). Furthermore, the densities of populations for these simulations indicate that several highly populated intermediate conformations were sampled. In summary, the current analysis establishes the relevance of aMD simulations for investigating nucleotide-dependent conformational transitions. Subsequent discussion will therefore concentrate on further analysis of these aMD simulations of wild-type Ras that sampled the GDP-to-GTP (forward) and GTP-to-GDP (backward) transitions.

8Clustering of trajectory conformers was used to visualize the dominant conformations sampled by each simulation (Figure 2 and Table 1). The most populated cluster in the wild-type GDP-to-GTP transition (black in Figure 2E–H and Figure S4) corresponds to the lower of the three basins in Figure 1G and overlaps with the dominant conformation sampled during cMD simulations of the same system (Figure 1E). Its overall structure is intermediate between GDP and GTP states (RMSD of 1.3 Å from both GDP and GTP representatives, PDB codes 4q21 and 1qra respectively). Members of this cluster have an intermediate α2 orientation and a closed active site loop2-SI and loop4-SII that more closely resembles the GTP configuration. The second most populated cluster (yellow in Figure 2E–H) has more distinctive GTP like characteristics (RMSD of 1 Å from GTP and 1.6 Å from GDP representatives) including a closed loop2-SI and loop4-SII active site and a reoriented GTP-like α2 helix. Interestingly, the PC projection (Figure 2G) and RMSD values (minimum value 0.4 Å) indicate that members of cluster 2 closely resemble the crystallographic GTP-bound A59G structure (PDB code 1lf0). This structure has been suggested previously to be an intermediate and is characterized by a GTP-like SII conformation and a SI conformation that has undergone partial transition with the side-chain of Y32 adopting an orientation that is intermediate between that in wild-type GDP and GTP crystallographic structures,, (see Figure S2). The third cluster (green in Figure 2E–H) is again equidistant from GDP and GTP states (RMSD of 1.3 to 1.4 Å) and, similar to cluster 1, is characterized by an intermediate α2 orientation. However, unlike cluster 1, loop2-SI and loop4-SII resemble the open GDP conformation. Clusters 4 and 5 are closer to the GDP conformation (with RMSDs of 1 and 0.3 Å) than to the GTP conformation (with RMSDs of 1.4 and 1.7 Å). However, cluster 5 has distinct conformations of loop3 and the C-terminal portion of α5.

9The GTP-to-GDP (backward) transition sampled a number of distinct conformational states (Figure 1D and Figure 2A–D). Clusters 5 and 2 have a similar loop2-SI conformation and resemble the starting GTP state (RMSD values of 0.9 and 0.8 Å, respectively). A major difference between clusters 2 and 5 is the orientation of loop3. Furthermore, conformations in cluster 2 project to a similar region of PC space to those in the second cluster of the forward GDP-to-GTP transition described above. The remaining clusters (black, green and pink in Figure 2A–D) have a strikingly open loop2-SI conformation, similar to that observed in crystal structures 1×cm, 1bkd, 1nvv, chain R of 1nvu, chain R of 1nvx, and chain B of 1×d2. Cluster 3 most closely resembles the GDP state both in terms of its PC projection (Figure 2C) and RMS distance (1.1 Å). Clusters 1 and 4 project close to the GTP-bound crystallographic structures in which Y32 has been replaced by a Cys-chromophore. However, RMSD measurements of over 1.4 Å indicate features distinct from Y32 mutants and wild-type GDP/GTP states. Taken together, these results suggest that in addition to several distinct conformations that were not observed previously, both the forward and backward transitions pass through a common intermediate similar to that captured in the A59G crystal structure.

10The temporal evolution of cluster membership in each trajectory (Figure 2D and 2H) indicates that the first 4 to 5 ns remain relatively close to the starting structure (with conformations classified as cluster 5 or 4). This is followed by multiple transitions between states with periods of sampling where one conformation predominates. For example, conformations from the 10 to 25 ns period of the GDP-to-GTP trajectory reside predominantly in the GTP-like cluster 2. Similarly, conformations from the 33 to 42 ns period of the GTP-to-GDP trajectory are classified as the more GDP-like cluster 3. These single cluster blocks are interspersed with periods of rapid interconversion between clusters, for example the 30 to 35 ns portion of the GDP-to-GTP trajectory (Figure 2H). Interestingly, the GDP-to-GTP trajectory evolves toward a GTP-like cluster 2 conformation and then transitions back to a intermediate cluster 3 conformation before returning back via clusters 4 and 5 (see Videos S1 and S2). This behavior is consistent with the presence of a higher barrier between cluster 2 and the main crystallographic GTP cluster, which we fail to cross, than exists between the intermediate conformations contained within our clusters. The high level of inteconversion between these clusters also suggests that they are energetically relatively close to one another.

11Analysis of the calculated structures indicates that certain side chain reorientations, diagnostic of GTP and GDP crystallographic states,,,, have only been partially realized in our simulations. We speculate that while the topology of the backbone provides the conformational blue print, specific side-chain interactions are required to stabilize the canonical GTP and GDP states. The apparent barrier between cluster 2 and the main crystallographic GTP cluster would therefore result from the energetic cost of reorganizing these side chains. A typical example is Y32, whose orientation in the simulated intermediates is neither fully solvent exposed nor hydrogen bonded with Y40 as in the main GTP and GDP crystallographic clusters respectively. Rather we note its similarity to that found in the crystal structures 1gnr and 1gnq.

12To examine whether the motions of one residue are related to the motions of another (distant) residue, the correlation of the displacements of all residue pairs were determined (Figure 3). As expected, the strongest positive correlations exist between covalently bonded residues and those residing within secondary structure elements (see Figures S5 and S6 for reference contact maps). Moving up the diagonal, the first area of notable correlations corresponds to helix α1 (residues 16 to 25) with loop2 (residues 26 to 40). The next area of significant correlation corresponds to the β2-loop3-β3 region (residues 38 to 57). The consistent appearance of correlated motions for these regions in each simulation and the cross-correlation (off-diagonal peak) with β1 (residues 2 to 10) highlights the subdomain-like structure and dynamics of these three N-terminal strands. The remaining strands, β4 to β6, display consistent positive cross-correlations with each other but not with strands β1 to β3. Note the off-diagonal peaks for β4 (residues 77 to 83) with β5 (residues 111 to 115) and β5 with β6 (residues 141 to 144). Moving further up the diagonal the next major correlations correspond to the SII region, encompassing loop4 and helix α2 (residues 58 to 74).

13Perhaps the most interesting feature of the plot is the pattern of correlation between α2 and α3-loop7 (residues 66 to 74 and 93 to 110). This feature is most evident in GTP-bound simulations and is largely absent in GDP-bound simulations. This pattern is particularly noteworthy as GDP-to-GTP aMD simulations exhibit these correlations only in portions of the trajectory that reside in a GTP like conformation (i.e. after the transition from GDP to GTP, see Figure 2G and 2H). It appears that the large rearrangement of helix α2 during the transition brings it into closer register with α3, hence facilitating the correlated motions of these regions. Furthermore, the correlations of these regions in the GTP-to-GDP aMD simulation reduce gradually as the conformation of α2 evolves toward a more GDP like state. These data thus show a novel GTP-dependent correlated motion in Ras that has functional implications (see below). Additional off-diagonal peaks include loop2 with loop10 (residues 26 to 37 and 145 to 151) and loop3 with α5 (residues 46 to 49 and 152 to 167). As discussed below loop3 and α5 are connected via several salt bridges whilst both loop2 and 10 directly interact with the bound nucleotide. These newly identified coupled motions suggest a dynamic linkage between the N-terminal nucleotide-binding subdomain and the C-terminal subdomain whose downstream residues are responsible for membrane binding.

14The dissection of the catalytic domain into two lobes or subdomains based on the correlated motions of the central β-strands is consistent with the localized nature of sequence variation between Ras isoforms. As previously noted, lobe 1 (residues 1–86) is strictly conserved in sequence and encompasses the P-loop and the switch regions; whilst lobe 2 (residues 87–171) contains amino acid variations that define functionally distinct H, N and K-ras isoforms. As isoform-specific properties include differences in nucleotide-state dependent membrane localization –, the segregation of both sequence variation and correlated motions implies that communication between lobes is likely to be of functional significance. The covalent connection between lobes is made by helix α2 of the SII region, which is the major dynamic element of the Ras structure. The current cross-correlation analysis indicates the existence of three additional non-covalent communication routes between lobes including loop3 to α5, loop2 to loop10 and α2 to α3-loop7. We speculate that residues at each of these sites may be important for nucleotide-dependent modulation of membrane attachment and lateral segregation by linking the switching apparatus to the membrane interaction apparatus. Indeed, alanine substitution of loop3 residues D47 and E49 produced a variant that is hyperactive in MAPK-signaling. These loop3 residues, together with their α5 salt bridge partners (R161/R164), have been shown to modulate the nucleotide-dependent membrane association of H-ras,. The other regions highlighted in the current study have thus far received little attention but likely warrant further investigation.

15In an effort to further understand the physical basis of the observed motions upon nucleotide exchange, we analyzed available structures with a simplified elastic-network normal mode method. The elastic network approach has the advantage that a single model, expressed in terms of Cα coordinates, leads to an objective expression of possible protein dynamics in terms of a superposition of collective normal mode coordinates. Consistent with previous studies on a range of systems, we note that only the “open” wild type GDP conformation yielded large overlap values between the crystal structure PCs and the low-frequency normal modes (overlap between NMA mode 1 and X-ray PC 1 of 0.57). Furthermore, the structural mobility predicted by NMA is very similar to that obtained from aMD simulations (overlap between GDP-to-GTP trajectory PC 1 and X-ray PC 1 of 0.69; Figure 4 and Table S2). This result implies that low-frequency global motions that are intrinsic to the open structure likely facilitate the observed conformational transitions to the closed GTP state. In other words, the fact that low frequency normal modes qualitatively capture the differences between available crystal structure conformations and have high overlap with the eigenvectors obtained from aMD simulations suggests that nucleotide-dependent dynamics is facilitated by the low frequency, global motions that are intrinsic to the structure. The aMD results discussed above indicate that the nature of the bound nucleotide attenuates these intrinsic motions.

16We have characterized the spontaneous transition between nucleotide-dependent conformational states of wild-type Ras with cMD, aMD and NMA. These functionally important transitions, achieved for the first time for wild-type Ras, are practically impossible to observe with cMD. Furthermore, NMA indicates that the dominant collective motion that occurs during these transitions is a low-frequency motion intrinsic to the structure.

17Mapping the reaction path sampled by aMD onto a PCA basis set derived from the distribution of crystallographic structures enabled identification of intermediate structures with unique switch orientations and/or distinct loop3, loop7 and α5 C-terminus conformations. Intriguingly, several of the highly populated intermediates have a close correspondence to known G59A and Y32C crystallographic conformers, both of which have been suggested to be intermediate structures,. The emergence of these conformations along with additional novel intermediates highlights the utility of aMD simulations to reliability sample conformational transitions. Furthermore, the current results imply that the G59A and Y32C variant conformations are accessible to wild-type Ras and that these mutations result in perturbations that localize the average structure at these intermediate positions. The functional relevance of these intermediates is reinforced by kinetic studies of G59A that indicated a reduced rate of nucleotide exchange; this has been linked to the need for tight nucleotide coordination during structural changes,. It will be interesting to see if our newly identified intermediates, some of which differ between forward and backward transitions, also have distinct kinetic behaviors related to nucleotide exchange and phosphate release.

18The pattern of correlated motions revealed by these simulations predicts novel nucleotide-dependent motions of potential significance in the signaling function of Ras. These include correlations of α2 with α3-loop7, loop2 with loop10 and loop3 with α5. Such dynamic linkages between the switching apparatus and the membrane interacting C-terminal region leads us to speculate that residues at each of these sites may be important for nucleotide-dependent modulation of membrane attachment. This is supported by recent experimental evidence for the role of loop3 residues D47 and E49 and α5 residues R161 and R164 in modulating the nucleotide-dependent membrane association of Ras,.

19Finally, low frequency normal modes qualitatively capture the differences between available crystal structure conformations and have high overlap with the eigenvectors obtained from aMD simulations. This result combined with aMD observations suggests that nucleotide-dependent dynamics is facilitated by low frequency, global motions that are intrinsic to the structure and that the nature of the bound nucleotide serves to attenuate these intrinsic low-frequency motions. Furthermore, the significant similarities of aMD, NMA and crystal structure PCA motions highlight the robustness of the observed motions.

20We believe that the current advanced simulation and analysis approach is equally applicable to a large number of structurally similar but functionally diverse P-loop NTPases such as kinesin and myosin. Such studies should uncover detailed dynamic behavior and help inform us about general principles and mechanisms underlying nucleotide-dependent conformational changes.

21All simulations were performed with the AMBER8 package and corresponding all-atom potential function ff99SB. Additional parameters for guanine nucleotides were taken from Meager at al.. All analysis was carried out using the Bio3D package,.

22Atomic models were prepared from three high-resolution crystal structures (PDB codes: 4Q21, 1QRA and 2Q21). These structures are representative of three distinct conformations highlighted by PCA; namely: wild-type GDP, wild-type GTP and mutant GDP. Each model was simulated with bound Mg2+·GDP and Mg2+·GTP. The latter was obtained for GDP structures by the addition and local optimization of the γ-phosphate onto the GDP of the corresponding crystal structure.

23All MD simulations were performed using periodic boundary conditions, TIP3P water and charge-neutralizing counter ions, with full particle-mesh Ewald electrostatics. Operational parameters included a 2fs time step and a 10Å cutoff for the truncation of VDW non-bonded interactions. Constant volume heating (to 300 K) was performed over 10ps, followed by constant temperature (300 K), constant pressure (1atm) equilibration for an additional 200ps. Finally, constant pressure constant temperature production dynamics was performed with both classical and accelerated MD implementations. The SHAKE algorithm was used to constrain all covalent bonds involving hydrogen atoms.

24Accelerated MD (aMD) extends the accessible time scale of conventional MD simulations by altering the underlying potential energy surface of the system under study. Acceleration stems from the addition of a non-negative boost potential that raises the energy within basins. Hence a trajectory propagated on this modified surface makes transitions from state to state with an accelerated rate. Furthermore, canonical ensemble averages of the system can be obtained by reweighing each point on the modified potential by the strength of the Boltzmann factor of the bias energy at that particular point. In the current study we apply the dual boost approach and corresponding potential developed previously,. Starting structures and standard operational parameters were identical to those used for cMD. The energy level, E, below which the boost is applied and tuning parameter, α, that modulates the depth and local roughness of basins in the modified potential were based on an earlier work.

25We employed the coarse-grained AD-ENM normal mode analysis approach developed by Zheng et al.. AD-ENM implements a single-parameter Hookean potential, which has previously been shown to yield low-frequency normal modes that are in good agreement with those obtained from more detailed, empirical, force fields. For further details see,

26PCA was employed to aid the interpretation of interconformer relationships. We utilized the previously reported PC basis set obtained from analysis of available Ras crystal structures. This basis set gives a clear separation of nucleotide-dependent conformational states. Projecting the Ras crystal structures and snapshots from MD trajectories into the sub-space defined by the largest PCs (along which the crystal structure variance is largest) results in a lower dimensional representation of the structural dataset (see Figure 2 for details). The resulting low-dimensional ‘conformer plots’, succinctly reveal the nature of conformational sampling during simulations. PCA was carried out on the individual trajectories using the same Cα atoms that were used in the analysis of the crystal structures. Conformer superposition was also based on the “core positions” obtained from crystal structure analysis.

27Distances between instantaneous trajectory conformations and the centroids of the main GTP and GDP clusters, reported as inserts in Figure 1, were calculated as the Euclidean distance between projected points in five dimensional PC space. Note that five PCs account for over 81% of the variance in the original distribution and produce a more succinct distance measure than the examination of average all-atom distances. This metric aids interpretation of an otherwise noisy signal as it is derived primarily from the concerted displacement of the switch regions comprising the secondary structure elements loop2 and loop4-α2 (residues 31 to 37 and 59 to 72).

28Structures from aMD simulations underwent average-linkage hierarchical clustering according to their pairwise RMSD distance matrix. Inspection of the resulting dendogram was used to partition structures into five dominant groups (ranked according to their populations). The closest structure to the average structure from each cluster, in terms of RMSD, was chosen as a representative for projection onto the PCA basis set described above.

29To identify protein segments with correlated atomic motions the cross-correlation coefficient, Cij, for the displacement of all Cα atom pairs, i and j, was calculatedwhere Δri is the displacement from the mean position of the i th atom determined from all configurations in the trajectory segment being analyzed (see and for further details).