Phase Separation on Bicontinuous Cubic Membranes: Symmetry Breaking, Re-entrant and Domain Facetting

We study the phase separation of binary lipid mixtures that form bicontinuous cubic phases. The competition between non-uniform Gaussian membrane curvature and line tension leads to a very rich phase diagram, where we observe symmetry breaking of the membrane morphologies and reentrant phenomena due to the formation of bridges between segregated domains. Upon increasing the line tension contribution, we also find facetting of lipid domains that we explain using a simple argument based on the symmetry of the underlying surface and topology.

Introduction -Lipid self-assembly can adopt an astonishing range of shapes and morphologies, from single bilayer structures to stacks and convoluted periodic structures [1]. Nature has, of course, exploited this polymorphism. A large number of organelles feature lipidbased structures, including synaptic vesicles, endoplasmic reticulum, and Golgi apparatus. At the same time, lipids are indispensable for detergency and foodstuffs industries [2], and membrane-based structures are increasingly exploited in biotechnological and biomedical applications, e.g. as efficient nanoporous scaffolds for tissue engineering [3] or for gene silencing with siRNA [4].
In this letter we will focus on one particular type of mesophases that lipid mixtures in water can adopt, the so-called bicontinuous cubic phases (BCP) [5][6][7], whereby the lipids form a triply periodic lipid bilayer that separates two percolating and non-intersecting water channels [8][9][10]. These phases have attracted attention due to their high surface area, continuity of the bilayer surface, and pore network. The amphiphilic nature of the lipids also allows other molecules to be embedded in them; for example, they have high propensity to enable membrane protein crystallization. Although the details remain unclear, it is thought that a combination of curvature induced phase separation on the cubic surface, a local destabilization of the cubic phase to a lamellar phase and a two dimensional reservoir of proteins provided by the cubic phase are responsible for the observed yield [11]. Here our interests are in the aforementioned curvature induced phase separation.
Both in the biological and synthetic systems, these lipid mesophases usually contain more than one lipid species. To the best of our knowledge, the distribution of different lipids across such a cubic surface, especially the possible demixing transitions under the influence of non-uniform curvature of the membrane structures, is still not well-understood. Most studies on lipid phase separation focus on much simpler membrane geometries, such as lipid vesicles and supported membranes [12][13][14]. From a biological perspective, lateral lipid organizations into domains and membrane curvatures are ubiquitous features, and are known to play an important role for the membrane functionalities [15,16]. From a materials perspective, understanding the distribution of species of interest on a BCP may be the first steps towards a systematic and rational functionalization of BCPs, where active species can be localized into targeted domains. Finally, our work provides a comprehensive phase diagram, with predictions of distinguishing features which we hope will stimulate experimental verifications.
This letter is organized as follows. We first show that if the two species do not interact but induce different bending rigidities, then a single type of curvature induced phase separation occurs at all non zero area fractions. Upon considering interactions between the species, we observe a multiplicity of new modalities for the phase separation, including the formation of bridges between previously disconnected lipid domains. Moreover, we observe facetting of domains for which we provide a simple explanation relying on symmetry and topology.
Segregation in absence of line tension -In this paper we consider a binary lipid mixture or, alternatively, a mixture of lipids and proteins that has formed a minimal surface S (with zero mean curvature everywhere), and ask what are the thermodynamically favoured repartitions of the species and how they depend on the bending rigidities and inter-species interactions. Here our focus is on the triply periodic surfaces which are known to be formed by lipid mixtures as well as mixtures of lipids and proteins in water [5][6][7]. We use standard notations P, D, and G for the primitive, Diamond, and Gyroid surfaces respectively. Since they are periodic, we characterize their properties per unit cell. To model a binary mixture on a curved surface, we use a straightforward extension of the Helfrich hamiltonian [17] to the case of a binary mixture on a minimal surface that reads [18][19][20][21]: where x denotes a point on S, dµ S (x) is the area measure on S at x, G(x) < 0 is the gaussian curvature at point x, g are the gaussian bending rigidities associated to the species A and B respectively. Our convention here is κ A g < κ B g , such that δκ < 0. It is also worth remarking that, since typically κ g < 0, B domains are softer than A domains. Such a model may represent a coexistence between L o (liquid ordered; A-rich) and L d (liquid disordered; Brich) domains, or alternatively, between lipid-rich and protein-rich domains.
To model the distribution σ A in a more tractable way, we first stress that any minimal surface with genus g embedded in a flat torus T 3 must contain 4(g − 1) zero Gaussian curvature points [22]. For the P, D and G surfaces, g = 3 and they have 8 zeros. Each zero is located at the centre of a hexagonal area we term as a patch (See Fig. 4 here and Fig. 1 in [18]). The unit cell of either of the P, D or G surfaces can thus be partitioned into 8 equivalent patches {Σ i } i=1..8 such that the unit cell sur- We characterise the repartition of the lipids on S by both the area fraction f i A of lipid A on each patch Σ i and the occupation number function σ i A (x) in it. For each patch Σ i , given σ i A , the entropy then reads A is a Lagrange multiplier that imposes the value of f i A . At low temperatures, the Fermi-Dirac distribution will reach a value close to unity for all points x of Σ i with an energy lower than λ i A . The lowest energy point p i in a patch Σ i is the symmetry point of the patch which has exactly zero gaussian curvature (cf. Fig. 4). Thus, at low T , the lipids A will fill the neighbourhood of p i until they reach a critical Fermi curve C F , where {x ∈ C F |δκG(x) = λ i A }, beyond which there is no more lipids of type A (cf. Fig. 4; a disconnected area occupied by lipid A is termed as a domain). Close to p i , one may use polar coordinates (ρ i , θ i ) and, as a crude approximation, the space is assumed euclidean and circularly symmetric near p i . This allows us to Taylor expand the function G about p i up to the second order so that the curvature energy reads H Σi Here C is a constant and µ(Σ i ) is the area of the hexagonal patch Σ i . Remarkably, in spite of the very crude approximations we have used, the predicted behaviour of the curvature energy H Σi is close to what we observed in simulations for the P-surface where the exponent is found to be α = 1.83 [18].
Next, upon minimising the total free energy F = C , it is easy to see that the ground state in repartition among the patches is always f i A = f A for all values of f A , corresponding to the 8 8 configuration in Fig. 4. Effect of the line tension -We have seen that, with only curvature, the A lipids are evenly distributed among the 8 available patches and formed dense domains in the neighbourhood of zero curvature points at low temperature. This begs the question of how this picture changes if the A − B interactions are not negligible i.e. if there are line tension effects arising with domain formation, which is a more realistic physical scenario. To answer this question, we now carry out computer simulations of the binary phase separation on the P-surface (the qualitative picture is the same for the D-and G-surfaces, as justified in [18]). There are several known approaches to model bicontinuous cubic membranes, from coarsegrained Molecular Dynamics simulations [23] to continuum field theoretical approaches [24][25][26]. Here we use Metropolis Monte Carlo simulations [27] to resolve the thermodynamics of the system.
In our approach, we explicitly discretize a piece of the P-surface contained in a cubic cell. This can be efficiently done with the help of the Weierstrass-Enneper (W-E) representation of minimal surfaces [22,28,29]. Upon discretization, the binary mixture can then be modelled as an Ising-like problem (see [18] for technical details). A spin variable s is associated to each site and takes either value 0 (for species B) or 1 (for species A). The curvature hamiltonian of Eq. (1) thus maps exactly onto a system of magnetic spins on a network N (S) with a node-dependent external magnetic field and reads: where δA i is the area of the tile i on the surface. In this language, at any finite f A , species A (spin variable s = 1) will occupy sites with the lowest value of δA i δκG i to min- imize the total energy, as we have analyzed with a different vocabulary in the previous section. To model the A − B interspecies interactions, we choose a short-range nearest neighbours interaction which directly translates into the line tension of the lipid domains: where J sets the magnitude of the exchange interactions, δL ij is the length of the edge shared by cells i and j, and (s i + s j − 2s i s j ) = 1 when s i = s j and 0 otherwise. Symmetry breaking -As it is evident, the hamiltonian in Eq. (3) is equivalent to an Ising model of ferromagnetism and therefore should lead to the same phenomenology: above a critical temperature T * (J), the system is paramagnetic and the two lipid species are mixed; while below T * , the system becomes ferromagnetic and a symmetry breaking favouring "lumping" of spins in spatial regions (segregation) occurs. There is, however, one crucial difference between the standard Ising model and our model. For the former, line tension effects always dominate demixing: domains of A lipids coalesce to minimize the overall interfacial energy. In our model, this coalescence mechanism competes with the curvatureinduced mechanism described in the previous section.
The first effect of line tension is to re-shuffle the (energy) ranking of configurations 8 k with k patches occupied by the A species by shifting down the low k configurations (because they have a lower interfacial cost) and up the high k ones (because they have a high interfacial cost). A first account of the competition between curvature and line tension consists in assuming that the total energy of a configuration 8 k at a given packing fraction as the sum of the free energy of individual patches of equal size. This summation approximation is valid when isolated domains are formed at the centre of the hexagonal patches, and one finds that increasing f A at fixed J /|δκ| always favours, eventually, higher k values, in agreement with Monte Carlo simulation results shown in the phase diagram in Fig. 5 The caveat is that this summation approach is only valid when the A species domains are disconnected. Above certain f A values, the lowest energy configurations are in fact those in which domains of lipid A span across multiple patches (see Fig. 5). These bridges between patches essentially make the domains interact negatively and in a non-pairwise fashion. The location of these bridges coincides with the lowest curvature energy regions at the patch boundary (c.f. Fig. 4). Taking these configurations into account, the phase diagram in Fig. 5(b) shows that the simple picture of Fig. 5(a) only holds for small J /|δκ| and f A . In fact, one observes re-entrant behaviours whereby a configuration 8 k previously unfavored in the disconnected regime, becomes refavored thermodynamically. We note that when bridges are formed in the 8 8 configuration (open crosses in Fig.  5(b)), the segregated and continuous phases are effectively inverted (lipid B domains are surrounded by A).
Domain facetting -Another distinguishing feature that appears with line tension is the facetting of the domains formed by the A lipids. This effect is shown in Fig. 3(b) where the domain almost draws a hexagon compared to Fig. 3(a) where the shape is more rounded, thus the term "facetting". To explain this, we recall that in general if the underlying manifold has an n-fold rotational symmetry, we expect the bounding curve that minimises the perimeter length of a domain with fixed area to be a regular n-gon whose sides are geodesics of the underlying manifold. Moreover, on an anisotropic curved surface, not all orientations of a regular n-gon are equivalent as they lead in principle to different total perimeter lengths. Thus, we interpret the bounding curve in Fig. 3  6-fold symmetry to be the curve that minimises both shape and orientation at the same time.
To test the above rationale, we estimate the geodesic curvature along the bounding curves of the two representative examples shown in Fig. 3(a) and (b) (cf. e.g. Ref. [30]). In these figures, the curvilinear coordinate l ∈ [0, 1] is the normalised arc length of each curve which enables the comparison of the geodesic curvature for curves with different total lengths in Fig. 3(c). In absence of line tension, the geodesic curvature g c is approximately constant around the boundary. With a large line tension, g c reaches very high values for l close to zero, but is much smaller than that without line tension as l approaches 1. This is consistent with the above explanation although it shows that the facetting is not perfect. It nevertheless sheds light on what happens as we approach the ideal facetting case: the geodesic vanishes almost everywhere except close to l 0 where it diverges. This divergence is representative of the wedge formed by the intersection of two geodesics of the 6-gon and whose angle γ can be estimated  [31,32]. Most synthetic BCPs, however, are formed using the lipid Monoolein, which is known to have a low bending rigidity with |κ G | ∼ κ m < 10k B T [33]. Here κ m is the (mean curvature) bending rigidity. Using these values and taking the typical lattice spacing of a BCP 10 − 100 nm [34][35][36], this leads to J /|δκ| in the range of O(0.1)-O(1) considered in this paper.
When the line tension effects can be neglected, the natural curvature of the surface alone is enough to a) in-duce segregation in all surfaces and b) the segregation is such that domains form in the same proportions on all available patches on the surface. We then confirmed this theoretical prediction by numerical calculations on the P-surface and looked at the effects of non-zero line tension. The latter gives rise to two important features. (i) Below the demixing critical temperature, it favours the formation of bigger domains in a fewer number of patches available on the surface that we characterize with a corresponding phase diagram. We also observe re-entrances in this patch-occupation space due to the formation of bridges between domains on neighbouring patches. Some of these morphologies should lead to distinguishing features (e.g. different x-ray scattering signatures due to the change in symmetry), and we hope this work will stimulate experimental works to verify our predictions. (ii) In the large line tension limit, we observed a facetting of the domains for which we provided a simple explanation and that we can relate to the curvature energy of a domain on a patch.
Predicting patterning on cubic membranes is the vital first step towards their systematic and rational functionalization. On one hand, the ability to localize molecular species by design into targeted domains can be beneficial for controlled release in drug delivery or of chemical substances [37,38], and for templating self-assembly [39] or phase separation [40,41] in the surrounding fluids. On the other hand, suppressing phase separation between lipid species or between lipids and proteins can be desirable in applications such as protein crystallization [11], where segregation at an incorrect stage can strongly hamper the efficiency of the applications.
There are also a number of avenues for future work. Firstly, here we have assumed that the BCP remains a minimal surface. A closer inspection based on the theory developed in [42] for domain-induced budding shows that the conclusions presented here can be qualitatively affected when κ m /|δκ| < 0.4 [18]. However, estimates of this ratio for a wide range of lipid bilayers and monolayers in the literature show that it is only rarely below 1 [43]. This suggests that the minimal surface assumption is very reasonable for realistic parameter values. Further work is however still needed to fully assess how membrane deformation, including budding instability, affects the phase diagram of multicomponent BCPs. Secondly, the present work tacitly assumes that the membrane domains are formed by lipids of the same species in the two leaflets (registration phase). Indeed, recent work on flat bilayers suggests that registered domains is the thermodynamically favoured phase for a wide range of lipid mixtures [44]. It would be interesting to relax this assumption to probe how curvature affects registration/anti-registration and how, in turn, registration/anti-registeration may affect the bilayer morphology. Thirdly, the system considered here provides an excellent setup to study how nonuniform curvature may affect the nature of the demixing phase transition.

