On Gaussian curvature and membrane fission

We propose a three-dimensional mathematical model to describe dynamical processes of membrane fission. The model is based on a phase field equation that includes the Gaussian curvature contribution to the bending energy. With the addition of the Gaussian curvature energy term numerical simulations agree with the predictions that tubular shapes can break down into multiple vesicles. A dispersion relation obtained with linear analysis predicts the wavelength of the instability and the number of formed vesicles. Finally, a membrane shape diagram is obtained for the different Gaussian and bending modulus, showing different shape regimes.

We propose a three-dimensional mathematical model to describe dynamical processes of membrane fission. The model is based on a phase field equation that includes the Gaussian curvature contribution to the bending energy. With the addition of the Gaussian curvature energy term numerical simulations agree with the predictions that tubular shapes can break down into multiple vesicles. A dispersion relation obtained with linear analysis predicts the wavelength of the instability and the number of formed vesicles. Finally, a membrane shape diagram is obtained for the different Gaussian and bending modulus, showing different shape regimes.
Vesicle formation is a fundamental process in many biological systems, e.g., the Golgi apparatus 1,2 , the synaptic system 3,4 , or enveloped viruses 5,6 . The Golgi apparatus constantly releases transport vesicles filled with proteins that are carried to other parts of the cell. In the synaptic nerve terminals, vesicles are filled with neurotransmitters and released by exocytosis 3 . Moreover, many viruses are enveloped by a lipid membrane which mediates the fusion of the virus with the host cell membrane; some examples are HIV-1, herpesviruses, the Ebola virus 5 , and coronaviruses 7 like SARS-CoV-2.
There is a fair amount of research on cellular membrane deformation and morphology 8,9 . However, there is little understanding when topological transitions occur and Gaussian curvature plays a role [8][9][10] . The Gauss-Bonnet theorem states that the integral of the Gaussian curvature over a surface is proportional to the surface Euler characteristic 11 . This assures that the Gaussian curvature term is topologically invariant. This leads to the term being ignored for homogeneous systems. However, it is fundamental for topological transitions like fusion or fission. Lipid bilayers exhibit different stable configurations, depending on the values of the Gaussian and bending energetic moduli. There are no direct experimental measurements of the Gaussian modulus, although a method has been proposed recently 12 . Indirect measurements give a negative value of about κ ′ ≈ −15K B T 13 , and molecular dynamics simulations give similar results 14 . The negative sign implies that the energetic term of the Gaussian curvature favors fission, since fission increases the Euler characteristic. The most common fission event is the formation of a vesicle.
In all cases what seems a necessary requirement for fission is a large membrane curvature on the area where a vesicle is to be generated, which can be modeled with a spontaneous curvature. The final fission of the vesicle is often mediated by very specific proteins, although this is not always the case. Large spontaneous curvature can suffice to produce fission [15][16][17] . This can be accomplished by interactions among membrane-bound proteins 15 or by a osmotically induced pearling instability 16 .
There has been extensive research for protein mediated fission. The dynamin superfamily 8,[18][19][20] or the ESCRT machinery 5,21,22 are two different sets of proteins that mediate in budding and fission. For example, the dynamin Drp1 is considered a major component in mitochondrial division; other dynamics are considered to either help or be necessary for the fission process. Previous numerical work 23,24 on budding or fission at mesoscopic scales are based on the evolution of the membrane up to the instant prior to fission. Pearling of tubes has also been studied with the inclusion of the Gaussian energy 25,26 , and some regimes where the Gaussian modulus makes the tube topologically stable are studied in those works. In some instances, on the other hand, the neck which connects the vesicles being formed may be stabilized by lateral segregation of membrane components [27][28][29][30][31] .
We have developed a three-dimensional model to study the dynamical evolution of a membrane, including not only the mean curvature energy term but also the Gaussian energy term. This is done by a phase-field methodology, which has been used to study a variety of systems based on the Helfrich theory for cellular membranes [32][33][34][35][36][37][38][39][40][41] . Our description is mesoscopic considering a two dimensional diffuse interfase, in contrast to the microscopic view, which describes the rearrangement of bilayers 13 . The phase field approach allows to study not only the www.nature.com/scientificreports/ equilibrium shapes, but also the dynamics of the formation of vesicles, in the same spirit as the Landau-Ginzburg framework.
Numerical integration results show that without the inclusion of the Gaussian curvature the system exhibits a pearling instability, which has been already observed in many simulations 33,42,43 and experiments 44,45 . Pearling happens in membranes due to the spontaneous curvature, and the Gaussian curvature may lead to the fission of the pearls. Our phase field model explores the fission events where the Gaussian curvature is relevant, and the results compare extremely well with the ones observed in experiments.

