Nucleation of symmetric domains in the coupled leaflets of a bilayer

We study the kinetics governing the attainment of inter-leaflet domain symmetry in a phase-separating amphiphilic bilayer."Indirect"inter-leaflet coupling via hydrophobic mismatch can induce an instability towards a metastable pattern of locally asymmetric domains upon quenching from high temperature. This necessitates a nucleation step to form the conventional symmetric pattern of domains, which are favoured by a"direct"inter-leaflet coupling. We model the energetics for a symmetric domain to nucleate from the metastable state, and find that an interplay between hydrophobic mismatch and thickness stretching/compression causes the effective hydrophobic mismatch, and thus line tension, to depend on domain size. This leads to strong departure from classical nucleation theory. We speculate on implications for cell membrane rafts or clusters, whose size may be of similar magnitude to estimated critical radii for domain symmetry.


I. INTRODUCTION
A phase-separating lipid (or other amphiphilic) bilayer may access competing equilibrium and metastable phase coexistences, due to the presence of two leaflets subject to competing inter-leaflet couplings (Figs. 1, 2a) [1,2].A "direct" inter-leaflet coupling [3][4][5][6][7] promotes registered (R) bilayer phases, in which both leaflets are locally dominated by the same species.An "indirect" coupling from hydrophobic tail length mismatch favours uniform bilayer thickness [8][9][10][11] and thus promotes antiregistered (AR) phases, in which the leaflets are dominated by different species so that the bilayer is locally compositionally asymmetric.Depending on the choice of overall leaflet compositions, two (e.g., R-R or AR-AR) or three bilayer phases may coexist.For example, asymmetric overall leaflet compositions can lead to two approximately symmetric bilayer phases and a highly asymmetric one (R-R-AR), first observed and explained in [6].Phase equilibria of coupled bilayer leaflets were subsequently studied using phenomenological free energies [3,5,12].However, a full description of the indirect coupling described above requires a model that microscopically incorporates hydrophobic mismatch [1].
Using such a model, we have shown how coexistence of antiregistered phases can be kinetically preferred due to the effect of hydrophobic mismatch [1] so that, before equilibrating, a quenched bilayer must escape a (typically metastable) locally asymmetric state.The understanding of this novel statistical thermodynamics will allow greater control over artificial membranes, and the underlying interactions are expected to play a role in transmembrane organisation of rafts or clusters in vivo, with possible relevance to signalling [13] and anaesthetic action [10].The predicted behaviour constitutes an example of Ostwald's "rule of stages" [14], by which a system FIG. 1. Partial phase diagram showing competing equilibrium R-R (black) and metastable AR-AR (red) phase coexistences [1].Spinodals enclose the regions of local stability.Cartoons of the dominant inter-leaflet arrangement in each bilayer phase are shown.The grey dotted line illustrates R-AR coexistence, which is briefly discussed in Section III A. Other phase coexistences not considered here are omitted [1].Parameters: ∆0 = 2 a, κ = 3 a −2 kBT , V = 0.6 kBT , J = 4 a −2 kBT , B = 0.48 a −2 kBT .
will pass through available metastable states on its way to equilibrium.Ostwald's rule is familiar (via different origins) in colloids, metallurgy and drug design.
Immediately after quenching a bilayer to a phaseseparating region of parameter space, any spinodal instabilities to which the uniform state is subject compete to determine the dominant initial demixing mode.A bilayer with roughly equimolar composition in each leaflet [15] can be subject both to an "R mode" with composition perturbations locally symmetric between leaflets, and a perpendicular AR mode with asymmetric perturbations (Fig. 1).If the AR mode is fastest-growing, spinodal de-composition to AR-AR coexistence occurs first, leading to local asymmetry throughout the bilayer.To equilibrate to R-R from a metastable AR-AR state, the bilayer must undergo nucleation of registered bilayer phases.Hence, three classes of kinetics arise: direct separation into equilibrium phases, equilibration via nucleation out of a metastable state, or trapping in a metastable state.For other overall compositions subject to competing instabilities, the competition of symmetric and asymmetric phases is qualitatively similar [2], though more complex.
In this paper, we focus on the nucleation energetics that determine whether equilibrium domain symmetry is reached from a metastable state.First, we introduce the model and discuss the interpretation of bilayer domain symmetry and asymmetry via phase diagrams with a composition axis for each leaflet, with reference to existing experiment and theory.We then identify the three classes of kinetics in simulation, guided by a linear instability analysis, and develop a theory for the nucleation of registered domains, which captures the interplay of bulk free energy with thickness mismatch occurring at the perimeter of a registered domain.Together with the linear stability analysis, the calculated nucleation energetics are consistent with the simulation results.We find that the effective hydrophobic mismatch between a domain and its surroundings is domain size-dependent, which causes strong departure from classical nucleation theory.

II. MODEL
A detailed description of the model appears in [1,2].We briefly recapitulate the model and its analysis in terms of phase diagrams and kinetics.A schematic of the lattice model and how it is coarse-grained is shown in Fig. 2.