FREE ENERGY OF A BINARY MIXTURE ON A MEMBRANE Free energy for single species membranes
For a membrane comprising of a single specie, the free energy of a given shape S can be described via the Helfrich hamiltonian [17]: where C ≡ (c 1 + c 2 )/2 is the local mean curvature, G = c 1 c 2 is the local gaussian curvature, c 1,2 are the local principal curvatures, C sp is the intrinsic mean curvature of the surface and κ m,g are the mean and gaussian bending rigidities respectively. Furthermore, the Gauss-Bonnet theorem states that [30]: where ∂S stands for the boundary of S, g c for the local geodesic curvature of the boundary [30] and χ(S) for the Euler characteristic of S. The key consequence is that, if the boundary and the topologies of the problem are fixed, minimizing Eq. (SM1) is equivalent to minimizing only the integral over the mean curvature C. In the absence of intrinsic curvature, the lowest energy solutions correspond to surfaces with exactly C = 0 at all points: such surfaces are called minimal surfaces [28]. Here we focus on the triply periodic surfaces which are known to be formed by lipid mixtures in water [5][6][7]. We use standard notations P, D, and G for the primitive, Diamond, and Gyroid surfaces respectively. Since they are periodic, we characterize their topology with their Euler characteristic per unit cell (e.g. in Eq. (SM2 )).

