In this work, a model for anisotropic interactions between proteins and cellular membranes is proposed for large-scale continuum simulations. The framework of the model is based on dynamic density functional theory, which provides a formalism to describe the lipid densities within the membrane as continuum fields while still maintaining the fidelity of the underlying molecular interactions. Within this framework, we extend recent results to include the anisotropic effects of protein–lipid interactions. As applications, we consider two membrane proteins of biological interest: a RAS–RAF complex tethered to the membrane and a membrane embedded G protein-coupled receptor. A strong qualitative and quantitative agreement is found between the numerical results and the corresponding molecular dynamics simulations. Combining the scope of continuum level simulations with the details from molecular level particle simulations enables research into protein–membrane behaviors at a more biologically relevant scale, which crucially can also be accessed via experiment.

Topics
Molecular dynamics, Dynamic density functional theory, Anisotropic interactions, Lipid bilayer, Proteins

Interactions in many-body systems are often assumed to be isotropic (i.e., dependent only on the distance); however, certain applications require a more detailed description. In particular, we explore the dynamics of membrane proteins by extending the work of Stanton et al.1 to account for the anisotropic interactions between proteins and the underlying cellular membrane. We demonstrate how our model reproduces the given density fields of the surrounding lipids and compare the expressive complexity gained by using anisotropic potentials for two types of complex proteins on the cellular membrane.

Considerable research has been dedicated to developing the theoretical description of liquids within the framework of statistical mechanics.2–5 In most applications of this research, the resulting derived quantities, such as correlation functions and interaction potentials, are expressed as radial functions assuming isotropic interactions, which is a natural choice for simple liquids. While those assumptions simplify the solutions to these mathematical models, they are not required, and such approximations may not be appropriate for fluid structures such as cellular membranes. To focus on the specific lipid–protein interactions, without adding further complexity and variables, these systems are currently maintained in a flat environment without membrane undulations or curvature. Furthermore, while certain investigations have explored the interactions between proteins and membrane curvature,6–8 an accurate generalization of deformation energies to multi-component membranes still remains an open question.9,10

Most biophysical phenomena of interest occur over much larger timescales than those associated with atomic processes. As such, continuum models provide a potential avenue for simulating biophysical systems due to their ability to access length scales on the order of microns and time scales on the order of seconds. Currently, these spatiotemporal scales are unachievable even with highly specialized supercomputers, such as Anton 3,11 or with particle-based molecular dynamics (MD) simulations at the coarse-grained level, where groups of atoms are represented with an interaction site.12 

Our systems of interest are cellular membranes, which comprise a lipid bilayer that contains many different types of membrane proteins. These membrane proteins are the most common targets for therapeutic drugs. It is estimated that ∼70% of all FDA approved drugs target the membrane proteins.13 The lipid bilayer itself is formed from hundreds of different types of lipids,14,15 each with subtly different properties. These different lipid types can affect protein function both directly through protein–lipid interactions and indirectly through changes in local membrane properties.16–18 For instance, studies show that if you increase the cholesterol concentration in the membrane, it can change the rate of activity of some membrane proteins.19 The interactions between the protein and the lipids are a complex, almost symbiotic relationship where the lipids around a protein can influence the protein behavior, yet the protein can also recruit and enrich certain lipids around itself, thereby defining its own immediate membrane environment. Recent computational work has demonstrated this influence of the proteins to reorganize local lipids around themselves, producing distinct spatial enrichment and/or depletion of the different lipids (a “lipid fingerprint”) that is quite distinct from one protein family to the next.20,21

One of the primary drawbacks of continuum models is that many of the details of the inter-molecular interactions are often lost. We attempt to overcome this deficiency by extending a dynamic density functional theory (DDFT) model developed in our previous study,1 in which a continuum description of the lipid density fields evolves based on correlation information trained on molecular dynamics simulations. This model enabled the radially symmetric replication of these lipid density patterns, but the distinct spatial lipid fingerprints around the membrane proteins require the implementation of an anisotropic model. These fingerprints have influences on the protein function and can define the propensity for protein aggregation. Furthermore, the anisotropic nature of the fingerprints also helps to drive specific interface interactions between proteins, a key feature for correct biological function and a potential target for therapeutics. Thus, being able to correctly recapture these complex lipid patterns is critical to establishing a more accurate protein behavior.

To focus on the specific lipid–protein interactions, without adding further complexity and variables, these systems are currently maintained in a flat environment without membrane undulations or curvature. Furthermore, while certain investigations have explored the interactions between proteins and membrane curvature,6–8 an accurate generalization of deformation energies to multi-component membranes still remains an open question.9,10

The remainder of this paper is organized as follows: The basic DDFT model is presented in Sec. II, along with its generalization to incorporate anisotropic lipid–protein interactions, as detailed in Sec. III. The methods of the underlying molecular dynamics simulations are described in Sec. IV. In Sec. V, we present two applications of this model describing the lipid–protein interactions for (1) a peripheral membrane protein complex formed with RAS, RAS-binding domain (RBD), and cysteine-rich domain (CRD) of RAF and (2) an integral membrane protein, β2AR, which is a member of a pharmacologically important family of G-protein coupled receptor (GPCR) proteins. Finally, concluding remarks and possible research questions that can be explored with the presented model are given in Sec. VI.

