A lattice model of ternary mixtures of lipids and cholesterol with tunable domain sizes

Much of our understanding of the physical mechanisms underlying the formation of raft domains in biological membranes is based on atomistic simulations of simple model systems, especially ternary mixtures consisting of saturated and unsaturated lipids, and cholesterol. To explore the properties of such systems at macroscopically large scales, we here present a simple ternary mixture lattice model, involving a small number of nearest neighbor interaction terms. Monte Carlo simulations of mixtures with different compositions show an excellent agreement with experimental and atomistic simulation observations across multiple scale, ranging from the local distributions of lipids to the macroscopic phase diagram of the system. The simplicity of the model allows us to identify the roles played by the different interactions between components, and the interplay between them. Importantly, by changing the value of one of the model parameters, we can tune the size of the liquid-ordered domains, thereby to simulate both Type II mixtures exhibiting macroscopic phase separation and Type I mixtures with nanoscopic domains. The Type II mixture simulation results fit well to the experimentally-determined phase diagram of mixtures containing saturated DPPC/unsaturated DOPC/Chol. When the tunable parameter is changed, we obtain the Type I version of DPPC/DOPC/Chol, i.e., a mixture not showing thermodynamic phase transitions but one that may be fitted to the same phase diagram if local measures are used to distinguish between the different states. Our model results suggest that short range packing is likely to be a key regulator of the stability and size distribution of biological rafts.