Free energy for a binary mixture on a minimal surface
If instead of a single species, the membrane comprises two different species, then different allowed mixture configurations may have different topologies and one cannot disregard anymore the gaussian curvature contribution to the energy. The simplest extension of Eq. (SM1) to a binary mixture would then read: where x denotes a point on S, dµ S (x) is the area measure on S at x, G(x) is the gaussian curvature at point x, f A is the imposed area fraction of species A (with f B = 1 − f A ), the field σ A (x) ∈ [0, 1] is the mean occupation number of species A at x and κ A,B g are the gaussian bending rigidities associated to the species A and B respectively. We then choose the free energyH S el (f A = 0) as reference so that, in practice, we look at the free en- is the starting point of our study.

General formulation
It can be shown that, locally, any minimal surface can be conformally mapped onto the complex plane via the Weierstrass-Enneper (W-E) representation [28]. More precisely, the W-E is a map from C to R 3 which, to a point (u, v) in an open subset of C, uniquely associates a point (x(u, v), y(u, v), z(u, v)) of R 3 that belongs to a minimal surface via: where g is a holomorphic function and f is a meromorphic function such that f g 2 is analytic.
The gaussian curvature G(x, y, z) at any point of a surface represented by Eqs. (SM5), (SM6) and (SM7) can be expressed as a function of w = u + iv via: (SM8) The negative sign is characteristic of minimal surfaces.
Since the mean curvature (c 1 + c 2 )/2 is zero everywhere, it implies that the gaussian curvature is always negative or zero. The W-E being a conformal map, it preserves the angles. Distances, however, are not conserved when mapping an infinitesimal segment from C to R 3 and are scaled by a factor Λ(w) given by Triply periodic Schwartz surfaces The three Schwartz surfaces P, D and G can be obtained by choosing: where θ B is the Bonnet angle such that θ B = 0 for the D surface, θ B = π/2 for the P surface and θ B = cotan(K(1/4)/K(3/4)) for the G surface. K is a complete elliptic integral of the first kind. Since e iθ B only changes the phase in the W-E representation, this means that the P, D and G surfaces are simply related by an isometry called the Bonnet transformation and share many of their physical properties.