A. Lattice model
A local bilayer patch is described as N lattice sites with top (t) and bottom (b)-leaflet lipids (Fig. 2b).The lipids' hydrophobic lengths ℓ t(b) i lead to the total bilayer thickness Model species S and U represent saturated and unsaturated lipids or, e.g., the liquid-ordered (L o ) and liquid-disordered (L d ) states in a ternary mixture [2].
The Hamiltonian is  ideal (i.e., preferred) hydrophobic tail lengths are ℓ t(b)i 0 = ℓ S0 for an S lipid at the top (bottom) of site i, or ℓ U0 for U , and each site is pairwise registered (R, SS or U U ) or antiregistered (AR, SU or U S).
V ≡ V 10 − 1 2 (V 00 + V 11 ) quantifies purely intra-leaflet interactions, such as those between headgroups.The "direct" coupling B promotes pairwise R between lipids, nominally by penalising tail structure mismatch (which we treat as implicit in tail length mismatch [16]) across the midplane.The particular mechanisms responsible for the direct coupling are not crucial to our model, however -for comparison with the literature we can simply estimate an effective strength of the conventional inter-leaflet mismatch energy γ [1], which is shown on Fig. 3 as well [4,5,7,[17][18][19].The hydrophobic "indirect" coupling J promotes pairwise AR, by penalising mismatch in the bilayer thickness profile.We also define J ≡ 4 J, which appears in the mean-field approximation of Eq. 1 used to derive the coarse-grained free energy.κ can be related to the area compression modulus κ A [1], and penalises variation from species-dependent ideal length.Weaker κ means the species can more easily adapt their tail length and structure to one another's presence.The mismatch parameter ∆ 0 ≡ ℓ S0 − ℓ U0 is cast as a length, but represents both tail length and structure mismatch; it couples to both the indirect and direct inter-leaflet couplings, J and B. Once fiducial values of the parameters are set, varying J alone approximates changing the mismatch in tail length but not structure (e.g., adding carbons to one species [9]), while varying B alone approximates varying the mismatch in tail structure, e.g., unsaturation.We arbitrarily choose ℓ S0 > ℓ U0 .The reference total thickness d 0 ≡ ℓ S0 + ℓ U0 is irrelevant in the absence of an external field acting on bilayer thickness.

B. Simulation protocol
We simulate a L 2 = N bilayer where L = 100 (script letters refer to the entire simulated bilayer, as opposed to a local bilayer patch).The Kinetic Monte Carlo scheme [2] resembles Kawasaki (spin-exchange) dynamics, governed by the Hamiltonian (Eq.1): lipids exchange between neighbouring lattice sites within their leaflet, thus mimicking diffusive evolution after a quench from a hightemperature, randomised initial configuration.Hydrodynamics are not included, but we do not expect this to significantly alter the conclusions; we address this further in the Discussion.
The overall leaflet compositions are conserved, since we do not consider flip-flop or exchange with the solvent, and are given by where is the total number of S lipids in the top (bottom) leaflet.R and AR lattice sites can interconvert (SU + U S ⇄ SS + U U ) to alter the "degree of registration" (microscopic transbilayer symmetry) where N SS is the total number of SS lattice sites, etc. λ can vary in the range Although the overall leaflet compositions are conserved, λ is not [2].Hence, the bilayer can phase-separate into either locally symmetric or asymmetric modes while the conserved overall leaflet compositions may, as here, be fully symmetric between leaflets.
We choose a mismatch parameter ∆ 0 = 2 a ∼ 1.6 nm, somewhat larger than the typical length mismatch in phospholipid mixtures, for which a (registered) phase thickness mismatch 2 nm [23,24] would imply ∆ 0 1 nm.Large ∆ 0 couples to both the indirect (J) and direct (B) couplings, increasing the energetic driving forces for both antiregistration and registration and making the competing phases clearer to interpret in simulation.The phenomenology is qualitatively similar upon reducing ∆ 0 [2], with the caveat that ∆ 0 affects the effective value of γ arising from a given B (Eq. 8).In Section V we mention the dependence of the nucleation theory on ∆ 0 .
We use V = 0.6 k B T in the mean-field theory, above the mean-field threshold V 0 ≡ 0.5 k B T for phase separation in the absence of other couplings, and use V = 0.9 k B T in simulation, where the corresponding threshold is V sim.0 = 0.88 k B T due to fluctuations [25].Although we cannot expect precise quantitative agreement between the mean-field theory and simulation, in this regime it appears the exact value of V is not crucial to the kinetics; for example, a value V = 0.9 k B T in the mean-field theory would yield a similar predicted landscape of relative R/AR growth rate ∆ω to that calculated in Fig. 3 [2].

III. REGISTERED AND ANTIREGISTERED PHASES
We now discuss bilayer domain symmetry and asymmetry in terms of competing bilayer phase coexistences, focusing on those relevant to the present work and briefly considering those relevant to bilayers of asymmetric overall leaflet compositions (Φ b = Φ t ).Upon coarse-graining the microscopic model Eq. 1 (Fig. 2c), a local bilayer patch is characterised by locally-averaged top and bottom leaflet compositions φ t(b) ≡ i φt(b) i /N .We derive a mean-field free-energy density f (φ t , φ b ) [1], which yields a phase diagram in (φ t , φ b ) space (Fig. 1) [1].The coexisting bilayer phases determine the local order parameter in each leaflet, given by the projection of a given tie-line endpoint onto the φ t or φ b axis.

A. Definition of registration and antiregistration
Some experiments, particularly with asymmetric overall leaflet compositions, may use separate fluorophores to image domain morphology in each leaflet [26][27][28].Interpreting these in (φ t , φ b ) space, and relating them to experiments that do not image the separate leaflets [24], requires care.We define a registered (R) bilayer phase as one in which both leaflets are dominated by the same species.Hence, in our model, an R phase is dominated by SS or U U lattice sites (Fig. 1) such that most lipids face one of their own species in the apposing leaflet.An antiregistered (AR) bilayer phase is one where the leaflets are dominated by opposite species.Hence, under our definition, "registration" is a property of a given homogeneous patch of the bilayer, describing approximate local compositional symmetry between the leaflets.
An alternative definition of registration is sometimes used, which we call "colocalised enrichment" [26,27].This describes a bilayer in which the regions of largest top-leaflet composition φ t (relative to the average in that leaflet ) spatially superimpose on the regions of largest bottom-leaflet composition φ b .Colocalised enrichment is therefore a property of domain morphology over the entire bilayer, not of an individual bilayer phase.It requires a tie-line of finite positive slope in (φ t , φ b ) space, so that: i) both leaflets contain domains of larger and smaller than average φ t(b) ; and ii) the domains of large φ t belong to the same bilayer phase as the domains of large φ b and are thus spatially colocalised with them.R-AR coexistence (Fig. 1 dotted line) can be accessed by a bilayer of asymmetric overall leaflet compositions.R-AR tie-lines have positive slope, and thus exhibit colocalised enrichment, although the AR phase is highly asymmetric in composition.Hence, the colocalised enriched domains in each leaflet reported in [27,28] (where the leaflets were separately imaged) are consistent with either an R-R or R-AR tie-line.Quantitative composition information would be required to unambiguously determine which was observed.
In the present work, we consider a bilayer of symmetric overall leaflet compositions for which R-R and AR-AR tie-lines compete.In this case, "registration" versus "colocalised enrichment" are practically equivalent.R-R has positive tie-line slope (colocalised enrichment), and both bilayer phases are compositionally symmetric (registered phases, under our definition).Conversely, AR-AR tie-lines are negatively-sloped so that enrichment in one leaflet colocalises with depletion in the other, and both bilayer phases are compositionally asymmetric (antiregistered phases, under our definition).
The literature is ambiguous.For example, full colocalised enrichment (which can reflect either an R-R or R-AR tie-line) is described in [27] as "registration".In contrast, in [10], "registration" is also defined as the presence of L o (or L d ) on both sides of the bilayer, which is more like our definition of R phases versus AR phases.This ambiguity is not surprising because often, as here, the two definitions outlined above are similar.The two concepts can clash for asymmetric overall leaflet compositions [28] where R-AR tie-lines can play a role.
In asymmetric supported bilayers in [27], the leaflets were imaged with separate fluorophores, revealing colocalised domains in both leaflets.The same was found in highly asymmetric vesicles in [28].On the face of it, these findings differ from some supported bilayers in [24] which exhibited compositionally asymmetric regions interpreted as domains in only one leaflet.These were inferred from variations in total bilayer thickness without imaging the separate leaflets.Height mismatch smaller than a known value for R-R was measured, indicating either R-AR or AR-AR-R coexistence.R-AR was then inferred by detecting that saturated lipids were predominantly in the top leaflet.The system was interpreted as having domains in only the top leaflet [24].
However, it is probable that the asymmetric bilayers in both [27,28] and [24] represent an R-AR tie-line.Thus, imaging the leaflets in [24] separately could have shown colocalised enrichment just as in [27,28].The top-leaflet gel domains in [24] may have apposed regions in the bottom leaflet weakly more "gel-like" than the average in the bottom leaflet.Similarly, the colocalised enrichment in [27,28] does not require that bottom [27] (or inner [28]) leaflet domains are truly L o -for R-AR they will, like the rest of the bottom/inner leaflet, be dominated by unsaturated lipids.It only requires that they are weakly more "L o -like" than the average in their leaflet, i.e., that the R-AR tie-line is tilted.Hence, the question of whether domains in one leaflet "induce domains in the other" becomes the question of whether the given tie-line is tilted enough for both leaflets to exhibit detectable domain formation when a different fluorophore is used in each.The degree of R-AR tie-line tilt will depend on the direct inter-leaflet coupling B and hence on molecular features.

