Machine learning-driven multiscale modeling reveals lipid-dependent dynamics of RAS signaling proteins.

RAS is a signaling protein associated with the cell membrane that is mutated in up to 30% of human cancers. RAS signaling has been proposed to be regulated by dynamic heterogeneity of the cell membrane. Investigating such a mechanism requires near-atomistic detail at macroscopic temporal and spatial scales, which is not possible with conventional computational or experimental techniques. We demonstrate here a multiscale simulation infrastructure that uses machine learning to create a scale-bridging ensemble of over 100,000 simulations of active wild-type KRAS on a complex, asymmetric membrane. Initialized and validated with experimental data (including a new structure of active wild-type KRAS), these simulations represent a substantial advance in the ability to characterize RAS-membrane biology. We report distinctive patterns of local lipid composition that correlate with interfacially promiscuous RAS multimerization. These lipid fingerprints are coupled to RAS dynamics, predicted to influence effector binding, and therefore may be a mechanism for regulating cell signaling cascades.

RAS driven cancers are common (1), difficult to treat (2), and a major cause of death worldwide (3). KRAS, the isoform most frequently associated with disease, is mutated in nearly all pancreatic cancers and often in lung and colorectal cancers (4, 5). Only recently, with the development of covalent inhibitors of the G12C mutant (6), has direct targeting of RAS been successful, and more broadly applicable inhibitors are needed.

The RAS subfamily comprises peripheral membrane proteins with a conserved globular GTP-binding domain (G-domain) (7) that is tethered to the cell membrane by a prenylated ∼20-residue C-terminal domain called the hypervariable region (HVR) (8, 9). RAS proteins function as molecular switches whose active conformations, stabilized by GTP binding, interact with several protein effectors to control cell growth, proliferation, differentiation, and migration (10). Constitutive activation of oncogenic RAS perturbs several cellular signaling cascades, including the MAPK pathway, which RAS accesses via activation of RAF kinase at the plasma membrane (PM).

There is substantial interest in assessing the ability of RAS molecules to dimerize (11–13) or colocalize (14–17) at the membrane, because RAS-dependent RAF activation requires dimerization of RAF (18–20). Although wild-type KRAS4b, a common splice variant of KRAS (hereafter referred to as RAS), does not dimerize on two-component supported lipid bilayers (21), it preferentially colocalizes with anionic lipids in the liquid-disordered domains of giant unilamellar vesicles (22) and clusters on the scale of tens of nanometers in extracted PM sheets (14, 23, 24).

Preferential interaction of RAS with anionic lipids is mediated by 11 positively charged lysines in its HVR (25, 26). However, charge complementarity is insufficient to fully describe RAS nanoclustering, which is exquisitely sensitive to lipid composition (27–29). Even less is known about the influence of RAS–lipid coupling on RAS self-assembly and effector activation. While several feasible dimer interfaces have been reported (11–13, 15, 30–37), how RAS forms dimers, if at all, remains a major area of interest.

The fundamental challenge of investigating RAS activation events is that many of the proposed mechanisms involve time and lengths scales currently not accessible. For example, functional events in RAS dynamics that may preferentially depend on local depletion or enrichment of specific lipids are extremely difficult to observe directly, either in computational or biological experiments. Experimentally, we use single-particle tracking (SPT) to follow HaloTag-conjugated RAS via total internal reflection fluorescence microscopy in live HeLa cells. While the broad heterogeneity in lateral diffusion observed by SPT (Fig. 1) is indicative of multiple RAS subpopulations that have distinct patterns of interaction with lipids and other cellular components, it provides little insight on local lipid environments and their effect on RAS behavior. Similarly, we can use detailed molecular dynamics (MD) simulations to probe specific lipid environments, but systems large enough to support substantial fluctuations in lipid and protein concentrations cannot be practically simulated long enough to observe relevant fluctuations. Instead, we present a multiscale infrastructure that directly couples a macroscale continuum simulation capable of observing RAS clustering and long-range lipid rearrangement, with a massive ensemble of microscale MD simulations that provides detailed insights into local dynamics. Both scales are connected through a machine learning (ML) informed dynamic sampling process (38), which enables mapping of findings from the MD simulations onto the continuum simulation, resulting in microscale insights that are observable over macrospatial and temporal scales.