Dynamic density functional theory (DDFT) is a non-equilibrium extension of (classical) density functional theory in which the underlying free energy functional is constructed through a correlation expansion about a uniform density.22,23 Within this model, the density of each lipid species ρi(r, t) will evolve according to the generalized diffusion equation,
(1)
where Fρ(r,t) is the free energy functional in terms of the set of lipid densities ρ(r, t) = {ρi(r, t)}, β = 1/kBT is the inverse thermal energy, Di is the diffusivity of the ith lipid species, and δ denotes a variational derivative. Note that stochastic effects can also be added, as shown in Refs. 1 and 23–25. For a system of proteins embedded in (or tethered to) a lipid bilayer, the free energy can generically be expressed in the following form:
(2)
where D is the system domain and fℓℓ, fℓp, and fpp are the lipid–lipid, lipid–protein, and protein–protein energy densities, respectively. As the proteins are fairly isolated and are not highly collisional with each other, they are treated individually rather than in terms of a density field.
The primary input to the DDFT model is the underlying correlation information from the microscopic interactions between the lipids and proteins, which enter through the lipid–lipid, lipid–protein, and protein–protein energy densities. Within the Ramakrishnan–Yussouff (RY) approximation,26 which truncates this expansion to two-body correlations, we can express the lipid–lipid energy density as
(3)
Here, Δρi(r)=ρi(r)−ρ̄i is the density fluctuation about its mean ρ̄i⁠, cij(r, r′) is the two-body direct correlation function (DCF) between the ith and jth lipid species, and Λ is the de Broglie wavelength, which is squared under the assumption of two dimensions. The DCFs can be calculated using atomically resolved simulations such as MD. In the limit of a homogeneous and isotropic system, these functions will be invariant under spatial translation and rotation and simplify their functional dependencies as ρi(r)→ρ̄i (with ρ̄i being a constant) and cij(r, r′) → cij(r), where r = |r − r′|. The DCFs can then be related to radial distribution functions gij(r) between the lipid headgroups through the Ornstein–Zernike (OZ) relations,27 given by
(4)
where the total correlation function is simply defined as hij(r) ≡ gij(r) − 1 and the operation (∗) denotes a convolution, which is defined as
(5)
The OZ relations of use here, as radial distribution functions (RDFs), are readily accessible from atomistic simulations.
The membrane proteins are represented as one or more “beads” with position rk(t). Given P beads, the lipid–protein energy density can be expressed as
(6)
where uik(r) is the potential of mean force (PMF) between the ith lipid species and the kth protein bead. Note that some of these PMFs are potentially zero in cases where the protein interacts with only one leaflet. The calculation of the PMF is detailed further in Sec. III. Finally, the protein–protein energy density can be expressed as
(7)
where the continuous density fields have been replaced with a distribution of point beads, expressed here using Dirac delta functions δ(r), and each sum is now only over the bead indices k. The factor of a half accounts for double-counting with the sum. As in Ref. 1, we approximate these interactions with the Mie (or Kihara) potential, which is a generalization of the Lennard-Jones potential for non-spherical molecules.28,29 Furthermore, these Mie potentials were shifted to account for the finite size of the proteins as ukk′(r) → ukk′(r + r0), where r0 represents the exclusion radius of a given protein, calculated by determining the region in which the lipid–protein energy exceeded some cutoff value (i.e., the region in the MD simulations in which lipids were excluded due to the presence of the protein). In reality, the protein–protein interactions are far more complicated; however, a thorough treatment of these interactions is beyond the scope of this work.
In the original model,1 the PMF associated with the lipid–protein energy density [Eq. (6)] was calculated by assuming that each protein bead was spherical, thus resulting in radially symmetric interactions. To extend this model to non-spherical beads, we first introduce the generalized Ornstein–Zernike relation5,30 for a heterogeneous system of N components in equilibrium,
(8)
Under the assumption of homogeneity (but not isotropy), these correlation functions will be invariant under spatial translation and simplify their functional dependencies as ρi(r)→ρ̄i⁠, hik(r, r′) → hik(r − r′), and cik(r, r′) → cik(r − r′). We can further simplify Eq. (8) by translating the coordinates such that r′ = 0, without loss of generality, to yield the Ornstein–Zernike relation for homogeneous (but still anisotropic) systems,
(9)
While the N × P equations resulting from Eq. (9) can be used to generate DCFs from the total correlation functions calculated from simulations, an additional closure is necessary to further relate these correlation functions to the PMFs that are needed to construct the lipid–protein energy density [Eq. (6)]. An appropriate closure for this system that is also consistent with the RY approximation within the DDFT model is the hypernetted chain (HNC) equation, given by
(10)
where uik(r) is the interaction potential between the ith lipid species and the kth protein bead (see  Appendix A for the derivation). Note that the form of the HNC equation is the same for both isotropic functions f(r) and anisotropic functions f(r). Inverting Eq. (10), we obtain an expression for the lipid–protein PMF,
(11)
Once each PMF is calculated, the equations of motions for each protein bead can be determined using the lipid–protein and protein–protein contributions to the free energy. To incorporate missing degrees of freedom associated with fluctuations in the cytoplasm, the equations of motion can be further generalized in the form of Langevin equations, given by
(12)
(13)
where ξk is a vector Wiener process with zero mean and the delta-correlated variance,
(14)
The relaxation rate is obtained from the fluctuation–dissipation theorem as γk = kBT/mkDk, where Dk is the diffusivity of the protein, which can be calculated from atomistic simulations through the appropriate Green–Kubo relation.31,32
A Langevin equation for the rotation angle θk(t) was also used, which is given by
(15)
where Ik is the moment of inertia of the protein, and the torque is given by τk=k̂⋅(r×Fk)⁠. As with Eq. (12), the relaxation rate αk is connected to the Wiener process ηk(t) through the relation αk=kBT/IkDk⁠, with Dk being a rotational diffusion coefficient.

For both Eqs. (12) and (15), the drag coefficients {γk, αk} were found to significantly dominate the dynamics (which is consistent with the overall DDFT model), and thus, the acceleration terms r̈k and θ̈k were neglected in the numerical simulations.

Similarly to our previous iterations of this model, we have not assigned specific interaction potentials between the proteins as the focus is on the dynamic lipid rearrangements. However, there are repulsive terms to ensure that proteins do not overlap and that their anisotropic shape is maintained. Extensive details of this repulsive term and how it is implemented can be found in  Appendix B.

To test the robustness of the DDFT model, we utilized two example membrane protein systems (namely, the RAS-RBDCRD complex and β2AR). Each system was chosen as they offer differing challenges and considerations for the nature of the protein–lipid interaction patterns that have to be reproduced within the DDFT model.

The first membrane protein system we studied is the RAS protein, which is a peripherally membrane bound small GTPase that sits on the inner surface of the cell membrane. It consists of a globular domain (the G-domain) that is anchored to the membrane via a flexible, farnesylated peptide linker. The protein plays a key role in the RAS/RAF/MAPK signaling pathway, which is commonly dysregulated in cancer.33 RAS has distinct lipid fingerprints based on its orientation with respect to the membrane,34 and more pronounced lipid fingerprints are observed when RAS is bound to the RAS-binding domain (RBD) and cysteine-rich domain (CRD) of the RAF protein.35,36 Here, the RAS-RBDCRD complex is used as an example of a peripheral membrane protein [Fig. 1(a)], whose lipid fingerprint changes depending on the orientation of the protein on the surface of the membrane. The lipid dependence of the RAS-RBDCRD protein complex and the state definitions are adopted from previous large multiscale simulation studies of RAS-RBDCRD membrane dynamics.35,37

FIG. 1.
(a) From the cytoplasmic side and (b) through the membrane, views of a peripheral membrane protein complex (the RAS-RBDCRD complex) are shown. The RAS-RBDCRD complex is peripherally embedded into a typical mammalian plasma membrane model composed of eight different lipid types asymmetrically distributed across the leaflets. (c) From the cytoplasmic side and (d) through the membrane, views of an integral membrane protein (a G-protein coupled receptor, β2AR) are shown, embedded in the typical mammalian plasma membrane model. Only backbone beads of the proteins are shown in the surface representation and are colored in orange. Lipids are represented as licorice. Cholesterols are shown in yellow. Phospholipid tails are colored in black if fully saturated (DP), dark gray if mono-unsaturated (PO), light gray if di-unsaturated (DI), and white if polyunsaturated. Phospholipid headgroups are shown as spheres and are colored according to their types: PC in light blue, PE in blue, PS in green, and PAP6 in magenta.
View largeDownload slide

(a) From the cytoplasmic side and (b) through the membrane, views of a peripheral membrane protein complex (the RAS-RBDCRD complex) are shown. The RAS-RBDCRD complex is peripherally embedded into a typical mammalian plasma membrane model composed of eight different lipid types asymmetrically distributed across the leaflets. (c) From the cytoplasmic side and (d) through the membrane, views of an integral membrane protein (a G-protein coupled receptor, β2AR) are shown, embedded in the typical mammalian plasma membrane model. Only backbone beads of the proteins are shown in the surface representation and are colored in orange. Lipids are represented as licorice. Cholesterols are shown in yellow. Phospholipid tails are colored in black if fully saturated (DP), dark gray if mono-unsaturated (PO), light gray if di-unsaturated (DI), and white if polyunsaturated. Phospholipid headgroups are shown as spheres and are colored according to their types: PC in light blue, PE in blue, PS in green, and PAP6 in magenta.

FIG. 1.
(a) From the cytoplasmic side and (b) through the membrane, views of a peripheral membrane protein complex (the RAS-RBDCRD complex) are shown. The RAS-RBDCRD complex is peripherally embedded into a typical mammalian plasma membrane model composed of eight different lipid types asymmetrically distributed across the leaflets. (c) From the cytoplasmic side and (d) through the membrane, views of an integral membrane protein (a G-protein coupled receptor, β2AR) are shown, embedded in the typical mammalian plasma membrane model. Only backbone beads of the proteins are shown in the surface representation and are colored in orange. Lipids are represented as licorice. Cholesterols are shown in yellow. Phospholipid tails are colored in black if fully saturated (DP), dark gray if mono-unsaturated (PO), light gray if di-unsaturated (DI), and white if polyunsaturated. Phospholipid headgroups are shown as spheres and are colored according to their types: PC in light blue, PE in blue, PS in green, and PAP6 in magenta.
View largeDownload slide