Independence of the free energy on the member of the Bonnet family
In particular, having the W-E map in mind, the whole integral in Eq. (SM4) can be thought of as an integral in the complex plane. Using the fact that dµ S (x(w)) = Λ 2 (w)dudv, the Eq. (SM4) can be recast as: where A(S) denotes the atlas used in C to characterize S. It is worth noting that in the integrand of Eq. (SM12), the W-E functions f and g only appear via their complex modulus and therefore their contribution to the curvature energy would be unchanged by a phase factor. Thus, we concludethat the curvature energy of a binary mixture on a minimal surface is independent of which member of the Bonnet family is considered and, in particular, so is its ground state. In a similar fashion, a continuous model of the line tension contribution would read formally: where I(A−B) denotes the set of points belonging to the A−B interface on S and dµ l (x) the length measure on S. Again, the integral can be thought as an integral on C by virtue of the W-E map. Furthermore, the length measure on S can be expressed in term of the length measure in C via dµ l (x(w)) = Λ(w)|dw|. Eq. (SM13) can thus be rewritten as: |dw| Λ(w).
As before, the W-E functions f and g only contribute to the integral via their complex modulus and therefore H A−B is independent of the member of the Bonnet family under study. As a consequence, the phenomenology of symmetry breaking and bridging-induced reentrant behaviour depicted in the phase diagrams of the P-surface in Fig. 2 of the main text hold in fact for all three Bonnet surfaces. The details of the phase diagram may slightly differ, however, as the cubic cells of the G and D surfaces do not contain the same surface area as the P surface. This is outside the scope of this letter, and we will discuss these details in a separate publication.