Methods
The model. Phase-field approaches are suitable to model the dynamics of membranes that change their shape under certain conditions [32][33][34][35][36][37][38][39][40] . As the Gaussian curvature is an intrinsic property of the surfaces, no matter their dimensions or the metric relations that can be exerted within them 46 , it certainly has an influence on the way membranes can change their shape. However, Gaussian curvature has not been considered in phase field models because its contribution to the energy is a topological invariant. Nonetheless, it has to be included in the study of membrane dynamics because the whole curvature energy depends on it and not only on the usual bending rigidity modulus due to the mean curvature. Minimization of the entire curvature energy allows to describe, not only the shape changes that membranes must acquire, but some processes that involve a change of genus, such as fission and fusion.
Starting off from the expression for the bending energy due the mean curvature H and the Gaussian curvature K we have where κ and κ ′ are the bending modulus and Gaussian modulus, c 0 is the spontaneous curvature and the integral is calculated over the whole membrane surface Ŵ . In here we include the influence of the Gaussian curvature K in the dynamics of the system in order to model situations in which the genus of the membrane changes.
In terms of the principal curvatures of the surface, the free energy of Eq. (1) can be written as follows: where R i are the principal radii of curvature. Note that the mean curvature has dimensions of inverse length, and the Gaussian curvature has dimensions of inverse length squared.
A phase-field model of the Cahn-Hilliard type can be defined from Eq. (1), as in Ref. 32 , in which the authors express the free energy F of the system as an expansion of powers of a smooth scalar field φ : � ⊂ R 3 −→ R , that acts as an order parameter: Assuming that the system is isotropic and homogeneous, and given that the order parameter must have two stable phases, the energy density L can be written as L = (�[φ]) 2 , where the functional is defined as The parameter ε represents the width of the interface between the two phases. One of the stable phases, typically defined by φ = 1 , corresponds to the interior of the volume delimited by the membrane located at φ = 0 , whereas φ = −1 represents the outer environment.
It can be demonstrated that Eq. (3) with L = (�[φ]) 2 is equivalent to the expression of the bending energy of the surface in terms of the mean curvature H 32 . This would model the first half of Eq. (1), only the Gaussian term remains to be portrayed in a phase field approach.
The Gaussian curvature term can be defined in terms of the curvature tensor Q αβ of the surface as which in turn can be defined in terms of the gradients of the order parameter φ as 47 , where ∂ α = ∂/∂x α , and ∂ αβ = ∂ 2 /∂x α ∂x β . The free energy F = F S C + F G of the system is then represented by the spontaneous curvature model plus the energy density due to the Gaussian curvature, . Both terms determine the complete bending energy of the system. The time evolution of the phase-field is set according to the Cahn-Hilliard dynamics, since the volume is supposed to be locally conserved [32][33][34][35][36][37] . The variations of the free energy with respect to φ must then be subjected to diffusion, yielding a dynamic equation for the phase-field governed by 32 , We can express this last equation in terms of the order parameter and its spatial variations only (see "Appendix A"). The final result, after some algebra, is, The terms F K i represent the Gaussian curvature effect and are, explicitly, and The parameter σ [φ] is a Lagrange multiplier that depends on the field φ and assures area conservation 47 . One can determine σ [φ] by calculating the area S ∝ � |∇φ| 2 dV and demanding that dS/dt ≈ 0 . Using that Eq. (8) guarantees the conservation of the field φ , one obtains