(a) From the cytoplasmic side and (b) through the membrane, views of a peripheral membrane protein complex (the RAS-RBDCRD complex) are shown. The RAS-RBDCRD complex is peripherally embedded into a typical mammalian plasma membrane model composed of eight different lipid types asymmetrically distributed across the leaflets. (c) From the cytoplasmic side and (d) through the membrane, views of an integral membrane protein (a G-protein coupled receptor, β2AR) are shown, embedded in the typical mammalian plasma membrane model. Only backbone beads of the proteins are shown in the surface representation and are colored in orange. Lipids are represented as licorice. Cholesterols are shown in yellow. Phospholipid tails are colored in black if fully saturated (DP), dark gray if mono-unsaturated (PO), light gray if di-unsaturated (DI), and white if polyunsaturated. Phospholipid headgroups are shown as spheres and are colored according to their types: PC in light blue, PE in blue, PS in green, and PAP6 in magenta.

Close modal

In brief, in Ref. 35, ∼35 000 MD simulations were run using the Martini 212 coarse-grained (CG) force field. These simulations consisted of RAS and RAS-RBDCRD proteins on top of an eight-component plasma membrane mimetic38 and totaled nearly 100 ms of aggregated simulation time. Here, ∼7200 simulations with one RAS-RBDCRD protein complex were used and with frames saved every 2 ns, which corresponds to roughly 20 ms of aggregated simulation time. These frames were used to evaluate the state-dependent lipid fingerprints.

Each frame was separated into one of three defined protein orientational states (“za,” in which the CRD does not contact the membrane, or “ma” or “mb”) based on different rotation and tilt of the RAS G-domain relative to the membrane.35 The frames were aligned such that the RAS G-domain center of mass was positioned at the origin and then rotated so that the hydrophobic loops of the CRD were positioned on the +x axis. The aligned data were analyzed to define the average lipid enrichment in each of the different states.

The second membrane protein system we studied is the β2 adrenergic receptor (β2AR), which is a membrane of the G protein-coupled receptor (GPCR) family of proteins. It is a key cellular receptor that has many functions within the body, primarily as a response to binding adrenaline. In contrast to RAS, which sits on the surface of the membrane, GPCRs are typical integral membrane proteins [Fig. 1(b)] that completely span the membrane and interact with both leaflets of the lipid bilayer. Upon binding of a ligand on the outside of a cell, GPCRs undergo a conformational change from an inactive to an active state. This allows for binding of downstream signaling proteins inside the cell. The high-resolution structures of β2AR in active and inactive states are available within the Protein Data Bank: entries 4LDO39 and 2RH1,40 respectively. Thus, we have the opportunity to assess how the lipid fingerprints differ between the inactive and active conformational states of the protein.

The simulation data are taken from a previous paper investigating these conformational-dependent fingerprints in detail.41 In brief, for each protein state, 20 independent simulation setups were prepared and run using the Martini 2 CG model. Following a short equilibration, a 6 μs-long simulation was carried out. The first 1 μs part of each trajectory was discarded for further equilibration.

The aggregated 100 μs of data was aligned such that the center of mass of all seven transmembrane helices is located at the origin and rotated such that the top part of transmembrane helix 1 (residues 32–34) is aligned with the +y axis.

For both the peripheral RAS-RBDCRD complex and the integral β2AR, the bilayers in the MD simulations were maintained mostly flat, with large membrane undulations and curvature prevented (using methods described in Ref. 41). Due to the flat nature of the bilayer, lipid fingerprints of the protein systems were computed by time-averaging the lipid counts on a two-dimensional grid placed on the xy plane for each leaflet. Lipids were assigned to individual bins using the x and y positions of the first tail bead (R1 for cholesterol and C1A, D1A, and T1A for others). The density of each lipid is averaged over the distinct datasets for the different orientational (RAS-RBDCRD) or conformational (β2AR) states, and all density plots are shown at 10 × 10 nm2.

Having used the 2D lipid density data from the CG MD simulations to parameterize the anisotropic protein–lipid potentials within the DDFT model, relaxation continuum simulations of each protein system were carried out. The resulting 2D lipid densities were then generated from the continuum simulation and compared to the “ground truth” CG MD data to determine how robust the DDFT model is at reproducing the various protein–lipid molecular features. Specific details are discussed below, and comparisons between all lipid species are shown in Figs. 12–16 in  Appendix C.

For the peripheral membrane protein complex RAS-RBDCRD [Figs. 1(a) and 1(b)], the protein–lipid interactions are highly dependent on the position and orientation of the complex with respect to the membrane and can be broadly classified into three membrane states,35 each with clearly different lipid fingerprints (Fig. 2, bottom row). Many of the lipid fingerprint patterns calculated from the CG simulation data are distinctly anisotropic in nature, demonstrating the need for a similar anisotropic implementation within the continuum simulations. Using the state dependent anisotropic interaction potentials, the DDFT continuum simulations can capture the 2D density profiles with great precision (Fig. 2, top row). Lipids that display more generalized anisotropic features at the CG level (such as DPSM, Fig. 3) also demonstrate the same patterns in the continuum simulations. Encouragingly, specific lipids that show more specific structured features in CG simulations (such as the “rings” of DIPE density, Fig. 4) are also recapitulated at the continuum level. These comparisons reveal a remarkable qualitative agreement between many differing, nuanced types of lipid fingerprint features.

FIG. 2.
Comparison of the PAP6 lipid density produced by the MD simulations (left) and the DDFT model (right) for the three conformational states of RAS-RBDCRD: states ma (top), mb (middle), and za (bottom).
View largeDownload slide

Comparison of the PAP6 lipid density produced by the MD simulations (left) and the DDFT model (right) for the three conformational states of RAS-RBDCRD: states ma (top), mb (middle), and za (bottom).

FIG. 2.
Comparison of the PAP6 lipid density produced by the MD simulations (left) and the DDFT model (right) for the three conformational states of RAS-RBDCRD: states ma (top), mb (middle), and za (bottom).
View largeDownload slide

Comparison of the PAP6 lipid density produced by the MD simulations (left) and the DDFT model (right) for the three conformational states of RAS-RBDCRD: states ma (top), mb (middle), and za (bottom).

Close modal
FIG. 3.
Comparison between the DPSM lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.
View largeDownload slide

Comparison between the DPSM lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.

FIG. 3.
Comparison between the DPSM lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.
View largeDownload slide

Comparison between the DPSM lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.

Close modal
FIG. 4.
Comparison between the DIPE lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.
View largeDownload slide

Comparison between the DIPE lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.

FIG. 4.
Comparison between the DIPE lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.
View largeDownload slide

Comparison between the DIPE lipid density produced by the MD simulations (left) and the DDFT model (right) when RAS-RBDCRD is in the mb state.

Close modal

Furthermore, and importantly, only a few CPU hours of DDFT simulation wall time are needed to capture the lipid protein fingerprint, which is several orders of magnitude faster than comparable CG simulations. In addition, we can now add in multiple proteins and investigate larger systems, parameter sweeps, etc. (for example, Fig. 9). Looking at the interactions and relationships between fingerprints and proteins is much more practical with this method than with the much more computationally expensive MD simulations. In addition to facilitating more accurate studies at the continuum scale, using the DDFT simulations as a “pre-run” step can greatly reduce the required equilibrium time at the computationally more expensive CG simulation resolution, such as was done in using radially averaged continuum potentials in the Multiscale Machine-Learned Modeling Infrastructure (MuMMI),1,34,35,42 which can now be done with spatially resolved fingerprints.

The MD simulation of the integral transmembrane protein β2AR revealed very specific spatial patterns of lipid density around the protein in both the inner and outer leaflets of the membrane and in both the active and inactive states. These fingerprints can be unique for each family of proteins or even specific types of GPCRs21 and are distinct between the inner and outer leaflets. The DDFT continuum simulations do a remarkable job of reproducing various different features observed in the MD simulations.

In contrast to RAS-RBDCRD simulations where the proteins lay on the surface of the membrane, the β2AR proteins completely span both leaflets of the membrane and, as such, have an excluded volume within the membrane. The shape of these protein “voids” in the lipid densities is captured within both leaflets of the DDFT continuum simulations, reproducing the subtle differences in the active and inactive states of β2AR due to their conformational changes.

