Minimal Lattice Model of Lipid Membranes with Liquid-Ordered Domains

Mixtures of lipids and cholesterol are commonly used as model systems for studying the formation of liquid-ordered ($L_o$) domains in heterogeneous biological membranes. The simplest model system exhibiting coexistence between $L_o$ domains and a liquid-disordered ($L_d$) matrix is that of a binary mixture of saturated lipids like DPPC and cholesterol (Chol). DPPC/Chol mixtures have been investigated for decades both experimentally, theoretically, and recently also by means of atomistic simulations. Here, we present a minimal lattice model that captures the correct behavior of this mixture across multiple scales. On the macroscopic scales, we present simulation results of mixtures of thousands of lipids and Chol molecules which show excellent agreement with the phase diagram of the system. The simulations are conducted on timescales of hundreds of microseconds and show the morphologies and dynamics of the domains. On the molecular scales, the simulations reveal local structures similar to those recently seen in atomistic simulations, including the formation of gel-like nano-domains ($\sim$ 1-10 nm) within larger Chol-rich $L_o$ domains ($\sim$ 10-100 nm). The observed multi-scale behavior is related to the tendency of Chol to induce ordering of acyl chains on the one hand, and disrupt their packing with each other, on the other hand.

Introduction: The plasma cell membrane is characterized by a lateral heterogeneous distribution of its biomolecular components [1]. Specifically, it contains small (10-200 nm) domains, termed lipid rafts, enriched in saturated lipids (mainly sphingolipids), cholesterol (Chol), and often particular proteins [2,3]. While many aspects of this "raft hypothesis" (including even their very existence in cellular membranes) remain controversial [4], there are strong evidences that raft domains play a key role in cellular processes such as signal transduction [5], cell adhesion [6,7], and membrane trafficking [8,9]. Rafts are "liquid-ordered" (L o ) domains where the lipid hydrocarbon chains are more ordered and more tightly packed than the surrounding non-raft "liquiddisordered" (L d ) phase of the bilayer [10,11]. The L o domains have features resembling both the L d phase (lateral lipid mobility [12]) and the "gel" phase where the lipid acyl chains are fully extended, tightly packed, and spatially ordered [13].
Biophysical studies of L o domains are rarely performed on biological membranes that contain hundreds types of lipids. Instead, simple model systems have been investigated, typically ternary mixtures of saturated and unsaturated lipids and cholesterol [14][15][16][17][18][19][20][21][22][23][24]. An even simpler system exhibiting domain formation is a binary mixture of Chol and saturated lipids like DPPC. Specifically, L o domains appear in coexistence with an L d phase at intermediate Chol concentrations and temperatures slightly above the DPPC liquid-gel transition temperature, T m [25,26] (see phase diagram in fig. 1). The domains appear to be both small and transient, not macroscopically phase separated from the L d matrix [27]. Recently, the homogeneity of the L o phase itself has been called into question [28]. Neutron scattering data at high Chol concentration has been interpreted as if the L o phase contains highly ordered nano-domains residing in a somewhat less ordered environment [29]. A recent atomistic simulation study of DPPC/Chol mixtures has provided further support to this picture [30]. The simulations reveal the existence of Chol-free hexagonally-packed small clusters of acyl chains within liquid ordered domains that are rich in Chol, especially along the domain boundaries. These new findings revoke a renewed discussion on the phase behavior of DPPC/Chol mixtures.
In this letter, we present a lattice model addressing the multi-scale behavior of DPPC/Chol mixtures. It reproduces the local inhomogeneous structure of the domains on the one hand, and the phase diagram of the whole system on the other hand. The model features a considerably smaller number of parameters compared to classical lattice models of DPPC/Chol mixtures like the 10-state Potts-like Pink model [31] and the model of Ipsen et al. [25]. It is also different than Almeida lattice model [32] that uses a top-down thermodynamic approach that assigns free energies to each state (gel/liquid disordered/liquid ordered). A unique feature of the model presented herein is the presence of empty sites that represent small area voids. This facilitates investigations of the diffusive dynamics of the lipids in the different phases. The minimal nature of the model and its ability to provide multi-scale information on both structural and dynamic properties of the system helps us to unravel the main thermodynamic and molecular forces underlying the formation of liquid order domains. Specifically, we identify the exclusion of Chol from both the gel and the L d phases as the reason for its accumulation along the interfaces between them. The accumulated high concentration of Chol molecules triggers the ordering of the lipid chains in their vicinity.
Model and Simulations: DPPC lipids are modeled as dimers of two chains, while the Chol molecules are represented by single chains. Each site of the lattice can be in one of four states (s = 0, 1, 2, 3): It can be occupied by a DPPC chain in either a disordered (s = 1) or an ordered (s = 2) state, a cholesterol molecule (s = 3), or be empty (s = 0) representing a small area void. The introduction of empty sites (i) allows the Monte Carlo (MC) simulations to mimic the diffusive dynamics of the lipids, and (ii) enables setting correctly the area per lipid in the different states (see details below), rather than making it a state-dependent adjustable parameter (as in [25,31]). This is also the reason for modeling lipids as dimers rather than monomers. The disordered state of the DPPC chain (s = 1) is entropically favored by a free energy F 1 = −Ω 1 k B T (where k B is Boltzmann constant and T is the temperature) compared to the ordered state (s = 2), but two nearest neighbor chains in the ordered state have interaction free energy F 22 = −ǫ 22 < 0 representing the favorable packing of chains in this state. The competition between these terms drives a first order melting transition of the DPPC bilayer at T m = 314K [33]. Chol interaction free energy with a nearest neighbor ordered chain is equal to F 23 = −ǫ 23 < 0. We set ǫ 23 < ǫ 22 , which means that Chol prefers locations close to the ordered chains, but with interaction free energy that is weaker than their packing interactions with each other. The physical length of the lattice spacing, l, is determined by considering a pure DPPC system (without Chol). In the gel state, we expect the ordered lipid chains to cover almost the entire lattice with only a negligible fraction of empty lattice sites. Thus, the lattice spacing is set to l = 5.5Å, which means that the area per lipid (which occupies two sites) is a ≃ 52Å 2 [30], and the area per Chol is a C ≃ 26Å 2 [34]. In the L d state, the area per DPPC molecule increases to roughly a ≃ 60Å 2 [30]. This means that about 15% of the lattice sites in the L d phase must be empty. In order to accomplish this, we introduce an elasticity-like quadratic term with a minimum at a ratio of c r = (85/2)/15 ≃ 2.83 between the number of lipids with two disordered chains (N 11 ) and the number of voids (N 0 ). Taking all together, the energy assigned to each MC configuration is given by where the deltas are Kronecker deltas, and the summations are carried over all N s lattice sites in the first term and over all nearest neighbor pair sites in the second and third terms in Eq. (1). Defining an energy unit ǫ = 1, we set the packing energies in Eq. (1) to ǫ 22 = 1.3ǫ and ǫ 23 = 0.72ǫ. The entropy parameter in the first term in Eq. (1) is set to Ω 1 = 3.9, and the elastic modulus in the fourth term is e m = 85ǫ. These values are in reasonable consistency with values of related parameters used in classical models that inspired our work [25,[35][36][37].
The MC simulations are conducted on a triangular lattice of N s = 121 × 140 = 16940 sites (with periodic boundary conditions) that has an aspect ratio close to 1. They involve three types of move attempts -(i) displacements of randomly chosen Chol or a lipid chain (dimer rotation) to a nearest neighbor vacant site, (ii) exchange of a lipid having two ordered chains with a pair of nearest neighbor voids and vise versa (22 ↔ 00: density-changing), and (iii) flipping the state (ordered/disordered) of a randomly chosen lipid chain (1 ↔ 2: state-exchange). In a pure DPPC systems, we start by dividing the system into two halves, one in the gel and one in the liquid-disordered phase. The system is equilibrated until one of these phases takes over the other. As expected, the melting transition is first order and takes place at k B T m ≃ 0.90ǫ. Chol is added by taking an equilibrated pure system and replacing randomly chosen lipids with two Chol chains, after which the system is further simulated until it settles into the new equilibrium distribution.
Results: We simulated the system at different temperatures (T ) and different molar fractions of Chol (χ). The phase diagram is shown in the upper left corner of fig. 1, where the magenta dots indicate the locations in phase space of the simulated systems, and the blue dots the locations of systems from which snapshots (a)-(g) are taken. The phase boundaries (solid black lines) are taken from ref. [26], and we find our observations to be consistent with them. At low concentrations of Chol, the system is in the gel phase for T < T m [ fig. 1(b)] and in the L d phase for T > T m [1(a)]. In the former the lattice is almost fully covered by ordered chains (depicted in green), while in the latter the lattice is mostly covered by disordered chains (blue) with a visible fraction of empty sites (white). At high Chol (red) concentrations (χ > ∼ 0.25), the system is found in a single L o state, containing a relatively high fraction of ordered chains [1(g)]. Snapshots (d), (e) and (f), which are from the L d +L o coexistence region, clearly show the appearance of domains of ordered lipids within a sea of disordered chains. These domains are dynamic and have irregular morphologies, suggesting that the system is not macroscopically phase separated.
In order to quantify the properties of the different phases, we measured the area per lipid, a = 52Å sharp crossover from L d (high a, low M ) to L o (Low a, high M ) is observed for 0.12 < ∼ χ < ∼ 0.24, which is exactly the location of the L d + L o coexistence regime in phase space. Similar trends are also detected in fig. 2(c) which shows the average diffusion coefficient, D, of the lipids. In full agreement with atomistic simulation results [30], we find the value of D in the liquid disordered state to be two orders of magnitude larger than that of the gel where the mobility of the lipids is nearly non-existing. In the liquid ordered state, D is about an order of magnitude smaller (larger) than in the liquid disordered (gel) state, implying that this is a relatively viscous liquid.
The values of D in fig. 2(c) are quoted in physical time units that are related to the MC time by using the experimental value of D = 11.2 × 10 −8 cm 2 /s at T = 322.5K and χ = 0.1 (L d state) [30]. The simulations extend over ∼ 5 mses, i.e., 3-4 orders longer than the atomistic simulations in [30]. The Supplemental Material movie [38] shows 100 µs dynamics (with a time difference of 4.6 µs between consecutive frames) of liquid ordered domains in the L d + L o coexistence regime. The characteristic life-time of the domains in the movie is τ < ∼ 20µs. We performed simulations with different ratios of translation vs. state-exchange move attempts (see simulation details, above) and found that increasing the fraction of state exchange move attempts has an insignificant impact on the observed dynamics. This indicates that the association/dissociation of the clusters is mainly governed by the lateral diffusion of the lipids rather than by the change in their states.
Figs. 1(h1) and (h2) provide a magnified view of the regions marked by squares in snapshots (e) and (d), respectively. These regions contain about 500 lipids which are plotted with the same colors coding as in the larger snapshots. Each point is also marked by a colored border, indicating whether the site is part of the L d region (purple), L 0 region (yellow), or a gel-like cluster of hexagonallypacked ordered chains (black) [labeled as "hexag." in (h2)]. The algorithm with which the different regions are identified is described in the Supplemental Material [38]. Fig. 1(h3) presents a snapshot of similar size from an atomistic simulation study of DPPC/Chol mixtures [30].
[The labels "L d ", "L o ", and "hexag." appear in the original copy and were adopted by us in snapshot (h2).] The similarity between the results of the atomistic (h3) and lattice simulations (h1,h2) is striking. In both, we observe Chol-free gel-like nano-domains surrounded by regions with high density of Chol and ordered chains. The composition of the system, i.e., fraction of sites belonging to each of the three regions (averaged over time and over the entire lattice) are plotted in fig. 3(a)  the gel-like clusters exhibits a curious non-monotonic behavior: It increases from zero as the system enters the L d +L o regime and then drops back to zero as χ continues to grow.
To conclude, our findings on the inhomogeneity of the L o domains are consistent with previous experimental [29] and atomistic simulations [30] studies. They suggest that the gel-like regions should be treated separately and not be counted as part of the L 0 state. The composition of the latter is plotted in fig. 3(b). It contains ordered chains (φ 2 ) and Chol (φ 3 ), with a small amount of voids and disordered chains (φ 0+1 ).
Discussion and Summary: Formation of L o domains has been associated with three general mechanisms [27]: (i) curvature-composition coupling, (ii) microemulsion stabilization by lineactant hybrid lipids, and (iii) near-critical fluctuations. None of these mechanisms seems to be relevant to the system under investigation. The model involves no curvature energy term (first mechanism), shows no evidence for substantial presence of lipids with ordered and disordered chains along the interfaces between L d and L o regions (second mechanism). With only the first two terms in Eq. (1), the system Hamiltonian can be mapped onto the Ising model. However, we find no signs of criticality (third mechanism) near the edge of the L d + L o coexistence regime [point (f) in fig. 1], suggesting that with the addition of Chol [third term in Eq. (1)] and the elastic energy term (fourth term), the model no longer belongs to the Ising universality class.
Our results call for a new thermodynamic explanation of the formation of the L 0 domains and their internal inhomogeneous structure. Adding Chol molecules to the membrane in the L d phase (T > T m ) causes an increasing number of chains in their vicinity to flip from a disordered to an ordered state, with which Chol interacts more favorably. However, the attraction between the Chol and the ordered chain is weaker than the interaction between the ordered chains with each other (ǫ 23 < ǫ 22 ), and provides only partial compensation for the loss of chain configurational entropy. Therefore, as the concentration of Chol and ordered chains somewhat grows, the later partially segregate and form the gel-like clusters. The distribution of lipids and Chol in fig. 1(h1)-(h3) may be interpreted as a coexistence of regions in the gel and L d phases which, close to the first-order melting temperature [(T − T m )/T m < ∼ 0.02], are nearly-equivalent thermodynamically. Chol is excluded from the gel and little dissolves in the L d region. It accumulates at high concentration along the gel-L d interfaces and induces ordering of the lipid chains, creating a L o "buffer" region. This may be regarded as if the Chol molecules serve as lineactants between the ordered and disordered lipid chains. The composition of L o state does not change much in the coexistence region of the phase diagram [see fig. 3(b)], suggesting that it may be regarded as a metastable phase. Indeed, a single chain exchange from an ordered to a disordered state (while keeping unchanged the states of the others) results in increase in the system free energy for most of the chains in the L o phase (see data in the Supplemental Material [38]). This explains the prolonged lifetime of the domains. At higher values of χ, Chol "invades" the gel clusters and eliminates them [ fig. 3(a)], and the L o domains merge into a single L o phase.
To summarize, our model provides a dynamical picture of DPPC/Chol mixtures over an unprecedented wide range of length-and time-scales. Its minimal nature helps us identifying the disparity between the packing interactions of ordered chains with Chol and with each other, as the reason for their non-ideal mixing within the L o nano-domains. At larger spatial scales, there is no macroscopic L d /L o phase separation. The system arranges itself such that the Chol-rich L 0 regions separate the gel and Chol-poor L d regions. The basic thermodynamic mechanism proposed here may also apply to ternary model mixtures of saturated and unsaturated lipids with Chol (where, indeed, similar heterogeneous domains have also been recently identified [24]) and may be relevant to raft domains in biological membranes. This work was supported by the Israel Science Foundation (ISF), grant No. 991/17. TS thanks the Planning and Budgeting Committee of the Council for Higher Education (Israel) for supporting his post-doctoral fellowship.