Results
Numerical simulations. In principle, by solving Eqs. (9)(10), it is possible to model situations in which the topological genus of the membrane changes, as in vesicle formation. Thus, we performed three dimensional calculations using the same method of integration as in 32 . We used a finite-difference scheme for the spatial discretization and an Euler method for the temporal derivatives with the appropriate small time step of dt = 10 −5 to ensure enough accuracy and avoid artefacts 48 . For the Gaussian curvature term in Eq. (9) we implement the integration using the residue theorem over the positive complex plane. As the initial membrane, we choose a cylindrical shape of radius R and length L closed at the top and open at the bottom. Its base is in contact with the wall of the domain and we are using zero flux boundary conditions.
We start with the simplest phenomenon, the formation of two closed membranes from one, the Gaussian curvature controlling the fission of a single vesicle. The results are shown in Fig. 1 for the cases when (a) the Gaussian curvature term is not included in the free energy, and (b) solving the complete dynamical equation (9).
For the case when the Gaussian curvature is not considered ( κ ′ = 0 in Fig. 1a) a single neck forms, but there is no fission of the membrane. These results are consistent with previous works in which the Gaussian curvature is ignored 33 .
The addition of the Gaussian term in Eq. (9) makes the fission of the membrane energetically favorable and a vesicle is formed from the initial cylinder (Fig. 1b). We take −κ ′ ≈κ , as estimated in experiments and simulations alike 13,14 and in the figure we indicated the local free energy in units of κ . The size of this system only allows the cylinder to make a single vesicle. The dimensions of the vesicles can be predicted by a dispersion relation calculation shown below.
For longer cylinders multiple fission events are possible. In Fig. 2, we show the results for a cylinder of length 45. We obtained a sequence of single fission events from the tip of the cylinder downwards. After ∼ 270, 000 time iterations, the cylinder splits into five separate closed vesicles. Shape changes start simultaneously through the entire tube, but the speed of pearling and fission is not equal. The tip of the cylinder changes shape faster and splits first. In Fig. 2a we show the initial and final configurations of φ = 0 for a longitudinal medial section of the tube. Although Helfrich-Landau bending energy was originally derived for nearly-flat membranes it has been successfully used to explain membranes with large curvatures 49 .
The time evolution of the volume, the area and the energy contributions is shown in Fig. 2b. The bending energy (open circles) and the Gaussian energy (open triangles) terms are drawn in absolute units of κ . The changes in the Gaussian energy is half that of bending energy, which is in agreement with previous works 9 . The oscillating behavior is associated with the formation of vesicles in which the energy peaks are related to the necks narrowing and the following energy reduction is correlated with the formation of vesicles.  Fig. 3. A pearling instability and an ordered fission is still observed. The onset of the pearling instability in Fig. 3 also appears on previous works 33 , but here the pearls fission due to the Gaussian curvature contribution. We conclude that Gaussian curvature controls the formation of many small vesicles from one big elongated vesicle as seen in the experimental results 15,16 .
All simulations were carried out having a reservoir of volume and area at the base of the cylinder. Therefore, global conservation of area and volume is not necessary. However, conservation of area and volume for individual vesicles should hold. As the vesicles fission from the main membrane, they lose contact with any reservoir of volume and area. Thus, from the moment they split, the vesicles must maintain their area and volume. In Fig. 2b it can be noticed that the area and volume vary slightly, although these values for the individual vesicles remain stable.
In order to analyse the stability of neck formation, it is possible to study the effects of small perturbations around the flat interface. These perturbations are taken as plane waves of the form φ = φ 0 e iq·x−ωt , near φ = 0 , with small amplitudes, φ 0 ≪ 1 . Substituting these expressions into Eq. (9), and considering isotropic perturbations, one obtains the dispersion relation Here, δ represents terms of order O (φ 2 ) . The detailed derivation of the dispersion relation is given in "Appendix B". The main contribution to the instability comes from the first term in Eq. (11), in which the spontaneous curvature predominates. As the Gaussian curvature (due to κ ′ ) is not present in Eq. (11), it does not alter the region of unstable wavelengths. The Gaussian curvature acts by causing topological changes of the membrane right in the sites where the instability occurs. The periodicity of the instability is determined by a critical length, l c = π/q c , where q c > 0 given by the point in which the unstable branch ceases to be positive (ω(q) = 0) . This quantity represents a scale for the formation of the necks. The dispersion relation for typical values of the parameters is depicted in Fig. 4   www.nature.com/scientificreports/ Finally, we explore the regime in which the Gaussian bending modulus κ ′ determines whether vesicle formation happens or not. The results in Fig. 5 show the equilibrium configurations for a single continuous membrane with tubular shape. For negative κ ′ the transition from a single vesicle to to many small vesicles, occurs along the line where 8ε 4 κ + κ ′ = 0 (dotted line in Fig. 5). When κ ′ is positive there is the formation of a multiple selfconnected membrane, a continuous single membrane with multiple holes.

