Membranes form the boundaries of all living cells. Besides structural integrity, they provide a host of additional functions, chief among them the controlled exchange of material and information between a cell and its environment. Combining these tasks requires an interface that is both thermodynamically and mechanically stable while also being deformable on many length scales with biochemically available energies. Self-assembled lipid bilayers have filled this role from the earliest point of the phylogenetic tree we have any information on. 1 One might think that a minimalist lipid bilayer made from a single lipid species would suffice for this role, but this is not what biology has chosen to do. Even ignoring all the embedded proteins or attached filaments and carbohydrates, a membrane's lipidome is extraordinarily complex. For reasons still not well understood, a cell can easily contain more than a thousand different lipid types, 2 which are distributed unevenly between organelles or the two leaflets of any membrane and maybe even be laterally inhomogeneous within a single leaflet. In the late 1980s, first (and remarkably indirect) hints of small submicroscopic domains, which came to be known as "lipid rafts," sparked a lot of interest in lipid mixtures and their thermodynamics, [3][4][5][6][7] and while many physical and functional aspects attributed to rafts continue to be vigorously debated, [8][9][10] the main idea that lipids are often very non-ideally mixed seems to be widely accepted these days.
Of course, mixtures of hundreds of lipids cannot be studied systematically, and so idealized model systems that capture key aspects of the underlying biophysical situations are needed. Of those, an extremely rewarding type turned out to be ternary mixtures of a high melting lipid [such as dipalmitoylphosphatidylcholine (DPPC) or sphingomyelin], a low melting lipid [such as dioleoylphosphatidylcholine (DOPC)], and cholesterol. These mixtures are intriguing because, in a wide region of compositional phase space, they exhibit phase coexistence of two different fluid phases-a liquid ordered ("lo") and a liquid disordered ("l d ") one. [11][12][13][14][15][16] In particular, upon increasing the cholesterol content,
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp this coexistence region narrows and ultimately terminates in a critical point (of 2d-Ising universality class), giving rise to all the fascinating non-ideal mixing thermodynamics that comes with (near-) criticality. 17,18 It was subsequently shown that the same Ising-like criticality occurs naturally in certain biomembranes, 19 that it can be shifted by a wide variety of membrane-partitioning small molecules (such as anesthetics), 20 and that cells may adjust their lipidome to set the difference between this critical point and their growth temperature to some given value. 21 Taken together, these results show that cells deeply care about the non-ideal mixing thermodynamics of their lipidomes, but that simple ternary lol d systems could be an excellent model system to learn why that might be so. In an independent line of inquiry, it was discovered in the early 1970s that the two leaflets of the plasma membrane differ in lipid composition, [22][23][24][25][26] an observation that modern research has backed up and greatly refined. 27 Combined with the nontrivial mixing thermodynamics already mentioned, the question of which leaflet cholesterol would preferentially occupy, how phase separated domains (or transient fluctuations) in two opposing leaflets would interact, and how any of this would depend on the state of membrane curvature, we recognize a formidable physical problem at hand that intricately couples mechanics, thermodynamics, and geometry.
Computational models aiming to address this situation need to simultaneously represent an unusually wide spectrum of physical phenomena, but they must also be efficient so that membrane patches larger than near-critical correlation lengths (say, dozens of nanometers or many thousands of lipids) can be treated. A number of lattice based models have successfully captured non-ideal mixing behavior, [28][29][30][31][32][33][34][35] but they typically do not contain elasticity. Conversely, continuum models designed to capture curvature elasticity in some discretized way [36][37][38][39][40][41][42][43][44][45][46][47] tend to not also include the compositional degree of freedom, even though some notable exceptions that couple to lipid composition [48][49][50] or embedded proteins [51][52][53][54] exist. In addition, particle based models that properly deal with the subtleties of lol d phases, even if coarse grained (such as Martini [55][56][57][58][59][60][61] ), struggle to access the scales and equilibration times needed for such questions.
As always, no model is perfect, and so usefulness is a matter of balancing needs with constraints in view of the essential set of questions one wishes to tackle. Here, our aim is to expand a well-studied and highly coarse-grained particle-based (i.e., off-lattice) implicitsolvent lipid model so that it can capture the essential aspects of the lol d phase coexistence of lipid membranes. To this end, we tune lipid properties and lipid-lipid interactions in the so-called Cooke model, 62,63 specifically its recent extension that can suppress flip-flop and hence stabilize asymmetric membranes. 64 Among the key phenomena we wish to capture are fluid-fluid phase coexistence with a suitable order and lipid area difference, the characteristic dependence on cholesterol, a critical point, having phospholipids stick to their leaflets but permit cholesterol to flip, while maintaining the proper bending and stretching elastic behavior that had been previously well established. However, our objective is not to precisely replicate any particular experimentally studied ternary mixture. Although it is possible to design various kinds of lipids by adjusting a host of lipid properties, our goal is to examine the generic aspects of the physics of mixing an ordered lipid, a disordered lipid, and a third sterol-like component in as simple a way as possible. We deliberately overlook notable differences among these systems, such as the presence of a significant coexistence region between saturated lipids and cholesterol, resulting in a three-phase triangle in the ternary phase diagram.
The benefit of this minimalist level of description is the ability to explore larger systems for longer times. To illustrate this, consider that the coarse grained (and hence sped-up) time unit for the popular Martini model is about 64 ps, if mapped by matching lipid diffusion. 65 In contrast, for the Cooke model we will be working with, that time unit is 1.9 ns. 64 Consequently, the Cooke model traverses phase space ∼30 times faster than Martini, even before considering the further decrease in computational effort due to the significantly lower number of beads in a lipid and the fact that the Cooke model does not explicitly represent solvent particles. Overall, this makes our model around 100 times faster than Martini, allowing us to explore longer timescales and larger system sizes-at the expense of a lower resolution.
Our paper is organized as follows: after introducing biophysical observables of interest and how to measure them, we discuss the step-by-step expansion of the model-from single-lipid properties to binary mixtures and then to ternary ones. Following this, we discuss the resulting phase diagram and examine a number of relevant observables for the resulting mixtures, especially as a function of cholesterol mole fraction.
The Cooke model 62,63 is an implicit solvent coarse-grained framework for simulating lipid bilayers, typically via molecular dynamics. In its latest iteration, lipids are represented as linear sequences of four distinct beads (Fig. 1). 64 Cooke lipids possess various adjustable characteristics, including individual bead dimensions, lipid tail stiffness, and spontaneous curvature. The Lorentz-Berthelot mixing rules are used to determine effective repulsive pair distances σij from the bead diameters σii and σjj of beads of type i and j, respectively, The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
All lengths are measured in the coarse-grained length unit σ, which is close to a (phospho-)lipid bead size and which maps to ∼0.75 nm in SI units. 64 Cooke lipids are straightened by a harmonic spring with rest length 4σ between the two pairs of beads that are separated by a single third bead,
To the lowest order in angular deflection, this is equivalent to an angular bending potential with spring constant k bend σ 2 . The constant k bend is set to 10σ 2 ( is the coarse-grained unit of energy) for regular Cooke lipids, but here we will adjust it as a straightforward method for tuning tail flexibility and hence the emergent lipid order. We performed all simulations at a temperature k B Ts = 1.4, from which we can infer the coarse-grained (CG) unit of energy . Interaction strengths, especially the attraction between lipid tails of the same or of different types, are tuned using the model parameter wc. This parameter sets the range of attraction, and the larger it is, the more strongly lipids cohere. Tuning wc and the simulation temperature Ts hand-in-hand is the key to obtaining a fluid bilayer in this model. 62,63 After carefully examining the wc vs Ts phase diagram, the values wc = 1.6σ and k B Ts = 1.4 were chosen because the bilayer remains fluid for a broad region surrounding this state point. To distinguish between lipids of different types, we adjust the interaction strength wc,ij between the second and third beads of lipid types i and j, whereas the fourth (or final tail) bead always has wc = 1.6σ with respect to all other tail beads. The subscripts i, j ∈ {s, u, c}, where "s," "u," and "c" stand for saturated lipids, unsaturated lipids, and cholesterol, respectively.
The Cooke model has been successfully applied to a variety of questions in membrane biophysics, such as protein triggered vesiculation, 66,67 composition-curvature coupling, 68 peptide induced pore formation, 69 lipid nanoparticle design, 70 and several other examples. More closely related to the present topic, we have recently used this model to examine binary mixtures of cholesterol and generic phospholipids in a regime far from phase separation, 71 with the underlying question being: what physical factors determine the distribution of a flippable species between leaflets? The present paper aims to extend the scope of this work by enabling us to bring lol d coexistence into the picture, which requires us to examine the more complex phase behavior of ternary mixtures. The main tuning parameters at our disposal are lipid rigidity k bend (as a means to influence its tail order) and mutual interaction strengths wc,ij (which affect solubility preferences). We also slightly reduce the size of Cooke-lipids that mimic cholesterol to account for the smaller membrane footprint of this species. We provide further details on these modifications in Secs. III A and III B.
All simulations were performed with the ESPResSo package, 72 using a time step δt = 0.005τ, where τ = σ m is the CG unit of time. The system was thermalized by a standard Langevin thermostat 73 with a friction constant γ = 1mτ. Tensionless membranes were realized through constant pressure simulations on the basis of a Kolb/Dünweg barostat 74 with a piston mass Q = 0.01mσ 4 and a friction constant of γ V = 0.0002mσ 4 τ. All trajectories used in our study were between 80 000τ and 100 000τ long, with the first 20 000τ being discarded to allow the initial state to thermalize.
We will use a range of observables to characterize specific aspects of the lipid membrane, its constituent phases, and individual lipid species. To distinguish ordered from disordered phases, we resort to area per lipid, the hexatic order parameter, and the P 2 order parameter. To probe mechanics, we will measure the area expansion modulus KA, both for pure and for mixed systems. In addition, to probe the dynamic differences between coexisting phases, we measure lipid diffusion constants. In this section, we offer some details on how these observables are computed. Given that all measurements ultimately derive from simulation trajectories, which harbor inevitable time correlations, we use a standard blocking analysis 75 combined with a fit to an analytical prediction for the blocking plateau 76 for estimating the uncertainties. When a different error estimation method is applied, it has been described separately under the relevant subsection.
The area per lipid is calculated by dividing a membrane's relaxed surface area by the number of lipids per leaflet, which we always chose to be either 512 or 1024. This requires holding a membrane at zero tension, which we achieve via a semi-anisotropic NPT barostat. 74
The P 2 order parameter quantifies how persistent lipids align with the bilayer normal. Given a lipid's deflection from the bilayer normal (here: the head-tail director relative to the ẑ-axis), we calculate P 2 = P 2 (cos i), where P 2 (x) is the second Legendre polynomial, and where we average over both time and lipids. P 2 is bound between -1 2 and 1, with a value close to 1 indicating a high degree of alignment with the bilayer normal, whereas P 2 = 0 implies a completely isotropic distribution. (Negative P 2 values arise for preferential alignments away from the axis.)
The hexatic order parameter measures the degree of hexagonal packing in two-dimensional molecular structures, such as lipid bilayers. To compute it, the extent to which the neighbors of any particle display a local 6-fold rotational symmetry is analyzed. More specifically, for each lipid k, one calculates
where the sum extends over all N k neighbors of lipid k and where θ jk is the angle which the line from lipid k to one of its neighboring lipids j makes with respect to some arbitrarily chosen reference direction (say, the x axis). The magnitude of this complex number quantifies how well lipid packing creates a local 6-fold orientational symmetry. A value close to 1 indicates a high degree of hexagonal packing, while a value closer to 0 signifies complete disorder. We subsequently thermally average this magnitude, and if we care about the overall hexatic order of the membrane (or a subphase of it), we then also average over lipids.
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
The area expansion modulus, denoted as KA, quantifies the resistance of a lipid bilayer to changes in its surface area under isotropic biaxial strain. Specifically, it is defined as the constant of proportionality between area stress ΣA and area strain,
where A represents the membrane surface area under stress and A 0 denotes the relaxed surface area. Given the linear nature of the relationship, it is possible to determine KA by either applying controlled stress and measuring the resulting strain or, conversely, by imposing area strain and measuring the resulting stress. We have chosen the latter method because (1) it is easier to impose strain and (2) we can identify departures from linearity without risking instabilities. For added reliability, the bilayer rest area A 0 for each configuration was determined through constant (zero) tension simulations. We applied strains ranging from 1% to 5% for each composition and measured the resulting stress. The simulations ran for 100 000τ, with equilibrium being reached once the cholesterol flip-flop showed no net drift, typically after ∼20 000τ. The lateral stress profile was measured using the Irving-Kirkwood formalism, 77 extensively detailed by Hardy 78 for implementation in molecular dynamics simulations. We used blocking 75 to estimate the error bars on the total membrane tension. A linear fit, as described earlier, was performed on the stress-strain data, and the associated error bars for KA were derived from this linear fit.
The lipid diffusion constant D is a widely used observable for characterizing a membrane's in-plane dynamics. We specifically want to determine its value in the liquid-ordered (lo) and liquiddisordered (l d ) phases, aiming to unveil their dynamic differences and whether they qualitatively follow experimentally known trends.
The diffusion constant is commonly measured from the slope of the mean squared displacement (MSD) as a function of time lags. Unfortunately, estimating the associated uncertainty is much less straightforward since subsequent MSDs along a lipid's trajectory are highly correlated, and so the naïve fitting error of the slope vastly underestimates the actual uncertainty in D.
In response to this challenge, alternative strategies have been proposed. [79][80][81][82] These include optimizing the number of MSD values for fitting or explicitly incorporating these correlations into the fitting process. A notably more efficient procedure is to model the relationship between MSD and time lags using a stochastic linear model, which yields the diffusion constant (and its uncertainty) as one of the parameters in such an analysis.
Here, we employ one such method proposed by Bullerjahn et al. 82 to extract diffusion constants from molecular dynamics trajectories. It parameterizes the mean squared displacement of a finite time series of positions {X 0 , X 1 , . . . , X N-1 , XN}, after a time lag it using the following linear model:
where t denotes the length of the time interval between two consecutive observations and i = 1, 2, . . . , M ≤ N. Here, s 2 is the diffusional spread and a 2 is a static-noise parameter that accounts for any nondiffusional spread in the particle position, such as cage diffusion and correlated collisions. The estimators for a 2 and s 2 , along with their uncertainties, are found using the Generalized Least Squares (GLS) method given by a 2 GLS , s 2 GLS = arg min
where Resi = MSDia 2is 2 is the residual at step i and Γ is the covariance matrix of the MSDi values. It should be noted that the covariance matrix is evaluated at the optimal estimates a 2 GLS , s 2 GLS while a 2 and s 2 are varied. Then, by enforcing self-consistency, a 2 GLS and s 2 GLS (which are previously unknown) are found. The GLS technique uses the covariance matrix (calculated numerically from the discrete data points) to "undo" the time correlation among various MSDi values. This yields the best linear unbiased estimator for a 2 and s 2 , subsequently giving us an estimate of the diffusion constant. In addition to this, lower bounds for the variances of the GLS estimators can be inferred from the inverse of the Fisher information matrix associated with the likelihood function,
For d-dimensional motion, the associated self-diffusion constant D and its uncertainty δD can now be calculated from the estimated diffusional spread s 2 GLS and its variance var(s 2 GLS ) as follows:
The details of constructing this estimator from MD trajectories along with the underlying mathematical concepts and calculations are discussed in Bullerjahn et al. 82
A hidden Markov model (HMM) is a statistical tool used to analyze sequential data believed to arise from an underlying Markov process. The HMM assumes that the observed data, often represented as a sequence of symbols, are functions of the Markov states. The states themselves are unobservable ("hidden"), but they have a probabilistic relationship with the observed data. The HMM comprises two main parts: the state transition probabilities, which define how likely it is to move from one hidden state to another, and the emission probabilities, which specify the likelihood of generating each observed symbol from each hidden state.
Inspired by previous work in the field, 83,84 we used an HMM that was trained on data from a well-separated ternary mixture to
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp assess the degree of phase separation and understand the composition of the two phases in membranes of varying compositions. As training features for the model, we used the specific area of each lipid, distance between a lipid and its nearest neighbor, and the most frequently occurring lipid type in each lipid's immediate vicinity.
The resulting model categorized the lipids into two classes, representing the liquid-ordered and liquid-disordered phases. In Fig. 2, we illustrate the outcome of such an analysis by contrasting lipid composition and phase composition in two extreme regimes-the near-critical/high cholesterol case and a well phase-separated/low cholesterol scenario. This HMM analysis serves as a robust and dependable method for quantifying phase coexistence within a membrane. Well within the lol d coexistence region, the HMM finds clearly circumscribed phases that are extensive, stable, and persistent (see Fig. 2, bottom row). In contrast, near the critical point, the identified phases are small, scattered, and transient (see Fig. 2, top row). Of course, an HMM analysis does not obviate complications due to finite size effects, and the closer one gets to a critical point, the larger the system one needs to examine.
Beyond identifying phase coexistence, this approach also permits us to quantify the composition of the two identified phases. From this we can construct a comprehensive phase diagram, outlining in particular the binodals and the tie lines of the coexistence region. The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp will ultimately translate into lipid order, we capture this aspect via the spring constant k bend from Eq. ( 2), which suppresses lipid shape fluctuations by stretching the four beads into a straight alignment.
To examine the magnitude of this effect and its impact on phase behavior, we conducted simulations of one-component membranes consisting of lipids for which we tuned k bend within the range of 1σ 2 to 40σ 2 . Figure 3 shows the extent to which increasing k bend increases lipid hexagonal packing order [as measured by ψ 6 , see Eq. ( 3)] and decreases the area per lipid. Our data show a phase transition, indicating the emergence of an ordered gel phase for values of k bend somewhat larger than 20σ 2 . Experimentally, the area per lipid has been shown to be significantly smaller in gel phases than in fluid phases, for instance by about 20% in dimyristoylphosphatidylcholine (DMPC) 85,86 and dilauroylphosphatidylethanolamine (DLPE) 87 and by almost 30% in dihexadecylphosphatidylcholine (DHPC). 88 Differences of this magnitude are readily accessible in our model.
Finally, we also wish to represent cholesterol as one of the lipid species. While its chemical structure differs substantially from a phospholipid, our coarse-grained model does not have enough resolution to capture this. Instead, we will model an "effective" cholesterol molecule as one that has again the same tail-tail attraction (wc,cc = 1.6σ), as well as a tail bending constant k bend = kc = 10σ 2 . This is equal to the standard value used in Ref. 64 and lies between the saturated and unsaturated cases. Furthermore, we reduce the size of all cholesterol beads by 20%, thus giving CG cholesterol a 36% smaller area per lipid compared to phospholipids, in accordance with cholesterol and phospholipid area simulations performed by Leeb and Maibaum. 89 Finally, since the flip-flop rate of sterols vastly exceeds that of ordinary lipids, we allow our CG cholesterol to flip-flop freely between the leaflets without any energy penalization (which happens at a rate of about 10 -5 τ -1 ), unlike the saturated and unsaturated lipids, whose flipping we suppress using the "flip-fix" method discussed in Ref. 64.
We now have a (fairly simple) lipidome available from which we can pick components for subsequent mixing. Recall that our main target is ternary systems that exhibit lol d coexistence. These consist of a high-melting lipid (saturated, "s"), a low-melting lipid (unsaturated, "u"), and cholesterol ("c"). Our first task is hence to pick lipids representing "s," "u," and "c" and then examine how they mix-starting with the binary cases.
Let us begin with the "u"-"s" pair. A common choice in experiments for these is DOPC for the unsaturated lipid and DPPC for the saturated one. At either room or body temperature, DOPC is safely in the fluid phase, while DPPC is in the gel phase. In view of Fig. 3, we should hence pick the bending constants k bend for both lipids to the left and right of the transition that happens around k bend ≈ 20σ 2 . Since gel membranes equilibrate poorly and, upon further cooling, may develop additional phases with intricate order parameters, we do not wish to move too deeply into the gel region, and so we chose ks = 25σ 2 for the chain bending constant of the model saturated lipid. To then permit the corresponding disordered fluid-phase lipid to have an area per lipid that is between 20% and 30% larger than the gel area, we chose ku = 5σ 2 . Smaller values of k bend would create an even larger area contrast, but this was not needed. Furthermore, membranes begin to look too disordered when k bend drops below about 5σ 2 which, especially for such highly coarse-grained models, might risk bilayer integrity. All lipid stiffnesses are summarized on the left hand side of Fig. 4.
To fully specify the mixture, we also need to set an interaction constant wc,su between tail beads. Our initial hope that at the standard value of 1.6σ the difference in lipid rigidity alone FIG. 4. Parameters used to create ternary mixtures. The left column gives k bend for the three lipid types u, s, and c in units of σ 2 . The right matrix gives the interaction parameters wc for pairs of lipids in units of σ; for instance, wc,us = 1.4σ.
pubs.aip.org/aip/jcp would suffice to drive phase separation was not borne out. Highly coarse-grained lipids simply do not possess enough conformational tail entropy and, therefore, demixing requires an additional enthalpic assist. This is performed by reducing the tail-tail interaction parameter wc,su compared to its value for like-lipid interactions (wc,uu = wc,ss = 1.6σ). We found that demixing sets in near wc,su ≈ 1.45σ, but is initially too weak to create a wide enough coexistence region. A subsequent decrease widens the binodal, and we chose wc,su = 1.4σ. When mixing CG lipids and CG cholesterol, even at the same tail-tail interaction constant wc = 1.6σ, we observed a pronounced tendency for cholesterol to form aggregates. This is a direct consequence of its reduced size since the effect does not occur for CG cholesterols having the same size as normal lipid molecules. Since we wish cholesterol to fully mix with unsaturated lipids, we compensated by picking a tail-tail mixing parameter that was slightly larger than the standard one, wc,uc = 1.68σ. For the corresponding interaction with the saturated lipid, we used a larger mixing bias, wc,sc = 1.75σ, to account for the fact that cholesterol prefers saturated over unsaturated phases.
Since picking k bend,sat = ks = 25σ 2 places the saturated lipids into the gel phase, while a mixture of "s" and (sufficiently much) "c" will result in an lo phase, we expect a coexistence region across the sc-mixing axis. This is indeed what the simulations indicate, but it appears that the onset of coexistence occurs at such a small cholesterol fraction that we cannot reliably resolve it. (Identifying these structurally fairly close phases with an HMM also turned out to be very challenging.) As this is not a point in our phase diagram that we currently wish to characterize very precisely, we have forgone attempts to pin down the details of this coexistence. We emphasize, though, that this does not occur in binary mixtures of cholesterol with the unsaturated lipids, which always mix into a homogeneous fluid phase-as indeed they should.
All wc parameters are summarized in the interaction matrix on the right hand side of Fig. 4.
As an illustration of the resulting binary mixing physics, Fig. 5 shows the P 2 order parameter as a function of cholesterol concentration in binary "uc"-and "sc"-mixtures. The effects are strikingly FIG. 5. P 2 order parameters of the phospholipids, when mixed with cholesterol, vs cholesterol mole fraction ϕ. The error bars are of the order 10 -5 and, therefore, invisible on this scale. different: very little cholesterol is needed to melt the gel phase of pure "s"-lipids, leading to a strong reduction of P 2 . After a short coexistence region, the system enters a one-phase liquid ordered region in which P 2 stays mostly constant. In contrast, pure "u"-lipids are already in a fluid phase and have a lower value of P 2 . The addition of cholesterol then orders the initially disordered phase, but it takes fairly large amounts of cholesterol to induce relatively modest ordering.
We conducted a rough mapping of the phase diagram, shown in Fig. 6, with a special emphasis on the extent of the miscibility gap. Our ternary mixtures exhibited a clear separation into two distinct phases: one characterized by a higher degree of order, denser packing, and a larger abundance of "s"-type lipids; the other being more disordered, less tightly packed, and richer in "u"-type lipids. To identify these phases and determine their lipid composition, we employed the HMM analysis described in Sec. II D.
As we traverse the cholesterol axis through the two-phase region, a consistent trend emerges: with increasing cholesterol mole fraction, the tie lines become shorter, but their slope increases. This indicates that as the compositional difference between the ordered and disordered phases in terms of phospholipid mole fraction decreases, their difference in terms of cholesterol content increases (until, of course, toward the critical point, they all become identical). The overall weak tilt of the tie lines (especially at low cholesterol content), i.e., the relatively small enhancement of cholesterol in the ordered phase, is in line with experimental observations. 12,13,16 Another notable feature of the coexisting phases is that they exhibit anti-registration between the leaflets. To examine this phenomenon further, we divided the membrane into grid cells, as illustrated in Fig. 7, and checked the degree of anti-registration using the relative asymmetries of saturated and unsaturated lipids, ρ s and ρ u , defined (in each cell) as
Here, s± or u± stand for the number of saturated or unsaturated lipids within each grid cell in the + orleaflet. These ratios take the value 1 if all lipids of that respective type are in the upper leaflet and -1 if they are all in the lower leaflet; they are zero at equidistribution.
In domain-forming simulations that contain equal amounts of saturated and unsaturated lipids, we observed a pronounced negative correlation between ρ u and ρ s , with a Pearson coefficient smaller than -0.8, often approaching -1. Regions in one leaflet that are rich in saturated lipids face regions in the other leaflet that are rich in unsaturated lipids. This shows that our lol d domains are anti-registered.
Trans-bilayer domain coupling has been previously studied in experiments and simulations. While both registration and antiregistration have been observed, the registration appears to be the more common scenario-prominently visible in early lol d observations made in giant unilamellar vesicles (GUV). 90 Enoki et al. have confirmed these findings but discovered that in asymmetric GUVs, antiregistration is also possible. 91 Langmuir-Blodgett/Langmuir-Schaefer (LB/LS) techniques, but only if the membranes are separated sufficiently far from the substrate by a flexible polymer cushion. 92 Wagner et al. have proposed a simple model for phase registration that only involves a single order parameter, which only leads to registration, provided the coupling is strong enough. 93 The more refined theory by Williamson and Olmsted permits both registered and antiregistered phases. 94 Sharma et al. examine this in a simplified simulation model and show that inter-leaflet coupling energies can favor either of the two outcomes. 95 Weiner and Feigenson show in CG simulations using the Martini force field that upon increasing the cholesterol fraction to 35 mol. %, registration can give way to antiregistration, but also a sizable fraction of cholesterol molecules residing in the bilayer midplane. 96 Thallmair et al. also observe a strong effect of intercalated cholesterol but instead find it to support registration. 97 All in all, membrane inter-leaflet coupling remains a surprisingly subtle phenomenon, likely because the outcome depends on many factors, such as contributions to lipid(-domain) interactions due to enthalpy (e.g., van der Waals and charge), entropy (lipid fluctuations or whole leaflet undulations), tail interdigitation, hydrophobic mismatch, domain curvature and -tension, and other effects. Of note, registration generally changes the membrane thickness across the interface and hence comes with an additional energy penalty due to hydrophobic mismatch; this is not true for anti-registration. Since for a small enough system such line tension effects can overwhelm the phase contact contribution, which scales with contact area, finite size effects may favor anti-registration.
Our model does account for some of the more subtle domain wall effects, and our simulations might still be too small to evaluate the equilibrium phase. We also wish to emphasize that our number-4 tail beads, which contribute most strongly to inter-leaflet contact, always interact with the same interaction range wc = 1.6σ, irrespective of lipid type. This evidently removes a strong driving force toward domain registration, but the alternative of having the bead-4 interactions also reflect lipid type would strike us as trivially over-representing the enthalpic contribution to a far more subtle overall free energy balance happening at the leaflet interface.
At this point, we do not yet aim for a systematic investigation of the interleaflet domain physics. Our primary focus in this first paper is not a comprehensive mapping of the entire phase space or determining the precise location of the critical point, but rather incorporating the fascinating general physics of ternary lol d phase separation into a highly coarse-grained model that has been well-characterized in terms of its many elastic properties. This puts us in a position to study a large set of questions where these two pieces of physics interact, such as curvature induced sorting and domain formation, domain registration, its modification once membranes become asymmetric, and the effect of differential stress on all of this.
We have analyzed the diffusion process in two regimes of the phase diagram-close to the critical point and well into the phasecoexistence region-with the aim of acquiring a deeper understanding of the dynamic characteristics exhibited by the phases. Our investigation, which employed the analytical techniques outlined by Bullerjahn et al., 82 shows that lipids undergo normal diffusion, with a mean squared displacement (MSD) that is linear in time but a phase-dependent diffusion constant.
As anticipated, lipids in the l d phase have a larger diffusion constant than those in the lo phase. Specifically, in the case of a membrane containing only 10% cholesterol, which segregates into two compositionally very distinct phases, the diffusion constant for lipids in the disordered phase is ∼26% greater than that in the ordered phase. Experimental observations indicate substantially larger factors in the range of 3-20 for this speed-up, contingent upon different compositions and temperature conditions. [98][99][100] It is important to recall, though, that dynamic properties are not the primary target when developing most CG models (including the present one), and absolute numbers tend to be far off. Here we see that the relative difference is also not well captured, even though it is reassuring that the trend is correct.
As illustrated in Fig. 8, our measurement of the diffusion constant reveals that the constituent elements within the phases display dynamic behavior associated with the phases themselves rather than those specific to individual molecules. For instance, a saturated lipid in the disordered phase still diffuses faster than an unsaturated lipid in the ordered phase. This supports the notion of diffusion as a collective process dependent on phase and environment, not a mere lipid property. This has been confirmed in experiment. 101,102 Observing normal diffusion closer to the critical point is conceivably surprising, given that lipids move in and out of lo and l d phases, which themselves come and go out of existence, and which are associated with a characteristic correlation length ξ. One might have expected some cross-over behavior when the MSD is comparable to ξ 2 , or the time is comparable to ξ 2 D, but this is not what we observe. In fact, the analysis presented by Bullerjahn et al. 82 includes a metric that measures any statistically significant deviation from a linear MSD, but we get no indication of anomalous diffusion. Admittedly, the diffusion constants in the two phases are fairly close, even far away from the critical point, and approaching it would further blur any potential difference. If some anomalous component existed, it would be difficult to see in our case.
We ran simulations to examine how the area expansion modulus KA depends on cholesterol concentration and phase state-see
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp FIG. 9. Lateral stretching modulus K A as a function of cholesterol content in the membrane.
The black points correspond to the measurements that were obtained from a mixed phase inside the coexistence region (with an equal mole fraction of phospholipids), while the solid black line represents the best linear fit to these data. The orange and brown points represent compositions at the ends of the tie lines of the corresponding two-phase systems on the ordered and disordered binodals, respectively.
Fig. 9. The solid black points show the effective modulus for mixtures containing equal amounts of saturated and unsaturated lipids as a function of cholesterol concentration ϕ, up to ϕ = 50%. Recall that these compositions lie in the two-phase region (see the phase diagram in Fig. 6). Increasing ϕ is found to make the membranes more rigid; in fact, a simple linear fit to the data yields a remarkably low p-value of 0.007. This trend is plausible, considering that adding cholesterol increases the lipid order parameter and decreases the area per lipid (black points in Fig. 10), thereby rendering cholesterol-rich bilayers overall better packed and hence harder to stretch. It would seem reasonable to expect that the elastic behavior of these two-phase mixtures is a weighted average of the rigidities of the pure coexisting phases-say, a more rigid lo phase balancing a less rigid l d phase. While this appears to be true at elevated cholesterol content, the behavior is qualitatively different at small ϕ, as the triangle and star symbols in Fig. 9 show. While the pure ordered phase is always more rigid than the two-phase mixture, the pure disordered phase stops softening as ϕ decreases, until it is more rigid than the mixture (at ϕ = 20%) and ultimately even more rigid than the ordered phase (at ϕ = 10%).
It is not impossible for two-phase systems to be softer than their coexisting phases. Consider a van der Waals gas in its coexistence region; its bulk modulus is zero and hence smaller than that of either the gas or the liquid phase. This happens because, upon expanding the system, we simply evaporate liquid at no excess free energy cost. In the present case of coexisting lol d domains, it is conceivable that pulling the system melts part of the lo phase into an l d phase, but due to additional constraints on all concentrations the softening should be smaller than in the much simpler van der Waals case. Moreover, the magnitude of this effect ought to be ϕ-dependent, presumably stronger in regions where the phases differ more drastically-as we FIG. 10. The two plots display the area per lipid a (upper panel) and the P 2 order parameter (lower panel) over the same ϕ-range as Fig. 9, with each measurement corresponding to the respective point in Fig. 9. The open symbols show the onephase values extracted instead from the coexisting phases of the mixed systems. indeed find. We tried to test this by measuring the fraction of a two-phase system in the lo phase as a function of applied strain, but unfortunately our HMM classifier is not precise enough to support a reliable decrease upon pulling.
That (at ϕ = 10%) the disordered phase is more rigid than the ordered phase is even harder to comprehend, and we cannot rule out that it is just statistics. In fact, we found that KA values can change significantly after relatively minor adjustments to the compositions, and the HMM identification of the coexisting phase boundaries might again be too noisy to translate to a dependable trend on KA. This difficulty is further illustrated by the individual pure phase lipid areas and P 2 order parameters shown in Fig. 10: while they show the expected trend with cholesterol concentration ϕ, it would seem impossible that the lipid area of the mixture is larger than either of the two phases (see ϕ = 0.3 and 0.4) or that the mixture's P 2 value is smaller than that of either of the coexisting phases (see ϕ = 0.1), given that in both cases the values for the mixture would literally be calculated as averages over the two-phase simulation. These peculiarities arise because the pure phase values shown in Fig. 10 (and Fig. 9) are not derived from the compositions coexisting in a twophase simulation; instead, they come from two separate simulations of the two one-phase systems at the end of the tie lines, whose compositions are identified by the HMM analysis. However, the The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp HMM classification runs into obvious complications arising from lipids at the interface, and since order parameters depend sensitively on the composition near the transition, it is conceivable that the pure phase simulations do not capture the actually coexisting bulk composition of the two-phase simulations precisely enough. Indeed, re-calculating area and P 2 from the two-phase system (see the corresponding open symbols in Fig. 10) shows neither of the two peculiarities. This difficult-to-quantify systematic uncertainty hence appears to overwhelm the statistical errors (which are too small to be seen in Fig. 10) and, therefore, urges caution when interpreting these data. All the same, our analysis suggests that the elastic properties of mixed systems change significantly as we traverse the phase diagram, with additional subtleties arising once we enter the coexistence region. It is easy to imagine scenarios where this has important biological consequences, especially as we recall that the two leaflets are likely differentially stressed. Clearly, more work is needed to explore this in more detail.
The Cooke model is a highly coarse-grained implicit solvent lipid model that has been thoroughly characterized in terms of its elastic properties and which has recently been extended to permit the simulation of asymmetric membranes. In this study, we have refined its force field in two main ways. First, we introduced a variable lipid tail stiffness as a way to represent the notion of saturated and unsaturated lipids. Tuning it, we can change tail-or packing order, as well as area per lipid, over an experimentally relevant range and transition into a gel phase at sufficiently high order. Second, we created a minimalist version of a cholesterol-like molecule, which is smaller than regular lipids and flip-competent. By adjusting pair-interactions between these species, we can investigate binary and, more importantly, ternary mixtures, with a particular focus on "raft-like" compositions.
At the temperature we simulate, unsaturated lipids form a disordered fluid phase, while the saturated lipids are in a gel phase. Binary mixtures between these two ("phospho"-)lipids exhibit a wide miscibility gap. Adding cholesterol to unsaturated lipids weakly increases their fluid phase order, while very small amounts of cholesterol suffice to melt the gel phase of saturated lipids into an lo phase. Ternary systems of these lipids phase separate into an lo and an l d phase across a two-phase region whose width narrows as cholesterol content increases. The two phases become identical at a critical point somewhere between 50% and 60% cholesterol mole fraction, after which only one phase exists. The lol d tie lines separate phases that-especially at low cholesterol content-differ strongly in phospholipid make-up but contain fairly similar cholesterol fractions, while further addition of cholesterol reduces the phospholipid difference but increases the relative cholesterol content in the lo phase. All of this agrees qualitatively with published phase diagrams on such ternary systems, except that we are unable to identify a possible three-phase triangle in which a gel phase coexists with an lo and an l d phase, likely because the miscibility gap in binary mixtures of cholesterol and saturated lipids appears to be very narrow (in our case) and our limited system sizes and phase identification routines cannot discriminate the small difference between a gel and an lo phase.
We examined several basic thermal equilibrium, elastic, and dynamic properties of mixed systems to begin characterizing their physical behavior. In these first steps, we limited the number of variables and focused on mixtures with equal mole fractions of saturated and unsaturated lipids. At low cholesterol content, these phases separate into compositionally quite distinct but approximately equally sized lo and l d domains. In symmetric membranes, these domains anti-register across the two leaflets, which is unusual given that experimentally registration is more commonly observed. However, a variety of different drivers contribute to the question of domain (anti-)registration: phase contact across the bilayer midplane (which involves enthalpic and entropic contributions, as well as the possibility of tail interdigitation-all of which scales with domain area), hydrophobic mismatch at the domain border (which scales with its length and might induce lipid tilt), torques at the interface due to spontaneous curvature gradients (which also scales with border length but can induce membrane curvature), and even effects not contained in our model (such as a net trans-membrane electric field created by asymmetric dipole potentials arising in anti-registered domains). Given the different scalings with system size and the still moderate membrane dimensions we have treated (up to around 40 nm in CG-to-SI mapped units 64 ), it is not yet clear whether larger domains will also anti-register. Moreover, for unequal amounts of saturated and unsaturated lipids, the area fractions of lo and l d domains also become unequal, which makes anti-registration geometrically impossible but would still permit registration.
We measured the stretching modulus of mixed membranes for systems within the two-phase region as well as for the corresponding coexisting pure phases. The two-phase mixture was softer than either of the pure phases and, unexpectedly, the l d phase was slightly more rigid than the lo phase. Adding cholesterol, the rigidity linearly increased, which we attribute to a concomitant decrease in area per lipid and an increase in lipid orientational order.
As an example of a dynamical observable, we measured lipid diffusion constants. While the observed differences are significantly smaller than what is known to be true in experimental systems (which is unsurprising for coarse-grained models), all changes follow expected trends. In particular, we found that the diffusion constant of lipids is more strongly determined by the phase they reside in rather than by their own identity.
These findings help to tie down some of the basic physical behavior of our model and open the door to a biophysically highly relevant and significantly more challenging situation: membrane asymmetry. Eukaryotic plasma membranes have leaflets that strongly differ in several important physical observables, most conspicuously their lipidomes. The cytosolic leaflet is substantially enriched in unsaturated lipids [22][23][24][25][26][27] and (at least for human red blood cells) appears to have significantly more phospholipids than the exoplasmic leaflet. 103 Average leaflet compositions are very close to a mixture of saturated lipids, unsaturated lipids, and cholesterol, in roughly comparable proportions, and hence prone to lol d phase separation. However, the strong difference in leaflet saturation, combined with a possible differential stress 104 associated with both packing and preferential cholesterol partitioning, 71 can move both leaflets away from phase separation or even just critical fluctuations and, therefore, away from the phenomena widely believed to be responsible for membrane rafts. 10 If so, then the intricacies of The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp how two ternary mixtures interact across their leaflets, share cholesterol, and balance uneven stresses and curvature torques all become collective regulators of rafts and the signaling they enable. Moreover, the underlying dynamics of raft formation would then be coupled to the kinetics of membrane asymmetry maintenance, breakdown, and restoration-namely, lipid flip-flop and lateral diffusion, stress diffusion, larger scale membrane trafficking, and the operation of passive scramblases and active fl(i/o)ppases. Studying this in computational models requires them to be able to reach the relevant length-and time-scales, but also to capture the physics of stretching and bending elasticity, lipid spontaneous curvature and sorting, asymmetry and selective flip-flop, differential stress, and mixtures able to exhibit lol d phase coexistence and a critical point. Our goal in this paper was to get closer to such a model.