INTRODUCTION
Biological membranes are complex mixtures of many types of lipids, proteins and cholesterol (Chol) [1]. Though the lateral organization of such membranes is far from being fully elucidated, evidences indicate that it is highly heterogeneous [2]. Specifically, biological membranes contain small (10-200 nm) and dynamic domains, known as rafts, that are enriched in saturated lipids, cholesterol (Chol) and often particular proteins [2][3][4]. Rafts are liquid-ordered (L o ) domain "floating" in a sea of liquid-disordered (L d ) lipids [5,6]. The L d and L o phase are liquid crystalline ones, where the lipids diffuse in the membrane plane [7]. In the L o phase, the hydrocarbon chains are ordered, fully extended and tightly packed, similarly to the gel (S o ) phase [8]. Because of the enormous complexity of biological membranes, much of our understanding of the biophysical behavior and properties of raft domains is based on investigations of simple model membranes with few components, especially ternary mixtures of saturated and unsaturated lipids and Chol [9][10][11][12][13]. Studies of many different ternary mixtures reveal a phase diagram with regions of coexisting phases in the composition space, including coexistence between L o and L d phases which is believed to have relevance to raft domains in biological membranes [13][14][15][16]. Depending on the identity of the lipids and temperature, the coexistence of liquid phases may be of two types: In the ternary mixture of DPPC/DOPC/Chol, for instance, there are thermodynamic (macroscopic) phase separations between the liquid-disordered and liquid-ordered * ofarago@bgu.ac.il regions of the membranes over an wide range of temperature and composition [17]. When the doubly-unsaturated DOPC is replaced with POPC, which has only one unsaturated hydrophobic chain, no true phase separation is observed [17]. Feigenson [13] termed mixtures exhibiting macroscopic phase separation between liquid phases as Type II mixtures, and those where such phase coexistence is missing as Type I. Generally speaking, formation of small and transient domains may be associated with a vanishingly small line tension between the liquid phases. In Type II mixtures, this is expected in the one phase region of the phase diagram close to the critical demixing point [16]. In the two phase region, a comparatively higher line tension drives the macroscopically large segregation of the liquid-ordered phases. Additionally, there are also other mechanisms that may lead to the appearance of small domains, including in Type I mixtures not exhibiting macroscopic phase separation (see recent reviews in [14,18]), e.g., two dimensional microemulsion domains stabilized by line active molecules (lineactant) like hybrid lipids [19,20], or a modulated phase arising from a coupling between the local lipid composition and bilayer [21][22][23] or monolayer [24] curvatures.
Accessing the length and time scales relevant to lipid domains in complex membranes is computationally very challenging. Atomistic simulation is the putative method for characterizing the molecular organization in lipid membranes [25]. It provides a molecular scale picture which is not available experimentally. However, even for simple model membranes, all-atom large scale simulations are computationally too costly. Therefore, a coarse-grained (CG) description is employed to allow access to the larger scales. In the past 15 years, CG simulations involving the MARTINI force-field have be-come very popular in the membrane biophysics community [26]. In MARTINI simulations, several atoms are unified and represented by effective beads. One particular problem of MARTINI simulations which often leads to inconsistencies with atomistic simulations of similar system, is the chain entropy which is heavily involved in the phase transition but is not properly represented when using a smaller number of effective beads. Beyond CG MARTINI (and similar force fields) simulations of specific lipids, we have more (ultra) CG generic models designed for studying the general mechanisms and biophysical driving forces governing the thermodynamic behavior of lipid mixtures [27]. The highest degree of abstraction and computational efficiency is offered by lattice models [28][29][30][31][32][33][34][35][36][37]. Lattice models have been extensively used as a tool for studying phase transitions, and many predictions of these models have proven to be relevant to lipid mixtures [30,38]. Recently, we presented a lattice model for a binary mixture of the saturated lipid DPPC and Chol, that correctly captures the thermodynamic behavior of the system across multiple scales [39]. On the macroscopic scales, the lattice simulations reproduced the well-established phase diagram of DPPC/Chol mixtures, including regions where liquid-ordered domains of size ∼ 10-100 nm are observed in a liquid-disordered matrix. On the molecular scales, the simulations reveal the existence of gel-like nano clusters of size ∼ 1-10 nm within the larger Chol-rich liquid-ordered domains. Similar sub-structures have been seen in atomistic simulations of this system [40]. Fig. 1(i) shows the phase diagram of DPPC (saturated)/DOPC (doubly unsaturated)/Chol ternary mixture at T = 280 K. This canonical model system for studying the formation of liquid-ordered domains [16,17,41,42] has been identified as a Type II mixture. Depending on the mole fractions of the lipids and Chol, the mixture may be found in a single or coexisting phases. Other Type II mixtures also show similar phase diagrams, where the exact location of the phase separating lines depends on the identity of the lipids and the temperatures. Curiously, many Type I ternary mixtures also feature phase diagrams resembling Fig. 1 (i) [43,44], where the liquid phases are locally rather than macroscopically separated. The presence of small ordered domains in Type I mixtures is inferred from studies based on different experimental methods, e.g., NMR [16], small-angle neutron scattering [45], X-Ray diffraction [46], and fluorescence resonance energy transfer (FRET) [37]. These observations can be interpreted in one of two ways: (i) a coexistence between two phases with extremely small line tension, or (ii) a single phase with local heterogeneities. From a theoretical point of view, the difference between these two interpretations is marginal, but if the latter interpretation is adopted than the liquidliquid coexistence region in the phase diagram does not really exist. This controversy regarding the nature of a presumably local phase coexistence of liquid phases provides some of the motivation for the development of the model, presented herein, of ternary mixtures with tunable domain sizes. In the following, we present a minimal lattice model for a Type II mixture that fits very well the DPPC/DOPC/Chol phase diagram. The model contains a very small number of parameters, one of which directly controls the line tension between the disordered and ordered lipids. We demonstrate that by tuning the value of this parameter, we can change the nature of the mixture from Type II to Type I, i.e., convert macroscopically separated liquid phases into a mixture with traces of inhomogeneity at the local scale. Our study does not aim to provide a definite answer whether the liquid-ordered domains represent a distinct thermodynamic phase, but it explanis why they may be so interpreted.
In atomistic and CG simulations, the presence of liquid-ordered domains is often based on visual inspection of the lateral organization and chain conformations of the lipids and Chol. One particular problem in computer simulations is the fact that the liquid-ordered domains may themselves be heterogeneous and contain nanoscopic gel-like clusters. The presence of such Cholfree hexagonally-packed clusters of acyl chains in liquidordered regions that are rich in Chol have been observed in simulations [39,40] of binary DPPC/Chol mixtures, as well as in simulations of ternary DPPC/DOPC/Chol mixtures [47]. This local heterogeneous distribution of "domains within domains" is also consistent with experimental scattering data [48]. Another problem of detailed computer simulations are finite size effects. A simulation study of a ternary DPPC/DIPC/Chol mixture has pointed that the nature of the two phase region and domain sizes may be strongly influenced by the size of the simulated system [49]. This may explain result of recent atomistic simulations of DPPC/DOPC/Chol mixtures where phase coexistence between liquid crystalline L d and L o is studied in a system consisting of ∼ 1000 lipids [50]. That study, however, identified the system as a Type I mixture with domains of characteristic size of the order of 10 nm, which is in disagreement with clear experimental evidences indicating that this mixture is of Type II [17] featuring macroscopic phase separation. This, again, reinforces the importance of developing a computationally simple lattice model that allows simulations at large spatial and temporal scales. The attractiveness of lattice simulations lies not only in their computational simplicity, but also in their minimal nature which puts the focus on the factors that are essential for understanding the biophysics of lipid mixtures and the mechanisms governing their thermodynamic phase behavior.
The DPPC/DOPC/Chol lattice model presented below is an extension of our recent successful lattice model of binary DPPC/Chol mixtures [39]. As will be shown below, the model provides a multiscale picture of the ternary mixture over an unprecedented wide range of scales. At the macroscopic scales, the model yields very nice agreement with the experimentally-derived phase diagram of this mixture, featuring all single and coexisting phases. At the microscopic scales, the model unravels the inhomogeneous nano scale distribution of lipids and Chol. Specifically, the model agrees with atomistic simulations showing that the liquid-ordered domains contain small gel cores of highly aligned DPPC chains [47]. The simplicity of the model helps identifying that the origin of these gel-like clusters can be attributed to the effective packing interactions of the saturated chains with each other and with the Chol. We further demonstrate that the strength of the packing interaction between ordered DPPC and DOPC chains controls the size of the domains. Thus, we can artificially eliminate the thermodynamic L d + L o phase coexistence in DPPC/DOPC/Chol mixtures in favor of formation of smaller liquid-ordered domains. The resulting locally-inhomogeneous Type I mixture may be envisioned as either a single phase or as a coexistence of two liquid phases with negligible line tension which, quite surprisingly, also fit the DPPC/DOPC/Chol phase diagram.