As with the RAS-RBDCRD systems, the more generalized, global lipid density patterns are reproduced in the β2AR simulations. For instance, the POPE density that is decreased to the north/south of the protein and increased to the east/west of the protein (Fig. 5) is nicely reproduced.

FIG. 5.
Comparison of MD (left) and DDFT (right) simulations of POPE densities in the inner leaflet around the inactive state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) simulations of POPE densities in the inner leaflet around the inactive state β2AR.

FIG. 5.
Comparison of MD (left) and DDFT (right) simulations of POPE densities in the inner leaflet around the inactive state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) simulations of POPE densities in the inner leaflet around the inactive state β2AR.

Close modal

In addition, very different types of lipid density distributions are also captured in the anisotropic DDFT model. Cholesterol is known to tightly associate with very specific binding site regions at the membrane–protein interface of the GPCRs,40,43,44 such as β2AR. This can be observed in the CG lipid densities as defined, distinct “hot spots” of small regions of high density in the cholesterol plots. These specific cholesterol hot spots are recaptured extremely well in the DDFT densities (Fig. 6). Both the location and relative intensities of these hot spots are achieved. Thus, the DDFT model can capture both subtle, generalized patterns and very specific high-value densities associated with lipid binding to the protein.

FIG. 6.
Comparison of MD (left) and DDFT (right) cholesterol densities in the inner leaflet around the active state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) cholesterol densities in the inner leaflet around the active state β2AR.

FIG. 6.
Comparison of MD (left) and DDFT (right) cholesterol densities in the inner leaflet around the active state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) cholesterol densities in the inner leaflet around the active state β2AR.

Close modal

Finally, several interesting features of the bilayer density are reproduced in the DDFT model. The total lipid densities of each leaflet are maintained, demonstrated by the slight difference between the inner and outer leaflets (Fig. 7). This leaflet density difference is seen in the CG simulations and is expected due to the different properties of the lipids in the asymmetric bilayer model.38 In some regions around the protein, the densities of certain lipids (such as DIPE, Fig. 8) display somewhat hexagonal packing patterns. This is likely due to the order-induced packing of lipids around the protein and the associated protein-bound cholesterol. Again, with high-resolution simulations of the DDFT model, these specific detailed density patterns are captured.

FIG. 7.
Total lipid densities around the active state β2AR in the outer (top) and inner (bottom) leaflets compared between the MD (left) and DDFT (right) simulations. To emphasize how the total lipid densities surround the protein, the protein center of mass is shown as a black dot and the positions of the seven transmembrane helices are shown as colored circles.
View largeDownload slide

Total lipid densities around the active state β2AR in the outer (top) and inner (bottom) leaflets compared between the MD (left) and DDFT (right) simulations. To emphasize how the total lipid densities surround the protein, the protein center of mass is shown as a black dot and the positions of the seven transmembrane helices are shown as colored circles.

FIG. 7.
Total lipid densities around the active state β2AR in the outer (top) and inner (bottom) leaflets compared between the MD (left) and DDFT (right) simulations. To emphasize how the total lipid densities surround the protein, the protein center of mass is shown as a black dot and the positions of the seven transmembrane helices are shown as colored circles.
View largeDownload slide

Total lipid densities around the active state β2AR in the outer (top) and inner (bottom) leaflets compared between the MD (left) and DDFT (right) simulations. To emphasize how the total lipid densities surround the protein, the protein center of mass is shown as a black dot and the positions of the seven transmembrane helices are shown as colored circles.

Close modal
FIG. 8.
Comparison of MD (left) and DDFT (right) DIPE densities in the inner leaflet around the inactive state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) DIPE densities in the inner leaflet around the inactive state β2AR.

FIG. 8.
Comparison of MD (left) and DDFT (right) DIPE densities in the inner leaflet around the inactive state β2AR.
View largeDownload slide

Comparison of MD (left) and DDFT (right) DIPE densities in the inner leaflet around the inactive state β2AR.

Close modal

To quantify the improvement in lipid densities around a protein compared to our previous 1D model,1 we used a circularly averaged version of the 2D input densities to produce equivalent 1D isotropic potentials. These potentials were used to relax lipids around β2AR in the active state (Fig. 17). To determine the convergence of the relaxation and to define an error norm between the ground truth CG data and the two different DDFT models, a 10 × 10 nm2 square centered around the protein was extracted, interpolated to the same grid, and compared. The L2 norm of these differences was calculated and divided by the L2 norm for the CG data. This analysis resulted in a value of 0.16 for the anisotropic (2D) potentials, compared to a value of 0.42 for the isotropic (1D) potentials, thus a factor of 2.6× improvement by going from isotropic to anisotropic. This can qualitatively be observed in Fig. 17, where specific spatial regions of lipid enrichment are radially averaged, and thus, the signal is lost.

We have demonstrated the efficacy of modeling anisotropic interactions within the DDFT formalism. As an application, we applied this generalization to membrane proteins in which the naturally anisotropic interactions between proteins and cellular membranes are of biological interest.

The DDTF model presented has shown a very promising ability to reproduce a variety of different important features observed in the MD simulations. General, global density trends match the MD systems, as well as localized high-density hot spots that result from specific protein–lipid interactions. Furthermore, bilayer density features, such as order-induced ripples and distinctly shaped voids produced by the embedding of different protein conformations, are both faithfully captured. The model has demonstrated applicability to both peripheral and integral membrane proteins (each presenting their own challenges) and the robustness to reproduce density patterns that can be dependent on both protein conformation and orientation.

The anisotropic nature of the lipid environment around membrane proteins is a key factor that can contribute to the physiological function of the protein. These immediate environments have emphasized importance with regard to protein aggregation45 and higher-order protein–protein interactions. As such, the ability to accurately reconstitute both the macro- and micro-details of the lipid densities within the DDFT model allows for the investigation of these large-scale phenomena that require multiple proteins, increased system size, and a much longer sampling time. Figure 9, for example, displays the spatial distribution of proteins and cholesterol on a micron scale, allowing for the comparison of both far-reaching changes in lipid/protein densities and the associated changes in the fine detail of the lipid fingerprints that are maintained around the proteins. The computational burden to extensively sample the systems using particle-based CG MD simulations is often prohibitively high. Thus, a huge benefit of the DDFT model is the capability to explore this space quickly and with a vastly smaller computational overhead.

FIG. 9.
Micron scale simulation of 300 β2AR proteins (displayed as red points) showing the broad changes in global cholesterol density, with greater protein populations in regions of higher cholesterol (top left). A 300 nm zoomed-in region (top right) highlights the various protein orientations observed (red arrows). Finally, zooming in to a 10 nm region around a single protein illustrates the fine detail of spatial lipid arrangements (bottom row).
View largeDownload slide

Micron scale simulation of 300 β2AR proteins (displayed as red points) showing the broad changes in global cholesterol density, with greater protein populations in regions of higher cholesterol (top left). A 300 nm zoomed-in region (top right) highlights the various protein orientations observed (red arrows). Finally, zooming in to a 10 nm region around a single protein illustrates the fine detail of spatial lipid arrangements (bottom row).

FIG. 9.
Micron scale simulation of 300 β2AR proteins (displayed as red points) showing the broad changes in global cholesterol density, with greater protein populations in regions of higher cholesterol (top left). A 300 nm zoomed-in region (top right) highlights the various protein orientations observed (red arrows). Finally, zooming in to a 10 nm region around a single protein illustrates the fine detail of spatial lipid arrangements (bottom row).
View largeDownload slide

Micron scale simulation of 300 β2AR proteins (displayed as red points) showing the broad changes in global cholesterol density, with greater protein populations in regions of higher cholesterol (top left). A 300 nm zoomed-in region (top right) highlights the various protein orientations observed (red arrows). Finally, zooming in to a 10 nm region around a single protein illustrates the fine detail of spatial lipid arrangements (bottom row).

Close modal