B. Kinetics and competing phase coexistences
For most typical parameter choices, the R phases are lower in f (φ t , φ b ) than are the AR, hence R-R tie-lines are equilibrium and AR-AR metastable [1].As well as the coarse-grained bulk free energy f (φ t , φ b ), the underlying microscopic model yields free-energy costs for composition and thickness gradients, which describe the role of spatial structure and domain formation in the kinetic competition of metastable and equilibrium states.Linear stability analysis of the initial homogeneous state (Fig. 3) predicts growth rates of competing R versus AR instability modes and thus whether symmetric or asymmetric domains form first [1].The interplay between gradient and bulk free energies then governs the nucleation of symmetric domains (introduced in Section V), which determines the eventual fate of a bilayer that has initially become metastably AR-AR.The ability to derive both gradient and bulk free energies from defined microscopic interactions is a key advantage of the present model over a purely phenomenological approach.
In Fig. 3, R-R separation is equilibrium (except for a very small region Ba 2 /k B T 0.005).Below the "AR instability" line, AR-AR separation is also possible, but metastable.Below the "R instability" line, the initial state is not subject to R instability although R-R separation remains the equilibrium state.The red and blue colours show the difference in growth rates, ∆ω, between R and AR instability modes.For example, red signifies a faster-growing AR instability mode, so that AR-AR separation dominates initial demixing.
For physical parameter ranges [1], neither the R or AR mode is trivially dominant, so that moderate changes to lipid tail length mismatch (affecting the effective value of J) or tail structure mismatch (affecting B) can determine whether the initial instability after a quench leads to locally symmetric (R) or locally asymmetric (AR) domains.In addition, the long-lived AR-AR states simulated in [9,10] provide prima facie evidence that metastable trapping due to failure to nucleate R domains is possible for physical phospholipids.However, due to small simula-FIG.3. Red/blue colours and dashed/dotted lines: mean-field theory parameter map showing relative growth rates of AR versus R instability modes for a bilayer comprising equimolar mixed leaflets, from linear stability analysis of initial demixing [1], with ∆0 = 2 a, κ = 3 a −2 kBT , V = 0.6 kBT .The equilibrium state is R-R in all cases.Below the "AR instability" line AR-AR coexistence is possible but is metastable.Below the "R instability" line the homogeneous state is not unstable to the R mode, although R-R separation is still the equilibrium state.Overlaid dots from simulation (with V = 0.9 kBT ) show the average degree of registration λ at the end of the simulation time, λ = 0 (black, AR-AR) to λ = 1 (white, R-R).We identify three kinetic classes discussed in the text (circled 1, 2, 3, bold line marks approximate boundaries).Illustrative simulation snapshot sequences for each class are shown (L = 200 in the snapshots).Simulations are visualised with OVITO [29].
tion sizes in such studies, this could constitute a stable, not metastably trapped, state, an issue which we address in the Discussion.
Finally we note that, contrary to [30], AR-AR coexistence does not require an exactly equimolar (or equal area fractions) mixture in each leaflet.Firstly, an overall composition away from Φ t = Φ b = 0.5 can lie on one of a set of AR-AR tie-lines running parallel to the central one depicted in Fig. 1 [2].Secondly, even if the overall composition is outside any AR-AR tie-line, AR-AR-R coexistence can occur, in which the presence of some R phase is forced by the composition being far from equimolar, but the R phase coexists with two AR phases [2].In that situation, every region of the bilayer except the R phase is, as for AR-AR, locally asymmetric.Like AR-AR, for typical parameters, AR-AR-R is a metastable state, which could become trapped for strong hydrophobic mismatch or stabilised if domain size is limited (e.g., by simulation size).