METHODS
The lattice model of ternary mixtures presented herein is based on the model of DPPC/Chol binary mixture introduced in ref. [39]. Each chain is represented as a lattice point. Thus, the lipids with two acyl chains are represented as dimers while the Chol molecules occupy a single site. All Monte Carlo (MC) simulations are conducted on a triangular lattice of N s = 121 × 140 = 16940 (having an aspect ratio which is very close to unity) with periodic boundary conditions. The simulations are performed on a single planar lattice. This means that we do not explicitly introduce mechanisms speculated to be relevant for raft formation such as cross-monolayer interactions or curvature effects. A few lattice sites are left empty, thereby allowing the molecules to diffuse on the lattice which mimics the fluidity of the liquid membrane. The voids also allow the variations in the local density of the chains, which are expected in cases of coexistence between different phases. As in ref. [39], saturated DPPC chains may be in one of two states: ordered and disordered. Here, we introduce a second type of lipid with two unsaturated chains, say DOPC. These are always found in the disordered state because the system is simulated at temperatures higher than the low melting temperature of the DOPC. Thus, each lattice site obtains one of five possible states depending on the occupant: (i) a small area void (s = 0), (ii) a disordered (s = 1) or (iii) an ordered (s = 2) DPPC chain, (iv) a Chol (s = 3), (v) a DOPC chain (s = 4). In the lattice model, the detailed molecular picture of the forces between adjacent molecules is simplified and represented by nearest-neighbor interactions with effective parameters. The model Hamiltonian reads 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 other terms. The first term in Eq. (1) accounts for the fact that the disordered state of the DPPC chains is entropically favored by a free energy, −Ω 1 k B T (where k B is Boltzmann constant and T is the temperature), over the ordered state. The other terms represent effective packing attraction between an ordered DPPC chain with another ordered DPPC chain (second term), Chol (third term), and DOPC chains (fourth term). For the model parameters in the first three terms in Eq. (1), we use the same values as in ref. [39]: Ω 1 = 3.9, ǫ 22 = 1.3ǫ and ǫ 23 = 0.72ǫ, where the energy unit ǫ is such that the melting temperature of a pure DPPC membrane is T m = 314K = 0.9ǫ/k B [39,51]. The values of the model parameters agree to order of magnitude with experimental data for the enthalpy change at the melting transition of DPPC vesicles. 1 The fourth term account for the affin-1 The melting transition reflects a competition between the dis-ity between ordered DPPC and disordered DOPC chains, thereby directly influencing their degree of mixing and controlling the line tension between the liquid-ordered and liquid-disordered regions. We consider two different values ǫ 24 = 0.4ǫ and ǫ 24 = 0. The latter result in a Type II mixture that fits nicely to the DPPC/DOPC/Chol phase diagram, while the former yields a Type I mixture with small liquid-ordered domains that arguably may be fitted to the very same phase diagram. Note that although we could include interaction terms between other components (and, thus, refine the results presented below), we prefer not to introduce such interactions in order ordered and ordered states of the DPPC chains. The former is favored because of the entropy of the chains while the latter packs more efficiently with the surrounding chains. One should not expect ǫ 22 to be in perfect agreement with the meting enthalpy because this is an effective free energy parameter in a highly simplified model where the entropic and energetic influences of the implicit (missing) degrees of freedom are not entirely separable. Furthermore, there are entropic contributions that are not directly represented by the model parameters, e.g., the mixing entropy of the disordered lipids with the area voids. to keep the model minimal and highlight only the essential biophysical driving forces.
The MC simulations employ the following move types: (i) Displacements and rotations of randomly chosen lipids and Chol molecules (which are possible if there is a vacant neighbor site), and (ii) a change in the state of a DPPC chain (ordered to disordered and vice versa). Simulations are conducted at temperatures of 280K or 310K, at which the area per DOPC lipid is ≃ 0.67 nm 2 and ≃ 0.73 nm 2 , respectively [52]. At this temperature range, the areas per DPPC lipid and Chol are ≃ 0.55 nm 2 and ≃ 0.28 nm 2 , respectively [40,53]. We set the lattice spacing to l = 0.56 nm, which means that the linear size of the simulation cell is ≃ 70 nm. Depending on the composition of the simulated system, we set the number of voids such that the average density of the system matches the weighted average density of the constituting molecules. Systems are equilibrated until the memory of the initial configuration (either a random or a highly ordered distribution of the molecules on the lattice) is lost, after which we simulate the system for a period that is ∼ 10 times longer.
In the next section we present snapshots (see Figs. 1, 3, 4, and 6) where the sites of the lattice are drawn with distinct colors, denoting their belonging to liquiddisordered, liquid-ordered, and gel regions of the system. The partition of the sites into different categories is based in the following simple algorithm: The sites may be in one of five states: void, disordered DPPC chain, ordered DPPC chain, Chol, and DOPC chain. Each site is assigned with a score representing the degree of order of the associated state. Explicitly, liquid-ordered domains are rich in ordered DPPC and Chol, which receive positive scores of Sc i = 2 and Sc i = 1, respectively. Conversely, the liquid-disordered regions are populated with DOPC and disordered DPPC chains, which get negative scores of Sc i = −1, and Sc i = −0.5, respectively. Voids have a zero score. The grade of a site is given by where the sum runs over the six nearest neighbors of the site. Sites with non-negative (negative) grades, G i ≥ 0 (G i < 0), are associated with the liquid-ordered (liquiddisordered) regions. The maximum grade of G i = 14 is obtained for ordered DPPC chains surrounded by 6 other ordered DPPC chains. These sites are marked as gel.
After dividing the sites to the liquid-disordered, liquidordered, and gel regions, we scan the lattice one extra time and switch a liquid-disordered site to liquid-ordered if: (i) it is occupied by an ordered DPPC chain, and (ii) it has at least one neighbor belonging to a liquidordered or gel regions. The final iteration is needed in order to include ordered chains located at the periphery of the liquid-ordered domains that are surrounded by many disordered chains and incorrectly classified as part of the liquid-disordered regions.
For a given composition of components, we have generated many configurations starting at different initial distributions. We simulated the mixtures until they relax into states representing equilibrium and verified that different initial conditions evolve into similar phases: For instance, the snapshot in Fig. 1(b) shows a macroscopically large liquid-ordered domain in a liquid-disordered membrane (where some small liquid-ordered domains are also floating). This snapshot can be unambiguously interpreted as an indication that the equilibrium state of the mixture is a macroscopic L d + L o phase coexistence. All other independent MC runs with the same compositions also evolve into similar structures corresponding to macroscopically phase separation between the two liquid phases. Our goal in this study is not to plot the exact locations of the boundaries of the phase diagram. These are adopted from experiments [16], and we investigate whether the outcomes of the simulations are consistent with the experimental phase diagram. Indeed, a visual inspection of the equilibrated systems shows good consistency with the experimentally determined phase diagram. Our lattice simulations results also show reasonable agreement with atomistic simulations data for the compositions of the liquid-disordered and liquid-ordered regions (see Table I below).