We acknowledge Laboratory Directed Research and Development at the Lawrence Livermore National Laboratory (21-ERD-047) for funding. For computing time, we acknowledge Livermore Computing (LC) and Livermore Institutional Grand Challenge. This work was performed under the auspices of the U.S. Department of Energy (DOE) by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and under the auspices of the National Cancer Institute (NCI) by Frederick National Laboratory for Cancer Research (FNLCR) under Contract No. 75N91019D00024. Part of this work was supported by the NCI-DOE Collaboration established by the U.S. DOE and the NCI of the National Institutes of Health LLNL-JRNL-865310.

The authors have no conflicts to disclose.

T. Oppelstrup: Conceptualization (equal); Data curation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). L. G. Stanton: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). J. O. B. Tempkin: Data curation (equal); Formal analysis (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). T. N. Ozturk: Formal analysis (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). H. I. Ingólfsson: Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). T. S. Carpenter: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we derive the hypernetted chain (HNC) closure relation by formulating our system in terms of (classical) density functional theory. In equilibrium, the densities ρ(r) = {ρi(r)} will minimize the grand potential, given by
(A1)
where Fρ(r) is the free energy functional and {μi} is the set of chemical potentials acting as constant Lagrange multipliers to constrain particle numbers for each species. The free energy can be decomposed into ideal, excess, and external components as follows:
(A2)
The ideal (or non-interacting) component is given by
(A3)
where Λi=2πβℏ2/mi1/2 is the thermal de Broglie wavelength, with ℏ being the reduced Planck constant, mi being the ith particle mass, and β = 1/(kBT) being the inverse thermal energy. Furthermore, d is the dimension of the system, where we have d = 2 for the applications in this work.
The excess free energy can be expanded about the uniform states {ρ̄i} in terms of a correlation expansion,
(A4)
where the energy densities are given by
(A5)
(A6)
(A7)
Here, Δρi(r)=ρi(r)−ρ̄i is the deviation of the ith density from its mean ρ̄i⁠, and the kernels ci⋯(n)(r,r′,…) are the n-body direct correlation functions (DCFs) between n species. Note that the superscript of the two-body DCF cij(2)(r,r′) was previously omitted for brevity. In the same homogeneous limit, we expect the correlation functions to simplify their dependencies as cij(2)(r,r′)→cij(2)(r−r′)⁠, cijk(3)(r,r′,r″)→cijk(3)(r−r′,r−r″)⁠, and so on, which means that the integrals can be expressed as convolutions.
Finally, the presence of a single protein (bead) is modeled through the pairwise interactions with each lipid species,
(A8)
where ui(r) is the potential of mean force (PMF) between the protein and the ith lipid density. Note that the results will apply to lipid–protein PMFs associated with each protein (bead), and the remaining protein–protein interactions are irrelevant to the derivation of this closure. To further simplify the system, we can also truncate this expansion at two-body interactions, which is sometimes referred to as the Ramakrishnan–Yussouff (or hypernetted-chain) approximation. Not only is this approximation known to be fairly accurate for most liquids, but is also consistent with the approximations of the underlying DDFT model.
Upon the variational minimization of the grand potential with respect to ρi(r), we obtain
(A9)
In the absence of the external potential, we expect the density to be uniform and obtain an expression for the chemical potential,
(A10)
We can now express each density in terms of pair correlation functions between each species as ρi(r)=ρ̄igik(r)⁠, which also allows us to write Δρi(r)=ρ̄ihik(r)⁠. Equation (A9) can then be expressed as
(A11)
The sum on the right-hand side can be simplified using the Ornstein–Zernike relation for homogeneous (but still anisotropic) systems,
(A12)
which yields the algebraic expression
(A13)
Finally, rewriting this expression in terms of gik(r), we obtain the hypernetted-chain equation for anisotropic systems,
(A14)
Equations (A12) and (A14) together form a direct connection between each external potential ui(r) and the pair correlation functions gik(r).

Currently, our model does not capture the specifics of the protein–protein interactions. However, in future iterations, we plan to utilize umbrella sampling approaches46–48 to include these explicit interactions (although the MD simulation requirements for adequate sampling of the protein–protein interactions remain a challenge).

However, the model does have a protein–protein repulsive potential to ensure that the proteins maintain their correct cross-sectional volume within the membrane and that the proteins do not overlap. A protein perimeter is defined where the protein–lipid interactions have a sufficiently high repulsion. The interaction energy between a pair of proteins is taken as ϵ/r12, where r is the smallest separation between the perimeters of the proteins in the pair and ϵ = 1kBT. Due to the anisotropic nature of the proteins, this value of r (smallest distance between the perimeters of the protein) does indeed depend on the angle between the proteins.

In order to implement this angle-dependent distance for the repulsive potential, we calculate the known interaction perimeters of the protein. These perimeters, or protein “masks” (Fig. 10), are determined from an analysis of the lipid densities measured from the MD simulations. If the total lipid density drops significantly below the normal value of ∼1.8 lipids/nm2, we conclude that that there is a part of the protein occupying that region, and thus, another protein cannot overlap and also occupy that region.

FIG. 10.
The protein masks, which define the perimeter of the protein, are defined as the regions where the total lipid density in each leaflet drops significantly below the normal value of ∼1.8 lipids/nm2. These regions effectively represent the areas of the lipid bilayer that is occluded by the embedded membrane protein.
View largeDownload slide

The protein masks, which define the perimeter of the protein, are defined as the regions where the total lipid density in each leaflet drops significantly below the normal value of ∼1.8 lipids/nm2. These regions effectively represent the areas of the lipid bilayer that is occluded by the embedded membrane protein.

FIG. 10.
The protein masks, which define the perimeter of the protein, are defined as the regions where the total lipid density in each leaflet drops significantly below the normal value of ∼1.8 lipids/nm2. These regions effectively represent the areas of the lipid bilayer that is occluded by the embedded membrane protein.
View largeDownload slide

The protein masks, which define the perimeter of the protein, are defined as the regions where the total lipid density in each leaflet drops significantly below the normal value of ∼1.8 lipids/nm2. These regions effectively represent the areas of the lipid bilayer that is occluded by the embedded membrane protein.

Close modal

Using the different protein masks, defined for each protein state and their shape in both the inner and outer leaflets, we can implement the protein repulsion distances. Figure 11 illustrates that the way repulsion is calculated for a specific protein–protein orientation (however, the repulsion is indeed dependent upon and implemented for all combinations of protein–protein orientation angles). Based upon the angle between the two proteins, their closest distance (r) is computed. A repulsive potential is added that aims to separate the protein centers at that distance and with direction to the other center.

FIG. 11.
In the “contact distance” (black dashed line) indicated, the proteins have the same orientation. The center protein has masks for the inner and outer layers plotted in thick green and red lines, respectively. The thinner green and red lines illustrate where the inner and outer masks of the other protein would be in two different positions. The pink crosses show the centers of the potentials. Given that a protein will only interact with another protein via regions in the same leaflet, the contact distance is defined as the distance where either the green–green or red–red masks intersect.
View largeDownload slide

In the “contact distance” (black dashed line) indicated, the proteins have the same orientation. The center protein has masks for the inner and outer layers plotted in thick green and red lines, respectively. The thinner green and red lines illustrate where the inner and outer masks of the other protein would be in two different positions. The pink crosses show the centers of the potentials. Given that a protein will only interact with another protein via regions in the same leaflet, the contact distance is defined as the distance where either the green–green or red–red masks intersect.

FIG. 11.
In the “contact distance” (black dashed line) indicated, the proteins have the same orientation. The center protein has masks for the inner and outer layers plotted in thick green and red lines, respectively. The thinner green and red lines illustrate where the inner and outer masks of the other protein would be in two different positions. The pink crosses show the centers of the potentials. Given that a protein will only interact with another protein via regions in the same leaflet, the contact distance is defined as the distance where either the green–green or red–red masks intersect.
View largeDownload slide