Angle of intersection between two geodesics
The Euler characteristic of a compact domain D of A lipids that is not bridged to another domain on a neighbouring patch is 1. Applying the Gauss-Bonnet theorem (Eq. (SM2 )) to such a domain gives: Up to a factor the first term of the l.h.s of Eq. (SM15) is simply the curvature energy of the domain that we can denote H Σi el (f i ) in referring to notations introduced in the main text. By using the fact that a patch has a 6-fold rotation symmetry, we can split the boundary integral of the geodesic curvature into 6 equivalent parts. If each piece of the boundary is almost everywhere a geodesic, then the second term of the l.h.s of Eq. (SM15) has zero integrand everywhere except at points where the geodesics meet. The total value of the contour integral becomes simply a sum over intersection angles that we call θ. We thus have: Finally, the actual interior angle γ between two geodesics making the hexagonal-like facetted domain is in fact the complementary angle of θ and reads: . (SM17) Note that in the case where the curvature energy vanishes but the symmetry is still imposed, we retrieve the interior angle of a planar hexagon as expected.

Fundamental patch
All well behaved minimal surfaces admit a description in terms of a fundamental patch in R 3 that is repeated by using the symmetries of the surface. By the W-E representation, this fundamental patch is associated to a fundamental domain of the complex plane. For the P, D and G surfaces, the fundamental domain is the set of complex points with positive real part bounded by the lines along the vectors (1 + i)/ √ 2 and 1 and by the circle of radius √ 2 whose center is located at the point −(1 + i)/ √ 2. The top left of Fig. 4 shows the fundamental domain in C. The triangular tessellation is obtained by using the Surface Evolver package [19] and the images plus the management of the network structure have been performed with the Mathematica software [20]. For the sake of illustration, Fig. 4 shows a coarse tessellation of the fundamental domain. The tessellations we used in the paper are typically 100 times finer. By using Eqs. (SM5)-(SM7) and (SM10) and (SM11), we get a three dimensional realization of the fundamental patch that is represented in the top right of Fig. 4. Then, following Ref. [29], we can generate first a full hexagonal patch of the P-surface (Σ) by replicating and stitching together 12 fundamental patches as seen in the bottom left of Fig.  4. The full cubic cell representation of the surface S is then obtained by combining 8 such hexagonal patches with the right symmetry operations as illustrated on the bottom right of Fig. 4. A similar procedure, albeit with different arrangements of the fundamental patches, can also be carried out for the D-and G-surfaces.