Type II DPPC/DOPC/Chol mixtures
Setting ǫ 24 = 0 in Eq. (1) for the DOPC interactions while keeping the values of the other model parameters as in our previous work on DPPC/Chol binary mixture [39], yields a ternary DPPC/DOPC/Chol type II mixture model that reproduces the anticipated phase diagram. Our simulations at T = 280K are summarized in Fig. 1. Fig. 1(a)-(h) show representative equilibrated snapshots from simulations with different compositions, the location of which in the ternary phase diagram are marked in red in Fig. 1(i). The purple points in Fig. 1(i) indicate additional simulated compositions whose snapshots are not shown. In the snapshots, we color sites of the liquid-disordered, liquid-ordered, and gel regions with purple, yellow, and black, respectively The blue lines in the phase diagram mark the experimentally-determined phase boundaries of DPPC/DOPC/Chol mixtures, as adapted from ref. [16]. We note that our aim herein is not to determine the phase lines but demonstrate that the model yield results that fit them. Indeed, we see in Fig. 1 that the snapshot configurations fit very well into the different regions in the phase diagram. DPPC/DOPC/Chol is a Type II mixture and, indeed, we observe in the two-and three-phase regions that the different phases are macroscopically segregated. At the microscopic scale, we notice that the different phases are locally inhomogeneous. This is best seen in fig 1(h)  ing in a liquid-disordered sea. This arrangement, which is thought to be a simplified model of lipid rafts in biological membranes, clearly differs from the L d + L o coexistence appearing in snapshot (b).
Another example of a local inhomgeneity is seen in snapshot (g) of the L o phase, where small gel clusters appear in the liquid-ordered region. Conversely, in the gel phase in snapshot (e) we see local liquid-ordered impurities. These inhomogeneities do not constitute a different phase but rather represent local variations in composition. The sites marked by black color in the snapshots are the sites occupied by ordered DPPC chains with six ordered DPPC nearest neighbors. Since the gel phase in (e) also contains small amounts of DOPC and Chol, not all of its sites, only most of them, can satisfy this criterion. In both phases, we also find mesoscopic liquiddisordered islands containing mostly DOPC lipids. The presence of islands of disordered lipids within the ordered phases has also been observed experimentally [10]. The fluorescence microscopy micrographs appear as uniform L o or S o phases only when no DOPC is present. This experimental observation may indicate that the unsaturated lipid has extremely low miscibility in the ordered phases. Consistently with snapshots (e)-(g), the disordered islands did not coalesce into a single region on the experimental time scales [54]. It has been speculated that this follows from the slow dynamics of the lipids when the mixture has a very high density of ordered chains. In our simulations, we indeed noticed that the system tends to jam at this stage.
Small gel clusters residing within the L o phase have been also detected inside liquid-ordered domains in binary DPPC/Chol mixtures, suggesting that the thermodynamic origin of this observation is similar in both systems. This can be understood considering that the first three terms in Eq. (1) constitute the binary mixture Hamiltonian in ref. [39], and that we set to zero the coefficient in the fourth term accounting for the DOPC interactions. Thus, the picture inside the L o phase remain quite similar, as demonstrated in Fig. 2 showing smaller portions of different simulated systems. Fig. 2(a) is a magnification of the region marked by a red square in snapshot (b) of Fig. 1. Depending on the state of the lattice site, it is marked by the following color: white (s = 0, void), light blue (s = 1, disordered DPPC chain), green (s = 2, ordered DPPC chain), red (s = 3, Chol), and grey (s = 4, DOPC chain). Additionally, the boarders of the sites are marked with the same color coding as in the snapshots in Fig. 1 in order to indicate in which region (liquid-disordered -purple, liquid-ordered -yellow, and gel -black) they reside. Fig. 2(b) shows a magnified view from our previous simulations of DPPC/Chol mixtures [39], featuring similar local distributions of gel nano domains inside the liquid-ordered region. Fig. 2(c) is taken from our simulations of Type I ternary mix- tures, to be presented in the following section [see red box in Fig. 4(b) below]. In Type I mixtures, the domains are smaller but, nevertheless, exhibit, similar local hetrogeneities. The thermodynamic origin of these arrangements is the strong packing attraction between the ordered DPPC which drives their clustering into small gel domains. Chol, which prefers the vicinity of ordered DPPC chains (and, therefore, induces ordering of DPPC) is expelled from these clusters because its attraction to them is weaker then their mutual attraction. It accumulates in the liquid-ordered region, where it acts as a sort of a lineactant between the gel and the liquid-disordered region of the system.
The location of a second order critical point on the binodal curve separating the one phase and the two-phase L d + L o regions is roughly near the intersection of bin-odal curve with the line corresponding to DOPC mole fraction, φ DOPC = 0.4 [along which points (h) and (b) reside] [16,17,41,42]. The incomplete segregation of the system in (h) into small liquid-ordered domains and a liquid-disordered matrix is viewed as a precursor to the phase transition. According to this picture, the density fluctuations in the one phase region are expected to intensify in the proximity of the critical point, leading to a broader distribution of domain sizes and the emergence of a percolation liquid-ordered cluster. Evidences supporting this scenario may be found in Fig. 3 showing a typical configuration DPPC/DOPC/Chol at mole fraction composition of 0.23:0.38:0.39. The distribution of liquid-ordered domain sizes is indeed broader here compared to Fig. 1(h) which, presumably, is located further from the critical point. Moreover, one can notice that the liquid-ordered domains in Fig. 3 almost percolate across the lattice.

