Pattern formation by curvature-inducing proteins on spherical membranes

Spatial organisation is a hallmark of all living cells, and recreating it in model systems is a necessary step in the creation of synthetic cells. It is therefore of both fundamental and practical interest to better understand the basic mechanisms underlying spatial organisation in cells. In this work, we use a continuum model of membrane and protein dynamics to study the behaviour of curvature-inducing proteins on membranes of spherical shape, such as living cells or lipid vesicles. We show that the interplay between curvature energy, entropic forces, and the geometric constraints on the membrane can result in the formation of patterns of highly-curved/protein-rich and weakly-curved/protein-poor domains on the membrane. The spontaneous formation of such patterns can be triggered either by an increase in the average density of curvature-inducing proteins, or by a relaxation of the geometric constraints on the membrane imposed by the membrane tension or by the tethering of the membrane to a rigid cell wall or cortex. These parameters can also be tuned to select the size and number of the protein-rich domains that arise upon pattern formation. The very general mechanism presented here could be related to protein self-organisation in many biological processes, ranging from (proto)cell division to the formation of membrane rafts.


I. INTRODUCTION
Spatial organisation into inhomogeneous patterns is an essential feature of living organisms, from the macroscale to the cellular level. In the later case, organisation of the plasma membrane and the cytoplasm into specialised domains is more commonly referred to as cell polarity. [1,2] This spatial organisation of the cell is necessary in order to coordinate important processes such as cell division, differentiation, or directed cell migration.
As early as in 1952, Turing realised [3] that very simple systems that are initially in a spatially homogeneous state can spontaneously self-organise into spatially inhomogeneous patterns. However, it is generally believed [1,2] that the generation of polarity in cells is the result of a tightly-controlled orchestration involving complex signalling networks and active processes such as the reorganisation of the cellular cytoskeleton. Nevertheless, active systems such as the cytoskeleton have been shown to undergo simple pattern formation, [4] and there also exist cells for which polarisation is presumably not generated by the cytoskeleton. [5][6][7][8][9][10][11] The underlying mechanisms in these systems are however not well understood.
Very recently, [12] a system was identified in which cell polarisation appears to be controlled by a relatively simple pattern-formation mechanism. In the coccal bacterium Staphylococcus aureus, essential proteins involved in lipid metabolism were seen to distribute in inhomogeneous spatial patterns, that could be explained by a model that considers the dynamics of curvature-inducing proteins on a spherical membrane. However, the model first introduced in Ref. 12 is very general, and we expect that it might be able to describe the formation of protein patterns on the surface of other types of cells, as well as in model systems consisting of lipid vesicles and proteins. In this work, we will explore in full generality and detail the predictions of such a model.
The basic idea behind the model is presented in figure 1. A closed, initially spherical membrane contains proteins that impose a spontaneous curvature C p on the membrane (in general, the proteins might be attached to the membrane from the cytoplasmic or the exoplasmic sides, or they might be transmembrane proteins embedded in the membrane). [13,14] If the proteins did not induce any curvature, a random, homogeneous distribution of proteins would be favoured by thermal fluctuations, that is, entropic forces (in the absence of direct attractive protein-protein interactions). However, if the curvature induced by the proteins is large enough, bending contributions to the free energy of the system can lead to an effective attraction between proteins and to the formation of spatially inhomogeneous patterns in protein distribution and membrane curvature. The details of membrane-mediated protein-protein interactions have been thoroughly studied in the past. [15][16][17][18] Furthermore, we will consider the possibility of geometric constraints on the membrane, such as the tethering of the membrane to a rigid cell wall/cortex or the existence of a membrane area reservoir at non-zero tension. Interestingly, it was recently shown that solid particles such as proteins can sense the local membrane curvature imposed by geometric constraints on the membrane. [19] repulsion attraction C p FIG. 1: The proteins (yellow) impose a spontaneous curvature C p on the membrane (blue). Depending on the interplay between curvature, entropic, and membrane tethering/tension forces, proteins might repel, resulting in a spatially homogeneous spherical membrane (turquoise), or they might attract, leading to the spontaneous formation of inhomogeneous patterns of membrane curvature and protein density.
Here, we have found that, in realistic situations, spontaneous pattern formation can be induced either by an increase in the surface density of curvature-inducing proteins, or by a decrease in the strength of the geometric constraints on the membrane. Furthermore, these two parameters can also control the size and number of protein-rich (highly curved) and protein-poor (weakly curved) domains. These mechanisms could be exploited by cells in order to trigger spatial organisation of the plasma membrane on demand, and could in principle be replicated in artificial model systems.
The paper is organised as follows. In section II, we present the continuum model for the energetics and dynamics of the system, and examine the linear stability of the dynamical equations for the shape of the membrane and the protein density distribution. In section III, we explore spontaneous pattern formation in the system as a function of all relevant parameters. Finally, in section IV we discuss the applicability and consequences of our results in real biological or biomimetic systems.