IV. KINETIC PHASE DIAGRAM
Fig. 3 shows the results of simulations of phasetransition kinetics for varying indirect coupling J and direct coupling B, to model varying lipid tail length mismatch and structural mismatch.The overlaid greyscale dots signify whether registered or antiregistered domains dominate at the end of the simulated time (t = 10 6 Monte Carlo Steps) by measuring the degree of registration λ, with λ ≈ 0 (black) corresponding to AR-AR coexistence and λ ≈ 1 (white) to R-R.Each dot on Fig. 3 is an average of four independent trajectories.This simulated "kinetic phase diagram" agrees semiquantitively with the theoretical linear stability analysis of initial demixing (red/blue on Fig. 3): inside the blue region, where R-R coexistence should be accessed directly, the simulation exhibits full registration.If instead the AR mode is fastest (red), we find two possibilities -the bilayer may reach R-R coexistence by nucleation out of the AR-AR state, or remain metastably trapped in AR-AR.Fig. 4a shows the successful formation of registered nuclei, which grow from the boundaries of antiregistered domains.In Fig. 4b, a different random quench in the simulation leads to failure to nucleate despite unchanged parameters, illustrating the stochastic nature of the nucleation process.
Having identified the expected three classes of kinetics, we next model the energetics of nucleating registered domains, to clarify the fates of the metastable AR-AR state (the second and third kinetic classes identified in Fig. 3).

V. NUCLEATION THEORY
The metastability of AR-AR implies that arbitrarily small composition fluctuations decay.To reach R-R, registered domains must be formed by composition fluctuations sufficient to generate a nucleus of a registered phase large enough for the area-dependent payoff in bulk free energy to outweigh the penalty for hydrophobic mismatch at the edge.We now study the energetics of this nucleation.For simplicity we consider the R nucleus to be compositionally uniform (thus ignoring mixing effects at the domain boundary), and use the response of the thickness profile to mismatch at the nucleus' boundary to calculate the nucleation energetics.Fig. 4a exhibits typical nucleation of registered domains.Our goal is to model the role of hydrophobic mismatch around the edge of R nuclei in determining whether they successfully grow (Fig. 4a) or decay (Fig. 4b).For the reference (zero free energy) state, we take a registered bilayer at equilibrium with no transmidplane mismatch and domains coarse enough for edges to make a vanishing energy contribution.We then assume that a dominant AR instability mode has led to an initially antiregistered state, which incurs everywhere an energy cost, relative to the registered reference state, due to the direct coupling B. We then introduce a circular R nucleus of radius R, thus removing the direct coupling energy within the domain's area but introducing a thickness mismatch around its perimeter (Fig. 5).
In the continuum limit of the lattice model, the free energy is where a is the lattice spacing and the J term captures the thickness gradient penalty from the corresponding neighbour interaction (Eq. 1) in the limit of small lattice spacing.We neglect contributions from the Isinglike interaction V .It is thought that, as hydrophobic mismatch increases, it becomes the dominant contribution to the line tension of (registered) domains [23,31] (this does not imply that thickness mismatch uniquely determines phase-transition temperature [32]).However, in Appendix B we briefly discuss an upper bound for the influence of V on the following nucleation energetics.The leaflets' compositions are reflected in their ideal thicknesses ℓ t,b 0 (r).The actual thickness profiles ℓ t,b (r) adopted by the leaflets minimise the free energy G cont .The actual thickness difference and total thickness profiles are denoted ∆(r) ≡ ℓ t (r) − ℓ b (r) and d(r) ≡ ℓ t (r) + ℓ b (r).The calculations underlying the following results are detailed in Appendix A.