Fig. 1.

Experimental input and computational approach. (A) Diffusion mapping of single molecules of KRAS4b tethered to or within 100 nm of the PM in a 16 × 16-μm2 region of a live HeLa cell accumulated over 10 s. (B) Crystal structures of wild-type KRAS in active (green and blue; GppNHp-bound) and inactive (gray; GppCH2p-bound configurations). (C) Average macro model lipid composition. (D–F) The Multiscale Machine-learned Modeling Infrastructure (MuMMI). (D) Representative snapshots of each of the different lipid distributions in the inner leaflet of a 0.3 × 0.3-μm2 region of the full 1 × 1-μm2 macro simulation; color saturation indicates local lipid density. (E) Schematic illustrating latent space encoding of lipid composition in 30 × 30-nm2 membrane patches. From the candidate patches (blue and green), those that are most dissimilar (green) to existing (white) CG simulations are selected and used to spawn new CG simulations. (F) Representative CG simulation systems (water not shown). (G) Improvement of macro model parameter inputs from feedback iteration. (H) Distribution of CG simulation duration and (Inset) number of RAS per patch.

More specifically, to enable simulations of active wild-type RAS, we solve its crystal structure bound to the GTP analog GppNHp at a resolution of 2.5 Å (Fig. 1 and ). We choose to focus the effort on wild-type GTP-loaded KRAS4b because GTP hydrolysis will not occur on the time scale of the simulations, and because available structures of wild-type or oncogenic KRAS4b bound to the RAS binding domain (RBD) of RAF1 are similar. This present work will serve as the foundation for understanding activation of RAF by GTP-loaded RAS. We mimic the composition of a biologically relevant PM (39) by employing an asymmetric eight-lipid mixture (40), here called the average RAS lipid composition (ARC) and is composed of cholesterol, phosphatidylcholines (PC), phosphatidylethanolamines (PE), phosphatidylserine (PS), phosphatidylinositol bisphosphate (PIP2), and sphingomyelin (SM), with varying acyl chain length and saturation (Fig. 1 and ).

Multiscale Simulation Framework

As discussed above, the behavior of RAS is mediated by the local lipid environment, which is not readily accessible experimentally. In principle, computational models can provide the necessary details for small neighborhoods around RAS, but they require a priori knowledge of which protein confirmations and local lipid environments are of interest. Traditional computational studies typically choose an average lipid composition and focus on providing many replicates and long-running trajectories. However, in vivo, proteins experience a wide range of different lipid configurations and, especially considering effects such as charged lipids, a protein is in fact unlikely to experience an average lipid environment. Instead, we expect the local lipid composition to be heavily influenced by the proteins and the codependence of the lipid and proteins is one of the key properties of interest. Manually selecting protein conformations and sets of lipid compositions can provide valuable insights, but requires a prior knowledge of what amplitudes and along what dimensions one should adjust the simulations. Unfortunately, running an MD simulation large enough and long enough to naturally express a much wider range of lipid environments and protein conformations is practically not feasible. Here, we introduce MuMMI, the Multiscale Machine-learned Modeling Infrastructure (Fig. 1 and ), which combines the relatively low computational cost of MD simulations in local neighborhoods around proteins with a diverse sampling of environments equivalent to running an otherwise infeasible MD simulation. As discussed in more details below, the resulting system efficiently samples local environments from a macroscale model, such that the resulting ensemble of MD simulations can be used to study the entire macroscale model at what amounts to MD resolution. For a full description of MuMMI, see , and for the description of the model generation and parameterization, see .