Discussion
In previous works 29 the Gaussian term has been worked out only in the case when the tube is stable, and pearled tubes are obtained, without pinching. Other works 51 focus on the microstructure of the membrane, which is basically a liquid crystal film, and suggest that topological defects could produce fission. By contrast, our model   www.nature.com/scientificreports/ results vesicle pearling and fission occur in sequence, starting at the tip, it is possible that new experiments with the appropriate proteins exhibits this behavior.
Here we have studied the fission of membranes (vesicle formation) by extending the bending free energy model through the introduction of the Gaussian curvature. The Gaussian energy term provides the pearling instability the ability to fission from the initial membrane and maintain the stability of the system. With this model, multiple fission events can be obtained in a single computation.
The dynamics of vesicle formation was studied numerically and an ordered fission of the tube, from the tip to the base of the tube, was obtained. The number of formed vesicles depends on the dimensions of the tubular domain and the value of the spontaneous curvature. This number could be predicted using the dispersion relation obtained from a linear analysis of the model.
Topological changes were explored taking into account the two bending modules, giving us a membrane shape diagram. We corroborate the existence of multiple vesicles for negative κ ′ for the values κ ′ < −8ε 4 κ . For Gaussian modulus κ ′ > −8ε 4 κ no topological transition occurs. For positive κ ′ the result is a multiple selfconnected membrane.
To summarize, numerical calculations based on this model describe fission of membrane tubes into multiples vesicles. Pearling and vesicle fission occur in sequence from the tip of the tube to the base. Using the bending modulus and the Gaussian modulus we can obtain a membrane shape diagram, and we explored the different regimes. The model can represent systems where appropriate proteins are required for fission by means of the Gaussian energy.

Appendix
Appendix A: Calculation of the complete dynamic equation for φ. Variations of F S C in Eq. (8) are explicitly, where is the chemical potential and σ [φ] is the surface tension coefficient, which is implemented as a Lagrange multiplier that ensures local area conservation.
It is necessary to express Eq. (1) in terms of φ so we might be able to establish the variations of the energy due to the Gaussian curvature similar to Eq. (12). These variations must also be subject to diffusion so we can obtain the dynamic equation that dictates the evolution of the surface. Variations of the free energy due to the Gaussian curvature can be written as We write this last equation in terms of the curvature tensor Q αβ , that is, in terms of the gradients of φ , and then calculate each of the terms of Eq. (14) separately.
First, variations of K with respect to φ can be written as: The second term of Eq. (14) is and the third term is: By adding up Eqs. (15), (16) and (17), and grouping similar terms we have: It is possible to express the variations of the free energy due to the Gaussian curvature in terms of K itself: . www.nature.com/scientificreports/ It is notable that the second term has a singularity when φ = ±1 (that is, in the bulk) and it also approaches to zero near the interface φ = 0 . Thus, the main contribution to the curvature energy of the surface comes from the first term of Eq. (19), which depends on the variation of the Gaussian curvature with respect to the order parameter φ.
It is more appropriate to write the dynamical equations in terms of φ . Thus, Eq. (19) is equivalent to The dynamic equation that dictates the evolution of the surface is which is defined only in terms of φ and its gradients and κ ′ = 3 √ 2εκ ′ /4.

Appendix B: Dispersion relation.
In order to analyse the stability of the membrane, one can study the effect of small perturbations of the flat interface. In particular, one considers that perturbations take the form of plane waves: φ = φ 0 e iq·x−ωt , where q = q α and x = x α are the wave vector and Cartesian coordinates ( α = 1, 2, 3 ), around the membrane φ = 0 . We assume that the amplitude φ 0 ≪ 1 is small. Taking into account Eq. (9) of the main text and substituting the approximation in each term, one obtains since The next term on the right hand side of Eq. (9) is, On the other hand, one obtains, The last equality is to be expected since in the plane wave approximation the functional associated with the Gaussian curvature is Now, if one substitutes Eqs. B1-B4 in Eq. (9), one obtains Assuming an isotropic perturbation q α = q with small wave number, and using φ := δ for the terms of the order O (φ 2 ) , one finally obtains the dispersion relation From this result, one observes that the spontaneous curvature is responsible for the onset of the instability and that the mean and Gaussian curvature do not affect the region of the unstable wavelengths.