Monte Carlo simulations
As emphasized in Eq. (SM12), the curvature energy can be recast in terms of a sum over points on a euclidean (complex) plane of a curvature field that multiplies a scaling field. Moreover, the discretized surface S on the bottom right of Fig. 4 is made of a network of cells N (S) which are either an original version or a replica of a cell in the fundamental patch. Thus, the whole set of values of the curvature and scaling fields on the whole network is determined solely by that of the sub-network of cells in the fundamental patch. The particular topology of the P-surface (of genus 3 in a cubic cell) is then accounted for by the topology of the network i.e. by assigning the right neighbours to each cell. If we add a species field s into the picture such that s i = 1 if cell i ∈ N (S) contains species A and s i = 0 otherwise, then the whole problem becomes that of paramagnetic spins on a network subject to an effective node-dependent magnetic field whose magnitude is G(w)Λ 2 (w)∆(w) and where ∆(w) denotes the euclidean area of the triangular unit at point w in the complex plane. Since the effective magnetic field is non uniform and non trivial, there is no simple explicit analytical expression for the thermodynamically favoured composition morphologies. By splitting the system into 8 equivalent patches, we could however suggest, as discussed in the main text, what would happen at low enough temperatures. In particular, a first approximation scheme neglecting the effect of curvature on the area measure and Taylor expanding the curvature field about its zero point p i suggested that the curvature energy If we set the Ising parameter J to zero, we can then probe the low energy curvature energy as a function of The red triangles are MC data points for the energy H Σ i el (f i A ) of an A domain on a patch Σi. The lines correspond to the best fit to these data points with a scaling law behaviour The solid red line corresponds to the best fit exponent value α ≈ 1.83, while the green dashed line is the best fit to the data with imposed α = 2 which shows quite good agreement with the data.
domain size for a single patch. As we see in Fig. 5, the proposition that the curvature energy goes as a power low is in very good agreement with MC data. In addition, the approximation that α = 2 is in remarkably quite good agreement with the data.

VALIDITY OF THE MINIMAL SURFACE ASSUMPTION
In the current study, it is assumed that the underlying surface remains minimal during segregation of the lipids/species. However, in general domain formation can induce local membrane deformation which in turn can destabilise the entire morphology of the membrane itself. In this section we will have a closer look at how phase separation induced bud formation may affect the conclusions drawn in the manuscript.

Budding on flat multicomponent membranes
Based on the work by Lipowsky [42] on flat membranes, domain formation may lead to a budding phenomenon driven by the line tension between the two coexisting demixed phases. Budding occurs when the line tension energy cost overcomes the bending energy penalty.
Consider a membrane domain of area A = πL 2 = 2πR 2 (1 − cos θ), where R is the radius of curvature of the deformed membrane domain, θ is the contact angle of the domain with respect to the horizontal plane, and L is the domain radius if there is no deformation to the membrane. We will now consider the competition between two energy terms: a) line tension energy that increases with the perimeter of the domain, 2πJR sin θ; and (ii) bending energy which depends on the domain area and curvature, 2κ m A/R 2 . Here J is the line tension and κ m is the mean curvature bending rigidity. We have also assumed that there is no mean spontaneous curvature. Fig. 6 shows the total energy (normalised by the line tension energy for a flat domain) as a function of the reduced membrane mean curvature L/R that plays the role of an order parameter. L/R = 0 corresponds to a flat membrane, i.e. no deformation. L/R = 2 corresponds to a complete bud formation. As we see, there are three regimes depending on the membrane domain size L: (i) For 0 < L < L * = 4κ m /J, budding is unfavourable and L/R = 0 is the global minimum configuration; (ii) For L * ≤ L < L o = 8κ m /J, budding is favourable but there is an energy barrier for its formation; (iii) Finally, only for L ≥ L o that the energy barrier for bud formation disappears. Here the flat membrane geometry is completely unstable.