MuMMI employs a macro model (Fig. 1) to explore interactions between lipids and proteins over large length and time scales. The macro model represents a lipid bilayer as a set of density fields describing the density of each lipid type as a function of space and evolves them using dynamic density functional theory. Proteins, RAS in our case, are represented as single, rotationally averaged beads, which interact with the lipids and other proteins through defined interaction potentials evolved using MD equations of motion. In this work, the macro model consists of a 1 × 1-µm2 PM, eight-lipid type (ARC) bilayer that is resolved at roughly nanometer resolution and with 300 RAS molecules on the inner leaflet. The RAS proteins exist in one of two distinct states representing different orientations with respect to the membrane, which they can transition between. During the MuMMI multiscale simulation, the macro model ran for 150 µs. The theory behind the macro model is described in Stanton et al. (41) and summarized in . The macro model parameterization is described in , while the MD simulations used to derive the initial parameters and the initial RAS state definition are described in , respectively.

While not providing sufficient resolution to probe molecular processes, the macro model faithfully models the overall lipid dynamics and naturally creates local fluctuations of the average membrane composition in which certain lipid types are enriched or depleted, driven by membrane mechanics and protein–membrane interactions. In particular, the macro model explores a wide range of lipid environments beneath the RAS proteins and protein stoichiometries mimicking the diverse configurations encountered in a cell PM (). More specifically, MuMMI focuses on local, here 30 × 30 nm2, subregions (hereafter, patches) and their corresponding proteins at both the macro- and at near-atomistic scales. One potential solution to explore the phase space of patches would be to randomly sample from the macro model, creating the corresponding coarse-grained (CG) Martini (42) MD simulations (Fig. 1 ) and ultimately assembling a rough approximation of the expected behavior at the macroscale. However, the patch distribution is highly nonuniform with some configurations occurring orders-of-magnitude more frequently than others. Nevertheless, compared to the relevant biological time scales even rare events, occurring at a frequency of 150 μs in the macro model, should not be considered outliers. In fact, rare lipid compositions might trigger signaling events that drastically amplify their overall impact and it is these cases we are most interested in studying. Therefore, MuMMI implements a ML-driven dynamic-importance (DynIm) sampling (38) to select a maximally diverse subset of patches to evaluate with CG simulations.

MuMMI uses a ML-driven selection process, which first encodes local lipid concentrations (patches) into a 15-dimensional latent space using a variational autoencoder that learns to identify similarity among macro configurations and then uses a novelty-sampling scheme to select the patches explored by the macro model that are most dissimilar to the ones already simulated at the CG scale (Fig. 1 and ). The ML model was trained prior to running the multiscale simulation using 302,000 patches generated from a separate run of the macro model; it was used for inference to dynamically select novel patches generated during the MuMMI multiscale simulation. (Computational limitations prohibited retraining the autoencoder on-the-fly; therefore, the pretrained ML model was used throughout.) Selected patches are used to spawn new CG simulations creating a maximally diverse set of CG simulations expected to explore a much wider range of RAS behavior. Simultaneously, MuMMI maintains sufficient information about the nonselected patches to accurately reconstruct the true, unbiased distribution of patches in latent space (). Note that the reconstructed distribution has a significantly wider support than the one that could have been recovered from a traditional random sampling () of comparable size. In other words, by deliberately including rare patches in a correctly weighted manner, our distribution is equivalent to a substantially longer macro simulation and thus is able to explore more lipid environments and protein configurations than otherwise possible. Ultimately, the weighted distribution allows computation of accurate, population-wide averages of any property of interest computed from the CG simulations. Furthermore, as more computing resources become available, the sampling of patches with corresponding CG data becomes denser, ultimately converging to a state at which, for every patch created by the macro model, MuMMI has either the corresponding CG simulation or one for a patch so similar as to be considered equivalent. At this point, the ensemble of CG simulations, together with the macro model data, jointly represent a simulation at macrolevel time and length scales but effectively at MD precision.