A. Energetics
We will adopt a continuum elastic model of a closed membrane, which might represent a model vesicle or a biological cell, and study the stability of spherical shapes to perturbations in the presence of curvature-inducing proteins that decorate the membrane. The shape of a quasi-spherical membrane can be written in spherical coordinates as R(θ, φ) = R[1 + u(θ, φ)]r, where R is the radius of the unperturbed sphere, u(θ, φ) is a scalar function that describes the deviations from the sphere, andr is the radial unit vector, see figure 2. The distribution of proteins on the membrane can be described in a similar way, with the surface number density ρ(θ, φ) = ρ 0 [1 + ψ(θ, φ)]. Here, ρ 0 is the average protein number density, i.e. ρ 0 = N/4πR 2 if N is the total number of proteins on the membrane, and the function ψ(θ, φ) represents the deviations from a homogeneous distribution of proteins.
We will assume that each protein covers a patch of membrane of area a 0 , and imposes a spontaneous curvature C p on the membrane, see figure 1. The bending free energy of the membrane can then be written within the spontaneous curvature model [20][21][22][23] as where κ is the bending rigidity of the membrane, and C is the local membrane curvature, where C 1 and C 2 are the two principal curvatures. The second term inside the integral represents the simplest possible coupling between protein density and local curvature. It can also be interpreted as a position-dependent spontaneous curvature C 0 (θ, φ) ≡ C p ρ(θ, φ)a 0 , which varies from C 0 = 0 in the absence of proteins, with ρ = 0, to C 0 = C p for full coverage of proteins, with ρ = 1/a 0 . The local membrane curvature C(θ, φ) can be written explicitly as a function of u(θ, φ), as described in ref 24.
Besides the bending contributions to the free energy, we need to take into account the entropic contributions due to the mixing and density fluctuations of the proteins. To lowest whereas the protein distribution is described by a scalar function ρ(θ, φ) represented by the colourcoding, e.g. yellow and blue could correspond to high and low protein density, respectively.
order, this contribution to the free energy can be incorporated as Here, χ and ξ are the compressibility and the correlation length of the protein density fluctuations, respectively. The first term in the integral penalises the creation of interfaces between high protein density and low protein density regions, whereas the second term penalises deviations from a homogeneous protein distribution.
We will also consider the effect of the tethering of the membrane to a cell wall or actomyosin cortex, by including a harmonic confinement potential of the form where k te is an effective spring constant per unit area, which in general may include contributions from specific interactions (i.e. proteins that directly link the membrane to the wall/cortex) as well as non-specific interactions such as steric repulsion, van der Waals attraction or electrostatic attraction/repulsion. Within this effective description, the cell wall/cortex is taken to be spherical and rigid (i.e. much more rigid than the membrane), and k te penalises deviations of the membrane position from the (optimal) equilibrium membranewall distance.
Lastly, we consider the possibility that the membrane is connected to a membrane area reservoir at constant membrane tension. A constant membrane tension is typical of biological cells, [25,26] and can be mimicked in model vesicle systems by the use of micropipette aspiration. The contribution of a membrane tension σ to the free energy is The total free energy can finally be written as the sum of these four contributions, with with the free energy density In addition, we will explicitly impose constraints on the volume enclosed by the membrane (representing osmotic balance), so that as well as on the total number of proteins N on the membrane, so that at all times.

B. Dynamics
The effective force exerted on the membrane in the radial direction will be balanced by a frictional force, leading to a dynamical equation for the shape of the membrane as a function where L u is a transport coefficient corresponding to the membrane mobility.
On the other hand, the dynamical equation describing the diffusion of the proteins on the membrane can be written in the form of a continuity equation with a current density J = −L ψ ∇µ, where L ψ is another transport coefficient and µ = δF/δψ(θ, φ) is the chemical potential. Putting all together, the dynamical equation for the protein density becomes The Laplacian operator on a sphere can be written as This operator is diagonal in the basis of spherical harmonics Y m (θ, φ). In particular, it

