A Rationale for Mesoscopic Domain Formation in Biomembranes

Cell plasma membranes display a dramatically rich structural complexity characterized by functional sub-wavelength domains with specific lipid and protein composition. Under favorable experimental conditions, patterned morphologies can also be observed in vitro on model systems such as supported membranes or lipid vesicles. Lipid mixtures separating in liquid-ordered and liquid-disordered phases below a demixing temperature play a pivotal role in this context. Protein-protein and protein-lipid interactions also contribute to membrane shaping by promoting small domains or clusters. Such phase separations displaying characteristic length-scales falling in-between the nanoscopic, molecular scale on the one hand and the macroscopic scale on the other hand, are named mesophases in soft condensed matter physics. In this review, we propose a classification of the diverse mechanisms leading to mesophase separation in biomembranes. We distinguish between mechanisms relying upon equilibrium thermodynamics and those involving out-of-equilibrium mechanisms, notably active membrane recycling. In equilibrium, we especially focus on the many mechanisms that dwell on an up-down symmetry breaking between the upper and lower bilayer leaflets. Symmetry breaking is an ubiquitous mechanism in condensed matter physics at the heart of several important phenomena. In the present case, it can be either spontaneous (domain buckling) or explicit, i.e., due to an external cause (global or local vesicle bending properties). Whenever possible, theoretical predictions and simulation results are confronted to experiments on model systems or living cells, which enables us to identify the most realistic mechanisms from a biological perspective.