At the microscale, CG MD simulations are used to capture interactions of proteins and lipids. Martini CG simulations provide over two orders-of-magnitude speedup compared to all-atom simulations and are well suited to capture overall lipid–lipid, protein–lipid, and protein–protein interactions, especially at longer length and time scales (42–45). The reduced degrees-of-freedom come with a price, and model limitations, such as lack of directional hydrogen bonds, need to be considered when interpreting the results, especially for protein–protein interactions (46, 47) (). In MuMMI, each selected macro model patch is converted into a CG simulation containing ∼140K particles. Membranes are built with subgrid resolution and statespecific placement of one or more RAS molecules sampled from a library (Fig. 1 and ), resulting in a large sampling of lipid compositions and RAS configurations (). These CG simulations are performed on GPUs using ddcMD (48) () and analyzed in situ (). The on-the-fly analysis of each CG simulation is weighted based on its representative patch in macro model example () and used to continuously refine the macro model protein–lipid parameters for each protein state (Fig. 1 and ), incrementally improving the accuracy of the macro model with improved sampling.

MuMMI represents an emerging new capability that delivers a macroscale model at experimentally relevant time and length scales with effective microscale resolution. The simulations presented in this study represent a massive ensemble comprising 120K, 1- to 4-μs-long, unique but correlated CG simulations, simulated over several weeks using one of the world’s largest supercomputers, Sierra, and consuming ∼5.6 million GPU hours (). These CG simulations capture and quantify RAS–lipid dynamics across local compositions observed in the much larger macro simulation. A new massively parallel workflow (49) was developed to orchestrate MuMMI () and MuMMI is capable of utilizing some of the largest supercomputers in the world; its functionality is also equally applicable to modest and more-commonly available compute resources.

RAS–Lipid Colocalization

Lipid fingerprints are spatially resolved lipid concentration signatures around proteins (50) and can be clearly observed for RAS in the PM both at the macro and CG resolution (see below). The macro simulation reveals lipid redistribution around RAS (Fig. 2), indicative of discrete lipid fingerprints, and a propensity for RAS to self-associate (Fig. 2). The latter is consistent with experimental RAS clustering as seen in atomic force microscopy experiments with simple lipid mixtures (22, 51).

Fig. 2.

Lipid-dependence of RAS colocalization. (A) Representative macro model inner leaflet lipid densities around RAS, shown separately (small boxes) and together (large central box). Color saturation indicates local lipid density. (B) Population ratio of RAS multimer sizes observed in the macro simulation vs. a random uniform distribution. (C) Average number of RAS in macro model regions (radius 5 nm) along the primary embedding dimension (PED) from function preserving projection analysis (). Vertical lines denote thresholds used to define high-RAS (HRC), initial average (ARC), and low-RAS (LRC) lipid compositions. (D) Distributions of inner leaflet lipid concentrations for all patches with RAS (black), HRC (red), and LRC (blue). The ARC is represented by a dashed vertical line. (E) Surface plasmon resonance sensorgrams of RAS adhesion to liposomes with the LRC, ARC, and HRC lipid compositions. Each subplot contains multiple traces representing distinct RAS concentrations (twofold dilutions, 60 to 0.1 μM).

To explore the connection between the local lipid environment and RAS clustering we use function preserving projections (52), which uses ML to identify two principal directions that are linear combinations of local lipid concentrations and maximally correlate these concentrations with the number of RAS in a patch (Fig. 2 and ). A subselection within the high/low RAS regions of the embedding is used to find unimodal lipid compositions that are either RAS depleted (low RAS colocalization, LRC) or enriched (high RAS colocalization, HRC) as classified in Fig. 2 and . In particular, the HRC is enriched in PIP2 and depleted in cholesterol and polyunsaturated lipids, whereas the LRC exhibits opposite trends. These findings on our in silico membrane corroborate RAS’s preferential association with PIP2 on experimental artificial membranes (53, 54). However, a discrepancy still exists between reconstituted systems and the observed enrichment of PS and phosphatidic acids in RAS nanoclusters on intact-cell PM sheets (14, 24). While some lipid concentrations, such as 15.5% PIP2 in the HRC, are extreme for the global PM, our macroscale model indicates their presence in spatially localized regions due to concentration fluctuations and protein-mediated enrichment. Surface plasmon resonance confirms that, compared to the initial average lipid composition (ARC), RAS binding is weaker to the low LRC and stronger to the high HRC compositions (Fig. 2), consistent with specific lipid concentrations influencing RAS recruitment and colocalization.