C. Linear stability analysis
To leading order in u and ψ, and taking into account the constraints (7-8) on the enclosed volume and total number of proteins on the membrane, we can write equations (9) and (11) as and We can write the solutions u(θ, φ, t) and ψ(θ, φ, t) as a sum of spherical harmonics, which provide a complete set of orthogonal functions on the sphere, so that where u m and ψ m are the amplitudes of the corresponding modes, and we have = 0, 1, 2... and |m| ≤ . However, the constraints (7-8) on the enclosed volume and total number of proteins on the membrane imply that the zero-amplitudes u 00 and ψ 00 cannot be varied independently. Explicitly imposing these constraints results in expressions for u 00 and ψ 00 as a function of the squared amplitudes of all modes Equations (17) and (18) imply that u 00 and ψ 00 are a function of the higher-order amplitudes, and furthermore, that they are of quadratic order (they are equal to a sum of u 2 m and u m ψ m terms). For this reason, the u 00 and ψ 00 terms are negligible to linear order, and we can rewrite (16) as Inserting (19) into (14) and (15), we can rewrite the dynamical equations as separate equations for each of the ≥ 1 modes. Introducing a rescaled time variable as well as dimensionless parameters the equations become and The solutions to (22) and (23) will have the form Inserting these solutions back into (22) and (23), and setting the determinant of the coefficients to zero, we can obtain an equation for the growth rates λ of the characteristic modes of the system, which reads and where we have defined the parameter The two characteristic modes of the system given by the solutions to (25) can finally be written as Because b in (25)  The physical significance of the five dimensionless parameters is the following. The parameter W represents the protein-induced spontaneous curvature, and increases both with the average density ρ 0 of proteins on the membrane and with the characteristic spontaneous curvature C p of these proteins. The parameter K represents the strength of the confinement of the membrane by its interaction with the rigid cell wall/cortex. The parameter T represents the magnitude of the membrane tension. The parameter P compares the correlation length of the protein density fluctuations to the size of the cell or vesicle. Given that correlation lengths are typically of the order of nanometers whereas cell or vesicle sizes are of the order of micrometers, P will generally be small, and will decrease or increase with increasing or decreasing cell/vesicle size, respectively. Finally, the parameter M compares the typical timescale of the changes in membrane shape (R 2 /L u κ) to that of changes in protein distribution (χR 2 /L ψ ). Importantly, we note that all five dimensionless parameters are always positive.