Type I version of DPPC/DOPC/Chol mixtures
In the previous section we saw that in Type II mixtures, like DPPC/DOPC/Chol, small liquid-ordered domains are expected in the one phase region, i.e., above the binodal line but not too far from the critical point. But what about Type I mixtures that do not exhibit a true thermodynamic coexistence? In this section we demonstrate that such mixtures may still feature small liquid-ordered domains that float in the sea of a liquiddisordered membrane. Explicitly, we consider the same system as in the previous section and only change the value of the model parameter ǫ 24 from 0 to 0.4, thereby allowing better mixing of the ordered DPPC and disordered DOPC chains. As we will now see, for this value of ǫ 24 , the mixture may arguably be regarded as a Type I version of the Type II DPPC/DOPC/Chol mixture studied in the previous section: It lacks thermodynamic phase transition on the one hand, but fits the DPPC/DOPC/Chol phase diagram if local measures are used to distinguish between the one phase and two phase regions. Fig. 4 summarizes our simulation results at T = 280K for ǫ 24 = 0.4ǫ. The color coding used here is similar to the one in Fig. 1, and the phase diagram in (i) is identical to Fig. 1(i), i.e., drawn with the blue phase boundaries taken from [16]. One clearly sees that the change in the value of ǫ 24 eliminates any signature of macroscopic phase separation, leaving only evidences of local phase separation, e.g., in snapshots (h), (b), (c), and (d). These snapshots indicate that the mixture is now of Type I. The question to be addressed now is whether the simulation snapshots in Figs. 4(a-h) still fit to the phase diagram in Fig. 4(i), which is the phase diagram of the Type II DPPC/DOPC/Chol ternary mixture. We first note that snapshots (a), (g), and (e) are consistent with the expectations to observe single L d , L o , and S o phases at the respective compositions. Identifying a two phase region when the system displays no macroscopic phase separation is, of course, harder. Snapshot (f), presumably showing coexisting L o and S o regions, exemplifies why the interpretation of the results is ambiguous. Recall that the gel regions are hexagonally packed arrangements of ordered DPPC chains. Therefore, we associate the gel phase in the snapshots with the lattice sites containing DPPC ordered chains with six nearest neighbors of the same kind. Both the liquid-ordered and gel phases appearing in (g) and (e), respectively, are highly populated with ordered DPPC chains. They differ in the fractions of chains obeying the six nearest neighbor criterion, which is high in (e) and low in (g). Looking at snapshot (f), one can interpret the observed distribution as a coexistence of regions of L o and S o phases, which are separated locally rather than macroscopically because of the smallness of the line tension between them. It is, however, impossible to rule out that the sequence of snapshots (e)-(f)-(g) simply shows a gradual increase in the fraction of sites classified as gel. In this scenario, snapshot (f) is interpreted as a single phase which is locally inhomogeneous with regions that look more gel-like and regions resembling the liquid-ordered phase.
The main question pertains to the coexistence region of two liquid phases, L d + L o . As noted above, the justification for drawing the binodal curve in the phase diagram (i) is based on an ambiguous interpretation of scattering data. The problem is nicely exemplified by snapshots (h) and (b) which, supposedly, reside on different sides of the binodal but display visually similar distributions of liquid-disordered and -ordered domains. This is not a result of an incorrect classification of the lattice sites to liquid-disordered and -ordered domains. The algorithm employed herein works well and yields very good consistency with data from ref. [50] reporting on atomistic simulations of a Type I version of DPPC/DOPC/Chol mixtures. Table I lists the average compositions of the liquid-disordered and liquid-ordered domains, measured for mixtures with overall DPPC/DOPC/Chol mole fraction of 0.35:0.35:0.30 (located in the L d + L o region of the phase space). The agreement between the lattice and atomistic simulations is more than fair, even though that the domain classification algorithms used here and in ref. [50] are different. A visual inspection of different snapshots indicates that the domain sizes in both works are also similar. We note here that it is not clear why the atomistic simulations fail to exhibit the experimentally-observed macroscopic phase separation of DPPC/DOPC/Chol mixtures. It may be due to a problem with the force fields or an artifact of the small system size. We nevertheless use them as a benchmark for our simulations of the Type I version of DPPC/DOPC/Chol mixtures.
The problem lies not in the accuracy of the algorithm for classifying sites to liquid-disordered, liquid-ordered, and gel regions, but in its resolution. The discussion above on the differences between snapshots (e), (f), and (g) has taught us that the algorithm, which is based on the inspection of the local vicinity of each lattice site, is too coarse and may not be sufficient for distinguishing between a single inhomogeneous phase and local coex-istence of two phases. This observation suggests that in order to resolve the differences between configurations (h) and (b), we must inspect more closely (i.e., define additional measures to quantify) the distribution of lipids at the nanometric small scales. Experimentally, coexistence of different phases may be inferred from NMR spectra which is sensitive to their symmetry and orderliness [16,41]. In our model, we can gain such insight from the statistics of the ordered DPPC chains that are classified by the algorithm as either in the liquid-ordered or gel regions. (Ordered DPPC chains barely exist in the liquid-disordered region.) Explicitly, we check all these ordered DPPC chains and compute of the probability P o of these chains to have n o = 0, 1, 2, 3, 4, 5, 6 other ordered DPPC neighbors. Our analysis is summarized in   Fig. 5 suggests that the most distinguishing feature between configurations in the one phase and the L d + L o two phase regions is the probability of having six ordered neighbors, which is quite negligible in the former and somewhat more substantial in the latter. This suggests that the presence of small gel clusters in the ordered domains is a defining property of the L o phase. In the L d + S o region, P o is a monotonically increasing function of n o . This is a characteristic feature of the gel phase. In the three phase region, the probability function exhibits features of both the L d + L o and L d + S o probabilities. As in the former, P o is not monotonically increasing; but similarly to the latter, there is a significant probability to detect an ordered chain with n o = 6 neighbors, indicating that the gel phase is also present.
Type II mixture at physiological temperature At T = 310K, most of the DPPC chains are found in the disordered state, which is understandable considering the proximity to their melting temperature T m = 314K. At this temperature, a single L d phase is expected throughout the entire phase diagram except, perhaps, for mixtures with a very high concentrations of DPPC [16]. Nevertheless, this phase may still contain small liquid-ordered domains, as demonstrated in Fig. 6, showing equilibrated configurations at compositions similar to points (b), (c), and (d) in Fig. 1 [denoted in Fig. 6 by (b2), (c2), (d2), respectively]. This is an important observation if one considers the ternary mixture model as a framework for understanding the formation of raft domains in biological membranes. As the snapshots in Fig. 6 demonstrate, our minimal model is capable of producing liquid-ordered domains of sizes < ∼ 10 nm even at the physiological temperature. These are half to one order of magnitude smaller than rafts in biological membrane, which is a well known limitation of ternary mixture model systems. This is why additional factors and mechanisms, e.g., curvature composition coupling [21][22][23]55], as well as the presence of specific proteins and the cell cytoskeleton [56,57], are believed to play a role in stabilizing larger domains. However, our results suggest that the same effect of increasing the domain sizes may be also accomplished by tuning differently the model parameters. After all, biological membranes are not ternary mixtures. They are composed of many lipid species with different melting temperatures and packing interactions, and our study suggests that even small changes in these parameters may affect dramatically the size distribution of the liquid-ordered domains.