RAS Diffusion and Preferential Interaction with PIP2

Peripheral membrane proteins display dynamic diffusional properties as they interact with the hierarchical architecture of the PM (55). SPT of RAS proteins identified three discrete populations with diffusion coefficients represented here by values acquired from HeLa cells: one dominant fast mobile component (0.96 ± 0.03 μm2/s), an intermediate component (0.26 ± 0.03 μm2/s), and a slow component (0.05 ± 0.01 μm2/s) (56). Defining the disposition of RAS in these discrete diffusional populations is challenging, and subdiffraction microscopic techniques do not provide sufficient resolution to accurately measure the size of these RAS complexes. Coupling conformational fluctuations of protein structure with dynamic diffusional properties of RAS on the membrane can only be bridged using simulations.

It is well documented that experimentally measured RAS has a preference for liquid-disordered membrane patches enriched in anionic lipids. A variety of different experimental techniques, using different lipid mimetics (supported lipid bilayers, bilayers, nanodiscs, or liposomes) composed of simple two-lipid mixtures (21, 25) to those composed of five or more lipids (36, 57), all indicate that KRAS has a preference for PS (22, 57) and phosphatidylinositol. Lysine residues within the HVR of RAS (27) are known to interact with PS in mammalian PMs and, as a consequence, impact the propensity of RAS to form nanoclusters. Previous MD simulations have identified that RAS can exist in multiple membrane-bound orientational states (58–61) while bound to two component lipid mixtures. These simulations are limited in their complexity of the lipid mixtures represented, even though the complexity of the lipid composition in the PM is well established (62). In addition, the length of the simulations is not long enough to explore even the shortest lifetime of RAS molecules identified in mammalian cells (56). Therefore, to represent more biologically relevant models, extended time scales and increased complexity of the lipid bilayer are needed in such simulations.

In this study, the diffusion properties of the lipids and RAS on ARC were evaluated in simulation as well as in experiments using fluorescence lifetime correlation spectroscopy (FLCS). All lipid types show a broad distribution of diffusion coefficients among the 88,392 unique CG simulations with one RAS molecule (Fig. 3 and ). This variation reflects the dependence of lipid diffusion on membrane composition. For example, a 10% change in cholesterol content slows the diffusion of all lipids except PIP2 by a factor of 2 (). For a more complete discussion of the impacts of membrane composition, see . The ARC diffusion rates calculated from FLCS experiments are shown in Fig. 3. All lipids except PIP2 diffuse at similar average rates in both FLCS and in the ensemble of simulations, whereas RAS and PIP2 diffuse two to three times slower (Fig. 3 ). Note, absolute diffusion numbers can vary significantly depending on experimental method used (63), as well as simulation force field, size, and thermostat (64, 65) () and should only be considered as semiqualitative. Here, all simulation diffusion numbers are scaled by the Martini canonical factor of 4 (42). Interestingly, in the simulation a direct interaction with RAS slows the diffusion of PIP2 (Fig. 3). In CG simulations with multiple RAS molecules, increased aggregation of RAS is associated with local enrichment of PIP2 and depletion of PAPS (1-palmitoyl-2-arachidonoyl-sn-glycero-3-phosphatidylserine) (Fig. 3), emphasizing the potential importance of charge density. RAS aggregation also slows the diffusion of both RAS and PIP2 (Fig. 3), while the diffusion of other lipids remain unaltered.

Fig. 3.

Diffusion of RAS and lipids. (A and B) Distributions of lateral diffusion rates, D, for lipids and RAS in (A) CG simulations and (B) FLCS on the ARC. (C) Values of D for PIP2 conditioned on interaction with monomeric RAS in the CG simulations. (D) Proportion of each lipid type in RAS’s first solvation shell (within 0.55 nm), normalized by its molar fraction in the initial average lipid mixture (ARC), as a function of RAS aggregate size. (E) Values of D as a function of RAS aggregate size. (F) Mean square displacement (MSD) curves of JF646-labeled RAS from SPT on supported lipid bilayers with the low (LRC), average (ARC), and high (HRC) RAS lipid compositions.