III. RESULTS
A positive value of the mode amplitude λ + implies that fluctuations of this mode will grow instead of decaying, and therefore modes with λ + > 0 are unstable. If, by small changes in one of the system parameters W , K, T , P , or M , one of the modes λ + switches from having a negative value to having a positive value, the system will exhibit spontaneous pattern formation. In the following, we will explore the conditions under which spontaneous pattern formation occurs.
First of all, we note that, as described above, the = 0 mode cannot vary independently as it is fixed by the constraints on the enclosed volume and total number of proteins, see (17)(18).
Furthermore, by substituting = 1 in (29), we find the mode amplitudes −K and −2M (1 + 2P ), which can never be positive, implying that the = 1 mode can never become unstable.
It can, however, become marginally stable in the particular case of K = 0, i.e. in the absence of tethering to the cell wall. This reflects the fact that = 1 deformations of the membrane shape are equivalent to spatial translations, and that the curvature energy of the membrane is invariant to such translations. The presence of the cell wall, however, breaks translational invariance. All things considered, instabilities and therefore spontaneous pattern formation can occur only for higher modes ≥ 2, which we will discuss below.
The larger solution λ + of (25) will be positive, with λ + > 0, if and only if c < 0. Using the definition of c in (27), this condition can be rewritten as which serves as a definition of W , the critical value of the parameter W above which mode becomes unstable. Going back to the definition of W in (28), the inequality (30) implies that an increase in the average density ρ 0 of curvature-inducing proteins beyond a critical density will trigger an instability with spherical harmonic mode in both the shape and protein distribution of the membrane. Furthermore, the critical protein density that is needed to trigger an instability decreases with increasing protein spontaneous curvature C p . Importantly, we note that the critical value W is independent of the parameter M , and therefore depends only on three parameters, P , K, and T . In fact, the parameter M drops out of all relevant equations in the following, so that pattern formation in the system turns out to be governed by only four dimensionless parameters: W , K, T , and P . This is a consequence of the fact that M is a mobility parameter that related the timescale of changes in membrane shape to that of changes in protein distribution, and as such it only affects the dynamics of the system.
Alternatively, the instability condition c < 0 can be written as which define K and T , the critical values of K and T , respectively, below which mode becomes unstable. Going back to the definitions of K and T in (21), the inequalities (31) and (32) respectively imply that the shape and protein distribution instability can also be triggered by a decrease in the tethering strength of the membrane to the cell wall/cortex, or by a decrease in the membrane tension. Once again, we note that the critical values K and T are independent of the parameter M .
As outlined in the previous two paragraphs, the parameters that could presumably be actively controlled by a biological cell or tuned in experiments with model vesicles are W , i.e.
the density of proteins on the cell surface, K, i.e. the tethering strength of the membrane to the cell wall/cortex, and T , the membrane tension. The parameter P , on the other hand, represents the correlation length of the protein density fluctuations, i.e. the typical distance at which proteins can sense each other, and will in general be fixed for a given system. It therefore makes sense to explore the behaviour of the system when W , K, and T are varied for a fixed value of P .
Using (30), in figure 3 we have plotted the lines W = W (K) for ≥ 2, using T = 0 (i.e. negligible membrane tension) and three different values of P , namely P = 0.1, 0.02, and 0.005. For a vesicle/cell of radius 1 µm, these values of P would correspond to correlation lengths of ξ = 320 nm, 140 nm, and 70 nm, respectively. In the region of low W and high K, depicted in grey, the spherical state with a homogeneous protein distribution is stable.
As W is increased from low values, the system will hit the instability of the first unstable mode, with a given value of which will depend on the value of K. Alternatively, if K is decreased from high values, the system will also hit the instability of the first unstable mode with a given which will depend on the value of W . The higher the value of K, the higher the value of of the first unstable mode as W is increased. Similarly, the higher the value of W , the higher the value of of the first unstable mode as K is decreased.
There are important differences in the way in which W and K act to trigger pattern formation. Independently of the value of K, and even for K = 0, a sufficiently high W will always lead to pattern formation. On the other hand, a decrease in K can only lead to pattern formation if W is above the critical value W (K = 0). Furthermore, we note that of K lead to the instabilities of higher-order modes with > 2. As P is decreased, as in (b), the first unstable mode at vanishing K is now = 3: the mode = 2 is not the first unstable mode for any value of K. When P is decreased even further, as in (c), = 4 becomes the first unstable mode at vanishing K, and neither = 2 nor = 3 are the first unstable modes for any value of K. This trend continues as P is decreased further, with progressively higher order modes becoming the first unstable mode at vanishing K. Moreover, we note that, as P is decreased, the critical value of W above which pattern formation occurs moves closer and closer to W = 1.
In figure 3 we have explored the stability behaviour of the system as a function of W and K, for fixed vanishing tension T = 0. Considering a fixed non-zero tension T > 0 leads to the same qualitative behaviour of the system as a function of W and K. Furthermore, the behaviour of the system as a function of W and T for fixed K is qualitatively identical to that as a function of W and K for fixed T , leading to instability lines analogous to those in figure 3. We thus omit these results for the sake of brevity.
As just described, in order to characterise the system, it is particularly important to identify the first unstable mode when W is increased, that is, the mode with smallest W for given values of P , K, and T , which we will denote as * W . The critical value of W above which the first unstable mode becomes unstable is then W * ≡ W * W ≡ min (W ). The boundaries between the regions in the three-dimensional (P, K, T ) parameter space in which modes and + 1 are the first unstable mode for increasing W can be obtained from the inducing proteins) and high K (i.e. strong confinement of the membrane due to tethering to the cell wall/cortex), the spherical homogeneous state is stable (grey region). As W is increased or K is decreased, the system will hit an instability with a given value of . If W or K keep increasing or decreasing, respectively, they will hit the instabilities of further modes. The higher the value of K or W , the higher the value of of the first unstable mode as W is increased or K is decreased.
The parameter P represents the correlation length of the protein density fluctuations.
condition W = W +1 , which can be written explicitly using (30) as In figure 4, we have used equation (33) to explore pattern formation in (a) the (P, K, T = 0) plane and (b) the (P, K = 0, T ) plane. For any point in (P, K, T ) space, we can obtain the critical value W * above which pattern formation occurs, using (30). This information is also colour-coded in figure 4. Several important observations can be made: (i) Once again, we see that K and T have qualitatively similar effects in pattern formation. (ii) Both an increase in K or T , as well as a decrease in P lead to increasingly higher-order modes being the first unstable mode. (iii) In most regions of the parameter space, the critical value W * above which pattern formation occurs is very close to 1. The only exception is the region of P ≈ 1 and large K or T , in which W * can be much larger than one.
A particularly important case, with regards to its experimental relevance, is that of a Alternatively, we could ask ourselves what is the first unstable mode * K when K is decreased for given values of P , W , and T , or equivalently, the mode with largest K for given P , W , and T . The critical value of K below which the first unstable mode becomes unstable is then K * ≡ K * K ≡ max (K ). The boundaries between the regions in the threedimensional (P, W, T ) parameter space in which modes and + 1 are the first unstable mode for decreasing K can be obtained from the condition K = K +1 , which can be written explicitly using (31). The resulting (P, W ) stability diagram for the particular case of T = 0 is shown in figure 5(a). In the same way, we can find the first unstable mode * T when T is decreased for given values of P , W , and K, with a critical value given by T * ≡ T * T ≡ max (T ), and boundaries in the (P, W, K) parameter space given by T = T +1 , which can be written explicitly using (32). The resulting (P, W ) stability diagram for the particular case of K = 0 is shown in figure 5(b).
Once again, we find that K (the strength of the tethering of the membrane to the cell wall/cortex) and T (the membrane tension) behave in a qualitatively similar way. As expected from figure 3, an instability can only occur for decreasing K (or T ) if W is sufficiently high. This minimum value of W required for pattern formation approaches W = 1 for small P . Indeed, figure 5 illustrates very clearly a striking feature of the system: for low values of P (which are the most typical given that the correlation length of protein density fluctuations ξ is normally much smaller than the membrane radius R), values of W only slightly above 1 can lead to the instability of modes with very high when K or T are decreased. This is evidenced by the high density of boundary lines in the region of low P 1, W 1.