A. Antiregistered background
The nucleus appears within a phase-separated AR-AR background.In practice, R domains typically form at AR-AR domain boundaries (Fig. 4).The effects of this on nucleation are briefly discussed in Appendix B, but for the following calculations we are ignoring the Ising contribution V (which provides the only source of line tension at AR-AR boundaries).Hence, the pattern and orientation of the initial AR domains is unimportant, and we assume a spatially uniform state with S lipids in the top leaflet, U in the bottom.For simplicity we assume strong compositional segregation so that the phases are approximately pure, with ideal leaflet thicknesses given by Minimisation of Eq. 5 yields for the difference in actual leaflet thicknesses.The species in each leaflet adapt their tail lengths (hence degree of tail ordering [16]) to one another's presence, balancing tail stretching with the direct coupling energy.Relaxing the assumption of strong compositional segregation would FIG. 5. Schematic structure of a registered SS nucleus introduced to an antiregistered background.The leaflet thicknesses inside and outside the nucleus are calculated from Eq. 15.The midplane (dashed) shifts discontinuously at r = R to maintain smooth outer contours (cf.[9,30]).S (light) and U (dark) lipids illustrating composition are superimposed.Primes indicate that leaflet thicknesses outside the nucleus differ from those inside.
not qualitatively affect the physics.The physical content of Eq. 7 is simply that transbilayer-mismatched lipids in an AR phase cause a free-energy cost; in strong compositional segregation this energy density is (see Eq. A5) Since no mismatches in total ideal thickness exist, d AR is uniform:

B. Introducing a registered nucleus
We consider the energy change induced by a registered nucleus of the S species (Fig. 5), with larger ideal thickness than its surroundings (the energetics are identical for a nucleus of the shorter species U ): The top and bottom leaflets now have the same composition (ideal thickness) within the nucleus, removing the inter-leaflet mismatch within the region r ≤ R and leading to an area-dependent free energy payoff via the B term of Eq. 5.The actual thickness difference profile is now The discontinuity arises from the discontinuous composition at r = R and the fact that, in contrast to the total thickness d(r), no energy in the Hamiltonian penalises variation in the thickness difference ∆(r).
If the composition interface were not sharp one would consider the coupling between composition and thickness gradients.Moreover, one could also allow for sliding into antiregistration at a registered domain boundary, which smears out thickness mismatch [2,9,30].These would not qualitatively alter the fact that the registered nucleus' boundary experiences an energy cost from thickness mismatch, but could help reduce it.
Calculating the total thickness profile after the nucleus is introduced yields Here, I n and K n are nth order modified Bessel functions of the first and second kind respectively.Their spatial dependence is controlled by a decay length which quantifies the competition between hydrophobic mismatch and stretching.Eq. 12 shows how the mismatch in ideal total thickness at r = R distorts the thickness profile both inside and outside the registered nucleus.The decay length ξ controls the lateral distance over which d nuc.(r) is perturbed (Fig. 6a).
For large ξ, the R nucleus does not reach its ideal thickness even at its centre, r/a = 0.This is an important point; for small nuclei R ∼ ξ, the effective thickness mismatch between the nucleus and its surroundings is less than its value ∆ 0 in the limit of large domain size (Fig. 6b).The thickness at the centre of the nucleus is given by Hence, small domains have smaller hydrophobic mismatch so experience a smaller effective line tension.This dependence of domain height on domain size could be observed by atomic force microscopy.Since J and κ will typically be of the same order of magnitude [1] such that ξ ∼ a, it may only be significant for very small domains.However, increasing thickness mismatch of registered domains as they grow has been reported in molecular simulation [33], as predicted here.
After introducing the nucleus, the individual thicknesses ℓ t, nuc.(r) and ℓ b, nuc.(r) of the top and bottom leaflets are given by Assuming no empty space inside the bilayer, the leaflet thicknesses fully determine the outer contours of the bilayer once the midplane position is specified.Noting that ∆ nuc.(r) (thus ℓ t, nuc.(r) and ℓ b, nuc.(r)) is discontinuous at the nucleus boundary r = R, we expect this discontinuity to occur in the midplane position to maintain smooth outer contours and minimise hydrophobic exposure [30].This can also be seen in molecular simulations containing AR domains [9].The resultant bilayer structure is shown in Fig. 5. Outside the R nucleus, the leaflets adapt their thicknesses (thus tail structure) to one another's presence (Eq. 11), i.e. the S species are shorter than ℓ S0 and the U longer than ℓ U0 .Inside the nucleus, the registered region is able to achieve its ideal thickness only once the nucleus grows large enough.