Introduction
The plasma membrane, a complex mixture of lipids and proteins, forms a selective barrier for eukaryotic cells [1,2], yet its role goes far beyond a simple frontier delimiting the cell interior and exterior. In the original 1972 model by Singer and Nicolson [3], the plasma membrane was seen as a more or less homogeneous mixture in which the proteins represent about 50% of the total mass. Since then, this model has known regular improvements, leading to the deciphering of an increasing organizational complexity, notably at the sub-micron level, and a growing understanding of the biophysical and biochemical mechanisms at play [4][5][6][7][8][9][10]. There are evidences now that the membrane plays a crucial, active role in a large amount of biological functions [1] such as viral and bacterial infection, immune response, cell adhesion, transport of solutes or signaling, to name a few, and its organization is directly related to its biological functions. The membrane is made up of a lipid bilayer (mainly phospholipids, sphingolipids and cholesterol), the arrangement of which minimizes contact between water and the hydrophobic tails of these amphiphilic molecules [11]. The above-mentioned lipid mixture has been shown to present in eukaryotic cells two main distinct phases, termed liquid-disordered (Ld phase, cholesterol-poor) and liquid-ordered (Lo phase, cholesterol-rich) [7,[12][13][14] that undergo phase separation [15][16][17][18] for a specific composition range and below T ∼ 20-30 • C depending on the nature of the lipid mixture, as it was already understood in the 1980's [19]. The minimal requirement for such liquid phase coexistence seems to be a ternary mixture of low-and high-melting temperature lipids, and cholesterol [20].
The plasma membrane also contains various inclusions made of peripheral or integral proteins [1,12]. Thanks to recent in vivo and in vitro experimental developments such as single particle tracking (SPT) [21][22][23], fluorescence (or Förster) resonance energy transfer (FRET) [20,24], atomic force microscopy (AFM) [25][26][27][28][29], super-resolution microscopy techniques (STED, PALM, STORM, SIM, see the review [30]) or small-angle neutron scattering (SANS) [31][32][33], it has been observed and it is now widely agreed that cell membrane components are heterogeneously distributed, and are generically organized into functional lipid and protein sub-micrometric or nanoscopic domains (nanodomains for short). In particular, a common view pictures lipid nanodomains as lipid "rafts" [5,10,12,17,34], a consensual definition of which was formulated in 2006 [35]. They are dynamic supramolecular assemblies of nanometric length-scale (10 to 200 nm) enriched in sphingolipids and cholesterol (the above Lo phase) within a "sea" of Ld phase. Specific proteins are then supposed to be targeted to them to perform their biological functions. However, the concept of "raft" remains controversial, especially in living cells [14,17,34,[36][37][38][39][40][41][42][43][44]. In particular the mechanisms for membrane rafts being so small are still debated. Indeed, creating many small domains instead of a single large one generates a longer boundary between Lo and Ld phases, resulting in an increased interfacial energy. Even their existence as just being the consequence of a lipidic meso-phase separation, independently of proteins interacting with lipids and other proteins, is debated. Nevertheless, at present, it can safely be concluded that most membrane proteins, as well as lipids, segregate into different domains via specific mechanisms. All these domains have proven to be key players in the above-mentioned biological functions, but a full and consensual understanding of their biophysical and biochemical origin is still lacking.
The aim of this review is precisely to expose the different models that have been proposed from a physical perspective and to deliver a rationale to better understand the mesophase formation and stability on length-scales of biological interest. In the biophysical literature, the observed domains below the diffraction limit are often named "nanodomains" because their size ranges from tens to hundreds of nanometers. In what follows, (nano)domains will refer either to protein domains (also known as "protein clusters" [30]), or to lipid domains with a composition different from the membrane bulk, or to lipid domains with a phase state different from the membrane bulk (e.g., Lo domains in a Ld phase), or (and most likely) to any combination of these three situations. Indeed, most continuous theories that will be under consideration in this review are sufficiently general to embrace such a large spectrum of physical realities [15,45,46]. When needed, we will specify which system is considered in particular. Apart from cell biology, alike membrane patterned morphologies can even be observed in different contexts of biology, such as pollen grains, fungal spores or insect eggs [47].
Two major classes of processes stand out and will be addressed here: mechanisms considering the membrane in thermodynamic equilibrium, and active or non-equilibrium processes, which integrate the dynamical interactions of the membrane with the cytosol, especially active membrane recycling.
In thermodynamic equilibrium, different models, many of them being shown in this review to be based upon the symmetry-breaking principle, provide explanations for small domains formation. In this respect, two main classes of mechanisms can be distinguished: spontaneous buckling and bending-related phenomena. Lipowsky indeed described in 1993 that membrane domains can arise due to a spontaneous symmetry breaking mechanism [48]. When the system is quenched below the phase-separation temperature, domains coarsen, which increases their perimeter and then the interfacial energy. They eventually bend undergoing buckling in the third dimension (perpendicular to the membrane plane), so that the interface energy cost between the two phases is lowered. We will discuss into detail below that they are then trapped in a metastable state.
A second category of equilibrium models (and experiments) show that domain formation can rely on the coupling between curvature and local composition. The global frame of this type of mechanisms was underlined by Seul and Andelman in 1995 [45], bringing into play a competition between short-range attraction (leading to a finite line tension at interfaces between different phases) and longer-range repulsion, for example generated by electrostatic (e.g., dipole-dipole) interactions between lipids or proteins, and/or by their coupling to the membrane curvature as in the present case. A typical patterning length-scale emerges set by the relative strengths and ranges of attraction and repulsion [46], and can be tuned by varying control parameters such as temperature, membrane tension or constituent concentrations. Generally speaking, local curvature occurs whenever the up-down symmetry is locally broken due to the presence of different lipid compositions in the two leaflets and/or non-symmetric proteins, or different types of solutes that can adsorb on the two sides of the membrane. By "broken up-down symmetry", we mean that for a variety of possible reasons that will be discussed further in this Review, the upper and lower leaflets of the bilayer, or equivalently the inner and outer leaflets of a vesicle or cell, do not play the same role.
A third type of membrane symmetry breaking in equilibrium is related to the presence of lipid species with different numbers of unsaturations, leading to thicker domains (often enriched with cholesterol, such as Lo phases) with a higher bending modulus inducing hydrophobic mismatch. This results in total demixing in a plane bilayer below the phase-separation temperature. But, on a curved surface like a vesicle, the larger these domains are, the more difficult to bend they are. The minimization of the bending energy cost then favors the division of large thick domains into smaller ones that better accommodate the membrane to the global spherical shape [49].
In this review, we will not consider in depth the potential role of so-called "linactants" (or line-active molecules), the two-dimensional equivalent of surfactants [11]. Such molecules (proteins or hybrid lipids having one saturated and one unsaturated hydrocarbon chain) localize at the interface between coexisting phases and lower the interface line tension. They potentially stabilize patterned structures and micro-emulsions, which are actually strongly related to the concept of rafts [12]. The interested reader can refer to the Reviews [7,14,18,50] and to our Discussion Section.
On the other hand, non-equilibrium models describe how active membrane recycling enters in competition with domain growth and then modulates their size in the stationary regime. Perpetual recycling of the cell plasma membrane due to material exchange with the cytosol, due for example to exocytosis and endocytosis of membrane patches, continually mixes the plasma membrane components and breaks too large membrane domains.
This review is organized as follows. After an introduction, Section 2 puts into perspective the principal mechanisms in equilibrium accounting for mesophase separation in biphasic planar membranes or vesicles. The mechanisms deal with either separation of lipidic phases such as Lo and Ld ones, or condensation in domains of interacting proteins embedded in lipid bilayer. We make the distinction between the so-called "weak segregation limit" in the vicinity of the miscibility critical point where the field theory at Gaussian order is often used. In principle, field theories at Gaussian order are only valid above the critical temperature and far from it [15]. However, applying them close to the critical point amounts to use mean-field theories that provide reasonable orders of magnitude of critical temperatures and critical exponents that are very useful as a first step. To go beyond these simple estimates, far more involved renormalization techniques must be used, beyond the scope of this review. The "strong segregation limit", well below the critical temperature, where effective models using the notion of line tension are more efficient. We stress how up-down symmetry breaking is then involved in the destabilization of the macrophase in favor of the mesophase. In combination with experiments, the theoretical tools to tackle these questions go from all-atom and coarse-grained molecular dynamics [51] and Monte Carlo simulations [52], to analytical modeling, notably using a continuous, field-theoretic description of the membrane [45]. Section 3 addresses the characterization of out-of-equilibrium membrane steady-states. Using tools from out-of-equilibrium statistical mechanics, either numerical or theoretical, the reviewed studies predict domain-size distributions that will have to be confronted to experiments in the future. This will enable biologists and biophysicists to identify the mechanisms that are actually at play in living cells. Finally, a discussion and conclusion Section discusses the different models from a cell biology perspective and proposes possible routes to be followed in future studies. Table 1 gives an overview of the principal notations used throughout the review. Notably, the natural energy scale will turn out to be the thermal energy k B T 4.3 × 10 −21 J at physiological temperature. Recycling time (out-of-equilibrium membranes) 3.1 [56,57]

In Thermodynamic Equilibrium
Two-dimensional patterning has been known for decades in condensed matter physics and its origin can generically be accounted for by a "competing-interactions model", as described by Seul and Andelman in their seminal paper [45]. In their article, mesophases or "modulated phases" are described in a variety of contexts, such as semiconductor surfaces, superconductor films, liquid crystal films, ferrofluids, polymer mixtures, diblock copolymers, Langmuir films, and biphasic biomembranes.
The first theoretical works on the formation of curved meso-structures in membranes have been done by Leibler [58] and Leibler and Andelman [59] in the 1980's by adopting a field-theoretic approach. They showed how coupling the membrane composition to its mean curvature can lead to modulated phases. At the mean-field level the system made of two types of molecules A and B is described as a function of the thermodynamic variable φ which is the area fraction of molecules of type A, and therefore 1 − φ is the one of molecules of type B.
If the coupling between composition and membrane height fluctuations is neglected, a second-order macro-phase transition occurs between a disordered phase where lipids A and B are mixed and an ordered phase where two phases, the A-rich and B-rich ones (alternatively Lo and Ld phases), coexist. The critical area fraction is φ c = 1/2 and the critical temperature, noted T c , is related to the Flory interaction parameter χ ∝ (α A − α B ) 2 where α A,B are the electronic polarisabilities of molecules A and B in this context. Close to T c , the grand potential per unit area can be expanded in powers of ψ = φ − 1/2 [15] where f (ψ) is the free energy and µ the chemical potential (both per unit area). The coefficients in the expansion, m = α(T − T c ) (with α > 0) and s are proportional to T − T c , whereas u is almost independent on T and h = µ − µ c − γ(T − T c ). In particular the equation of state is given by ∂w/∂ψ = ∂ f /∂ψ − µ = 0, and the coexistence curve in the (ψ, T) plane can be drawn close to T c . Note that if the area fraction is different form the critical value 1/2, the temperature at which the coexistence occurs is smaller and noted T d (φ) < T c [15]. It is called the "phase-separation" or "demixing" temperature.
To take into account the spatially varying composition φ(r), one constructs the phenomenological Landau-Ginzburg Hamiltonian where a term in b 2 (∇φ) 2 is introduced, characterizing the energy cost associated with local variations of the concentration. The coefficient b is related to the pair interaction between molecules A and B. Below T d , the system phase-separates into two phases and an interface sets up. Minimizing this Hamiltonian (with m negative, and s = 0 which can always be fixed by shifting φ) leads to the classical interfacial tension [11,60], in this 2D case a line tension λ [18,48], with λ ∝ b∆φ 2 /ξ OZ ∝ |m| 3/2 √ b/u where ∆φ = |m|/u is the difference in composition of the two phases and ξ OZ = b/|m| is the Ornstein-Zernike correlation length of composition fluctuations in the bulk and far from the critical point [15,54].
However, the lipids do not evolve on a stiff planar surface but on a fluctuating elastic membrane. In the case of an homogeneous lipid mixture, Helfrich proposed in 1973 that the elastic energy of a membrane or vesicle is given by the integral [2,11,55] over the whole membrane surface A. Here H = (c 1 + c 2 )/2 is the local mean curvature, K = c 1 c 2 the Gaussian curvature, and C sp the spontaneous mean curvature. It measures the more or less pronounced membrane tendency to bend spontaneously, upward or downward. The bending modulus κ typically falls in the 10 to 100 k B T interval [12,14,61]. In this respect, thicker Lo phases are known to have higher bending rigidities that Ld phases. For example, measured ratios are κ Lo /κ Ld 5 in reference [62] and 4 in reference [63], a range of values confirmed in more recent studies measuring ratios up to 10 [32,33]. Membrane stability requires that −2 < κ G /κ < 0 (see, e.g., [64]). The membrane surface tension is denoted by σ (see Figure 1a). Here σ appears as a Lagrange multiplier controlling the membrane area. There are several alternative definitions of the surface tension that coincide in the high tension limit [65]. It is imposed by external constraints, and cannot exceed the so-called "lysis" tension, on the order of 10 −2 N/m for usual lipids such as 1,2-Dioleoyl-sn-glycero-3-phosphocholine (DOPC).
In the case where the membrane adopts a globally planar configuration, it can be described by a height function h(r) measuring the distance to a reference plane, as in Figure 1a. Then H ∇ 2 h/2 when the fluctuations of h remain small. The coupling between local composition and membrane height fluctuations was introduced by Leibler and Andelman in the following effective field-theoretic Hamiltonian (by skipping the irrelevant constant and ignoring the Gaussian curvature at this stage): where S is now the projected area in the (xOy) plane. The constant Λ couples the composition φ(r) and the membrane mean curvature ∇ 2 h/2, and is connected to the difference C 1 between the spontaneous curvatures of lipids A and B, through Λ = −κC 1 (here C sp = C 1 φ). A value |Λ| 5 pN has been proposed by fitting experimental data [66]. It is quite consistent with the typical values κ ∼ 10 k B T and C 1 ∼ 0.1 nm −1 discussed below. A value as large as Λ ∼ 100 pN can reasonably be reached for more rigid Lo phases and stronger spontaneous curvatures.
(a) (b) Figure 1. (a) Membrane simulation snapshot. The blue lipids are 1,2-dipalmitoylphosphatidylcholine (DPPC), the yellow ones di-C16:2-C18:2 PC (DIPC) and cholesterol appears in green. The pink thick line is the 1D interface delimiting the liquid-ordered (Lo) phase (cholesterol rich, φ(r) 1) and the liquid-disordered (Ld) phase (cholesterol poor, φ(r) 0). This interface has an energy cost per unit length, homogeneous to a force λ, named the line tension. Although fluctuating, the membrane is globally planar and parallel to the plane (xOy). The height of the membrane above this reference plane (xOy) is measured by the height function z = h(r). The membrane is taut with surface tension σ as illustrated by the four gray arrows parallel to the plane (xOy). Adapted from an original membrane image generated by the MARTINI force field, reproduced with the courtesy of Matthieu Chavent. (b) In the cluster-phase scenario, proteins are described as individual objects embedded in a continuous fluctuating 2D lipid mattress, also represented by a height function. They gather because of short-range attractive forces but long-range repulsion between clusters limits their growth. In the present case, each individual protein (in purple) locally imposes a spontaneous curvature to the membrane that is represented as an elastic sheet (in orange). In a taut membrane, an effective long-range repulsion between proteins ensues (see text). The units are arbitrary.
By integrating over the field h(r) (which in the present Gaussian approximation can be done exactly), the effective Landau-Ginzburg Hamiltonian becomes where the stiffness b has been "renormalized" to b = b − Λ 2 /σ, which can now be negative for either low enough surface tension σ, or low b or low enough line tension λ, or high enough coupling Λ. It signals the onset of a curvature instability of the membrane. If b < Λ 2 /σ, a "pattern" spontaneously emerges with the characteristic length-scale (see also the beginning of Section 2.1.1 below): Note that this field-theoretic approach does not predict whether the emerging patterns look like roundish "bubbles" of an A-rich phase in a B-rich phase, or more elongated domains, or labyrinthine structures, or even alternate stripes [45]. We shall come back to this point below.
Cluster phases are another instance of patterned structure extensively studied in soft condensed matter physics because they also appear in a wide variety of contexts [67]. For example, if charged colloidal particles are in suspension in water, they feel a long-range mutual repulsion in low ionic force solutions. If small polymers are added to the solution, colloids also experience a short-range attraction usually named a depletion interaction [68] (in addition to the classical steric repulsion). Above a critical colloid concentration, small colloid clusters emerge in equilibrium, the typical size of which is also set by the attraction-to-repulsion ratio and the concentration [67]. The underlying mechanism is as follows: short-range attraction tends to favor phase separation, but long-range repulsion between domains limits their growth. The macro-phase separation where a single large colloid assembly would be formed is not reached because the long-range repulsion destabilizes too large clusters. The ideas that we will present below amount to extend this mechanism to protein assemblies. Whereas the preceding field-theoretic approach applies to both lipid and protein nanodomains because it is based on local concentrations, the cluster-phase scenario is more specific to protein clustering because it deals with assemblies of individual objects in a continuous solvent (proteins in the 2D lipid mattress), as illustrated in Figure 1b [69].

Weak-Segregation Limit in the Vicinity of a Critical Point
In this section, we assume that the lipid mixture is close to the miscibility critical point [15,16,18], the so-called weak-segregation limit. The temperature T is therefore close to (just above or just below) the critical temperature T c . In this limit, density fluctuations play an important role and one can construct a field theory using the Landau-Ginzburg formalism that includes weak density gradients and therefore wide interfaces between fluctuating lipid domains, which do not have a precise shape. On the contrary, in the strong segregation limit, where T is well below the critical temperature T c , the two separated phases form well defined lipidic domains with sharp boundaries, and one can compute the free energy of chosen patterns, controlled by the interfaces between domains. This situation will be considered in the Section 2.2.

Curvature-Composition Coupling in Planar Membranes
Spontaneous curvature C sp measures the more or less pronounced membrane tendency to bend spontaneously. Several mechanisms can lead to local spontaneous curvature C sp = 0 of lipid bilayers. The most evident one in cells is the marked difference of composition between the inner (cytoplasmic) and the outer (exoplasmic) leaflets of the plasma membrane [1]. Some lipids have a global conical shape (to say it briefly, either "big" head-group and "small" tail(s), or "small" head-group and "big" tails) and locally concentrating them in one of the leaflets leads to non-vanishing spontaneous curvature [12,70]. More precisely: (a) The difference in lipid composition between both leaflets is important in cellular membranes, and it is maintained by the active cell metabolism. It can lead to bilayer spontaneous curvature if both leaflets conspire in this direction, because the bilayer curvature results form the difference in the spontaneous curvature of the monolayers [12]. The spontaneous curvatures of the main lipids found in plasma membranes are listed in [71] and they can be as large as 0.3 nm −1 for cholesterol or 1,2-Dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE). This is global on the whole membrane, but it can be accentuated locally due to the membrane lateral heterogeneity. For example, it has been shown on the basis of coarse-grained molecular dynamics simulations that mean curvatures of about 0.1 nm −1 can be attained in asymmetric membranes containing separated Lo and Ld phases on one leaflet and pure unsaturated lipid on the other leaflet [72]. (b) The difference in the aqueous solution composition on the two sides of the membrane is maintained by the cell [1]. As explained by Lipowsky in 2013, a difference of solute concentrations, including ions and small molecules, generically leads to spontaneous curvature when they adsorb onto the membrane surface, for purely entropic causes [73]. The membrane "bends away from the exterior compartment if the concentration c ex in this compartment exceeds the concentration c in in the interior compartment". For a single solute with different concentrations across the membrane, the spontaneous curvature is given by where is the membrane thickness, Γ max is the maximal surface density by adsorption and K d is the equilibrium constant of adsorption. Putting realistic numbers in this relation (in particular κ ≈ 20k B T and ≈ 4 nm [74]), one gets C sp ≈ 1/20 nm −1 . The membrane itself is supposed to be up-down symmetric here. The adsorption of biopolymers is also examined in this work [73]. However, it is demonstrated to lead to smaller spontaneous curvature in realistic regimes of parameters. (c) The area difference between both leaflets can also lead to global spontaneous curvature.
For instance, an area difference of ∼ 10% leads a to spontaneous curvature C sp ∼ 10 −2 nm −1 [71]. This is the keystone of the area-difference-elasticity (ADE) model that has been developed to explain the rich shape variability of homogeneous lipid vesicles, in particular in function of their reduced volume v [12,75].
Locally curving a membrane under tension leads to an increase of the Helfrich energy due to the induced area excess. We shall see below that this energy penalty is lower when splitting a large curved domain into smaller ones and is then favorable, in spite of the increased interfacial energy penalty. In this respect, the field theory given above was very basic as compared to the actual developments of the three last decades, an overview of which is now proposed. In this section, we focus on theoretical results obtained above the critical temperature where Gaussian approximations are useful. Presumably they remain valid in the vicinity of the critical point. Many results are similar to what has been observed in the context of micro-emulsions, made of a bicontinuous mixture of oil, water and amphiphile where one order parameter φ is now the local amphiphile concentration and the other one is the local concentration difference between oil and water [76,77].
To be more quantitative, the membrane structure is well characterized by the structure factor (or power spectrum) S(q), which is the Fourier transform of the spatial correlation function of lipid composition [15,18]. At the level of the linear response, S(q) characterizes the amplitude of the response to an external perturbation of the local composition Fourier transformφ(q) at a wave-vector q. It is therefore the generalization of the susceptibility of the system at a given wave-vector q.
At the Gaussian level where w(φ) = mφ 2 (with m > 0) one can generically write the effective Landau-Ginzburg hamiltonian asH The classical structure factor for an inhomogeneous mixture has the form S(q) = 1/(α(T − T c ) + bq 2 ), a Lorentzian of half-width ξ OZ . It diverges at the transition T = T c at q = 0, the signature of a macrophase separation. By plugging in the coupling to the shape fluctuations, it is simple to show from Equation (5), that S(q) has a second maximum at q c = |b |σ 2 /(Λ 2 κ) when b < 0. By introducing the Helfrich correlation length ξ = √ κ/σ and using Λ = −κC 1 , one can rewrite Hence this secondary maximum occurs as soon as C 1 is larger than C * 1 , i.e., for relatively large difference in spontaneous curvatures. Indeed, we recall that C 1 is the difference between the spontaneous curvatures of both phases, in other words C sp (φ) = C 0 + C 1 φ at the lowest order in φ.
This phase is a homogeneous (or liquid) but structured one, also called "Structured Disordered" (SD) phase, where correlations between lipids at a typical distance 1/q c exist, which reflects a tendency towards order. The binary mixture is more susceptible to form transient structures of this size when it is perturbed. This regime of structured liquid on a length on the order of ξ is consistent with the nanodomains of size ranging from 10 to 100 nm, for realistic cell membrane elastic parameter values [78,79].
Within this picture, this secondary maximum diverges when m = |b |q 2 c /2, i.e., for T * = T c + |b |q 2 c /(2α) > T c [80]. Note that T * depends on C 1 through q c . It is a monotonically increasing function of C 1 . Hence a modulated phase appears for temperatures T ≤ T * where the system undergoes mesophase separation [78,79,81] which has been observed in simulations [72,82]. The occurrence of stripe phases, square phases and hexagonal phases has been predicted using the single mode approximation and phase diagrams have been constructed for two and even three component fluids [81].
Schick and coworkers argued that in plasma membranes both modulated phases (for which S(q) diverges at q c ) and structured liquid (as soon as C 1 > C * 1 ) can be observed [78,79,83]. This would be the origin of lipid rafts, the typical size of which is 2π/q c 100 nm. Note that this domain size is equal or larger than the Helfrich correlation length, ξ, according to Equation (9), that is fixed by the balance between bending and surface energy. The classical interpretation for the criterion C 1 > C * 1 is the following: creating a bent domain is possible if the gain in bending energy ≈ ξ 2 κC 2 1 , is larger than the cost in demixing, ≈ b(∇φ) 2 ξ 2 ≈ b. The last relation comes from the fact that on domains of size ≈ ξ, |∇φ| ≈ ξ −1 .
In references [84][85][86], it is argued that a more correct way to introduce the coupling between composition and curvature is to write explicitly in the Helfrich Hamiltonian that the spontaneous curvature C sp is a function of φ. The full bending energy term in Equation (4) is then replaced by . This amounts to set Λ = −κC 1 as already stated and more importantly to add a new composition-dependent surface tension term κC 2 1 φ 2 (r)/2 in the Hamiltonian Equation (4). The inverse of the structure factor becomes S −1 (q) = α(T − T c ) + bq 2 + κC 2 1 /[1 + (qξ) 2 ] which leads to an additional effective interaction mediated by the elastic membrane between lipids of the same kind, where J 0 is the Bessel function of 0th order and K 0 is the modified one. Hence the interaction is repulsive due to higher membrane deformation (and therefore elastic cost) between similar lipids when they are closer. This interaction of bending nature has range ξ and an amplitude set by κC 2 1 and is screened at large distances r ξ due to the surface tension cost. It also modifies the critical temperature through T c = T c − κC 2 1 /α, because the parameter m in Equation (2) becomes m + κC 2 1 , as well as the value of q c according to [85] for C 1 > C * 1 . Hence in this case, for C 1 > 2C * 1 , i.e., for large spontaneous curvature differences, the domain size can be smaller than ξ. More importantly, this new term forbids the formation of mesophases by excluding the divergence of the structure factor at finite q, and only structured disordered phases can occur.
Using this curvature-induced mechanism, MacKintosh considered the alternative effect a quadratic confinement term Kh 2 (r) in the Hamiltonian altogether with a vanishing surface tension. This also leads to modulated phases [84]. However, the characteristic wave vector is different since the term σq 2 is replaced by a constant K in the Fourier space. The second maximum of the structure factor then obeys q 2 c ≥ κC 2 1 /(2b). This confinement can also be seen as a model for the cell wall pinning by proteins and/or tethering to the actin cortex [87,88] (the membrane tension was restored in [88], and the analytical treatment was performed in spherical geometry there). The Leibler-Andelman mechanism described above does not properly take into account the bilayer structure of the membrane. It can be fully relevant when some biomolecules adsorb only onto the external leaflet as in the model system by Shimobayashi et al. [86] where the ganglioside GM1 is added in the solution. However, although the average lipid composition of real biomembranes is in general asymmetric [79,83], both cases where the local lipid compositions of the two leaflets are correlated (registration) or anti-correlated (anti-registration) are in principle possible [85,89].
Hence a more detailed description of the bilayer requires the introduction of two composition fields, which correspond either directly to the compositions in one of the two types of lipids, relative to the average ones,φ 1,2 in the two monolayer (1) and (2), ψ 1,2 = φ 1,2 −φ 1,2 , or to the linear combination of them ψ ± = (ψ 1 ± ψ 2 )/2. This idea has first been proposed by McKintosh [84], and studied recently by Shlomovitz and Schick [79] for two different lipid types in the two leaflets and by Gueguen et al. [85] for two different average compositions of the same couple of lipids in the two leaflets. The phase diagrams become richer since the numbers of "masses" in the field-theory, m = α(T − T c ), is no more one but three, associated with the terms in m + ψ 2 + , m − ψ 2 − , and m 0 ψ + ψ − in the Hamiltonian (assuming that the coefficient b is the same in the two leaflets). In both cases, in addition to the modulated phase and the macrophase separation regions, the liquid phase (i.e., no order at long range) is now divided in 3 sub-regions: the true liquid one, and the two structured disordered regions corresponding to the two fields ψ + and ψ − (Figure 2a).

Page 7 of 13
obtain a saturaollowing Hirose isordered (SD) in the thermoh correlations in height fluctualation length on /2J) 1 /4 √ C 1 for nce between the large number of n between them zed by J, which tches (and thus l perimeter. (related to the e for m 0 = 0, can be checked ith previous res introduced by natural Hamilne introduces a q. (46), and forr at q * . , i.e. for differivergence at fihe formation of ensional struc-= ξq for ℓ + ≡ 0.3 ξ (note that − ) and for vari-We classify the sponds to the 4 C * 1 (red curve q = 0 and de-This behaviour ited by the blue lack line correseparation, i.e. stem follows a xhibits a diveres (green curve ̸ = 0. The SD− ̸ = 0 such that dered phase for stence of q * + ̸ = 0 e blue and green ), We have checked numerically that for m 0 > 0, (ξq * ) 2 ∝ Hence, as shown in fig. 2

RESULTS
The phase diagram of our model, as obtained below, depends directly on three intensive parameters: the temperature, T; the strength, G, of the coupling between curvature fluctuations and the composition of the inner leaf; and L, the strength of the coupling between compositions of the inner and outer leaves. This latter coupling has been estimated by May (42), utilizing theoretical arguments of Collins (43), to be such that 0.1 < L/nk B T < 1.0. We choose the lower limit and can display the diagram in the plane of temperature, T, and dimensionless coupling, G/(b f g) 1/2 . For convenience, we will discuss the system in which the average values of the composition differences, f and j, are zero. Because the relative amounts of SM and PC in the outer leaf of the plasma membrane are comparable, as are the relative amounts of PS and PE in the inner leaf, this is a particularly interesting situation to examine. The phase diagram is shown in Fig. 2 a.
There are four distinct phases: Fluid, liquid-ordered, liquid-disordered, and modulated.
For any coupling and at sufficiently high temperature, there is first a fluid phase. It is characterized by the vanishing of the ensemble average of all Fourier components of the order parameters. As the temperature is lowered, for small values of the coupling G, a critical point is encountered FIGURE 2 Phase diagrams of the system for three different couplings between the leaves. Each diagram exhibits four phases: For any coupling G, there is at sufficiently high temperatures a fluid phase; at lower temperatures and small G, the system separates into two fluid phases, liquid-Model Raft in Both Leaves of a Bilayer 1409 The interpretation is different in both cases. Shlomovitz and Schick fix some values to the species in each leaflet and using the Flory mean-field model, they construct a phase diagram in the (C 1 /C * 1 , T) plane, as displayed in Figure 2b. In particular, they study the case where both leaflets have different critical temperatures, that we denote by T in c and T out c . For example, in asymmetric lipid bilayers made of sphingomyelin, phosphatidylcholine in one leaflet, and phosphatidylserine and phosphatidylethanolamine in the other one, T in c ≈ 200 K and T out c ≈ 300 K. Furthermore, they assume that the coupling constant |Λ| = κ|C 1 | is on the order of 100 pN which implies that if the membrane bilayer is in the fluid phase, this phase is a micro-emulsion. They show that if T > T in c a structured liquid can appear on each leaflet: on the inner leaflet only for intermediate values of C 1 /C * 1 (region SDin in the figure), and on both leaflets for larger C 1 /C * 1 values. For even larger values, a mesophase (or modulated phase) appears. In particular in region SDin, the two leaflets are decoupled such that the upper leaflet still behaves as a ordinary liquid. The rafts, however, are identified by those authors as the transient curved domains that may appear in the region denoted SD in the figure, where both leaflets are in the structured liquid regime.
On the contrary, Gueguen et al. [85] construct a phase diagram at fixed temperature in the (m 0 , C 1 /C * 1 ) plane. In principle many parameters indeed depend on T in an ill-defined way, notably the masses but also the surface tension [65] and anticipating how the phase diagram behaves with T is quite challenging. Since the same mixture of lipids is assumed to be present on each leaflet, the critical temperature T c is the same. However, due to the difference in concentration, the demixing temperatures T d can be different. The parameter m 0 is related to these different average concentrations in the two bilayers. According to the Flory theory, one has For large C 1 a structured liquid regime appears with transient curved domains (phase denoted SD-) and if m 0 is increased one observes the occurrence of a structured liquid with both curved, and flat and thicker, transient domains (phase SD+). If Lo lipid rafts are assumed to be flat, the transient thick domains appearing in the SD+ region of the phase diagram are good candidates for being rafts.
Another mechanism has been proposed by Schmid and coworkers to explain the occurrence of small flat nanodomains [14,90,91]. In this model, a composition-curvature coupling similar to the one described above is introduced for each leaflet composing the bilayer. An elastic cost is introduced when the bilayer has a width different from the equilibrium one, similarly to the hydrophobic mismatch penalty close to proteins inserted in the membrane [92]. When the two leaflets bend towards the exterior or the interior of the membrane, this elastic cost forbids the formation of large domains leading to small nanodomains instead of a macrophase separation. Schmid et al. argued that these domains are candidates for lipid rafts, due to their small size of a few nanometers, fixed by a balance between the bending modulus and the area compressibility [14].

Bending Modulus-Composition Coupling in Planar Membranes
Apart from a coupling between the spontaneous curvature and the composition, it has been proposed that for symmetric membranes a coupling through a composition-dependent bending modulus, κ(φ), should occur. For instance Lo phases are thicker than Ld ones, which leads to κ Lo > κ Ld as already discussed. The simplest way to deal with this fact is to interpolate linearly between both values as κ(φ) κ 0 + κ 1 (φ − φ 0 ). Hence a new term (κ 1 /2)φ(∇ 2 h) 2 enters the Hamiltonian of the field theory, beyond the Gaussian order. In addition to the above coupling in φ(∇ 2 h) 2 , additional terms can be introduced in φ(∇h) 2 or higher-order ones where φ is replaced by φ 2 or (∇φ) 2 .
On the other hand several works have studied the interaction between two rigid inclusions in a fluctuating membrane [92][93][94][95] and shown that it is attractive and proportional to k B T which means that it is a membrane fluctuation-mediated interaction (see the so-called "Casimir" forces in Section 2.3). It is natural to anticipate that the more rigid Lo domains will also interact in a similar way. By performing a cumulant expansion, Dean and Manghi [96] indeed showed that the resulting fluctuation-induced interaction between lipids is a Casimir-like attractive interaction V(r) ∝ −k B T/r 4 at intermediate distances. However, the interaction at short distances V(r) ∝ −k B Tκ 1 w (φ 0 )/(κ 2 0 r 2 ), where w(φ) has been defined in Equation (2), is larger than the van der Waals attractive interactions and can be either repulsive or attractive depending on the sign of κ 1 w (φ 0 ). In the case where κ 1 w (φ 0 ) < 0, mesophases can occur due to this repulsive fluctuation-induced interaction for large values of κ 1 and small values of w (φ 0 ). Note that it has been proposed using coarse-grained simulations that κ(φ), where φ is the area fraction of short amphiphiles, is not linear, as assumed above. It decreases when φ increases for φ < 0.6 but increases slowly for higher values of φ [97,98]. These numerical works show that −4 ≤ κ 1 /(k B T) ≤ 2 which suggests that the condition κ 1 w (φ 0 ) < 0 can be reached in real systems.

Vesicles
We begin with an homogeneous vesicle. All the above theories where developed for planar membranes for which the average curvature vanishes, d 2 r ∇ 2 h = 0, a constraint that is relaxed in spherical geometry. It is important to note that the vesicle can be studied either at constant volume V or at constant pressure jump between the interior and the exterior of the vesicle, ∆p = p in − p out . This pressure jump can be an osmotic one due to solutes to which the membrane is impermeable, or from a theoretical perspective, a Lagrange multiplier to enforce the volume constraint. This pressure jump fixes therefore the projected area A s = 4πR 2 of the sphere having the same volume, related to R by V = 4πR 3 /3, and plays exactly the same role as the frame tension for planar membranes [65]. At the Gaussian order and for large surface tensions, σ κ/R 2 , both are connected through the Laplace law ∆p = 2σ/R.
In particular, this expression can be interpreted as follows: the limit σ → ∞, where thermal fluctuations are suppressed, is equivalent to the high pressure limit ∆p → ∞, itself equivalent to the limit where the reduced volume v ≡ 6 √ πV/A 3/2 goes to 1. In contrast, for low surface tensions, σ κ/R 2 , the vesicle is no more quasi-spherical and a local equilibrium shape condition has been derived by Zhong-can and Helfrich [99] which generalizes the Laplace law. Note however, that in the limit σ κ/R 2 , high order terms must be taken into account in the Hamiltonian, which renormalize the surface tension [65]. The field-theoretic framework presented above can be easily extended to the case of quasi-spherical vesicles made of a lipid mixture. The membrane fluctuations of single-component vesicles have been described by the Helfrich hamiltonian expressed in the spherical harmonics basis (Y m ) [100][101][102][103]. Taniguchi et al. [104] investigated the equilibrium shapes of two-component vesicles using the same coupling as in Equation (4) when the temperature T is just below T c (weak segregation limit), and at fixed pressure jump ∆p. By performing a linear stability analysis and using the single-mode approximation, they obtained phase diagrams focusing on the role played by the domain boundary energy (parameter b) which favors fewer domains and the coupling energy proportional to Λ which favors high-mode states. They found a large number of domains for = 3, 4, 5... (the size of which is on the order of the vesicle one) for very small values of b, even for T > T c (see Figure 3a). This mechanism is exactly the same as the Leibler-Andelman one described in Section 2.1.1 for planar membranes.   (a) AT high T or low C 1 , the liquid phase is homogeneous; numbered regions correspond to modulated phases with the number of the most stable -mode ( = 1 to 5 here; = 1 corresponds to the macrophase separation). The superscript corresponds to the sign of C 1 . The bending modulus κ does not depend on the order parameter φ in the model. In this case, C 0 = ∆p = 0. Adapted from [104]. (b) Sketch of a quasi-spherical vesicle with two different types of lipids inducing either thicker (with a larger bending modulus κ 0 + κ 1 ) or curved patches with a local spontaneous curvature The case of symmetric bilayer with κ depending on φ has been studied numerically by Feigenson et al. for T < T c (see Section 2.2.2 below). In an attempt to unify both types of approaches, Gueguen et al. studied on an equal footing the cases where both the spontaneous curvature C sp and the bending modulus κ depend on the composition φ [85], for quasi-spherical vesicles at fixed surface tension σ. When the vesicle is made of two types of lipids of different sizes, C sp increases when the composition is locally non-symmetric between the two leaflets whereas κ increases when the longer lipids are in register (see Figure 3b). By again choosing a linear interpolation, C sp = C 0 + C 1 (ψ out − ψ int )/2 and κ = κ 0 + κ 1 (ψ out + ψ int )/2, phases diagrams in the plane (C 1 , κ 1 ) are constructed at Gaussian order (see Figure 3c). Even for C 1 = 0 (no curvature-induced mechanism), a structured liquid phase occurs (SD+) with transient highly correlated flat raft-like domains of larger bending rigidity, whereas no transient curved domains appear. They appear for high enough values of if m 0 = 0. If m 0 = 0, a more complex formula is given in reference [85] (see Figure 3c). Again, the structured disordered phase occurs when the elastic cost associated with difference κ 1 in bending moduli is larger than the interfacial energy ∝ b.
This regime is experimentally accessible given the measured values of κ Lo /κ Ld already discussed. This effect cannot be seen in planar membranes. Indeed, κ * 1 is finite provided that C 0 is different from the average curvature 2/R imposed by the vesicle volume, whereas C 0 R = 2 is imposed in the large R limit corresponding to the planar case [85]. For non-vanishing C 1 , a large region with structured phases emerges in the phase diagram, as soon as C 1 C 0 and κ 1 3κ 0 (see Figure 3c).

Strong Segregation Limit
Well below the demixing temperature, appealing to field-theoretic techniques is more involved because it is necessary to take into account higher order φ 4 terms in the Hamiltonian. Alternative effective theories have been proposed in the literature, the most popular of which are now presented. They use the fact that in the strong segregation limit, the large line tension prevents the boundary B Lo−Ld between both phases to fluctuate much at short wavelengths. The phase-boundary energy can be written as its length |B Lo−Ld |, which is well defined, times the line tension, as we already explained it in the introduction of this Section 2 [45]: In the strong segregation limit, the typical values of λ are measured to be on the order of a few pN [25,33,[105][106][107].
This idea was first fully exploited in 1992 by Andelman and co-authors [108,109], in the case where the composition field φ and the curvature ∇ 2 h are coupled as in Equation (4), through the coupling constant Λ. In simplified geometries, i.e., striped planar membranes, and cylindrical or axially-symmetric vesicles, they could show that modulated mesophases are more stable than a macrophase separation in some favorable circumstances. As compared to the condition b < Λ 2 /σ for mesophase separation discussed below Equation (5), being far from the critical point leads to a different discussion. In the planar case, the emerging wavelength is i.e., in the low surface tension limit. In the opposite limit, d * diverges logarithmically when the membrane tension σ approaches the limiting tension signaling the onset of macro-phase separation.
In the particular vesicle geometries tackled in this preliminary work, the role of the pressure difference ∆p across the membrane was explored. At fixed membrane area, if ∆p = 0, then the volume is free to fluctuate and the vesicle adopts the most favorable geometry where mono-phasic domains are spherical caps with their prefered curvature. If ∆p = 0, numerical integration of the ordinary differential equation governing the vesicle shape shows how the vesicle is deformed by pressure. When the interior pressure is larger than the exterior one (∆p > 0 with our notations), the number of domains grows with ∆p in spite of the growing interfacial energy in Equation (16) because this allows the enclosed volume V to be larger as compared to a macrophase separation.
A rich literature has followed this pioneering work to account for mesophase separation in the context of coexisting Lo and Ld phase separated by a marked boundary.

Domain Buckling Induced by Line Tension-Spontaneous Symmetry Breaking
As illustrated in Figure 4, the first mechanism that we present is based upon spontaneous up-down symmetry breaking induced by the strong line tension λ. We first focus on planar membranes before dealing with closed vesicles.

1.
Planar membranes-We begin with the simplest form of this mechanism, as proposed in planar geometry by Lipowsky in 1992 [110]. We consider a single membrane Lo domain (denoted β-phase in this work) in a large planar Ld membrane (α-phase). Well below the demixing temperature, the boundary shape is close to a circle to minimize the interfacial energy. The total domain area is denoted by A β ≡ πL 2 (L is its radius in the membrane plane). Lipowsky first assumes that the surface tension is vanishingly small (σ = 0). If it buds in the third dimension, the domain adopts the shape of a spherical cap supported by a sphere of radius R, while the surrounding membrane remains flat ( Figure 4a). The interface is now a circle of radius N ≤ L. Adopting a mechanical approach where fluctuations are ignored, the total elastic energy E el of the domain is given by the sum of two antagonist contributions: the boundary line-energy E bound = 2πNλ that is proportional to the domain boundary length and tends to minimize it (by protruding in the third dimension); and the elastic Helfrich energy H Helfrich which disfavors bending. For a fixed domain area A β the value of the cap radius R is obtained by minimizing E el (R). A natural length-scale ξ I = κ/λ can be introduced, called the "invagination length". If κ β ∼ 100k B T for the Lo phase [106] and λ ≈ 1 pN far from the critical point, then ξ I ≈ 400 nm. When getting closer to the critical point, λ decreases as (T c − T) ν with ν a universal critical exponent equal to 1 in 2D biphasic systems in the 2D Ising universality class [15,110] and ξ I grows. We shall come back to these values later in the Discussion Section.
Lipowsky shows that if L < 4ξ I , then the optimal geometry is a flat domain (N = L); conversely, if L > 4ξ I , it is a complete sphere (N = 0), protruding upward or downward with equal probability. Differently said, this simple model without surface tension proposes that above a critical line tension spontaneously breaking the up-down symmetry is energetically favorable because it reduces to zero the interfacial energy cost in spite of the increased bending energy (Figure 4b).
The transition is first order. When C sp = 0, the symmetry between the two sides of the membrane is explicitly broken. If the value of L/ξ I is sufficiently small, the minimization of energy E el (R) predicts intermediate equilibrium values of N, between 0 and L. This situation for which the bending energy and the interfacial energy balance each other is called incomplete budding, or "dimpling".

2.
Additional role of surface tension-The case where the surface tension is finite, σ > 0, has been explored in detail in [107]. As bending stiffness, membrane tension applied in the membrane plane favors flat domains and comes in opposition to interfacial energy minimization. In this case also, and without necessarily appealing to spontaneous curvature, incomplete budding occurs above a critical line tension, through the spontaneous symmetry breaking principle ( Figure 4c). The transition from flat to dimpled domains is now continuous whereas it was discontinuous without tension. More quantitatively, it is proven in this work that the critical line tension is given by λ c 8κ β /L in the limiting case where the domain area A β = πL 2 κ β /σ. Here κ β is the domain stiffness, which can be different from the surrounding membrane one, κ α . Coming back to the notations used in the paragraph just above, this condition reads L = 8ξ I at the critical point, which is twice the transition value found when σ = 0. This means that in the interval 4ξ I ≤ L ≤ 8ξ I , budding is energetically favored when σ = 0, but becomes less stable than the flat geometry as soon as σ is positive, even if small. Furthermore, just above this critical value, the contact angle at the domain boundary scales as | | ∝ √ λ/λ c − 1. The Lo domain continuously but rapidly deviates from the flat state. By up-down symmetry, the domain is equally likely to bud upward ( > 0) or downward ( < 0). When C sp = 0, this symmetry is again explicitly broken, and phase diagrams can also be inferred [107].

3.
Vesicles-Jülicher and Lipowsky addressed the same question in the case of biphasic vesicles with spherical topology [111]. As above, different situations exist, but the up-down symmetry (more precisely the exterior/interior symmetry in this case) is explicitly broken on a vesicle. As stressed in the field-theoretic approaches presented in Section 2.1.3, a new ingredient can come into play here, namely the conservation of the volume V enclosed by the vesicle, or equivalently the pressure jump across the membrane, ∆p, which can be controlled through the osmotic pressure difference. The control parameter is, e.g., the reduced volume v ≡ 6 √ πV/A 3/2 ≤ 1, measuring the deviation to a sphere (for which v = 1). If two domains coexist as above, describing the membrane through an elastic continuum theory, the minimization of the total energy provides the so-called "shape equations", from which the equilibrium vesicle shape under the relevant constraints is derived. In particular, it depends on the relative area fractions and on the different parameters (bending moduli κ, saddle-splay moduli κ G , spontaneous curvature C sp ), which can in principle be different in the two phases. Indeed, even though it is not a pre-requisite, the difference between the bending moduli of the two phases now likely plays a role, contrary to the planar case, because both phases are bent in this geometry.
A rich phase diagram can be computed by minimizing the membrane energy, still neglecting thermal fluctuations. In this case as well, budding can be incomplete or complete, a closed vesicle then being connected to the main vesicle through an infinitesimal "neck". However, a strong volume constraint v 1 (or equivalently large σ), where the shape is quasi-spherical, can act against the budding process but does not, in general, suppress it. The reader can refer to [111] for further details. These results have been confirmed by numerical coarse-grained modeling (4-bead lipids and explicit solvent) based on dissipative particle dynamics, where both area and volume are conserved [112].

4.
Experiments-Fluorescence microscopy experiments [105,106,113] have later validated this theoretical approach on free-floating giant unilamellar vesicles (GUV) made of ternary mixtures of saturated lipids, unsaturated lipids and cholesterol, well below the demixing temperature, which display separated Lo and Ld phases ( Figure 4d). In reference [114], the reduced volume v of GUVs made of a DPPC/DOPC/cholesterol mixture is controlled by varying the osmotic pressure. If one starts from a spherical vesicle, domains bud (inward or outward according to the experimental conditions) when the enclosed volume decreases. Following these original studies, a series of papers studied the experimental counterpart of these theoretically predicted circular, budded Lo domains and established phase diagrams [63,107,115]. When the cholesterol concentration was increased above ≈ 35%, a reversal phenomenon was observed, now with Ld domains in a Lo continuous background. The domain sizes were typically observed to be in the micron scale. We have previously explained that if σ = 0 [105], then the critical radius L above which domains buckles is 8ξ I ∼ few µm with the above value of ξ I ≈ 400 nm. Experiments and theory are compatible. Even though in a less evident manner, AFM experiments also suggest that budding exists in planar geometry [25], as predicted by theoretical approaches in the relevant regimes of parameters.

5.
Elastic interaction between budded domains-In these experiments, it is also observed that domains sometimes coalesce [105,115] but that this process is very slow and does not follow the usual laws of coarsening [116]. The reason is that budded domains repel each other when they come in close proximity because they deform the elastic membrane, in an enhanced way if they are very close. This repulsion has even be very well quantified experimentally [63,107,115] and shown to be compatible with theoretical predictions. A supposedly metastable configuration is then observed with long but finite lifetime. After several hours, all Lo domains eventually coalesce and one ends with a complete macro-phase separation. Note also that coarsening is not always trapped and that the existence of normal coarsening has been correlated to a vesicle reduced volume v very close to 1 [115]. Indeed, budding requires excess area that is only available if the vesicle is at least slightly deflated.
As a matter of fact, the complete proposed scenario is as follows: after quenching below the demixing temperature and once domain have nucleated, normal coarsening is initiated, with small but growing nanoscopic domains. Being small, these domains are flat as demonstrated above [63]. When their size reaches the critical value, all these domains suddenly buckle and coarsening is then trapped in the metastable state [107,115]. The lateral organization of domains observed on phase-separated Sphingomyelin(SM)/DOPC/cholesterol vesicles in [117] has been attributed to this inter-domain repulsion, and the force between domain measured. Strong slowing-down of domain coarsening observed in DPPC/DOPC/cholesterol GUVs [118] was also attributed to budding, even though the inter-bud repulsion was not explicitly appealed to in this work. In contrast, when budding is avoided on sufficiently taut vesicles, no slowing-down is observed with respect to the expected dynamical exponent [119].
We have presented here the membrane buckling mechanism due to the high line tension λ at the domain boundary well below the phase-transition temperature. It does not require any spontaneous curvature C sp of the Lo phase nor any difference of bending rigidities κ between both phases. However, cell plasma membranes are believed to be relatively close to the phase transition temperature, for which the value of λ is lower than the pN order of magnitude considered above, and the lower limit of the domain radius L required to have buckling is thus significantly larger than 1 µm. Below this radius, domains coarsen and L continuously grows with time. It is thus unlikely that the main driving mechanism leading to the observed nanodomains in cells is buckling. We shall come back to this issue in the discussion section.

Competing Interactions: Phase-Dependent Bending Modulus
Bending elastic moduli of lipid bilayers κ are known to depend on the lipid species and on the phase state (Lo or Ld) as already discussed [120]. When peptides or proteins are included in the bilayer, the effective value of κ also significantly depends on their concentration [61,121,122]. In the vesicle geometry, the bending elastic energy necessarily plays a role since the membrane is always globally curved. The way the different membrane constituents spatially organize in the membrane plane can be more or less favorable with respect to this bending constraint.
In two contemporary series of papers, Lipowsky and collaborators on the one hand, and Feigenson and collaborators on the other hand, have investigated in detail the role of the different bending elastic moduli between two separated phases, e.g., the Lo and Ld ones, and how they can lead to mesophase separation. These investigations combine experiments on GUVs, analytical considerations and numerical simulations ( Figure 5).
relatively rigid domains cannot bend easily to shorten the domain boundaries, the increase in line tension energy cannot be overcompensated by a reduction in bending energy. One might expect that the bending energy of the majority b phase could be reduced by having more a domains, but this reduction in b bending energy is not sufficient in the absence of spontaneous curvature. In Sec. 3.3, we will show that the introduction of spontaneous curvature can stabilize such morphologies with multiple a domains.

Effect of volume constraint
The volume of lipid vesicles is determined by osmotically active particles in the aqueous solution. Indeed, these vesicles adapt their volumes in such a way that the osmotic pressure inside the vesicles is balanced by the osmotic pressure in the exterior solution. In Fig. 4(a) and 4(b), we display morphology diagrams of vesicles with relatively large bending rigidity ratios as a function of membrane composition and reduced volume v as defined in eqn (27).
For v ¼ 1, the vesicles attain spherical shapes and exhibit the morphology I1, which minimizes the line tension energy. The reduced volume v represents a useful control parameter since it can be easily changed by osmotic deflation or inflation of the vesicles. If one starts with a spherical vesicle and decreases the reduced volume v for small values of the composition variable c (a) , the vesicle undergoes the morphological transition I1 / II1, see Fig. 4(a) and 4(b). During this transition, both the rigid a domain and the flexible b domain change their topologies: the large a domain splits up into two smaller a domains, which form two weakly curved caps, whereas the flexible b domain changes its topology from a cap to a narrow belt along the ''equator'' of the vesicle, see Fig. 3(a). For larger values of the composition variable c (a) , the vesicle attains additional morphologies such as I3 and I4 with more than two flexible b domains. It is interesting to note that, in the absence of a volume constraint, the corresponding parameter regions of the morphology diagram are mainly occupied by II1, compare Fig. 2(a).
The multi-domain morphologies I3 and I4 have already been observed for lipid vesicles, see Fig. 1  to the previously measured tielines in the liquid-liquid coexistence region.
Modulated phases occur at similar r-values along a given tieline The first tieline examined (tieline 1) is located slightly above the bottom phase boundary of the Ld þ Lo region (Fig. 2). Six starting compositions with defined DSPC/ CHOL ratios were chosen along tieline 1 to investigate the modulated phase compositional window (compositions T1A-T1F; Table 1). In all of the chosen compositions, the majority of GUVs appeared uniform at r < 15%, where domains are nanoscopic (22). Modulated phase patterns started to occur at r ¼ 15% for composition T1A, located on the far right of tieline 1 (Fig. 3 A). At r ¼ 20%, patterned domains were observed in compositions T1A-T1D (Figs. 2 and 3 A). For compositions T1E and T1F, this modulated phase window began at r ¼~25%. The fraction of patterned GUVs in compositions T1E and T1F, where the predominant phase had changed to Ld-rich, was smaller than that found for compositions T1A-T1D (Fig. 3 A), where the predominant phase was Lo-rich. In T1F, only~25% of the GUVs analyzed displayed modulated morphologies; the other GUVs were either uniform or displayed large round domains (Fig. 3 A, Table S1). For all of the compositions examined on tieline 1, the modulated phase window ended at r ¼~30-35% (Fig. 3 A), with higher r-values showing rounded macroscopic domains. A typical progression of GUV morphologies from the onset to the end of a modulated phase window is shown for composition T1B in Fig. S1.
We observed that the occurrence and behavior of modulated phases followed rules that govern phase separation along thermodynamic tielines. On the far right of tieline 1, Ld is the minor phase, forming thin stripes, honeycomb, or stripe-like patterns that resembled two-dimensional (2D) bubbles on the predominantly Lo phase surface of the GUVs (Fig. 2, I-L). As the Ld phase fraction increases, either by moving toward the left side of the tieline or by increasing r, the thin stripes start to coarsen (Fig. 2, C and G), typically giving rise at higher r-values to larger dispersed Lo domains with irregular edges (Fig. 2, E and  H; Fig. S1 F). The mole fractions of Ld and Lo on patterned GUVs follow the lever arm rule (see Fig. 6).    Lipowsky and his collaborators proposed two different approaches to clarify the role of two different bending moduli [53,64]. Buckling as discussed in Section 2.2.1 is not at play here because the line tension λ is below its critical value λ c . The Lo and Ld phases can have a non-zero spontaneous curvature C sp , but this is not a prerequisite for the competing-interactions mechanism under consideration here. As above, the reduced vesicle volume v can be either controlled by the osmotic pressure difference across the membrane, or free to fluctuate, thus assuming in both cases the membrane to be permeable to water [2]. In 2009, Gutlederer et al. [64] studied this problem by appealing to a numerical solution of the mechanical energy minimization problem, again ignoring thermal fluctuations and thus configurational entropy. A morphology diagram was obtained, with or without volume constraint. The configuration minimizing the interfacial energy at the phase boundary is the macrophase separation, with a single domain of each phase. However, when decreasing the line tension λ or increasing the κ Lo to κ Ld ratio or the Lo area fraction φ Lo ∈ [0, 1], it was found that macro-domains split into two or more smaller ones in order to deal with the bending energy extra cost. An additional volume constraint through the parameter v leads to additional morphological transitions of multi-domain vesicles, in particular oblate-to-prolate ones.
To study the role of thermal fluctuations in this context, Hu, Weikl and Lipowsky later performed Monte Carlo simulations on a model of dynamically triangulated biphasic vesicles below the demixing temperature T d [53] (Figure 5a), where each triangle is either Lo or Ld. Adjacent triangles interact through an Ising-like model. By construction, Monte Carlo simulations fully take thermal fluctuations into account. In thermodynamic equilibrium, the obtained morphology diagram is even richer than the previous one, with macro-domains splitting in up to 6 smaller ones, so that most curvature is absorbed by the more bendable Ld phase. Here also, constraining the reduced volume v can have noticeable effects. For v very close to 1, the vesicles are forced to adopt nearly spherical shapes with the same curvature for both phases. They display a macrophase separation to minimize the interfacial energy. As before, when decreasing v, the vesicles undergo morphological transitions to reduce the bending energy in spite of the increased interfacial cost. Some morphologies predicted by this approach are successfully compared to experimental images [124,125]. These morphologies have also been explored numerically by solving the two-fluid Stokes equation with a thermal noise in a two-dimensional curved surface [126]. It is observed that several shapes can be in competition and that the final (meta-)stable shape may depend on the system history.
The first paper of the series co-authored by Feigenson on the same topic [127], in 2011, examined the experimental phase-behavior of 4-component DSPC/DOPC/POPC/cholesterol lipid mixtures by fluorescence microscopy (POPC is 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine and DSPC is 1,2-distearoyl-sn-glycero-3-phosphocholine). A rich variety of patterning morphologies was observed on GUVs, comparable to what is shown in Figure 5c, such as roundish small domains, honeycomb arrangements, stripes, or labyrinths. The main finding of this work is that mesophases are observed provided that a fourth lipid is added to the mixture. The DOPC to DOPC+POPC molar ratio, denoted by ρ, is to be in a relatively narrow window to avoid either a macrophase separation or an homogeneous mixture (at least homogeneous as it appears above the Rayleigh resolution limit). This patterning propensity of this 4-component lipid mixture was explored further by the same group in Refs. [31,49,123,128], using both Monte Carlo simulations and performing a more systematic exploration of the parameter space. The reason why a fourth lipid species is required to observe mesophase separation in this context presumably comes from the lower line tension λ between the Lo and Ld phases. Line tension has been shown to increase with ρ because the thickness mismatch between Lo and Ld phases has been shown to grows linearly with ρ thanks to SANS experiments [31]. Thickness mismatch goes from 0.6 nm when ρ = 0 to 1.1 nm when ρ = 100 %. The line tension is then estimated to be a fraction of pN in the range of values of ρ where modulated phases are observed. When it grows beyond a limiting value, in the pN range, the interfacial energy is too large and the macrophase is more stable than the patterned mesophase in the GUVs of interest. In contrast, in the small 60 nm-vesicles of reference [31], a large domain presumably remains unstable even when ρ grows because its bending energy remains large as compared to several smaller domains. A moderate line tension also forbids buckling as above.
The simulations [123] are based upon a model comparable to the one discussed just above, with the vesicle discretized in small triangles. In this case however, the concentration field is defined on the lattice vertices, contrary to Lipowsky's model where it is defined on the triangular faces. The system is also simulated below T d (Figure 5b). The initial vesicle radius is taken to be R = 25 µm and the volume is not constrained, but non-concave shapes are forbidden. Even though the bending moduli κ are ∼ 10 times higher than measured ones (up to 5000 k B T for κ Lo ), the model reproduces well the variety of morphologies observed in experiments. The line tension is also relatively weak, λ ∼ 10 −2 pN, which indicates in principle the proximity of a critical point [18]. Moderately increasing the line tension makes domains larger. Values of λ in the pN range again lead to macro-phase separation instead of meso-patterning. This work also demonstrates that increasing both the vesicle radius R and the line tension λ while keeping R/λ constant keeps morphology unchanged. In a subsequent experimental work [49] complementing Ref. [127], the parameter space is explored following tie-lines in the coexistence region of the phase space (a tie-line is such that the composition of each phase, Lo or Ld, is fixed along it [17]). The ranges of the DOPC to DOPC+POPC ratio ρ for which a modulated phase can be observed have been characterized thoroughly along the tie-lines. The agreement between experimental and numerical morphologies is confirmed in this study, at least on a qualitative level. This series is continued by a fourth paper [128] where a new ingredient is added in the simulations, namely dipole-dipole repulsion between lipids with a 1 to 3 nm range depending on the solution ionic strength. This repulsion was not included so far because the membrane was coarse-grained at a larger, micrometric length-scale. The dipole density is assumed to be about twice as large in the Lo phase as in the Ld one. This model belongs to the competing-interactions class because the repulsion is antagonist to the lipid-lipid attraction concretized by the finite line tension below T d [45]. Small domains of a size > 10 nm are observed in the simulations. As anticipated, the higher the dipole density in the Lo phase, the smaller the domains. Vesicle modeling embracing both length-scales (nanoscopic and micrometric) remains to be done because this study was limited to patches of 60 nm due to the necessary low level of coarse-graining (see the Section 4.3 for further discussion on this topic). This work also studied how the line tension is renormalized in the Ising universality class when changing the coarse-graining level, but the coupling to membrane fluctuations is still lacking in this respect.
Coming back to the situation without dipole-dipole repulsion, maybe the main finding of this work [128] is that "background curvature is necessary for the stability of patterns", thus pointing out the explicit up-down symmetry breaking mechanism. Considering membrane patches with the same simulations parameters except an increasing radius of curvature R, it is observed that first domains broaden and eventually disappear to the profit of a macro-phase separation (when R 50 µm), as seen by inspection of Figure 6. Indeed, as in Lipowsky and coworkers' studies cited above, bending the most rigid Lo phase requires elastic energy. Missing spontaneous curvature, this rigid Lo "carapace" must be broken up into smaller fragments in order to conform to the curved spherical shape. The interstice between the Lo domains is filled by the more bendable Ld phase, which minimizes the elastic energy. This theoretical conclusion seems to be confirmed by AFM experiments on a 4-component lipid mixture containing both POPC and DOPC, as above (DSPC has been replaced by brain-sphyngomyelin (bSM)) [29]. This systematic study points out the absence of "isolated disk-like domains" in flat geometry (where R → ∞) resulting from phase-dependent bending moduli, even at the 10 nm length-scale. Figure 6. From right to left, the vesicle radius of curvature (written above each snapshot) is increased while the total area of the patch is held fixed. Above a critical curvature, there is a transition from a Lo/Ld mesophase to a macrophase separation. Reproduced from [128], with permission from the American Physical Society, Copyright 2014.
To confirm the decisive role played by line tension in the coexistence region below the demixing temperature on GUVs, a very recent work by Feigenson and coworkers [33]

explored a large variety of 4-component mixtures of the (high-T m PC or sphingomyelin)/low-T m PC(s)/cholesterol type.
In particular, the bSM/DOPC/POPC/cholesterol mixture of [29], as just mentioned above, is also included in this study and does display mesophases on GUVs, thus confirming our previous statement. This work also refines further the parameter range, notably the value of ρ, where mesophase separation is visible by optical microscopy (with respect to either homogeneous phases or macrophase separation). The line tension of the Lo/Ld interface was measured using the flickering spectroscopy method and confirmed to dwell in the 1 pN range. More precisely, patterned morphologies were observed for moderate values of the line tension 0.3 λ 1.1 pN. The possibility that thermal fluctuations break up domains below 0.3 pN is evoked. Above ≈ 1.1 pN, macrophase separation in large round domains is found to be favored.

Competing Interactions: Spontaneous Curvature Induced by Membrane Inclusions
The approaches in Section 2.2 were focussed on the competition between line tension on the one hand, and bending elastic energy, through the dependence of the bending rigidity κ (and possibly the saddle-splay modulus κ G ) upon the local membrane local composition φ on the other hand. As already discussed in Section 2.1, another local parameter likely depends upon the local membrane phase state or local composition, namely the spontaneous curvature C sp . As seen in Section 2.1.1, several mechanisms can lead to local spontaneous curvature C sp = 0 of lipid bilayers. As far as membrane inclusions are concerned, their molecular conformation, for example, can also confer them a conical shape that locally bends the membrane and thus plays a similar role. We shall mainly insist on this second source of spontaneous curvature in the present section. This section is not included in the previous strong segregation limit one because it deals with membrane inclusions such as proteins instead of lipid phases. However, the effective theories proposed in this context often consider protein domains with little fluctuating boundaries, i.e., large line tension, and connections will be made with continuous theories as above.
The basic mechanism is that the membrane can deal with the curvature constraint by splitting the macro-phase into small domains, such that regions with a strong spontaneous curvature C sp are dispersed throughout the membrane. This will be shown to lower the overall bending energy.

Inter-Protein Short-Range Forces
Before dealing with the role of spontaneous curvature in the stabilization of mesophases, we need to briefly review the origins, ranges and intensities of inter-protein attractive forces that lead to their phase-separation in biomembranes, and the interplay with lipidic phase-separation. Indeed, the fact that they are inserted in (or adsorbed on) the membrane is the source of effective interactions (also called "potential of mean force") propagated by the elastic and fluctuating membrane. They can be either of mechanical origin, or of entropic nature, resulting from the important thermal fluctuations of the membrane. Many works have calculated these forces propagated by the membrane, for example the pioneering works of references [93,[129][130][131][132] and many of those cited below, using the Helfrich free energy to describe the membrane [55] and different kinds of boundary conditions imposed by the inclusions.

1.
Electrostatic, van der Waals and hydrogen-bond interactions-Polar and charged amino acids at their surface can interact when two proteins come in close proximity. The Debye length λ D ∼ 1 nm in water at physiological salt condition [2] sets the typical range above which these interactions are screened. Inside the apolar hydrophobic membrane region where the dielectric constant is weaker, the range can be somewhat larger, of a few nanometers [133,134]. The range of van der Waals and hydrogen-bond interactions is also nanometric.

2.
Hydrophobic mismatch-Integral proteins have transmembrane domains consisting of alpha helices with hydrophobic amino-acid side chains, buried inside the hydrophobic core of the lipid membrane. Protein and membrane hydrophobic core thicknesses do not necessarily match.
Since exposure of hydrophobic residues to the aqueous solvent is energetically unfavorable, the membrane must be deformed in the case of significant mismatch [12,135]. If two (or more) proteins are in proximity, the overall energy penalty depends on their distance d. As above, an effective force ensues (Figure 7b). It is attractive when both mismatches have the same sign and repulsive in the converse case. The energies at play go from a fraction of k B T to several k B T, depending on the degree of hydrophobic mismatch, and the range of these forces is few nanometers [136][137][138][139]. It has been suggested that hydrophobic mismatch forces are not pairwise additive [140].

3.
Casimir interaction-This attractive interaction, of entropic origin, is named by extension of the Casimir interaction in quantum physics (the attraction between conducting plates mediated by quantum fluctuations in the electromagnetic field). Here it results from the transverse thermal fluctuations of the elastic membrane. The number of vibrational degrees of liberty of a membrane in which two (or more) inclusions are embedded depends on their mutual distance d.
The potential of mean force thus depends on d, and has been shown to behave as ∼ −k B T/d 4 in the case of vanishing membrane tension σ [131,141]. The calculation can be extended to the case σ > 0, where the interaction energy decays much faster with d, as ∼ −k B T/d 8 when d ξ, and as ∼ −k B T exp(−d/ξ) when d ξ [142,143].

4.
Depletion (or excluded-volume) forces-Attractive depletion forces (Figure 7a) are well characterized in soft condensed matter when large particles evolve among smaller ones, and play a role in physical biology (see [2] for example). In the present case, they are due to the 2D osmotic pressure laterally exerted by the surrounding lipids on large transmembrane proteins (larger than lipids). It should be far less pronounced for peptides. Roughly speaking, when two proteins are far away, the lateral osmotic pressure is isotropic and no net force ensues. When the relative distance becomes on the order of the lipid lateral size (<1 nm), the interval between the two inclusions tends to be depleted in lipids, and the pressure is not isotropic anymore. This tends to bring proteins closer when they are about a nanometer away [144,145]. The ensuing binding energy is on the k B T range, even though the actual value depends on the model details.

5.
Lipid wetting-Some lipids are known to have a preferential affinity for given proteins species [39,[146][147][148], in particular but not exclusively because they better match their hydrophobic length. Even above the phase-transition temperature, the protein can nucleates a small "halo" of such lipids, the range of which is on the order of magnitude of the composition correlation length ξ OZ (see Section 2). This mechanism known as "wetting" [132,149,150] is reminiscent of the "lipid annulus" or "lipid shell" concepts that have become popular in the biophysical literature a dozen of years ago [40]. When two proteins approach close enough for their halos to overlap, they tend to assemble because it reduces the net interfacial energy. An effective attractive force ensues (Figure 7c). This nucleation mechanism can also promote the formation of a lipid halo of a thermodynamic phase that would be unstable in absence of the inclusion. A similar mechanism has been demonstrated to emerge in a very illustrative way [151]. In all cases, the range is set by the correlation length ξ OZ .
This force is enhanced near a miscibility critical point because the composition correlation length ξ OZ grows significantly. Exactly at the critical point, a long-range, power-law decrease of the potential of mean force at large inter-inclusion distance d has been predicted by a conformal field theory approach, with exponent −1/4, and confirmed by Monte Carlo simulations of the Ising model [152]. Coarse-grained molecular dynamics simulations on a model membrane and a phenomenological Ginzburg-Landau theory have explored the same mechanism in the case of peripheral proteins adsorbed onto the bilayer and interacting preferentially with one lipid species (among two). They drawn similar conclusions [147]. The binding energy at close range for two identical particles is also found in the k B T range.
Note that this mechanism is specific to the protein species and the lipids with which it preferentially interacts because the halos must be miscible if the interaction is attractive. In the case where they are immiscible, the force can even become repulsive instead [147,152]. Small alterations in lipid chemical structure can thus lead to dramatic changes in the membrane organization. This mechanism has been evidenced in model membranes [146]. The majority of these forces mediated by the lipid membrane are attractive, notably the ubiquitous Casimir and depletion ones, which explains why proteins are found in small clusters in bio-membranes. However, the variability of protein structures and the combination with possibly repulsive forces likely tune the strength of the attraction in function of the interacting proteins. It will be used in point Section 2.3.5 below. Since the intensity of each force is on the order of k B T, the typical binding energies att between membrane proteins are a few times larger than the thermal energy k B T, and are sufficient to drive protein macro-phase separation if no additional forces counter-balance them at long range. This order of magnitude has been confirmed by refined numerical simulations [137,147,153,154]. High-precision AFM experiments have been able to confirm the orders of magnitudes of range and intensity of these forces in the case of ATP-synthase c-rings [155]. The range was found to be ≈ 3 nm and the binding energy att ≈ 3.5k B T.

The Cluster Phase Scenario
As introduced at the beginning of this Section 2, cluster phases are patterned structures resulting from a competition between strong short-range attraction and longer-range, weaker repulsion between proteins [67]. Thus this mechanism also belongs to the competing-interactions model class [45]. It is reminiscent of the micellization theory [11]: domains, named clusters here, have a non-extensive free energy because of the long-range repulsion.
To make the argument more precise, we write a functional form of the free energy F(k) of a cluster containing k proteins, supposed to be identical in a first step [156]: with f 0 , ρ 0 , χ 0 > 0 and α > 1. The first term of F(k) is its extensive part, where f 0 is the energy gained by a protein when it joins the cluster. It is proportional to the short-range protein-protein binding energy att , as discussed above. The second term accounts for the proteins at the cluster boundary, that are energetically disfavored as compared to those in the cluster bulk. The coefficient ρ 0 is thus proportional to the cluster line tension λ cl , itself related to att . Finally, the last term takes long-range repulsion into account in an effective way, by growing faster than the cluster size k. In the extreme case where each protein pair inside the cluster experiences the same repulsive energy rep , the total repulsive energy is rep k 2 /2, thus α = 2 [157]; χ 0 measures the intensity of the repulsion. This expression is the so-called "droplet approximation" with an additional non-extensive term (i.e., a term growing faster than k). Note that in Equation (19), F(k) has been written as a function of (k − 1) so that F(1) = 0, as desired. This is a technical choice that has no incidence on the statistical mechanics of clusters for which k * 1. Replacing (k − 1) by k and subtracting a constant to ensure that F(1) = 0 leads to the same conclusions. This term makes too large clusters unstable because they have an unfavorable energy cost, and leads to the observed mesophase separation [11], as demonstrated via Monte Carlo simulations [69], analytical calculations based on the principles of statistical mechanics in the limit where interactions between clusters are negligible (i.e., in the dilute limit) [156], or more phenomenological arguments [158,159].
As a result, above a protein threshold concentration φ c , the cluster size distribution is bimodal. A protein dilute "gas phase", constituted of protein monomers and rare dimers or trimers, coexists with large multimers, or clusters, with a typical protein number k * depending on the balance between attraction and repulsion. Just above φ c , one finds that k * ∝ (ρ 0 /χ 0 ) (α−1/2) −1 [156]. This scaling can be simply obtained by balancing the second and third terms of F(k) in Equation (19), accounting respectively for short-range attraction, via the line tension, and longer-range repulsion. The typical size k * of multimers, excluding monomers, is expected to grow only slowly with φ above φ c [156]. This is consistent with experimental observations on membrane protein clusters, both in vitro [160] and in vivo [28,161,162]. In contrast, in a wide range of parameters, it is also found that the average cluster size k , including monomers, is proportional to the protein concentration: k ∝ φ above φ c [69,156]. This is also consistent with soft matter experiments [67].
Note that in this scenario, the clusters are seen as dynamical entities, with a well-defined boundary, but that however, permanently exchange isolated molecules with the surrounding membrane (the above "gas phase") [43,69], because the binding energies at play are on the same order of magnitude as the thermal energy k B T. Experimental observations have confirmed this prediction [23], showing with the help of dual color fluorescence microscopy and particle tracking how a single protein enters a cluster, dwells in it for a duration of several seconds, and eventually leaves it to diffuse again freely in the cell membrane. This allows a rapid and efficient remodeling of clusters, in function of the evolution of inter-protein interaction parameters, themselves depending on protein conformational changes, e.g., upon receptor activation. For example, in Ref [163], it has been measured in living cells that the spontaneous curvature of some G-Protein Coupled Receptors (GPCR) significantly changes upon ligand binding. The close lipidic environment can also depend on the protein conformation [38]. Dynamical re-targeting to different clusters after activation can then be envisaged.

Spontaneous Curvature Can Play the Same Role as a Long-Range Repulsion
When two up-down asymmetric inclusions locally bending the membrane are present at a given distance d, they feel a mutual interaction propagated by the elastic membrane ( Figure 7d).
As compared to the Casimir attraction evoked in Section 2.1.2, it comes from the mechanical deformation of the membrane and not from thermal fluctuations [93,143,164,165]. The interaction is repulsive at long range if both inclusions have the same orientation, and attractive in the converse case where they are "head-to-tail". The potential of mean force falls off as d −4 at large d in the tensionless case, and exponentially as exp(−d/ξ) for a taut membrane (in fact Bessel functions as in Equation (10), but with a different prefactor), with a range set by the Helfrich correlation length ξ = √ κ/σ. The intensity of the interaction is proportional to θ 2 σ, where θ is the cone half-aperture angle as illustrated in Figure 7d. Note that for a protein of radius a in the membrane plane, θ and C sp are simply related through θ = a C sp .
When an assembly of such identical proteins dwell in close proximity, they collectively bend the membrane and provoke an invagination as shown in Figure 7e. The situation becomes more involved because the energy cannot be written as a sum of pairwise additive interactions anymore. This question has been studied in details in Ref. [143], by calculating exactly the elastic deformation shape and its energy cost. The main finding is that the elastic energy of a k-cluster can be written in the form of Equation (19), with α ≈ 1.5. It follows that the spontaneous curvature imposed by proteins in a cluster is formally equivalent to a long-range repulsive force when expressing the free-energy F(k) as in Equation (19). As a long-range repulsion would do it, it makes large cluster unstable because the membrane deformation elastic energy grows faster than k and overcomes the gain in terms of line tension if k is too large.
A simple argument has been proposed in [143] to account for this non-extensivity. In the low-surface tension regime σ κ/R 2 , where R is the vesicle radius, Helfrich's free energy is dominated by the curvature term. The tension term can be treated as a perturbation. The spontaneous curvature C sp is imposed in the domain by the asymmetric inclusions. Minimizing Helfrich's free energy results, at order 0 in the perturbation, in a parabolic membrane invagination inside the cluster h(r) = C sp r 2 /2 if the cluster is centered at the origin. Note that in the present approximation of small gradients of h, a paraboloid of revolution is identical to a spherical cap. If r c is the cluster radius, Helfrich's free energy inside the cluster can be written The first term is a correction to the second one because σ is small. This term is proportional to r 4 c and thus makes the total membrane free energy non-extensive, as required. Since the cluster number of proteins satisfies k ∝ r 2 c , we infer that α = 2 and χ 0 ∝ C 2 sp σ ∝ θ 2 σ in Equation (19). When σ takes intermediate values more compatible with real cell tensions, the exponent α appears to be slightly lesser than 2. This effective repulsion term is always positive, even for weak asymmetry and weak tension.
A connection can be made with Section 2.2, where buckling of lipid domains was considered by following a continuous theory approach. The study was extended to the case C sp = 0 by several authors [48,107]. They naturally concluded that spontaneous curvature breaks the up-down symmetry and that budding then occurs in a preferred direction, regardless of the line tension value. This is the same mechanism as above, even though tackled with a different mathematical approach. They quantified how the spherical cap is deformed in this case when σ grows. A similar approach was proposed in parallel in Refs. [166,167] in the simplified case where all domains have the same size and form a regular two-dimensional array. It was concluded that the modulated phase is the most stable provided that the membrane tension is smaller than the same limiting value σ cr as given in Equation (17) (when both phases have the same curvature rigidity κ). It was also found that the number of domains grows proportionally to the concentration φ [166], and consequently that the typical domain size does not depend significantly on φ, as already predicted by the cluster phase approach above. In principle, the membrane elastic energy can be calculated in function of the domain size and the statistical mechanics of nanodomains can be tackled, as in the above cluster phase scenario. This remains to be done.
To our knowledge, the best available experimental test of this mechanism in a controlled way was performed by Shimobayashi and coworkers in 2016, by inserting GM1 gangliosides in the outer leaflet of GUVs [86]. GM1 gangliosides are known to insert preferentially in the Lo phase and then to impose a strong local spontaneous curvature to the membrane [168]. Thus Lo domains are endowed with a local spontaneous curvature in this system, whereas Ld domains are not. First the ternary DOPC/DPPC/cholesterol GUVs were prepared without GM1. The temperature was set below T d . To avoid buckling and the ensuing trapping in a metastable state (see Section 2.2.1), the experimental protocol ensured that the vesicles were close to spherical. Consequently, the Lo and Ld phases were fully separated. Then GM1 was added in the solution and got inserted in the outer leaflet. After several hours, the macrophase separation disappeared to the benefit of a mesophase in almost all vesicles, thus proving experimentally that large domains with a spontaneous curvature are unstable. An interesting intermediary stripe morphology was observed before turning to roundish small curved domains, which were stable for several days. Higher concentration of GM1, thus higher spontaneous curvature lead to smaller domains, as expected.
In a cell plasma membrane, there is no reason why proteins should preferentially curve the membrane in one direction, inward or outward.
If one assumes that there are roughly as many inclusions imposing an inward spontaneous curvature as there are bending it outward, they conceivably cancel each other out and lead to a globally vanishing membrane spontaneous curvature. This simplistic argument has been refuted in reference [143]. It has been shown by using Helfrich free energy that mixing differently oriented inclusions in a same nanodomain is energetically unfavorable because it imposes geometric frustration on the height function h(r). Consequently, inward and outward bending inclusions phase-separate and the cluster phase mechanism then applies to each orientation, leading to distinct inward bent clusters and outward bent ones. Such two differently oriented clusters feel an attraction at long range but a strong repulsion at short range (when their separation becomes shorter than ξ) [143].
2.3.4. Sources of (Local) Spontaneous Curvature C sp = 0 We have listed at the beginning of Section 2.1.1 the variety of up-down-symmetry breaking mechanisms involving lipids. Proteins also come into play to locally curve the membrane in collaboration with lipids [71], as follows: 1.
The transmembrane part of an integral protein has no reason to be up-down symmetric, not least because the cytosolic and extracellular protein regions do not have the same biological function. This is either apparent in the molecular shape of transmembrane proteins or can be inferred from their behavior in biophysical experiments [71,163,[169][170][171][172][173][174]. However, it seems difficult to infer the spontaneous curvature from the sole molecular shape displaying up-down symmetry breaking, for the reasons that we discuss now.

2.
Peripheral proteins naturally break the up-down symmetry [70,175,176], to a degree that depends in particular on the depth of penetration of the hydrophobic domain of the protein into the bilayer [71]. Numerical evidence can be found for example in Ref. [147], where the small shoulder on the interaction potential at intermediate range indicates a weak repulsion. More generally, anchored molecules can play the same role, as it was non-ambigously demonstrated in reference [86] on experimental proofs. 3.
The coupling between lipid composition and protein wetting by lipids is also a potential source of local curvature if a protein recruits different lipids in the two leaflets, themselves promoting markedly differential local curvature of the two leaflets.
Of course, local spontaneous curvature likely results from the interplay of these different mechanisms, notably in a nanodomain where several lipid and protein species cohabit.

Diversity of Membrane Proteins and Biological Specialization of Clusters
The above theoretical developments deal with a single protein species. In really, a cell plasma membrane usually contains hundreds of different populations that ensure as many biological functions. Clusters are thought to be specialized on the biological level, by gathering one or a few protein species out of the hundreds ones present in the membrane in order to perform a determined biological function. In spite of their structural similarity, several membrane proteins have even been observed to segregate in separate clusters [161,[177][178][179][180][181][182], thus enforcing the idea of an efficient, fine-tuned sorting of membrane proteins in function of their mutual affinities. From a statistical mechanics perspective, the more natural way to account for these observations is to suppose that fine-tuning of protein-protein short-range attractions enables their targeting to different, specialized clusters. In this respect, a mechanism for specialization of clusters as been proposed as follows [183,184]. We have seen above that even though some short-range forces are generic to (almost) all membrane proteins, such as Casimir or depletion forces, most of them depend on the protein specificities (molecular shape, width of the hydrophobic core, surface charge density, physico-chemical affinity for specific lipids). Consequently, a very wide range of protein-specific modulation of the short-range affinity att can be confidently envisaged [179,185,186].
The concept of protein "family" can be introduced, a family regrouping the protein species that tend to co-localize in the same specialized clusters (e.g., a receptor, a G-protein and an effector) [9,187]. We denote by q the number of such families, by att,= the binding energy (or affinity) between proteins of the same family and by att, = the binding energy between proteins of two different family. We introduce the difference of affinities ∆ = | att, = − att,= | that indicates the degree of preference that proteins have for partners of their own family. Appealing to a mean-field Flory-Huggins theory [188], one defines a new modified droplet approximation [183]. For example, if q = 2, where k 1 and k 2 are respectively the numbers of family-1 and family-2 proteins in the cluster, k = k 1 + k 2 , x = k 1 /k and ν is the average number of neighbors of protein in the cluster bulk (ν 6 in 2D). The last terms account for the mixing entropy. The Flory-Huggins theory states that there exists a critical value ∆ c below which the two families are mixed in the clusters. In contrast, if ∆ > ∆ c , a majority color emerges in each cluster, as desired. This argument can be extended to any value of q and it has been demonstrated that the critical value becomes exact in the large-q limit in spite of the mean-field approximation: This value has been confirmed by Monte Carlo simulations up to q = 5 [184]. The main interest of this result is that the logarithm grows very slowly with q. Even if q = 10 3 or even 10 4 different protein families cohabit in the same cell plasma membrane, ∆ c does not exceed 3 k B T. We have seen in the previous paragraphs that such a modulation of the short-range affinity is easily attainable thanks to the different forces propagated by the membrane.
We have thus shown that a "smart" mechanism for nano-sorting and cluster specialization is at hand without appealing to complex active processes. Of course, future studies will have to deal with the variability of the parameters att,= and att, = with the proteins species at play, whereas we have considered so far that each of them can take only one value.

A Unifying Rationale: Up-Down Symmetry Breaking
In Sections 2.2 and 2.3, we have examined into detail the three main mechanisms leading to the splitting up of the macroscopic Lo phase below T d in order to balance the elastic energy and the interface length. Note that the configurational entropy is also increased when splitting the macrophase, which should in principle be included in more refined theories.
Without any source of explicit up-down symmetry breaking, a strong line tension λ well below the phase-separation temperature is sufficient to promote spontaneous symmetry breaking leading to two equally probable buckled configurations of the Lo domains, upward or downward, the precise shape of which depends on the applied surface tension σ. In the planar topology, the equilibrium configuration remains the complete macrophase separation, but metastable configurations with long lifetimes can emerge. They are constituted of budded domains that feel mutual repulsion mediated by the elastic membrane. In the spherical topology of GUVs, configurations with several domains instead of a single macrophase have been observed, both experimentally and numerically.
We have also examined two mechanisms where the up-down symmetry is explicitly broken, theoretically predicted and experimentally observed. On the one hand, even in the case of an up-down symmetric membrane, the symmetry can be explicitly broken because of the finite curvature imposed in the spherical topology. Like a "carapace" forced to be curved, the Lo phase is broken into smaller pieces when it is bent because it is more rigid than the surrounding Ld phase. An equilibrium mesophase separation emerges.
On the other hand, membrane inclusions (such as proteins) or asymmetric lipid composition of both leaflets can lead to local spontaneous membrane curvature C sp (r) = 0 coupled to the local composition φ(r). In this case, the elastic energy of the domain grows faster than its area (it is non-extensive), which sets a maximum domain size above which domains become thermodynamically unstable. In this case also, a patterned mesophase is stable in equilibrium, also called a cluster phase in the case of membrane proteins.
These three mechanisms can be compared to the ones presented in Section 2.1.1 in the weak segregation limit, even though quantitative comparison is difficult because the technical details are not always easy to put in direct correspondence between the different models. No spontaneous symmetry-breaking mechanism was identified in the weak segregation limit where the line tension is too weak to induce buckling. In contrast, differences in rigidities κ or in spontaneous curvatures C sp have been shown to lead to rich phase diagrams in both limits, especially in the spherical geometry. As compared to the strong segregation case, the weak segregation one proposes, in addition to mesophases, the existence of SD phases where density fluctuations are transient with a typical size corresponding to a maximum of the structure factor. Density fluctuations are known to exist in any biphasic system, with a typical size set by the Ornstein-Zernike correlation length ξ OZ in absence of coupling. Because shape and composition are coupled, SD-phase fluctuations are different since their size is proportional to the Helfrich correlation length ξ, characteristic of shape fluctuations. It is the natural length since at distances larger than ξ, the membrane mechanics is only controlled by the surface tension σ and bending effects become negligible. This is illustrated by Equations (9) or (11) for coupling to spontaneous curvature and detailed in [85] for coupling to bending rigidity. Presumably, their typical size ensues from the same physical mechanism as in the mesophases, as follows. In the vicinity of a miscibility critical point, critical phenomena theory predicts that the domain size distribution decays with a power law cut off at a size set by ξ OZ [52]. But too large domains, even though transient, should become energetically unfavorable in the coupled case for the same reason as in the strong segregation limit, because of their increased elastic energy induced by either κ 1 or C 1 (or both). They split into smaller domains, the typical size of which is now set by q −1 c . In the vicinity of the miscibility critical point,we have seen that κ 1 does not play any role in the planar case at Gaussian order, whereas it does in the spherical case. This is consistent with the findings in the strong segregation limit, where the role played by κ 1 is directly related to the fact that R is finite (see Figure 6). More quantitatively, if only the dependence of κ with φ is considered in this spherical geometry, and not the spontaneous curvature, then one must focus on the C 1 = 0 vertical line in the phase diagram of Figure 3c. When κ 1 grows above κ * 1 as given in Equation (15), the liquid phase becomes a structured disordered one where density fluctuations correspond to thicker transient domains. Increasing further κ 1 leads to even more structured mesophases, as in the strong segregation case.
As far as spontaneous curvature is concerned, a temperature T * (C 1 ) has been defined in Section 2.1.1. Below T * and for large enough C 1 , mesophases are favored by spontaneous curvature as in the strong segregation limit. Above T * , the structured disordered phase now corresponds to small transient domains of the spontaneously curved phase. If only spontaneous curvature is taken into account and no dependence of κ with φ is considered, one must focus on the κ 1 = 0 horizontal line in Figure 3c. When C 1 grows above a certain threshold on the order of C * 1 , the liquid phase becomes a structured disordered one. In the planar geometry, Figure 2b also reveals that the macrophase below T * is destabilized by spontaneous curvature, leading again to a modulated mesophase. The same observation can be made in Figure 2a

Active and Out-of-Equilibrium Processes
Alternatively to the approaches in thermodynamic equilibrium discussed so far, some works have discussed for a dozen years [5,39] the possibility that out-from-equilibrium, active processes might explain the observed domain finite size. Generally speaking, the mechanisms at play appeal either to the active remodeling of the cortical cytoskeleton (which is coupled to the membrane [6]) or to the rapid traffic of membrane components. Traffic occurs to and from the membrane, either by exchange of monomers or by endocytosis and exocytosis of membrane patches through vesicles of ∼ 100 nm in diameter [189], with a surface of 10 4 to 10 5 nm 2 . This active exchange of membrane material is often called "membrane recycling". The relationship between domain size and endocytosis rates has been explored experimentally in at least one biological system [190]: the data suggest that blocking the endocytic pathway leads to an increased number of large clusters of E-cadherins (containing more than ∼ 100 copies), supporting the idea that the growth of large domains is impeded by recycling. It was concluded by the authors of reference [190] that "endocytosis targets large E-cadherins clusters, while not perturbing the mechanisms of clustering at scales smaller than 100" copies.
Experiments also indicate that the typical amount of plasma membrane endocytosed per hour is equivalent to the total surface of the cell plasma membrane of a fibroblast [191]. Faster recycling times for membrane components have even been measured, on the order of 10 min or even shorter (see [192] and references therein). As illustrated in Figure 8, recycling enters in competition with the growth of domains, as driven either by domain diffusion followed by their coalescence, or by Ostwald ripening [116,193] (see below). Here the domain coalescence is still driven by a finite line tension below T d which tends to minimize the interface length between Lo and Ld regions. Note that as in Abney and Scalettar's work [194], Gheber and Edidin had already anticipated such a mechanism in 1999, above T c however [195]. In their model, an additional meshgrid of obstacles hinders diffusion, mimicking the presence of the cortical actin cytoskeleton, and ensures a sufficient lifetime of domains by preventing the too fast diffusion of membrane molecules. Note also that if the transport between the membrane and the cytosol is passive, without any energy consumption, the membrane is actually in equilibrium, in contact with a molecule reservoir. One recovers the equilibrium case discussed in the previous sections, but now described in the grand-canonical ensemble [196].
Active membrane recycling can fracture large assemblies or domains, the growth of which (below the demixing temperature) is thus arrested at a finite typical radius L * in the steady state [56,57,197,198]. As in the equilibrium case, these models have first been designed to account for the finite size of Lo lipid domains, however, at this stage they are sufficiently general to be applicable to protein clusters as well. Both cases will be treated in parallel in this section. discussed in further detail below; 0 when there is no recycling. The kernals k n;m and k 0 n;m control, respectively, the rates of domain scission in which one domain of n m monomers breaks into two, of size n and m, and domain fusion, in which two domains containing n and m monomers fuse to form a single domain of size n m [see Fig. 1(a)]. A similar approach has been rather successful in describing the kinetics of wormlike micelles [20].
We assume that two domains fuse whenever they diffuse into contact. Thus k 0 D=s, with D a characteristic diffusion constant for the domains and s the area of the smallest (monomeric) domains. The characteristic, microscopic ''diffusion'' time scale D 1=k 0 s=D 10 ÿ5 s for biological membranes. We propose to neglect any size dependence of the diffusion constant of the domains, which is a fair approximation due to the nonintuitive nature of hydrodynamics in quasi-two-dimensional systems [21]. This simplifies the analytic analysis and should only weakly affect our results, which rely primarily on the fact that the recycling rate is much slower than 1= D .
The scission kernal k n;m is now determined by the principle of detailed balance, an approach that is appropriate provided the longest intradomain relaxation time following a fusion or scission event is shorter than the domain collision time. This should be satisfied on average for fluid domains [22] where the rate of subsequent interdomain events should then converge to that found at equilibrium for each (pair of) domain(s), in spite of the fact that the size distribution may be far from equilibrium. By inspection of Eq.
where k 0 n;m k 0 1= D is treated as a constant in what follows. This kinetic scheme, without recycling, will result in the usual slow growth of domains (coarsening) [23]. In order to analyze the effect of lipid recycling we first consider the ''monomer deposition/raft removal'' (MDRR) recycling scheme in which raft lipids and proteins are brought to the membrane at random as ''monomer'' sized units with a rate j on , possibly on vesicles that consist of predominantly nonraft membrane, and entire rafts are lost from the membrane at random with a rate j off , irrespective of their size [see Fig. 1(b)]. Thus n appearing in Eq. (3) takes the form n j on n;1 ÿ j off c n ; (5) with i;j the usual Kronecker delta. It is easily shown that j on =j off at steady state. In general Eq. (3) requires numerical solution (see Fig. 2 and [24]), although an asymptotic solution is possible in the most interesting regime 1, as discussed below.
It can be seen from Fig. 2 that the domain size distribution is broad; indeed there is significant contribution to the total area fraction from domains with n & n 2 . The distribution depends only weakly on when is itself large because all domain scission events are then rare. As a result of this we will see that full numerical solution of the master equation may not always be necessary.
To investigate this ''scissionless'' large regime further, we proceed by neglecting all scission terms in Eq.
where we will later have to check our solutions for selfconsistency, which will translate to establishing a lower bound for . The equations _ c n 0 for n 1; 2; . . . ; can then be used to build up c n recursively. For MDRR recycling given by 10 100

Models
A basic but illuminating scaling argument was proposed by Foret in 2005 [56]. When quenching a homogeneous lipid mixture below the demixing temperature, it phase-separates. Under certain realistic conditions, the typical domain size coarsens as L(t) ∝ (D a t) 1/3 [116], where D is the two-dimensional diffusion coefficient, and a is the typical lipid (or protein) lateral size, in the nanometer range. Additionally, traffic to and from the membrane sets a typical recycling timescale τ, as discussed above. Recycling prevents equilibration beyond the timescale τ, and the coarsening stops at a typical domain radius L * ∝ (D a τ) 1/3 . Assuming D 0.1 µm 2 /s for cell membrane constituents [74], a smaller than few nm for a typical membrane lipid or protein, τ ≤ 1 hour and ignoring numerical prefactors in the expression of L * , one gets L * 1 µm. This rough estimate is encouraging, just one order of magnitude larger than typical nanodomain sizes. This argument has been completed and confirmed by several theoretical studies in the last decade. They dwell on a variety of theoretical techniques, such as linear stability analysis in the Fourier space [45,56,198,200], or numerical and/or analytical resolutions of stochastic equations belonging to two principal categories: (i) Master Smoluchowski's coagulation equation [57,196,199,201]. For example, Turner et al. [57] studied the coagulation equation where c n is the area fraction of domains or clusters containing n monomers and the constants k n,m and k n,m are scission and fusion rates, respectively. Finally, ζ(n) controls recycling. Both cases in Figure 8a,b correspond to ζ(n) = j on δ n,1 − j off c n because only monomers are injected in the membrane with a rate j on while domains are recycled to the cytosol with a rate j off independent of their size n. (ii) Non-linear reaction-diffusion equations [56,197,198,200,202], which can also be seen as Cahn-Hilliard equations [15,116] suitably modified to take recycling into effect. For instance, in Ref [56], the Cahn-Hilliard equation is used. The order parameter φ(r, t) again measures the local concentration of one phase (say the Lo phase),φ is its equilibrium value, τ is the recycling rate as discussed above, and A is a kinetic constant proportional to the lipid diffusion coefficient. The function w(φ) is the Landau-Ginzburg expansion of a membrane patch free-energy, at the relevant order in powers of the order parameter φ, as already evoked in Section 2. Recycling is embedded in the first term of the r.h.s, while the second term describes the diffusive transport of the conserved order parameter φ.
The different models can be further categorized as follows (see also references [196,203] for alternative discussions): 1.
The off-rates (from the membrane to the cytosol) can be size-dependent [190] or not [57,199].
In the former case, it means for example that endocytosis is able to extract patches from the membrane with a limited size set by the endocytosed vesicle typical size [201]. A "recycling correlation length" can also be introduced in the modified Cahn-Hilliard equations, mimicking the spatial range of recycling processes, i.e., the typical size of membrane patches recycled through vesicle traffic [197,202]. In the models of Refs. [56,57,196,198,200], only monomers are locally extracted from the membrane.

2.
The on-rates (from the cytosol to the membrane) are also size-dependent. Several models only inject monomers or tiny domains in the membrane [56,57,190,196,[198][199][200][201] because they do not assume any pre-order in the exocytosed patches or because they assume direct exchange of monomers from the cytosol to the membrane, e.g., for peripheral proteins. Indeed, Foret argues that the traffic should be modeled differently for peripheral and transmembrane proteins [196], because the former are preferentially exchanged as monomers between the cytosol and the membrane, while the latter preferentially escape and join the membrane by endo-and exocytosis, respectively. Another approach assumes that domains with a characteristic size are directly injected in the membrane [197,202].
To finish with, an alternative viewpoint has considered intermittently active proteins, which undergo conformational transitions driven by external energy sources of any form (such as light for photoreceptors, or ATP for other active proteins) [204,205]. The model assumes that the protein affinity for the Lo or Ld phases depends on its conformation, for example because of hydrophobic mismatch energy. Otherwise macroscopic domains are then split into finite-size ones. Indeed, large protein assemblies are destabilized because protein changing conformation continually leave them, as it is the case with recycling above. Accordingly, the typical domain size is demonstrated to decrease with enhanced protein activity.

Results and Prospects
These models agree on the existence of finite-size domains in the steady state instead of the macrophase that would prevail in equilibrium. However, they can lead to quantitatively different predictions concerning the steady-state distributions p(L) of domain sizes, depending on the recycling scheme details. As illustrated in Figure 9, the different approaches can predict 1.
Or large, power-law distributions p(L) ∝ L −δ , with a large domain-size cut-off above a limiting upper size L max , above which domains become extremely rare [57,190,201].
(a) (b) When going from curve "a" to curve "f", the number of monomers in the membrane increases because more and more monomers are injected in the out-of-equilibrium membrane from the cytosol. Above a critical value, multimers nucleate and the distribution becomes bimodal as in curves "d", "e" and "f". The limit between monomodal and bimodal distributions is curve "c". Adapted from [196]. (b) Power-law distributions p(n) ∝ n −3/2 with a cut-off n max (vertical dashed line) above which p(k) decreases exponentially. Here the coalescence rate between two clusters is independent of their size. Clusters smaller than a size n v are recycled entirely as a whole, while larger clusters are fragmented and lose an area n v during a recycling event. The "whole cluster recycling" limit corresponds to n v → ∞ and the "monomer recycling" limit to n v = 1. Adapted from [201]. For both figures, see the cited references for more details on the recycling dynamics.
All models predicting a bimodal distribution inject only monomers in the membrane from the cytosol and are based on the Ostwald ripening dynamics inside the membrane (exchange of monomers between clusters). These two prerequisites will have to be confirmed with a solid justification in future studies. In contrast, the way off-rates depend on the domain size does not seem determinant in setting the shape of p(L).
Note also that Smoluchowski's coagulation-equations where clusters of any size can coagulate or split into smaller ones have been known for long to lead to power-law distributions with a cut-off in quite generic circumstances [206]. The fact that clusters of any size can coagulate or not depend on whether they feel mutual repulsion or not. Mutual repulsion can be attributed to membrane spontaneous curvature inside clusters [199]. The fact that clusters can split into smaller multimers depends on the line tension λ. If λ is too large, creating two smaller new clusters is quite expensive in terms of interfacial energy and Ostwald ripening is favored. Hence the shape of the distribution p(L) certainly strongly depends on the way the different membrane constituents interact at short and longer range.
To conclude this section on active out-of-equilibrium processes, it is worth emphasizing that experiments are absolutely required to clarify the rich but somewhat involved landscape of the theoretical predictions that we have described. To our knowledge, the experimental approach proposed in Ref. [190] is the unique one to bring insight onto the effect of membrane recycling rates on membrane domain sizes. Experimental are thus wanted to assert (or not) the relevance of the sometimes contradictory assumptions on which the different models discussed above rest.

A Variety of Mechanisms in Equilibrium in the Strong Segregation Limit
We have focussed in this review on biphasic fluctuating membranes that phase-separate in two phases at sufficiently low temperature. We have named these phases Lo and Ld even though the formalism at play can in principle be extended to a wider context, in particular to protein clusters in an homogenous lipid mattress. We have classified the various mechanisms involved in the stabilization of mesophases and we have argued that many of them dwell on the broken symmetry between the inner and outer (or upper and lower) leaflets.
Even in the case of simplified lipid bilayer mixtures, where only three or four lipid species coexist, the phenomenology is already dramatically rich, as exemplified by Feigenson and collaborators' studies. Here we discuss the different mechanisms presented from the perspective of their relevance in the cell biology context. Which mechanism(s) is/are likely to be at play in a real cell to promote the observed nanodomains having a 10 to 200 nm diameter, commonly called "rafts" in the cell biology community?
Buckling induced by a strong line tension (see Section 2.2.1), as first deciphered by Lipowsky and collaborators in the 1990's [110,111], has been shown to lead to micrometric domains in the metastable state where buckled domains feel mutual repulsion. We have seen that the line tension required to promote buckling of Lo domains is λ > 8κ Lo /L for domains of radius L in taut membranes. For κ Lo ≥ 100 k B T as measured experimentally, buckling of domains of radius < 100 nm require λ > 30 pN, an order of magnitude larger than measured line tensions. The situation is even worse for 10-nm nanodomains, the admitted lower bound of raft sizes. These values seem to exclude buckling as a generic realistic mechanism in the cell biology context. In addition, measured critical temperatures T c lie below the physiological temperature, in the T c = 20 to 30 • C range, even in giant plasma membrane vesicles [207]. It ensues that Lo and Ld phases are probably not separated in living cells. However, this last argument must be taken with care because many membrane constituents are missing in model membranes. They might shift the actual critical temperature to higher values.
In the phenomenon studied by Lipowski's and Feigenson's groups (see Section 2.2.2), the more rigid Lo phase "carapace" must be split into smaller domains to accommodate it to the global spherical shape of the GUV [64,127]. With the typical bending moduli of Lo phases given above, and GUVs with diameter R > 10 µm, the domains are observable by classical optical microscopy and are again on the micrometer scale. Getting much smaller nanodomains requires to increase the background mean curvature 1/R. This can be achieved on highly curved cell substructures such as microvili, small cellular vesicles, dendrites or axons, but it is unlikely on flatter membrane regions where nanodomains are also observed. The same remark as above concerning the critical temperature holds.
In contrast, domains ensuing from the spontaneous curvature C sp = 0 induced by some of their constituents can readily lead to nanodomains for realistic parameter values because they make too large domains unstable, as discussed in Section 2.3. We have explained how a large variety of mechanisms can promote membrane spontaneous curvature by playing on the asymmetry between both leaflets. Without excluding that the two previous mechanisms are also more marginally at play, spontaneous curvature is probably the most efficient mechanism in the strong segregation limit.
Linactants are also known to reduce the line tension λ and to promote metastable micro-emulsions by localizing at the interface between coexisting phases, as surfactants do in many soft condensed matter phenomena [11,18]. This role has been hypothesized for cholesterol or asymmetric saturated/monounsaturated lipids, like POPC and 1-stearoyl-2-oleoyl-sn-glycero-3-phosphocholine (SOPC) , although this last possibility has been contested (see references [20,31,33,208] and references therein). From a practical point of view, reducing the line tension λ with the help of linactants is comparable to driving the system closer to a critical point where λ vanishes [14].

Nanodomains and Critical Density Fluctuations
The previous section dealt with the strong-segregation limit. In Section 2.3.5, we have underlined that the density fluctuations in the weak-segregation limit are quite specific when membrane shape and composition are coupled. In absence of coupling, biphasic systems close to criticality belong to the Ising universality class and the characteristic fluctuation length-scale is set by the Ornstein-Zernike correlation length ξ OZ which diverges at the critical point and remains large in its vicinity [12,15,18]. Above T c , the domain-size distribution decays with a power-law up to a maximum size set by ξ OZ [209]. By contrast, when coupling is switched on, the nanodomain size in the Structured Disordered phases is set by the Helfrich correlation length ξ which has no reason to be comparable to ξ OZ , and the domain-size distribution is peaked at ξ where the structure factor has a maximum.
Since critical density fluctuations are good candidates for the nanometric rafts predicted by biochemical techniques [7,9,14], they have been probed by a large variety of experimental approaches. Indeed ∼ 10 nm nanodomains in ternary lipid mixtures can be observed by FRET and electron spin resonance [20,24], SANS [31] or AFM [27,29]. However, the domain size distribution cannot be determined through these experimental techniques and only conjectures can be made concerning their precise nature. It cannot easily be ascertained whether they are just the manifestation of density fluctuations in the vicinity of a critical point or correspond to the structured disordered phases evoked above.
To clarify this issue, Veatch and collaborators observed by fluorescence microscopy giant plasma membrane vesicles (GPMVs) directly isolated from cells [207]. As compared to the model three-or four-component membranes discussed so far, GPMVs retain many of the original plasma membrane constituents, lipids and proteins, even though some up-down asymmetry is likely lost during the GPMV formation [18]. They also separate into Lo and Ld phases at low temperature. Even though the resolution is not as good as in AFM experiments, fluctuations above T c could be analyzed at scales ≥ 150 nm, i.e., above the Rayleigh limit. Their observations were qualitatively consistent with all GPMVs having near-critical compositions, suggesting that plasma membranes have compositions that are naturally designed to reside near a miscibility critical point in vivo. Thermodynamic quantities such as the correlation length above T c , the line tension below T c , the compressibility and the diffusion coefficient could be inferred close to a critical point. In spite of important error bars, the authors concluded that all of them indicated that the system belongs to the 2D Ising universality class, i.e., that of an ordinary binary mixture close to its critical point. A similar study was conducted on GUVs made of a ternary lipid mixture close to its critical point and the authors were led to similar conclusions [210]. However, in such studies subject to Rayleigh limit, the presence of a peak in the structure factor at short wave-lengths 150 nm, as predicted by coupled theories, cannot be excluded. Connell and collaborators circumvented this limitation in 2013 by examining very accurately the bilayer composition fluctuations with an AFM [27]. Lo and Ld phases can easily be distinguished because they do not have the same thicknesses. Importantly, the authors observed that "The region of critically fluctuating 10-100 nm nanodomains has been found to extend a considerable distance above T c to temperatures within the biological range, and seem to be an ideal candidate for the actual structure of lipid rafts in cell membranes." Domains of size ≥ 10 nm could be observed at temperatures up to 7 • C higher than T c . They also concluded that supported DOPC/SM/cholesterol bilayers belong to the 2D Ising universality class because the measured structure factors have the expected shape and the critical exponent associated with the correlation length is close to the 2D Ising-class value, ν = 1. However, the membrane observed by AFM are supported, which suppresses wave-length fluctuations above some cut-off depending on the strength of the substrate-membrane attraction, and potentially inhibits the effects of shape-composition coupling. These experiments have not definitely ruled out the existence of structured disordered phases either.
To finish with, a potentially important structural difference between domains below the demixing temperature T d and density fluctuations close to the critical point has recently been suggested through coarse-grained molecular dynamics [211]. Well below T d , there is a strong overlap between Lo domains in opposed leaflets (at least in planar geometry), which we have called registered domains above. In contrast, the overlap is negligible in the case of density fluctuations (anti-registration). This difference could be fully relevant from a biological perspective.

Needed Theoretical Clarifications
When observing 2D patterned morphologies, it appears that in some circumstances, domains are more or less roundish, like bubbles, whereas they become more elongated, leading to interdigitated, labyrinthine structures, in other situations (see, e.g., Figure 5c). They can even adopt a relatively regular striped geometry in some cases. This issue was already examined in the articles of Andelman and collaborators [45,59] because shape instabilities have been studied for long, back to the 1970's, for example in ferrofluid films in a magnetic field or polar phospholipid Langmuir films at suitably chosen temperatures and pH. Fluorescence microscopy imaging on vesicles made of a SM/DOPC/cholesterol mixture below the phase-separation temperature have shown that there is a transition from stripes to bubbles when the reduced volume v grows and approaches 1, presumably because of increasing surface tension [117], as later confirmed by field-theoretic arguments based on the Hamiltonian (4) and the single-mode and weak-segregation limit approximations in planar geometry [66], and much more recently by new experiments [209]. By using the same kind of arguments, a transition from bubbles to stripes is also predicted when increasing the concentration of the minority phase at fixed thermodynamic conditions or lowering the temperature [59,127,166,167]. However, as pointed out in Ref. [66], a consensual understanding of the relative stabilities of the different morphologies remains necessary. Experiments with controlled surface tension are also still lacking to guide theoretical developments. In the even more complex cell biology contexts, experimental data tend to confirm that increasing the protein concentration leads to more elongated protein clusters [162] (see Figures 5 and 5S of this reference). More systematic experimental studies in connection with theoretical works are also needed.
The understanding of the role of electrostatics in shaping the membrane, even though already considered in the 1990's [212], also remains lacunar. As already indicated in Section 2.2.2, reference [128], followed by reference [33], studied numerically how dipole-dipole or charge-charge repulsion between lipids, having a nanometric range because of screening, competes with the shorter-range attraction between lipids of the same phase. As expected [45], competing interactions are found in this work to lead to small domains, of size ∼ 10 nm, provided that line-tension is not too large (in which case macrophase separation occurs instead). As discussed in the review by Schmid [14], the uncertainty of how to finely model membrane leads to contradictory predictions concerning the resulting domain size. In addition, we also remarked above that models embracing both nanoscopic and micrometric length-scales remain to be developed. Very few works have tackled them on an equal footing so far, probably because the level of coarse-graining cannot be the same for the different scales which makes numerical studies difficult. A multi-scale approach is needed. If both electrostatic lipid-lipid repulsion at nanoscopic scale and a mechanism as described above leading to larger domains are present, does one expect a hierarchy of domains, with nanoscopic domains embedded in microscopic ones? Such a scenario has been proposed on analytical grounds in [134] where, again appealing to the Gaussian approximation of a 3-field theory, a region of the parameter space was found where nanometric "raft" domains dwell in larger, micrometric "nonraft" ones, themselves distinct from the "background" components. Such hierarchy has also been suggested by different experimental techniques [213][214][215], but its connection with theoretical findings will have to be confirmed in future studied.
Furthermore, we have evoked that some studies presented in this Review include Gaussian curvature and the saddle-splay modulus κ G in membrane modeling [53,63,64]. It is generally assumed that κ G ≈ −κ to ensure membrane stability [64]. However, in most analytical and numerical works, the corresponding second term in Equation (3) is neglected, by arguing that at fixed topology, its integral over the whole membrane surface is an invariant due to the Gauss-Bonnet theorem. But if κ G turns out to depend significantly on the lipid phase, this argument fails and a more rigorous approach is needed. Jülicher and Lipowsky have tackled this question in 1996 [111], using the fact that when κ G,Lo = κ G,Ld , the integral over the whole surface can be simplified, still owing to the Gauss-Bonnet theorem. It is proportional to κ G,Lo − κ G,Ld and can be written as the integral of the geodesic curvature along the interface between both phases, see Equation (24) of reference [111]. It can potentially favor a longer interface and thus a mesophase [216]. However, Jülicher and Lipowsky's study was only dedicated to the shape of vesicles undergoing macrophase separation, and not to the stability of mesophases. To our knowledge, a systematic study of the interplay between the different contributions to the membrane energy remains to be done to ascertain the precise role played by the Gaussian curvature and clarify when it can be neglected or not.

Needed New Experiments
In Section 3, we have reviewed the different mechanisms proposed in the literature to take active membrane recycling into account. At this stage, there is no definitive reason to favor mechanisms in equilibrium with respect to active ones. However, nanodomains can be observed in membrane fragments peeled off from cells long before their observation. This tends to prove that active processes are not absolutely required to maintain nanodomains. For example, retina cell fragments from dead human and murine corpses have been observed by AFM. They presented very well characterized rhodopsin nanodomains [28]. This does not preclude that active processes can come in addition to equilibrium arguments and reorganize the mesophases. Berger and her coworkers have recently combined such out-from-equilibrium arguments and a competing-interaction model as discussed in Section 2.3 to explore how recycling modifies the equilibrium domain sizes [199] (see Figure 8b above). The main finding of this analytical work is that when taking recycling into account, the typical cluster size at steady state increases very slowly, logarithmically with the recycling rate. Using physically realistic model parameters, the predicted two-or three-fold increase of the typical cluster size due to recycling in living cells is likely experimentally measurable with the help of super-resolution microscopy, even though such experiments might be delicate to implement because the typical cluster size depends on many other parameters. In addition, in this model and contrary to others, an increase of the typical clusters size is predicted with increased recycling at fixed protein density, because only monomers are injected in the membrane and only monomers are exchanged between clusters (what we have called Ostwald ripening [193]). However, as already evoked in Section 3.2, cluster-or domain-size distribution shapes (Figure 9) should help to discriminate between the different propositions based on theoretical considerations. They will have to be confronted to high-precision experiments in the future, thus enabling biologists and biophysicists to identify the mechanisms that are actually at play in living cells.
Some experimental studies [24,49] have suggested a transition between two distinct modulated phases, one with a nanoscopic domain length-scale (possibly critical density fluctuations), and one with much larger domains. Does the apparently discontinuous character of the phase diagram just come from the diversity of the experimental techniques used to probe the different regimes of parameters? This point will have to be clarified in future studies. It is closely akin the problem raised in the previous section about the domains predicted at different scales by different theoretical approaches.
In this work, we have considered Lo domains and protein domains on an equal footing, because they share many common characteristics from a theoretical perspective. Real rafts in cells are assumed to gather together specific proteins in a Lo phase nanodomain. In the specialized literature, deciphering the respective roles of lipids and proteins in the stabilization of rafts is sometimes controversial. Does lipidic mesophase-separation precede targeting of protein to rafts or do proteins gather and bring their specific lipid annuli with them, thus enriching nanodomains in some lipid species? This somewhat sterile debate seems useless to us because a nanodomain is stabilized by the aggregation of all its constituents, lipids and proteins, that feel mutual attraction at short range, as discussed in detail in this Review. However, the important question is to understand what sets their finite, nanoscopic size. Imbedding membrane proteins to three-or four-component lipid GUVs will probably constitute a step further towards the understanding of how lipids and proteins collaborate in real cells.

Concluding Remarks
Some major pieces of the jigsaw puzzle have been identified, but the subtle way in which they assemble is still the object of passionate debates. We have argued that the best candidates in equilibrium to be the so-called rafts in real cells are: (i) nanodomains induced by the curvature of some of their constituents, in the strong-segregation limit, as discussed in depth in Section 2.3; or (ii) composition fluctuations in the vicinity of a miscibility critical point, as described in Sections 2.1 and 4.2. Both mechanisms likely coexist in real cells. However, a somewhat confusing situation prevails in the biophysics literature, in particular because the physical concepts are not always used in their strict domain of validity. As pointed out by Veatch and Keller, "it can be confusing to use the language of phase separation when describing small features [composition fluctuations] because they differ in many ways from large domains [meso-or macro-phases]" [17]. In addition, for both kinds of domains, their size distributions are potentially modified by active, out-of-equilibrium phenomena (see Section 4.4).
A way to disentangle further this long-standing issue might be to characterize better domain lifetimes. Nanodomains must have sufficiently long lifetimes to allow their constituents to interact [20]. Measured lifetimes of composition fluctuations are on the order of 1 s, due to their rapid destruction by diffusion [207]. Is it sufficient to enable nanodomains to play their biological role of micro-reactor? To our mind, this is one of the major issues that ought to be challenged in the forthcoming years.