A. Estimation and control of model parameters in real systems
We have shown above that pattern formation in a spherical membrane containing curvature-inducing proteins is controlled by the four dimensionless parameters W , K, T formation occurs always stays in the proximity of W * ≈ 1, except in the extreme case of very high K (or T ) and P ≈ 1 simultaneously, see figure 4. Going back to the definition of W in terms of the dimensionful system parameters in (28), we see that the requirement W ≥ W * ≈ 1 implies that κχρ 2 0 a 2 0 C 2 p 1. Here, κ is the bending rigidity of the membrane, χ is the compressibility of the protein density fluctuations, ρ 0 is the average density of curvature-inducing proteins on the membrane, a 0 is the lateral area of a single protein, and C p is the protein spontaneous curvature. Furthermore, in the limit of low protein density ρ 0 a 0 1, the compressibility χ can be approximated by χ ≈ 1/(k B T ρ 0 ), [27] where k B is Boltzmann's constant and T is the temperature. In this limit, the requirement for pattern formation thus becomes (κ/k B T )ρ 0 a 2 0 C 2 p 1. In general, the protein spontaneous curvature C p will be of the order of the (inverse) characteristic length of the protein √ a 0 , so that we can take C 2 p a 0 ≈ 1. We finally conclude that pattern formation typically occurs for average protein densities satisfying Here, ρ 0 a 0 is simply the dimensionless area fraction of membrane covered by the protein.
Typical values of the bending rigidity of membranes range from 10 to 100 k B T , leading to a critical protein coverage of the order of ρ 0 a 0 ∼ 0.01-0.1. Importantly, the range obtained selfconsistently validates the low protein density assumption made above. Furthermore, such coverages are within the range achievable both in biological cells as well as in model vesicles.
In this picture, a biological cell could up-or down-regulate the expression of the curvatureinducing protein in order to switch between patterned and non-patterned conformations, see Let us now turn to the dimensionless parameters K ≡ k te R 4 /κ and T ≡ σR 2 /κ, which both act as geometric constraints on the membrane: K represents the confinement of the membrane due to its interaction/tethering to the cell wall or cortex, whereas T represents the membrane tension, which acts to minimise the cell membrane area. It is interesting to note that, while model membranes such as Giant Unilamellar Vesicles show clear shape fluctuations due to thermal excitation of bending modes, [28] eukaryotic cells or bacteria do not show such fluctuations. The latter is an indication that, in such systems, membrane confinement and tension must overpower bending, and consequently that in these systems The typical tension of cellular membranes, on the other hand, has been extensively measured for different cell types, and can range from σ = 3 pN/µm for epithelial cells up to about σ = 300 pN/µm for keratocytes. [25,26] Using the range σ = 3-300 pN/µm for the membrane tension, together with the estimates κ = 10-100 k B T for the bending rigidity of the membrane and R = 1-10 µm for a typical cell radius, we obtain values of T ranging from 7 to 7 · 10 5 . Cells can actively regulate their own tension in order to maintain homeostasis. [25] In this way, cells could switch between patterned and non-patterned conformations by actively decreasing or increasing the tension of their plasma membrane, see figure 5(b).
Furthermore, triggering of pattern formation via a decrease in membrane tension could be explored in experiments using model Giant Unilamellar Vesicles aspirated by micropipettes, which allows direct experimental control over the membrane tension.
Let us finally examine the dimensionless parameter P , defined as P ≡ ξ 2 /R 2 , where ξ is the correlation length of the protein density fluctuations and R is the radius of the cell or vesicle. The correlation length ξ is a measure of the distance at which proteins or protein clusters can sense each other, typically via membrane-mediated interactions in the absence of other long-ranged interactions. Previous work [15][16][17][18] has shown that the typical length scale of membrane-mediated interactions is the size of the curvature-inducing element itself, so that we can use an estimate of ξ = 10-20 nm. On the other hand, the radius of cells or cellular compartments, as well as of model vesicles, can range between R = 100 nm and 10 µm. With this, we find a range of P ∼ 10 −6 -10 −2 . In this range of values with P 1, as described above, pattern formation is tighly controlled by the number of curvature-inducing proteins, with an instability occurring as soon as W 1. Moreover, the value of P directly controls the -order of the first unstable during pattern formation, and as a consequence controls the typical size of the protein-rich, highly-curved domains. The consequences of this fact are discussed in the following section.
B. Biological relevance