To further evaluate how the lipid composition affects RAS diffusion, we use SPT to study RAS on supported lipid bilayers with low (LRC), average (ARC) and high (HRC) RAS colocalization compositions. Mean square displacement plots from 50,000 to 100,000 tracks obtained on each bilayer show that RAS undergoes anomalous diffusion that is slower and more confined on HRC compared to LRC and ARC (Fig. 3). A more rigorous hidden Markov modeling (HMM) analysis of the SPT data on HRC and LRC identifies three interconverting states in which RAS diffusion is slow (∼0.1 μm2/s), medium (∼0.5 μm2/s), and fast (∼4 μm2/s) (). It also reveals that the lipid-dependence of RAS diffusion (Fig. 3) results from a decrease in slow and medium diffusion rates and an increased population of slowly diffusing RAS in the HRC (). One interpretation of these data is that the high PIP2 concentration in HRC directly correlates to large RAS aggregate that diffuses much slower compared to less PIP2 and a smaller RAS cluster on LRC, although other factors can also slow diffusion.

Taken together, simulation and experiments identify a lipid composition (HRC) that recruits more RAS (Fig. 2 C–) and favors RAS multimerization or confinement (Fig. 3), slowing its diffusion (Fig. 3). Lipid-driven RAS clustering provides an attractive hypothesis as an initial organizational step in nucleating a molecular signaling platform for efficient recruitment of RAS effectors.

RAS–RAS Association

RAS multimers are hypothesized to function as signaling hotspots, and defining RAS–RAS interaction is critical to understanding the role of colocalization in signaling competence (14). Indeed, specific RAS–RAS interfaces have been identified in previously published works (11, 13, 15, 30, 31). However, analysis of 10,939 MuMMI simulations with RAS–RAS contacts reveals a broad distribution of interfacial arrangements (Fig. 4), including published interfaces (13, 31, 32, 66) (Fig. 4). The resulting modest (approximately twofold) interfacial preferences are inconsistent with the existence of a distinct dimer interface (Fig. 4). As such, distinct lipid fingerprints may facilitate initial RAS association through preorientation of the protein. Due to these multiple potential interfaces and heterogeneity seen in RAS–RAS association (Fig. 4), it is problematic to deconvolute a nonrotationally averaged lipid fingerprint for RAS multimers. Nevertheless, the local lipid environments of RAS aggregates are enriched/depleted in different lipid species, including enriched in PIP2 and depleted of ordered lipids in comparison to RAS monomers (Figs. 3 and 4 , and ). This lipid–RAS codependence is consistent with a signaling mechanism in which lipid–RAS interactions concentrate RAS prior to RAF binding, and in which subsequent RAF dimerization is facilitated by this RAF-independent RAS colocalization.

Fig. 4.

RAS–RAS interactions. (A) Distributions of radii of gyration based on RAS aggregation number. (B) Sampling of the protein–protein interfaces projected as a density map of two angles. The first angle (1stRASΘ, in red) is between the vectors center of mass (COM) of the first RAS G-domain to COM of the second RAS G-domain (1stRAS GCOM→2ndRAS GCOM) and 1stRAS GCOM→1stRAS T35. The second angle (2ndRASΘ, in blue) is between the vectors 1stRAS GCOM→2ndRAS GCOM and 2ndRAS GCOM→2ndRAS T35. Illustrative examples of the interfaces defined by the angles are shown. White dashed regions indicate various interfaces presented in the literature (13, 31, 32, 66). Data are only plotted for 1stRAS GCOM→2ndRAS GCOM < 6.0 nm from simulations where the RAS proteins start apart (). (C) PIP2 remodeling based on RAS–RAS G-domain separation, dG-G, shown with G-domain centers of mass on the x axis. Data show PIP2 density (red heatmaps), PIP2 density integrated over −4 nm < y < 4 nm (black lines), and a model reflecting translation of static PIP2 density distributions (blue dashed lines). Differences between the integrated density and the static translation model indicates regions of enrichment and depletion of PIP2 lipids during RAS–RAS association. (D) Preferential binding coefficients showing the enrichment of PIP2 among lipids in contact with RAS, dPIP2, as a function of the number of PIP2 per CG patch, shown separately for RAS dimers and monomers. The larger dPIP2 values for dimers indicate that dimers formation is favored by higher PIP2 concentration.