In the “contact distance” (black dashed line) indicated, the proteins have the same orientation. The center protein has masks for the inner and outer layers plotted in thick green and red lines, respectively. The thinner green and red lines illustrate where the inner and outer masks of the other protein would be in two different positions. The pink crosses show the centers of the potentials. Given that a protein will only interact with another protein via regions in the same leaflet, the contact distance is defined as the distance where either the green–green or red–red masks intersect.

Close modal

The comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in the ma, mb, and za states as resolved in the MD and DDFT simulations is shown in Figs. 12–14, respectively. The comparison of all the lipid densities around the active and inactive states of β2AR as resolved in the MD and DDFT simulations is shown in Figs. 15 and 16, respectively. The comparison of all the inner leaflet lipid densities around the active state of β2AR as resolved in the MD simulations and DDFT models containing either anisotropic (DDTF-2D) or radially symmetric (DDFT-1D) potentials is shown in Fig. 17.

FIG. 12.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state ma as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state ma as resolved in the MD and DDFT simulations.

FIG. 12.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state ma as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state ma as resolved in the MD and DDFT simulations.

Close modal
FIG. 13.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state mb as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state mb as resolved in the MD and DDFT simulations.

FIG. 13.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state mb as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state mb as resolved in the MD and DDFT simulations.

Close modal
FIG. 14.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state za as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state za as resolved in the MD and DDFT simulations.

FIG. 14.
Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state za as resolved in the MD and DDFT simulations.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the RAS-RBDCRD complex in state za as resolved in the MD and DDFT simulations.

Close modal
FIG. 15.
Comparison of all the lipid densities around the active state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.
View largeDownload slide

Comparison of all the lipid densities around the active state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.

FIG. 15.
Comparison of all the lipid densities around the active state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.
View largeDownload slide

Comparison of all the lipid densities around the active state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.

Close modal
FIG. 16.
Comparison of all the lipid densities around the inactive state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.
View largeDownload slide

Comparison of all the lipid densities around the inactive state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.

FIG. 16.
Comparison of all the lipid densities around the inactive state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.
View largeDownload slide

Comparison of all the lipid densities around the inactive state of β2AR as resolved in the MD and DDFT simulations. The data are presented separately for the inner and outer leaflets, as well as for the total density in either leaflet.

Close modal
FIG. 17.
Comparison of all the inner leaflet lipid densities around the active state of β2AR as resolved in the MD simulations and DDFT models containing either anisotropic (DDTF-2D) or radially symmetric (DDFT-1D) potentials.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the active state of β2AR as resolved in the MD simulations and DDFT models containing either anisotropic (DDTF-2D) or radially symmetric (DDFT-1D) potentials.

FIG. 17.
Comparison of all the inner leaflet lipid densities around the active state of β2AR as resolved in the MD simulations and DDFT models containing either anisotropic (DDTF-2D) or radially symmetric (DDFT-1D) potentials.
View largeDownload slide

Comparison of all the inner leaflet lipid densities around the active state of β2AR as resolved in the MD simulations and DDFT models containing either anisotropic (DDTF-2D) or radially symmetric (DDFT-1D) potentials.

Close modal
1.
L. G.
Stanton
,
T.
Oppelstrup
,
T. S.
Carpenter
,
H. I.
Ingólfsson
,
M. P.
Surh
,
F. C.
Lightstone
, and
J. N.
Glosli
, “
Dynamic density functional theory of multicomponent cellular membranes
,”
Phys. Rev. Res.
5
(
1
),
013080
(
2023
).
https://doi.org/10.1103/physrevresearch.5.013080
Google Scholar
Crossref
Search ADS
 
2.
J. K.
Percus
 and
G. J.
Yevick
, “
Analysis of classical statistical mechanics by means of collective coordinates
,”
Phys. Rev.
110
(
1
),
1
(
1958
).
https://doi.org/10.1103/physrev.110.1
Google Scholar
Crossref
Search ADS
 
3.
J. M. J.
van Leeuwen
,
J.
Groeneveld
, and
J.
de Boer
, “
New method for the calculation of the pair correlation function. I
,”
Physica
25
(
7–12
),
792
–
808
(
1959
).
https://doi.org/10.1016/0031-8914(59)90004-7
Google Scholar
Crossref
Search ADS
 
4.
J. A.
Barker
 and
D.
Henderson
, “
What is ‘liquid?’ Understanding the states of matter
,”
Rev. Mod. Phys.
48
(
4
),
587
(
1976
).
https://doi.org/10.1103/revmodphys.48.587
Google Scholar
Crossref
Search ADS
 
5.
J.-P.
Hansen
 and
I. R.
McDonald
,
Theory of Simple Liquids: With Applications to Soft Matter
(
Academic Press
,
2013
).
Google Scholar
 
6.
J. B.
Fournier
, “
Nontopological saddle-splay and curvature instabilities from anisotropic membrane inclusions
,”
Phys. Rev. Lett.
76
(
23
),
4436
(
1996
).
https://doi.org/10.1103/physrevlett.76.4436
Google Scholar
Crossref
Search ADS
PubMed
 
7.
M.
Fošnaric
,
K.
Bohinc
,
D. R.
Gauger
,
A.
Iglic
,
V.
Kralj-Iglic
, and
S.
May
, “
The influence of anisotropic membrane inclusions on curvature elastic properties of lipid membranes
,”
J. Chem. Inf. Model.
45
(
6
),
1652
–
1661
(
2005
).
https://doi.org/10.1021/ci050171t
Google Scholar
Crossref
Search ADS
 
8.
M.
Fošnarič
,
A.
Iglič
, and
S.
May
, “
Influence of rigid inclusions on the bending elasticity of a lipid membrane
,”
Phys. Rev. E
74
(
5
),
051503
(
2006
).
https://doi.org/10.1103/physreve.74.051503
Google Scholar
Crossref
Search ADS
 
9.
W.
Helfrich
, “
Elastic properties of lipid bilayers: Theory and possible experiments
,”
Z. Naturforsch. C
28
(
11–12
),
693
–
703
(
1973
).
https://doi.org/10.1515/znc-1973-11-1209
Google Scholar
Crossref
Search ADS
PubMed
 
10.
A. J.
Sodt
,
R. M.
Venable
,
E.
Lyman
, and
R. W.
Pastor
, “
Nonadditive compositional curvature energetics of lipid bilayers
,”
Phys. Rev. Lett.
117
(
13
),
138104
(
2016
).
https://doi.org/10.1103/physrevlett.117.138104
Google Scholar
Crossref
Search ADS
PubMed
 
11.
D. E.
Shaw
,
P. J.
Adams
,
A.
Azaria
,
J. A.
Bank
,
B.
Batson
,
A.
Bell
,
M.
Bergdorf
,
J.
Bhatt
,
J. A.
Butts
,
T.
Correia
 et al, “
Anton 3: Twenty microseconds of molecular dynamics simulation before lunch
,” in
Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis
(
Association for Computing Machinery
,
2021
), pp.
1
–
11
.
Google Scholar
Crossref
Search ADS
 
12.
S. J.
Marrink
,
H. J.
Risselada
,
S.
Yefimov
,
D. P.
Tieleman
, and
A. H.
de Vries
, “
The MARTINI force field: Coarse grained model for biomolecular simulations
,”
J. Phys. Chem. B
111
(
27
),
7812
–
7824
(
2007
).
https://doi.org/10.1021/jp071097f
Google Scholar
Crossref
Search ADS
 
13.
J. P.
Overington
,
B.
Al-Lazikani
, and
A. L.
Hopkins
, “
How many drug targets are there?
,”
Nat. Rev. Drug Discovery
5
(
12
),
993
–
996
(
2006
).
https://doi.org/10.1038/nrd2199
Google Scholar
Crossref
Search ADS
 
14.
J. L.
Sampaio
,
M. J.
Gerl
,
C.
Klose
,
C. S.
Ejsing
,
H.
Beug
,
K.
Simons
, and
A.
Shevchenko
, “
Membrane lipidome of an epithelial cell line
,”
Proc. Natl. Acad. Sci. U. S. A.
108
(
5
),
1903
–
1907
(
2011
).
https://doi.org/10.1073/pnas.1019267108
Google Scholar
Crossref
Search ADS
PubMed
 