Cell division
Cell division requires polarisation of the cell, so that the spherical symmetry of the cell is broken, leading to two identifiable poles as well as an equatorial line. Spontaneous pattern formation via an instability due to the presence of curvature-inducing proteins, as described here, provides a simple mechanism for such a symmetry breaking. This occurs for the mode = 2 of the instability, see figure 6(a), which can be the first unstable mode as long as P > 1/18 0.056, as determined from (34) and displayed on figures 4 and 5. For a typical cell size of R = 1 µm, such values of P would correspond to a correlation length ξ for protein density fluctuations larger than 230 nm. This value appears too high for a typical protein, given that the correlation length is expected to be of the order of the protein size, i.e. a couple of tens of nanometers. It could, however, be a plausible value for protein clusters, composed of a few tens of proteins with lateral sizes of the order of 100 nm. If such clusters arose by a separate mechanism, a curvature-instability such as the one described here could lead to cell polarisation.
Several proteins, many of them related to cell division, are known to preferentially localise at the poles of bacterial membranes. [30][31][32][33][34] We note, however, that the generic mechanism proposed here is distinct and unrelated to the well-studied Min system, which serves to localise the FtsZ protein ring in rod-shaped bacteria. [10,11] The Min system involves both membrane bound as well as cytosolic components, and locates the bacterial equator via an oscillatory mechanism. The mechanism proposed here might however explain why FtsZ proteins can spontaneously self-assemble on vesicles, even in the absence of the Min system. [8,9] In eukaryotic cells, cell division is mediated by the cytoskeleton, in particular by the mitotic spindle and the cleavage furrow. The mechanism underlying the initial positioning of this cell division apparatus is however not well understood at the molecular level. [35,36] Even if the mechanism described here did not play a direct role in cell division, it provides a generic pathway for symmetry breaking and the initiation of division of spherical membranes into two equally sized daughters, using only a minimal number of ingredients.
As such, it could serve as a plausible mechanism for the division of protocells, as well as of synthetic cells in bottom-up synthetic biology. [37][38][39] 2. Large-scale protein organisation Going beyond = 2, the mechanism described here also predicts pattern formation with modes of intermediate , e.g. in the range = 5-10. Formation of such patterns would imply protein-rich, strongly-curved clusters with sizes on the order of 1/5 to 1/10 of the cell size, see figure 6(b). Such patterns have been observed in L-form [6,7] bacteria, as well as in coccal bacteria. [5] The patterns formed by PlsY and CdsA proteins (both essential to lipid metabolism) in Staphylococcus aureus, in particular, show a striking coupling between protein density and membrane curvature. [12] As seen in figures 4 and 5, and determined from (34), modes with ≤ 10 are expected for P > 1.4 · 10 −4 which, for a typical cell size of R = 1 µm, would correspond to correlation lengths ξ > 10 nm, well within the biologically plausible range.