C. Energy of a registered nucleus
Given the thickness profiles before and after the nucleus is introduced, the energy required to create a registered nucleus of size R is calculated as (see Appendix A) The negative term ∝ R 2 arises from the removal of inter-leaflet tail structure mismatch over the domain's area, and the positive term arises from hydrophobic mismatch.For large nucleus size R relative to the decay length ξ, the approximations . The hydrophobic mismatch term then becomes overall linear in R, acting like a standard line tension Γ = ∆ 2 0 κξ/8a 2 .Hence, for large nuclei or large stiffness κ, ∆G behaves as in 2D classical nucleation theory (CNT).For a ≈ 0.8 nm, κ ≈ 3 a −2 k B T , ξ ≈ a and ∆ 0 ≈ a at T = 300 K, we estimate Γ ≈ 2 pN between R and AR domains or Γ ≈ 8 pN between R-R domains of different species, quite close to existing estimates for phospholipids [34].Thus, although our simplified model does not FIG.7. Energy ∆G for a registered nucleus of radius R, calculated from Eq. 16, with ∆0 = 2 a and κ = 3 a −2 kBT .We vary the decay length ξ by varying J but hold the ratio of indirect and direct couplings constant at J/B = 2 (J/B = 8 line in Fig. 8).In 2D classical nucleation theory the ratio ∆G/R would yield a straight line of negative slope.Circles mark critical radii R * .Inset: behaviour of ∆G without normalisation by nucleus size.capture all details of hydrophobic mismatch and the concomitant bilayer deformation [30], the associated energy scale is well captured.Some illustrative nucleation energy curves are shown in Fig. 7.For CNT in 2D, the ratio ∆G/R would be a straight line of negative slope.Deviation from a straight line indicates deviation from CNT.If ξ is small (weak hydrophobic mismatch compared to stiffness), the behaviour is CNT-like over most of R; moreover the critical radius R * (defined as the value where ∆G(R) is maximised) occurs deep into the CNT-like regime.For larger ξ (stronger hydrophobic length mismatch relative to stiffness) the critical radius no longer occurs in the CNT-like regime, so non-CNT effects arising from nonlocal thickness deformation influence nucleation.

D. Nucleation barrier and critical radius landscapes
Fig. 7 (inset) implies that increasing J for a fixed ratio J/B and fixed κ (equivalently, increasing ξ) first increases, then reduces, the nucleation barrier.In the limit of large κ one would expect increasing J and B to simply increase the energy scale, hence the nucleation barrier.But when J is comparable to κ, the mechanism illustrated in Fig. 6b takes over; increasing ξ ≡ 2 Ja 2 /κ reduces the effective thickness mismatch for small nuclei, thus reducing the hydrophobic penalty for nucleation and reducing the nucleation barrier.The nucleation theory can be related to the simulations by plotting landscapes of the nucleation barrier ∆G(R * ) (Fig. 8) and critical radius R * (Fig. 9).A small barrier clearly facilitates nucleation, although the role of the critical radius is more subtle [35][36][37].One important factor is that for smaller critical radii, regions randomly enriched in registration from the initial quench are more likely to be near or above the critical radius, so can grow more easily.
The predicted energetics for nucleation are consistent with the classes of kinetics identified in Fig. 3.Where the R instability mode is fastest, nucleation energetics are irrelevant since the equilibrium phases are formed immediately.If the AR mode is fastest, a successful subsequent transition to equilibrium R-R coexistence generally occurs where the predicted nucleation barrier and critical radius are smaller.As we have used a relatively large value of the mismatch parameter ∆ 0 = 2 a, we briefly mention the dependence on it [38].The critical radii are independent of ∆ 0 , because ∆ 0 represents both tail length and structure mismatch so affects both the indirect and direct coupling strengths.Thus the only change to Fig. 9 would be in the effective value of γ shown on the secondary axis, which scales as ∆ 2 0 (Eq.8).On Fig. 8, in addition to the change in effective γ, the nucleation barrier ∆G(R * ) scales as ∆ 2 0 .