15.
G.
van Meer
,
D. R.
Voelker
, and
G. W.
Feigenson
, “
Membrane lipids: Where they are and how they behave
,”
Nat. Rev. Mol. Cell Biol.
9
(
2
),
112
–
124
(
2008
).
https://doi.org/10.1038/nrm2330
Google Scholar
Crossref
Search ADS
PubMed
 
16.
V.
Corradi
,
B. I.
Sejdiu
,
H.
Mesa-Galloso
,
H.
Abdizadeh
,
S. Y.
Noskov
,
S. J.
Marrink
, and
D. P.
Tieleman
, “
Emerging diversity in lipid–protein interactions
,”
Chem. Rev.
119
(
9
),
5775
–
5848
(
2019
).
https://doi.org/10.1021/acs.chemrev.8b00451
Google Scholar
Crossref
Search ADS
PubMed
 
17.
H. I.
Ingólfsson
,
M. N.
Melo
,
F. J.
van Eerden
,
C.
Arnarez
,
C. A.
Lopez
,
T. A.
Wassenaar
,
X.
Periole
,
A. H.
de Vries
,
D. P.
Tieleman
, and
S. J.
Marrink
, “
Lipid organization of the plasma membrane
,”
J. Am. Chem. Soc.
136
(
41
),
14554
–
14559
(
2014
).
https://doi.org/10.1021/ja507832e
Google Scholar
Crossref
Search ADS
 
18.
J. A.
Lundbaek
,
S. A.
Collingwood
,
H. I.
Ingólfsson
,
R.
Kapoor
, and
O. S.
Andersen
, “
Lipid bilayer regulation of membrane protein function: Gramicidin channels as molecular force probes
,”
J. R. Soc. Interface
7
(
44
),
373
–
395
(
2010
).
https://doi.org/10.1098/rsif.2009.0443
Google Scholar
Crossref
Search ADS
 
19.
M.
Gutierrez
,
K. S.
Mansfield
, and
N.
Malmstadt
, “
The functional activity of the human serotonin 5-HT1A receptor is controlled by lipid bilayer composition
,”
Biophys. J.
110
(
11
),
2486
–
2495
(
2016
).
https://doi.org/10.1016/j.bpj.2016.04.042
Google Scholar
Crossref
Search ADS
PubMed
 
20.
V.
Corradi
,
E.
Mendez-Villuendas
,
H. I.
Ingólfsson
,
R.-X.
Gu
,
I.
Siuda
,
M. N.
Melo
,
A.
Moussatova
,
L. J.
DeGagné
,
B. I.
Sejdiu
,
G.
Singh
 et al, “
Lipid–protein interactions are unique fingerprints for membrane proteins
,”
ACS Cent. Sci.
4
(
6
),
709
–
717
(
2018
).
https://doi.org/10.1021/acscentsci.8b00143
Google Scholar
Crossref
Search ADS
PubMed
 
21.
B. I.
Sejdiu
 and
D. P.
Tieleman
, “
Lipid–protein interactions are a unique property and defining feature of G protein-coupled receptors
,”
Biophys. J.
118
(
8
),
1887
–
1900
(
2020
).
https://doi.org/10.1016/j.bpj.2020.03.008
Google Scholar
Crossref
Search ADS
PubMed
 
22.
R.
Evans
, “
The nature of the liquid–vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids
,”
Adv. Phys.
28
(
2
),
143
–
200
(
1979
).
https://doi.org/10.1080/00018737900101365
Google Scholar
Crossref
Search ADS
 
23.
U. M. B.
Marconi
 and
P.
Tarazona
, “
Dynamic density functional theory of fluids
,”
J. Chem. Phys.
110
(
16
),
8032
–
8044
(
1999
).
https://doi.org/10.1063/1.478705
Google Scholar
Crossref
Search ADS
 
24.
D. S.
Dean
, “
Langevin equation for the density of a system of interacting Langevin processes
,”
J. Phys. A: Math. Gen.
29
(
24
),
L613
(
1996
).
https://doi.org/10.1088/0305-4470/29/24/001
Google Scholar
Crossref
Search ADS
 
25.
K.
Kawasaki
, “
Microscopic analyses of the dynamical density functional equation of dense fluids
,”
J. Stat. Phys.
93
,
527
–
546
(
1998
).
https://doi.org/10.1023/b:joss.0000033240.66359.6c
Google Scholar
Crossref
Search ADS
 
26.
T. V.
Ramakrishnan
 and
M.
Yussouff
, “
First-principles order-parameter theory of freezing
,”
Phys. Rev. B
19
(
5
),
2775
(
1979
).
https://doi.org/10.1103/physrevb.19.2775
Google Scholar
Crossref
Search ADS
 
27.
L. S.
Ornstein
 and
F.
Zernike
, “
Accidental deviations of density and opalescence at the critical point of a single substance
,”
Proc. Akad. Sci.
17
,
793
(
1914
).
Google Scholar
 
28.
G.
Mie
, “
Zur kinetischen Theorie der einatomigen Körper
,”
Ann. Phys.
316
(
8
),
657
–
697
(
1903
).
https://doi.org/10.1002/andp.19033160802
Google Scholar
Crossref
Search ADS
 
29.
T.
Kihara
, “
The second virial coefficient of non-spherical molecules
,”
J. Phys. Soc. Jpn.
6
(
5
),
289
–
296
(
1951
).
https://doi.org/10.1143/jpsj.6.289
Google Scholar
Crossref
Search ADS
 
30.
S.
Jorge
,
E.
Lomba
, and
J. L. F.
Abascal
, “
An inhomogeneous integral equation for the triplet structure of binary liquids
,”
J. Chem. Phys.
114
(
8
),
3562
–
3569
(
2001
).
https://doi.org/10.1063/1.1342813
Google Scholar
Crossref
Search ADS
 
31.
M. S.
Green
, “
Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids
,”
J. Chem. Phys.
22
(
3
),
398
–
413
(
1954
).
https://doi.org/10.1063/1.1740082
Google Scholar
Crossref
Search ADS
 
32.
R.
Kubo
, “
Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems
,”
J. Phys. Soc. Jpn.
12
(
6
),
570
–
586
(
1957
).
https://doi.org/10.1143/jpsj.12.570
Google Scholar
Crossref
Search ADS
 
33.
M.
Drosten
 and
M.
Barbacid
, “
Targeting the MAPK pathway in KRAS-driven tumors
,”
Cancer Cell
37
(
4
),
543
–
550
(
2020
).
https://doi.org/10.1016/j.ccell.2020.03.013
Google Scholar
Crossref
Search ADS
PubMed
 
34.
H. I.
Ingólfsson
,
C.
Neale
,
T. S.
Carpenter
,
R.
Shrestha
,
C. A.
López
,
T. H.
Tran
,
T.
Oppelstrup
,
H.
Bhatia
,
L. G.
Stanton
,
X.
Zhang
,
S.
Sundram
,
F.
Di Natale
,
A.
Agarwal
,
G.
Dharuman
,
S. I. L.
Kokkila Schumacher
,
T.
Turbyville
,
G.
Gulten
,
Q. N.
Van
,
D.
Goswami
,
F.
Jean-Francois
,
C.
Agamasu
,
D.
Chen
,
J. J.
Hettige
,
T.
Travers
,
S.
Sarkar
,
M. P.
Surh
,
Y.
Yang
,
A.
Moody
,
S.
Liu
,
B. C.
Van Essen
,
A. F.
Voter
,
A.
Ramanathan
,
N. W.
Hengartner
,
D. K.
Simanshu
,
A. G.
Stephen
,
P.-T.
Bremer
,
S.
Gnanakaran
,
J. N.
Glosli
,
F. C.
Lightstone
,
F.
McCormick
,
D. V.
Nissley
, and
F. H.
Streitz
, “
Machine learning-driven multiscale modeling reveals lipid-dependent dynamics of RAS signaling proteins
,”
Proc. Natl. Acad. Sci. U. S. A.
119
(
1
),
e2113297119
(
2022
).
https://doi.org/10.1073/pnas.2113297119
Google Scholar
Crossref
Search ADS
PubMed
 