Bud formation on a P-surface
Since the P-surface is a minimal surface with zero mean curvature everywhere, we will continue to assume that the spontaneous mean curvature is zero, as in the previous paragraph. The presence of non-uniform gaussian curvature on the P-surface may alter the numerical prefactors in L * and L o . However, since the domain formation occurs in the neighbourhood of very specific points on the P-surface i.e. 8 zero-curvature points each at the centre of a patch in the cubic cell, to first approximation, it is reasonable to assume that the above results from [42] to hold. As a result, we can approximate the two critical domain sizes as L * = 4κ m /J and L o = 8κ m /J. Correspondingly, when L < L * , we expect our assumption that the underlying surface remains minimal during phase separation to hold. Strong deviations to the results presented in the main text are only expected when L ≥ L * . Figure 7. Modification of the phase diagram of Fig. 2(a) in the main manuscript resulting from budding instability for κm/|δκ| = 1/4. The diagonally hashed region corresponds to a region where budding is preferable but there is an energy barrier for its formation, while the horizontally hashed region corresponds to a fully unstable bud formation.
To see which part of the phase diagram in Fig. 2 of the main manuscript is affected by the budding instability, we relate the area fraction f A occupied by the Alipids to the size L of the domains depending on which configuration 8 k they are in. For the purpose of this analysis, we will focus on cases where there are no bridge formations between lipid A domains. The total area occupied by the A-lipids on a cubic cell of the P-surface is f A µ S (S). If the A-lipids are partitioned in k patches among the 8 available, then the typical area per domain is f A µ S (S)/k ≈ πL 2 . It follows that requiring full stability against budding (L < L * ) is equivalent to requiring f A µ S (S)/k < π16(κ m /J) 2 . If the cubic cell bounding the P-surface has sides of length , then µ S (S) = 24 2 K(1/4)/K(3/4), where K(x) is the complete elliptic integral of the first kind [29]. Thus, the stability criterion becomes for configuration 8 k . As we see in Eq. (SM21), the set of stability lines depends on the ratio κ m /|δκ| which constitutes an additional parameter in our modelling. We have observed that for values κ m /|δκ| > 0.4, there is no intersection between any of the stability lines and the lipid repartition phases they correspond to within the parameter ranges of the present study, 0 ≤ J /|δκ| ≤ 2.
It is worth emphasizing that experimental and numerical estimates of the ratio κ m /|δκ| for various lipid systems [43] show that it is rarely below 1. As an example, consider the mixture of DOPC, sphingomyelin, and cholesterol forming coexisting L o and L d domains reported in references [31,32]. Here the difference in Gaussian bending moduli |δκ| 3×10 −19 J and the mean curvature bending modulus for the L o phase κ m 8×10 −19 J, which gives us κ m /|δκ| 8/3. Even if we use the mean curvature bending modulus for the L d phase κ m 2 × 10 −19 J, which may be more appropriate when bridges are formed such that the lipid B domains are now surrounded by A (see Fig. 2(b) and (c) in the main text), we still have κ m /|δκ| 2/3 > 0.4. Thus, this strongly suggests that our conclusions in the main text will hold even when taking into account segregation induced membrane deformation.
We can redo a similar calculation to that leading to Eq. (SM21) to determine the phase boundaries beyond which the membrane is completely unstable against bud formation. Not surprisingly we find that they satisfy a very similar equation: (SM22) To illustrate how the phase diagram is modified when budding instability is taken into account, we show the results for κ m /|δκ| = 1/4 in Fig. 7. The parameter regime susceptible to budding is for large area fraction of A lipids, f A , and large line tension, J /|δκ|. The hashed regions correspond to parameter regimes where budding is preferable. Energy barriers are present for budding to occur in the diagonally hashed region, while for the horizontally hashed region the P-surface is completely unstable against budding.