VI. DISCUSSION
In this paper we have studied phase-separation kinetics in a bilayer subject to competing symmetric and asymmetric instabilities.We have modelled the nucleation of symmetric (registered, R) phases where the initial spinodal instability has led to antiregistered AR-AR coexistence (local asymmetry everywhere).We made use of a microscopic lattice model that allows study of both phase equilibria and phase-transition kinetics, through coarse-grained theory and direct simulation.For realistic parameters, the competing symmetric and asymmetric instabilities can be comparable in strength, so that changes to molecule properties can tip the balance [1].
If the fastest-growing instability mode is R, spinodal decomposition directly into equilibrium R-R coexistence takes place (class 1 Fig. 3).If the AR mode grows fastest then a metastable AR-AR state forms, and nucleation is required to reach equilibrium (class 2 Fig. 3).Unsuccessful nucleation leaves the bilayer metastably trapped (class 3 Fig.3).
The key parameters for the initial demixing and subsequent nucleation energetics are the indirect inter-leaflet coupling J via hydrophobic mismatch, favouring antiregistration, and the direct coupling B which favours registration.Physically, the hydrophobic tail length mismatch or structure mismatch would affect the effective values of J or B respectively.Both J and B couple to the stiffness κ, which determines how easily mixed lipids can change their tail length and structure to adapt to one another's presence.Here we have focused on an equimolar mixture in both leaflets such that AR-AR coexistence competes with R-R; the kinetic considerations are qualitatively similar for other overall compositions subject to competing instability modes, in which other states including three-phase coexistence (AR-AR-R or R-R-AR) can enter [2].
Nucleation implies a critical radius, which a registered domain must exceed so that the penalty of thickness mismatch at the perimeter does not outweigh the bulk free energy gain.This does not mean domains automatically become registered beyond a certain size, as is sometimes implied [5,10,39].If nucleation is energetically prohibitive, the mere presence of large coarsening antiregistered domains does not guarantee registered domains will ever form, just as assembling a large volume of supercooled water does not guarantee a supercritical ice nucleus will form.That said, it is intriguing to con-sider whether some hydrodynamic or curvature-mediated mechanism (neither of which are included here) could dynamically assist very large antiregistered domains to overcome the nucleation barrier and become registered.
In the present simulations it appears that a nucleation barrier ∆G(R * ) 5 k B T inhibits nucleation on the simulated timescale for this simulation size.In principle the nucleation rate I (per unit area per unit time) is determined by ∆G(R * ) via I = I 0 exp (−∆G(R * )/k B T ), although the unknown kinetic prefactor I 0 severely limits the quantitative utility of such a picture.We also expect the critical radius R * to determine how large random regions of inter-leaflet symmetry during the initial quench must be in order for them to grow.The critical radius also captures a size-dependence of potential relevance to cell membrane rafts and clusters, discussed below.
By examining the effect of thickness mismatch at the boundary of a nucleating registered domain, we showed that the process departs strongly from classical nucleation theory due to deformation of the thickness profile, which reduces the effective line tension of small nuclei.For small nuclei, the deformation essentially spreads over the whole nucleus, reducing the effective hydrophobic mismatch.To our knowledge this behaviour has not been noted in previous theories which tended to focus on domains large enough to treat the boundary as a straight line [30,31].The prediction is supported by a recent molecular simulation of coarsening of a quenched bilayer [33], where registered liquid-ordered and liquiddisordered domains became respectively thicker and thinner as they grew through time.
The simulation method does not include hydrodynamics, which are expected to dominate domain coarsening beyond a lengthscale ∼ 10 −6 m [40] such that the purely diffusive dynamics simulated here would no longer apply.However, this lengthscale is far beyond both the size simulated here and the predicted critical radii (Fig. 9).Note that hydrodynamics cannot influence the free-energy landscape of the system, and thus cannot change the competing metastable and equilibrium states.
The phase equilibria of "macroscopic" domains (which coarsen until limited by bilayer size) are not influenced by edge energies and, for most reasonable parameters, we predict that R phases are lower in free energy so that R-R is the equilibrium state in bulk [1], as seen in fluorescence studies of large domains [41,42].At the opposite end of the size spectrum, experiments reveal pairwise antiregistration at the single-lipid level [11,43], as also reported in simulation [8] and predicted by our mean-field theory [1].In between these regimes, edge energies can influence the equilibrium state; for example, the AR-AR coexistence in [9,10] is probably metastable in the limit of large size but could be stabilised if the simulation box is too small to accommodate a supercritical registered domain.
This has important physical consequences, because cell membrane rafts or clusters are not macroscopic.Small domain size in vivo could be due to elastic repulsion [44], hybrid lipids [45], critical fluctuations [46], active recycling [47,48] or another mechanism.In either case, it is crucial to recognise that edge energies could influence the thermodynamic preference and even stabilise an AR-AR type of domain formation, in which one leaflet's local enrichment in longer species (relative to the background in that leaflet) colocalises with a relative depletion of such species in the other leaflet.
Although the cell membrane is maintained out of equilibrium, thermodynamic driving forces can be expected to play a role.The potential importance of this sizedependence is underlined by the fact that the estimated critical radii (Fig. 9) can be of the order of putative lipid raft sizes [49].The basic biophysical question is: if a cluster of longer lipids and proteins exists in one leaflet, does it colocalise a similar cluster in the opposite leaflet to maintain transbilayer structural similarity, or does it choose shorter lipids and proteins to maintain uniform thickness?Our work implies that finite-size effects, metastable states and phase-transition kinetics can be key in determining the answer.