DISCUSSION
We presented a simple lattice model for ternary mixtures of saturated and unsaturated lipids with Chol, involving only a minimal number of effective nearestneighbor interactions between the constituting components. The model succeeds in explaining experimental and atomistic simulation observations across multiple scale, ranging from the local distributions of lipids to the macroscopic phase diagram. The model Hamiltonian is constructed by adding a single interaction term, between unsaturated (DOPC) and saturated lipids, to a previously presented model of a binary mixture of saturated lipids (DPPC) and Chol. The minimal nature of the model allows us to identify the influences of the different effective packing interactions and the underling demixing mechanisms. The formation of the L o phase is obviously related to the affinity of Chol to the saturated lipids. This phase is inhomogeneous and contains gel-like nanoclusters whose origin is the particularly strong attraction of the saturated lipids to each other. The strength of the interaction between the unsaturated and saturated lipids, ǫ 24 , controls the nature of phase transitions of the ternary mixture and the sizes of the liquid-ordered domains. For ǫ 24 = 0, the simulated mixture exhibits a macroscopic phase separation and fits nicely to the phase diagram of a Type II DPPC/DOPC/Chol mixture. In this case, liquid-ordered domains are observed in the one-phase region near the critical point of the L d + L o coexistence binodal. For ǫ 24 = 0.4, we obtain the Type I version of a DPPC/DOPC/Chol mixture with nanoscopic domains. This system can be also fitted to the DPPC/DOPC/Chol phase diagram provided that local measures, such as the statistics of the ordered DPPC chains, are considered. In scattering experiments, nanometric features of this kind leave signatures in the scattering data, thus making the interpretation of configurations like in Fig. 4(b) ambiguous -either as coexistence of two liquid phases separated by a very small line tension or as a single phase with local inhomogeneities. Our aim here was not to decide which interpretation is correct, but to demonstrate that even if the phase transitions of the Type II DPPC/DOPC/Chol mixture are lost in the Type I version, local traces of phase coexistence still remain.
The motivation for studying lipid/Chol mixtures with a few components largely stems from the belief that they provide insight into the thermodynamic properties of cellular membranes. Because biological membranes are much complex physical systems, one should be cautious when coming to draw conclusions on their behavior from simpler model systems. On the one hand, the fact that biological membranes are mixtures of many different types of lipids suggests that it may be possible to construct a suitable model with many interaction terms that not only shows formation of liquid-ordered domains, but also reproduces their size distributions (in contrast to ternary mixtures where the domains are usually smaller than biological rafts). Several recent theoretical studies of multi-component mixtures have indeed demonstrated that strong affinity of certain components to each other may lead to demixing phase transitions that are more robust with respect to variations in temperature and the intermolecular interactions than mixtures with a few components [58][59][60]. On the other hand, it would be probably wrong to ignore other mechanisms beyond shortrange packing. This is particularly true for the inner cytoplasmic layer of cellular membranes that contains a low concentration of saturated lipids [61]. We can expect rafts to be influenced by the proximity to the cytoskeleton, the presence of proteins with affinities to certain lipids, and the curvature elasticity of the membrane (that depends on the local composition and induces coupling between the layers via the spontaneous curvature). Assessing the importance of these different mechanisms and understanding the interplay between them is going to remain an open question in the foreseeable future.