Lipid Dependence of RAS Orientation and Competence for Effector Binding

RAS is a closely tethered peripheral membrane protein, and the orientation of the G-domain controls access to its effector binding site through reversible membrane-based occlusion (67, 68). We define the orientation of the RAS G-domain with respect to the membrane surface (69). Briefly, a tilt angle defines the deflection of the long axis of helix 5 from the bilayer normal, and a rotation angle defines the direction in which that tilt occurs (Fig. 5). CG simulations with one RAS reveal statistical confinement of RAS to specific orientations; moreover, Markov-state modeling of RAS orientation reveals three kinetic states of RAS orientation in which tilting brings G-domain helices 3 to 5 (α-state) or β-sheets 1 to 3 and switch I (β and β′ states) toward the membrane (Fig. 5 and ).

Fig. 5.

RAS orientations. (A) Definition of tilt and rotation angles based on a reference orientation with the long axis of G-domain helix 5 perpendicular to the membrane surface. (B) Kinetic states of G-domain orientation. Dots represent Markov state modeling (MSM) microstates and color associates sampled orientations with the nearest microstate. Blue isocontours define the likelihood of membrane-based occlusion of RAF binding. (C) G-domain disposition in simulations with one RAS. (D–G) Representative configurations of α, β, and β′ states illustrating orientation-dependent competence for effector binding. (H) MSM state populations and transition rates. (I) State-specific distributions of the displacement of the H166 backbone bead from the bilayer center, dzH166. (J) G-domain tilt vs. dzH166.

The most populated RAS orientations are incompatible with effector binding (Fig. 5). We use the RBD of RAF (PDB ID code 4G0N) as the relevant effector and observe that α-states project the effector binding interface of RAS away from the membrane (Fig. 5 ), possibly positioning RAS to bind its effectors (59, 70). Conversely, effector binding is occluded in most of β′- and half of the β-state orientations (Fig. 5 ). Previously described exposed and occluded orientations (59, 70, 71) are included in the α and β/β′ states, respectively, and an experimentally driven model of RAS–RAF association straddles the dividing line between the α and β states (70). HMM analysis indicates that direct transitions between β′ and α are relatively rare (Fig. 5 and ) and interconversion between these states typically occurs via the β-state intermediate.

Consistent with the macro simulation, CG simulations reveal distinct RAS lipid fingerprints, which exhibit spatially complex patterns of lipid enrichment and depletion around the G-domain, HVR, and farnesyl anchor (). A variety of protein–lipid interactions contribute to the modulation of G-domain orientation by PM lipids, and distinct lipid fingerprints observed for the different RAS orientational states (Fig. 6) are specific enough that a convolutional neural network can predict RAS states with an accuracy of 82% based on inner leaflet spatial lipid distributions (Fig. 6 and ). Among the complex lipid signatures, regression analysis shows that enrichment/depletion of different lipids species affects the likelihood of transitions between states (Fig. 6 and ). For example, changes in the number of PIP2 lipids interacting with the G-domain emerge as a good predictor of transitions between β and β′ states, an orientational change that can modulate effector occlusion (Fig. 5).

Fig. 6.

RAS–lipid interactions. (A) Relative lipid densities around monomeric RAS. The center of mass of the G-domain is at the origin and the backbone bead of C-terminal HVR residue C185 is on the positive x axis. (B) The confusion matrix shows the accuracy of ML predictions of the RAS state based on local lipid compositions. Each row corresponds to an actual state (computed using tilt/rotation) and each column to a predicted state; off-diagonal elements are the errors made by the prediction. (C) Local lipid composition is predictive for transitions between states. Transition probabilities increase when noted lipid concentrations (+) increase and (−) decrease (see ). (D) RAF-occlusion vs. PIP2 content.