35.
H. I.
Ingólfsson
,
H.
Bhatia
,
F.
Aydin
,
T.
Oppelstrup
,
C. A.
López
,
L. G.
Stanton
,
T. S.
Carpenter
,
S.
Wong
,
F.
Di Natale
,
X.
Zhang
 et al, “
Machine learning-driven multiscale modeling: Bridging the scales with a next-generation simulation infrastructure
,”
J. Chem. Theory Comput.
19
(
9
),
2658
–
2675
(
2023
).
https://doi.org/10.1021/acs.jctc.2c01018
Google Scholar
Crossref
Search ADS
 
36.
R.
Shrestha
,
T. S.
Carpenter
,
Q. N.
Van
,
C.
Agamasu
,
M.
Tonelli
,
F.
Aydin
,
D.
Chen
,
G.
Gulten
,
J. N.
Glosli
,
C. A.
López
,
T.
Oppelstrup
,
C.
Neale
,
S.
Gnanakaran
,
W. K.
Gillette
,
H. I.
Ingólfsson
,
F. C.
Lightstone
,
A. G.
Stephen
,
F. H.
Streitz
,
D. V.
Nissley
, and
T. J.
Turbyville
, “
Membrane lipids drive formation of KRAS4b-RAF1 RBDCRD nanoclusters on the membrane
,”
Commun. Biol.
7
(
1
),
242
(
2024
).
https://doi.org/10.1038/s42003-024-05916-0
Google Scholar
Crossref
Search ADS
PubMed
 
37.
K.
Nguyen
,
C. A.
López
,
C.
Neale
,
Q. N.
Van
,
T. S.
Carpenter
,
F.
Di Natale
,
T.
Travers
,
T. H.
Tran
,
A. H.
Chan
,
H.
Bhatia
 et al, “
Exploring CRD mobility during RAS/RAF engagement at the membrane
,”
Biophys. J.
121
(
19
),
3630
–
3650
(
2022
).
https://doi.org/10.1016/j.bpj.2022.06.035
Google Scholar
Crossref
Search ADS
PubMed
 
38.
H. I.
Ingólfsson
,
H.
Bhatia
,
T.
Zeppelin
,
W. F. D.
Bennett
,
K. A.
Carpenter
,
P.-C.
Hsu
,
G.
Dharuman
,
P.-T.
Bremer
,
B.
Schiøtt
,
F. C.
Lightstone
, and
T. S.
Carpenter
, “
Capturing biologically complex tissue-specific membranes at different levels of compositional complexity
,”
J. Phys. Chem. B
124
(
36
),
7819
–
7829
(
2020
).
https://doi.org/10.1021/acs.jpcb.0c03368
Google Scholar
Crossref
Search ADS
 
39.
A. M.
Ring
,
A.
Manglik
,
A. C.
Kruse
,
M. D.
Enos
,
W. I.
Weis
,
K. C.
Garcia
, and
B. K.
Kobilka
, “
Adrenaline-activated structure of β2-adrenoceptor stabilized by an engineered nanobody
,”
Nature
502
(
7472
),
575
–
579
(
2013
).
https://doi.org/10.1038/nature12572
Google Scholar
Crossref
Search ADS
PubMed
 
40.
V.
Cherezov
,
D. M.
Rosenbaum
,
M. A.
Hanson
,
S. G. F.
Rasmussen
,
F. S.
Thian
,
T. S.
Kobilka
,
H.-J.
Choi
,
P.
Kuhn
,
W. I.
Weis
,
B. K.
Kobilka
, and
R. C.
Stevens
, “
High-resolution crystal structure of an engineered human β2-adrenergic G protein–coupled receptor
,”
Science
318
(
5854
),
1258
–
1265
(
2007
).
https://doi.org/10.1126/science.1150577
Google Scholar
Crossref
Search ADS
PubMed
 
41.
T.
Nur Ozturk
,
M.
König
,
T. S.
Carpenter
,
K. B.
Pedersen
,
T. A.
Wassenaar
,
H. I.
Ingólfsson
, and
S. J.
Marrink
, “
Building complex membranes with Martini 3
,” in
Methods in Enzymology
(
Elsevier
,
2024
).
Google Scholar
Crossref
Search ADS
 
42.
F.
Di Natale
,
H.
Bhatia
,
T. S.
Carpenter
,
C.
Neale
,
S. K.
Schumacher
,
T.
Oppelstrup
,
L.
Stanton
,
X.
Zhang
,
S.
Sundram
,
T. R. W.
Scogland
,
G.
Dharuman
,
M. P.
Surh
,
Y.
Yang
,
C.
Misale
,
L.
Schneidenbach
,
C.
Costa
,
C.
Kim
,
B.
D’Amora
,
S.
Gnanakaran
,
D. V.
Nissley
,
F.
Streitz
,
F. C.
Lightstone
,
P.-T.
Bremer
,
J. N.
Glosli
, and
H. I.
Ingólfsson
, “
A massively parallel infrastructure for adaptive multiscale simulations: Modeling RAS initiation pathway for cancer
,” in
The International Conference for High Performance Computing, Networking, Storage and Analysis
(
ACM
,
2019
), p.
57
.
Google Scholar
Crossref
Search ADS
 
43.
M.
Manna
,
M.
Niemelä
,
J.
Tynkkynen
,
M.
Javanainen
,
W.
Kulig
,
D. J.
Müller
,
T.
Rog
, and
I.
Vattulainen
, “
Mechanism of allosteric regulation of β2-adrenergic receptor by cholesterol
,”
eLife
5
,
e18432
(
2016
).
https://doi.org/10.7554/elife.18432
Google Scholar
Crossref
Search ADS
PubMed
 
44.
Y.
Wu
,
L.
Zeng
, and
S.
Zhao
, “
Ligands of adrenergic receptors: A structural point of view
,”
Biomolecules
11
(
7
),
936
(
2021
).
https://doi.org/10.3390/biom11070936
Google Scholar
Crossref
Search ADS
PubMed
 
45.
S.
Gahbauer
 and
R. A.
Böckmann
, “
Membrane-mediated oligomerization of G protein coupled receptors and its implications for GPCR function
,”
Front. Physiol.
7
,
226635
(
2016
).
https://doi.org/10.3389/fphys.2016.00494
Google Scholar
Crossref
Search ADS
 
46.
J. M.
Johnston
,
H.
Wang
,
D.
Provasi
, and
M.
Filizola
, “
Assessing the relative stability of dimer interfaces in G protein-coupled receptors
,”
PLoS Comput. Biol.
8
(
8
),
e1002649
(
2012
).
https://doi.org/10.1371/journal.pcbi.1002649
Google Scholar
Crossref
Search ADS
PubMed
 
47.
X.
Periole
,
A. M.
Knepp
,
T. P.
Sakmar
,
S. J.
Marrink
, and
T.
Huber
, “
Structural determinants of the supramolecular organization of G protein-coupled receptors in bilayers
,”
J. Am. Chem. Soc.
134
(
26
),
10959
–
10965
(
2012
).
https://doi.org/10.1021/ja303286e
Google Scholar
Crossref
Search ADS
 
48.
T. N.
Ozturk
,
N.
Bernhardt
,
N.
Schwartz
,
R.
Chadda
,
J. L.
Robertson
, and
J. D.
Faraldo-Gómez
, “
Mitigation of membrane morphology defects explain stability and orientational specificity of CLC dimers
,” bioRxiv:10.1101/2023.03.16.533024 (
2023
).
Google Scholar
 
© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
This work is licensed under a Creative Commons Attribution 4.0 International License .
Attribution