Nano-sized membrane rafts
The existence of protein-rich raft domains in the plasma membrane was controversial for some time, partly due to a conflation between the macroscale fluid-fluid phase separation observed in model lipid membranes with the observation of rafts in living cells. [40] Nevertheless, it is currently accepted that rafts are dynamic, fluctuating assemblies of proteins with sizes on the order of tens of nanometers. [40][41][42][43] The precise physical mechanism behind raft formation, however, is still a matter of debate. Currently proposed theories include that rafts are compositional fluctuations near the critical point of fluid-fluid phase separation in lipid membranes, [44] or that the actin cortex underlying the plasma membrane acts as a 'picket-fence' which inhibits the lateral diffusion of proteins and promotes the formation of nano-scale aggregates. [45] The model that we have presented here predicts that, under biologically reasonable parameters, curvature-inducing proteins can spontaneously self-organise into patterns that may be built from spherical harmonics with very high-order -modes, with 1. As a consequence, in such cases the typical size of the protein-rich domains (which goes as ∼ R/ ) will be much smaller than the cell size R, leading to domain sizes on the order of tens of nanometers for a micron-sized cell, see figure 6(c) for an example with = 100. It is therefore tempting to speculate that the mechanism presented here might also be connected to the existence of such nano-scale protein-rich rafts. Indeed, let us use the quantitative estimates of parameters obtained above, with typical values of membrane bending rigidity κ = 10 k B T , correlation length of protein density fluctuations of ξ = 10 nm, tethering strength of the membrane to the cell cortex k te ≈ 2.5 · 10 6 k B T /µm 4 , and membrane tension σ = 30 pN/µm.
For a small cell of radius 1 µm, we can calculate our dimensionless parameters as P = 10 −4 , K = 2.5 · 10 5 , and T = 7 · 10 2 . Using these values in equation (33), we expect the first unstable mode to be ≈ 50. This would correspond to a typical domain size R/ ≈ 20 nm.
For a larger cell of radius 10 µm, we calculate P = 10 −6 , K = 2.5·10 9 , and T = 7·10 4 . Using these values in (33), we find ≈ 500 for the first unstable mode, once again corresponding to a typical domain size R/ ≈ 20 nm. corresponding to the formation of nano-sized protein-rich membrane rafts. The mode = 100 can be triggered if P > 1.9 · 10 −8 .

C. Summary
To summarise, we have explored in detail pattern formation in spherical membranes that contain curvature-inducing proteins. Pattern formation arises from the interplay between membrane curvature energy, protein density fluctuations, and geometric constraints such as membrane tension and confinement forces due to the tethering of the membrane to the cell wall/cortex. We have shown that pattern formation in this system is controlled by just four dimensionless parameters, W , K, T , and P , defined in (21) and (28). These parameters represent the number of curvature-induced proteins on the membrane, the confinement of the membrane due to the cell wall/cortex, the membrane tension, and the correlation length of protein density fluctuations, respectively. In most circumstances, pattern formation is expected to occur as the result of an increase in the average surface density of proteins (i.e. the total number of proteins on the membrane surface), or of a relaxation of the geometric constraints on the membrane due to membrane tension or membrane tethering to the cell wall/cortex. The patterns that arise consist of protein-rich, highly-curved domains that alternate with protein-poor, weakly-curved domains. We hypothesise that spontaneous pattern formation as described here might be exploited by biological cells as a way to regulate their geometry in situations that require spatial organisation, symmetry breaking or polarisation of the cell, using only a minimal number of ingredients.