Significantly, occlusion of the RAF-binding interface by the membrane is relieved by PIP2 at molar concentrations above 2% (Fig. 6). The kinetic barrier between β and β′ states is related to a specific PIP2 binding site at the junction between the G-domain and the HVR, near H166 (Fig. 6). The exchange between metastable β and β′ states (Fig. 5 ), appears to be controlled by reversible membrane adhesion of this region of the HVR (Fig. 5), which must disengage to permit extensive G-domain tilting and transition from β to β′ (Fig. 5). In this work, we identify molecular details by which lipid–RAS interactions regulate occlusion, revealing a plausible strategy for inhibition of RAS function. Indeed, the identification of a small molecule that stabilizes KRAS in the β-state and prevents RAF binding validates this approach and motivates the search for more potent molecules (71).

Discussion

The unique capabilities of MuMMI to operate and provide insights at two vastly different temporal and spatial scales in a coordinated manner represent a fundamentally unique technology in computational biology. Furthermore, this study highlights the necessity of such a capability by demonstrating the importance of considering membrane proteins and lipids collectively, as a single unit with two highly interdependent components. Only through the combination of the macroscale sampling of diverse local lipid environments with the microscale dynamics at the molecular level can we observe these distinct RAS–lipid interaction patterns, which were subsequently used to design the supporting experiments.

The combination of simulation and experiment reveals that RAS not only remodels its local lipid environment, but that it is itself regulated by these distinct lipid patterns. For example, an increased concentration of GTP-loaded KRAS4b is associated with PIP2 enrichment and cholesterol depletion, while PIP2 drives functionally relevant orientational changes in RAS. The differential influence of lipid fingerprints on interactions of RAS with regulators, such as SOS (72), and effectors such as RAF, PI3K, and RALGDS, are unknown, but there is clear potential for lipid patterning to define recruitment sites for discrete effectors and thereby segregate RAS-mediated signal transduction pathways. For membrane proteins in general, there is mounting evidence from both simulation and experiment for these types of “symbiotic” relationships between proteins and their surrounding membrane environment (45, 73, 74), establishing another mechanism for the cell to impose fine control over molecular function.

The current work finds RAS–RAS interactions to be interfacially nonspecific and validates all previously proposed interfaces as part of a broad ensemble of possible interactions. This interfacial promiscuity suggests that RAS multimer formation is mediated by lipids rather than specific interfacial contacts, and leads us to conclude that nonenzymatically relevant mutations, some of which impact tumor growth (11), may impact the dynamics of RAS multimer formation through changes in protein–lipid interactions.

MuMMI uses cohesive multiscale coupling, which in this application reveals both the broad scope and fine details of membrane remodeling that underlies functionally relevant RAS–lipid dynamics. More importantly, MuMMI has been designed to be generalizable and, as such, opens new avenues for potential research and expands the range of questions that can be addressed using biomolecular simulation. With a preinvestment in macro model development or model tuning and significant computational resources, a MuMMI-style campaign could be run instead of nearly any simulation ensemble, resulting in a queryable broad ensemble over a large set of relevant conditions, instead of predefined sampling along a given reaction coordinate.

Materials and Methods

An overview of MuMMI is given in Multiscale Simulation Framework and a detailed description of MuMMI is provided in . The model generation and parameterization are described in . Simulation analysis is described in and experimental methods in . MuMMI is composed of numerus subcomponents, both freely available third-party applications and custom codes that all have been made available. Custom codes are: ddcMD (https://github.com/LLNL/ddcMD), ddcMDconverter (https://github.com/LLNL/ddcMDconverter), Maestro (https://github.com/LLNL/maestrowf), Flux (https://github.com/flux-framework), PyTarIdx (https://github.com/LLNL/pytaridx), DynIm (github.com/LLNL/dynim), MuMMI-core (that contains the workflow manager, https://github.com/mummi-framework/mummi-core), and MuMMI-ras (the remaining application components, https://github.com/mummi-framework/mummi-ras).