Appendix A: Calculations for nucleation theory
In the nucleation theory the response of the actual profiles ℓ t,b (r) to variations in the ideal thicknesses ℓ t,b 0 (r) is calculated as follows.Varying G cont (Eq.5) with respect to ℓ t and ℓ b gives a pair of differential equations that can be combined to yield independent equations for the thickness difference ∆(r) ≡ ℓ t In the above we have used the ideal thickness difference and total profiles, ∆ * 0 (r) ≡ ℓ t 0 (r) − ℓ b 0 (r) and d * 0 (r) ≡ ℓ t 0 (r) + ℓ b 0 (r).In the AR background the S and U species' lengths are not equal to their ideal values -the direct coupling B encourages equality of tail length (thus degree of ordering) across the bilayer.Eq.A1a is solved by Inserting Eq. 6 for the ideal thicknesses yields Eq. 7. The total thickness before the nucleus is introduced is uniform (Eq.9).After introducing the R nucleus, the difference in ideal leaflet thicknesses is zero within the R nucleus (∆ * 0 (r) = 0 for r ≤ R), which gives Eq. 11.The total thickness profile is found by solving Eq.A1b (which is expressible as a modified Bessel's equation of order zero).We require that the gradient vanishes at the centre of the nucleus (d d(r)/dr| r=0 = 0), and that the total thickness approaches its ideal value away from the domain (d(∞) = ℓ S0 + ℓ U0 = d 0 ).The profile d(r) and its gradient are required to be continuous at r = R in order to match the profiles inside and outside the nucleus.With ideal thicknesses given by Eq. 10, this yields Eq. 12.
To calculate the energy ∆G required to introduce the registered nucleus, we consider the initial energy G AR cont from inserting Eqs.7 and 9 into Eq.5, and the energy G nuc.
cont when the nucleus is introduced, from inserting Eqs. .(A5) ∆G out contains only positive contributions, resulting from deformation of the total thickness profile by the nucleus.∆G in contains similar deformation terms, but also a negative term resulting from the fact that there is now no inter-leaflet mismatch in the region r ≤ R. The sum ∆G in + ∆G out leads to Eq. 16.
Note that the negative term in ∆G in yields the estimated inter-leaflet mismatch energy density γ under our microscopic definition (Eq.8) [1].Near the strongsegregation regime this is a close approximation to the difference in bulk free-energy density between registered and antiregistered phases (which can be used to construct an alternative "macroscopic" definition of γ [1]), becoming exact in the strong-segregation limit assumed in Section V.However, in principle one could replace this term with the actual free-energy difference in cases where the correspondence between microscopic and macroscopic γ breaks down [1].

Appendix B: Effect of V on nucleation
In Section V we ignored the effect of V (headgroup interactions) on the nucleation energetics.This applies best where hydrophobic mismatch is the dominant source of line tension [31].Here we sketch an upper bound on the contribution of V to nucleation energetics of a registered domain.Fig. 10 shows an R domain nucleating at a boundary between AR domains (cf.Fig. 4).This relieves some AR-AR interfacial energy cost (like a particle in a Pickering emulsion).
Assuming the original AR-AR boundary is straight and the R nucleus is a circle, the length of AR-AR boundary removed by an R domain of radius R is equal to 2R.The length of R-AR boundary introduced by the nucleus is 2πR.It is extremely difficult to estimate the line tension caused by V , because i) the separated phases will not be fully pure and ii) the interface will relax its energy by smearing the composition change over some finite distance, and iii) this will be coupled to the thickness gradients.However, an upper bound can be obtained by assuming, as in the main text, that the phases are strongly-segregated and that the compositional interface is sharp.With these assumptions, the line tension of the AR-AR boundary is 2V /a, and the contribution from V to an R-AR boundary is V /a.Subtracting the AR-AR energy contribution and adding the R-AR contribution, the correction to Eq. 16 would be: Hence, the contribution of V is always to increase the nucleation energy (thus both the critical radius and nucleation barrier) of registered domains.For the parameters ranges covered by Fig. 3 (assuming V = 0.9 k B T as used in the simulations), we find that including the upper-bound correction from V increases the critical radius R * by at most ∼ 40 % for J 2 a −2 k B T , from the corresponding value on Fig. 9. (If ∆ 0 is reduced to 1 a, reflecting decreased tail hydrophobic and structural mismatch so that V is proportionally more important relative to the indirect and direct couplings, the increase can be ∼ 100 %).For ∆G(R * ), if we assume that the original R * is not altered much by including V , the associated increase in the nucleation barrier is estimated as (2V /a)(π − 2)R * .Alternatively, explicitly taking into account the change in R * resulting from including V , the increase in ∆G(R * ) in the parameter range covered by Fig. 3 is up to ∼ 100 % for J 2 a −2 k B T , from the corresponding value in Fig. 8.
We reiterate that Eq.B1 must overestimate the contribution V to nucleation energetics, possibly quite severely, so that the argument outlined here provides only a rough upper bound.
where φt(b) i = 1 if the top (bottom) of lattice site i contains an S lipid, φt(b) i = 0 if U .The species-dependent

FIG. 2 .
FIG. 2. (a) Mixed bilayer containing S and U (Lo or gel and L d -like) model species, illustrating the locally symmetric (R-R) and locally asymmetric (AR-AR) phase coexistences considered here.(b) Microscopic lattice model for coupled leaflets, which can be coarse-grained (CG) (c) to give the mean-field free-energy density f (φ t , φ b ) as a function of locally-averaged leaflet compositions, and analysed for kinetics of domain formation with the inclusion of gradient costs for domain boundaries [1].The lattice model can also be directly simulated.

FIG. 6 .
FIG. 6.(a) Calculated profile of total thickness d nuc.(r) relative to the AR background, when the nucleus is introduced.The ideal thickness relative to the surroundings is ∆0 = 2 a. Arrow: decreasing decay length ξ.(b) d nuc.(r) relative to the AR background and normalised by its ideal value inside the nucleus, ∆0.Arrow: increasing nucleus size.

2 ,
11 and 12 instead.The energy for introducing the nucleus is then ∆G = G nuc. cont − G AR cont .(A3) Since Eqs.11 and 12 split at the nucleus boundary r = R, we evaluate separate contributions to ∆G outside and inside the nucleus and find ∆G out = .(r) − d 0 )

FIG. 10 .
FIG.10.Schematic of an SS R nucleus of radius R forming at the boundary between AR-AR domains.