key: cord-0541361-oz903868 authors: Bottacchiari, Matteo; Gallo, Mirko; Bussoletti, Marco; Casciola, Carlo Massimo title: Topological transitions in fluid lipid vesicles: activation energy and force fields date: 2022-02-28 journal: nan DOI: nan sha: 410878ca02983e304311147075420004ec1df5ec doc_id: 541361 cord_uid: oz903868 Lipid bilayer membranes are the fundamental biological barriers that permit life. The bilayer dynamics largely participates in orchestrating cellular workings and is characterized by substantial stability together with extreme plasticity that allows controlled morphological/topological changes. Modeling and understanding the topological change of vesicle-like membrane at the scale of a full cell has proved an elusive aim. We propose and discuss here a continuum model able to encompass the fusion/fission transition of a bilayer membrane at the scale of a Large Unilamellar Vesicle and evaluate the minimal free energy path across the transition, inspired by the idea that fusion/fission-inducing proteins have evolved in Nature towards minimal work expenditure. The results predict the correct height for the energetic barrier and provide the force field that, by acting on the membrane, can induce the transition. They are found in excellent agreement, in terms of intensity, scale, and spatial localization with experimental data on typical protein systems at play during the transition. The model may pave the way for the development of more complete models of the process where the dynamics is coupled to the fluid internal and external environment and to other applications where the topological changes do either occur or should be prevented. Biological membranes are formed by a fluid lipid bilayer which can be described not only from a biochemical point of view but also from a mechanical one. The most important example of fluid lipid membrane is the plasma membrane, which defines the boundary between the inside and outside of cells. Topological transitions of fluid lipid membranes are involved in most of the fundamental processes of cell life, like endocytosis and exocytosis. An example of topological transformation is the merging of two membranes. This is the case of vesiclevesicle fusion or viral membrane fusion. Indeed, viruses enveloped by a lipid bilayer, such as HIV, Ebola virus, influenza, measles, rabies virus and SARS-CoV-2 can infect a target cell by fusion of their membrane with the cell plasma membrane. The aim of this fusion is to inject the viral genome into the cell [1, 2] . Viral infection of the * carlomassimo.casciola@uniroma1.it cell can also occur via endocytosis, in which the plasma membrane undergoes fission to internalize the virus via an endosome. Therefore, another important topological change is membrane fission, which is also fundamental for cell division and therefore for life [3, 4] . Topological transitions of lipid membranes are of great interest not only in biology and biophysics but also in medicine and in the pharmaceutical industry. Indeed, lipid-based nanoparticles are used for drug delivery, offering many advantages including biocompatibility, bioavailability, self-assembly and payload flexibility. For example, liposomes, namely lipid bilayer vesicles with a diameter in the range of 25 nm ÷ 1 µm, are able to carry both hydrophilic and hydrophobic drugs [5] . Micelles, closed lipid monolayers, are currently used in mRNA-vaccines against COVID-19 and many other lipid nanoparticle-mRNA applications are under clinical evaluation, e.g. for the treatment of cancer or genetic diseases [6] . Regardless of the specific application, all these nanoparticles are engineered to overcome many physiological barriers, some of which involve topological transitions of lipid membranes. For instance, a topological transformation is needed for the cellular uptake of drugs and others may be required in the cytosol, e.g., to enter the nucleus for DNA delivery or genome editing. Therefore, a better comprehension of the mechanisms underlying such barriers could be useful to enhance the delivery efficiency, especially moving toward precision medicine [7] . Fluid lipid membranes can be mechanically described using a continuum approach initially introduced in [8, 9] . Such a classical elastic perspective describes a membrane as a two-dimensional surface Γ with an energy density depending on its principal curvatures. An expansion of this density up to the second order in curvatures leads to the Canham-Helfrich Hamiltonian: Here, the first term on the right-hand side is the bending energy and the second is the Gaussian energy. M is the mean curvature of the surface, G its Gaussian curvature, m a spontaneous mean curvature that the membrane tends to adopt in absence of external forces, and k and k G are called bending rigidity and Gaussian curvature modulus, respectively. k can be experimentally measured in different ways [10] , whereas k G is more elusive due to the celebrated Gauss-Bonnet (GB) theorem, where χ(Γ) is the Euler characteristic of Γ and k g is the geodesic curvature of the surface boundary ∂Γ. In the vesicles case, since they are compact surfaces without boundary, the line integral vanishes, and χ(Γ) becomes equal to 2 − 2g, being g the genus of the surface. Therefore, the Gaussian energy term remains constant as long as no topological transitions occur, leading to the aforementioned elusive behavior of k G . A stability argument [11] shows that −2 < k G /k < 0 and in literature there is evidence [12] [13] [14] that k G ≈ −k. Hence, the Gaussian energy is expected to play a crucial role during topological transitions. Because of the scale invariance of the Canham-Helfrich free energy, for a given topology, vesicles shapes are dictated by their reduced volume v = V /(π D 3 ve /6), as well as by their reduced spontaneous curvature m 0 = mD ve , where D ve = A/π is the characteristic lenght of the vesicle under consideration, having area A and volume V . The Canham-Helfrich description is thought to hold for vesicles with a characteristic length D ve ≥ 40 l me [15] , being l me the lipid bilayer thickness, which is usually about 5 nm; otherwise, higher-order terms in the energy density could make a significant contribution. Therefore, Eq. (1) describes vesicles larger than the ones simulated by means of coarse-grained molecular dynamics (MD) and dissipative particle dynamics (DPD), which have been the most widely used techniques for in silico studies of topological transitions to date. Indeed, these computer simulations, which take into account the molecular details of lipid bilayers, allow monitoring in time morphological changes of small liposomes with a size below 50 nm [16] . For example, membrane adsorption of small solutes inducing budding and fission of nanoparticles has been simulated in [17] , membrane fusion and drug delivery with carbon nanotube porins in [18] , whereas the intermediate structures in which lipids are organized during fusion reactions have been investigated in a series of works including [19] [20] [21] [22] [23] . Theoretical description of these intermediates has been addressed, among others, in [24] [25] [26] [27] [28] and the same has been done for membrane fission in [29] . Membrane fusion and fission events have also been experimentally studied. For example, recently, controlled fission of cell-sized vesicles by low densities of membranebound proteins has been reported in [30] . Other examples of fission experiments can be found in [31, 32] , whereas, as regards fusion, merging of giant liposomes has been observed in [33] [34] [35] , the stalk intermediate in [36] , and activation energies for small liposomes fusion events have been measured in [37, 38] by means of kinetic analysis. In this work, we numerically study fission and fusion events of Helfrich-sized vesicles. In order to achieve our purpose, we use the phase-field framework [39, 40] , which allows us not to track the interface and naturally handle topological transitions. A phase-field version of the bending energy term alone has been initially introduced in [41] [42] [43] , leading to numerous applications like [44] [45] [46] [47] [48] [49] [50] [51] . Furthermore, in [52] , it has been pointed out that it is possible to retrieve topological information for phase-field models. Here, given its importance in topological transitions, we include the Gaussian energy term in the dynamics. This coupling has a rigorous mathematical derivation, which is reported in Section II. The new term numerically captures with high accuracy the energy jumps as expected by the GB theorem; see for comparison the recent paper [53] , where a different form of phase-field Gaussian energy was proposed to study pearling instabilities. Interestingly, Appendix A, the GB theorem can be derived from our new phase-field energy term. In this paper, by means of the string method [54] , a rare event technique, we compute a minimal energy pathway (MEP) between two spherical vesicles and a dumbbell-shaped one (see, e.g., [55] [56] [57] [58] for the string method applied to a phase-field approach). This configuration is of special interest for applications since most of the aforementioned experimental results refer to fusion/fission processes involving spherical vesicles. However, the approach we propose can be easily applied also to other vesicle configurations like, e.g., two spheres of different sizes, two non-spherical vesicles or to the limiting case of a sphere and a plane. The string method, coupled with the new phase-field model, provides the free energy barriers of the fusion and fission processes and the membranes configurations along the path. Moreover, since the phase-field extension of the Canham-Helfrich energy regularizes the behavior of the Gaussian term, we also compute the force fields needed to overcome these barriers. These forces are necessary to balance the elastic reaction arising from the bending and Gaussian energies, as well as from the membrane incompressibility. In the classical Canham-Helfrich approach, the computation of the bending forces would not have been so easy [59] and it would not have been clear how to include the Gaussian contribution during topological transitions. The computation of such forces will pave the way for exploring the effective mechanisms by which fusion and fission machinery work across the full scale of vesicles. The paper is structured as follows: Section II illustrates the main aspects of the underlying mathematics; this section is not essential for the discussion of the physical results. Section III details the obtained results, which are discussed in Section IV; final remarks are reported in Section V, whereas Section VI treats the more technical aspects (numerical scheme, string method, and the computation of the force fields). The classical Canham-Helfrich model relies on a sharp interface description of the membrane, which is treated as a (zero thickness) surface endowed with an energy density depending on the principal curvatures, Eq. (1). The model succeeds in describing many aspects of the vesicle dynamics but rules out the possibility of dealing with topological changes, unless unphysical surgical operations are conceived to cut and paste patches of the membrane. A viable alternative to the sharp interface description is to employ a smooth function defined on a domain Ω -the phase-field φ(x) -that discriminates between the inner and the outer environment of the vesicle assuming the limiting values ±1 in the two regions. The φ(x) = 0 level set represents the membrane midsurface Γ. The transition between the two limiting values takes place in a narrow region whose width is controlled by a small parameter and can be associated with the thickness of the lipid bilayer. The main advantage of describing the membrane with such a field lies in the fact that it enables topological modifications of the membrane, allowing to address the problem of vesicle fusion and fission. An energy E[φ] is associated with each field configuration and is such as to admit local minimizers of the form where d(·) is the signed distance function from the membrane midsurface Γ. We choose to define the signed distance such that n = ∇d computed on Γ is equal to the inward-pointing unit normal to the vesicle. Setting d * (x) = d(x)/ , we also require that lim d * →±∞ φ = ±1 and φ = 0 for d = 0. Therefore, ±1 are the values for the stable phases of the inside and outside bulk and the level set φ = 0 identifies the membrane midsurface. Physically, the energy functional should recover the Canham-Helfrich energy, (1), in the limit of small width-to-vesicleextension ratio. The proposed expression for the phase-field free energy functional is with and Here m is the spontaneous curvature of the membrane, taken to be positive if the membrane bulges towards the exterior, k is the bending rigidity and k G is the Gaussian curvature modulus. Henceforth, we will assume was already introduced in [41] to model the bending energy of the membrane while E G [φ] is the new term proposed here to account for the Gaussian energy. As anticipated, our purpose here is to show that, under the general ansatz (3) and in the sharp-interface limit ( /D ve = λ << 1), minimizing the phase-field free energy functional (4) is equivalent to minimizing the Canham-Helfrich free energy. Denoting with a prime the derivative done with respect to d * (x), a direct computation leads to where we have denoted with a bar the dimensionless lengths obtained by dividing by D ve . Therefore, in order to minimize E = E B + E G , as λ → 0, the leading-order term which has the solution Hence, is actually related to the width of the interface. Moreover, by repeating the computations done in [60] for the bending energy alone, it is possible to show that, also in the presence of the new Gaussian energy term, one finds f 1 (d * (x)) ≡ 0 (see the Supplemental Material [61] for the whole computation). Therefore, given that we are left with Denoting with k 1 and k 2 the principal curvatures, we have ∇ · n = −(k 1 + k 2 ) = −2M and n · ∇k i = k 2 i , with the result that (∇ · n) 2 + n · ∇(∇ · n) = 2k 1 k 2 = 2G. Now, noticing that for λ → 0 note [62] , and getting back to dimensional variables, the asymptotic behavior follows as i.e., the phase-field energy functional reproduces the Canham-Helfrich free energy in the sharp-interface limit ( /D ve << 1). It is worth noticing that the inclusion of the Gaussian energy, which is subdominant in λ, preserves the hyperbolic tangent form (12) of the leading order solution together with f 1 (d * (x)) ≡ 0, as for the more standard model with the bending energy alone [41] . Since f 1 (d * (x)) ≡ 0, the desired expression of the bending energy is retained at order λ −1 , and the accuracies O(λ) and O(λ 2 ) are guaranteed in Eqs. (14) , (15) , respectively. Furthermore, in our formulation, the phase-field Gaussian energy (7) has no singularities and actually depends at most on derivatives of order two, as it is possible to see by replacing ∇φ·∇∇ 2 φ with ∇ 2 |∇φ| 2 /2−H φ : H φ in (8) , where H φ is the Hessian matrix of the field. In many cases, it is important to fix the vesicle area and volume. Indeed, since lipids are insoluble in water, the number of membrane lipids is conserved. This fact, coupled with the observation that the membrane rupture tension is very small, implies that the vesicle area A cannot substantially vary at a fixed temperature. The volume V of the vesicle is instead determined by the osmotic conditions. In order to enforce the above constraints, in the phase-field context, one can use the functionals which respectively behave like the vesicle area and volume in the sharp-interface limit. Throughout the paper, an asterisk will denote the dimensionless quantities obtained using as the reference length and 8πk as the reference energy. The latter is the bending energy of an isolated sphere. A typical value of the bending rigidity is k = 20 k B T , with k B the Boltzmann constant and T the temperature. Figure 1 shows the Gaussian energy during a series of scissions of an unstable prolate shape into several spheres. The process is driven by the spontaneous curvature. The evolution equation is described in Section VI A together with the adopted numerical scheme. Details concerning the numerical validation are provided in Appendix B. Consistency of the present phase-field approach with the Gauss-Bonnet theorem is discussed in Appendix A. Here, it is only worth saying that the novel model is able to properly capture the Gaussian energy jumps due to topological transitions. FIG. 1. The phase-field Gaussian energy during a series of scissions of a prolate shape into several spheres. The energy jumps by −4πk for any division as prescribed by the Gauss-Bonnet theorem (k = −kG). The fission process occurs due to the presence of a spontaneous curvature m * ≈ 0.42. Time evolution is given by the Allen-Cahn gradient flow with M * = 8 (see Section VI A for more details on the dynamics, the adopted numerical scheme and dimensionless quantities). The inset shows the total energy E = EB + EG, which monotonically decreases in time, revealing the stability of the scheme. This z-axial symmetric simulation has been carried out in a [0, 36] × [0, 440] computational domain in the r * − z * plane with a 54 × 660 mesh, initial D * ve = 1/λ ≈ 109 and dt * = 4. There is no constraint on the area, which, at the end of the simulation, differs from the initial value by approximately 6.9%. Volume is conserved with a relative error smaller than 10 −7 with respect to its initial value. In the topological transition between two spherical vesicles and a dumbbell-shaped one, which are two stable states, the system goes through a sequence of configurations φ α (x) in the space of the phase-field, identifying a path which we parameterize by the normalized arclength α ∈ [0, 1]. An MEP for this transition is a curve on the energy landscape E[φ] connecting the two stable states φ α=0 (x) and φ α=1 (x), respectively, and such that it is everywhere tangent to the gradient of the potential (∂φ α /∂α ∝ δE[φ α ]/δφ), except at critical points [63] . An initial guess of the path is discretized in a string made up of N = 100 images corresponding to α i = (i − 1)/(N − 1). The initial guess is relaxed towards the MEP by means of the string method (see [54, 64, 65] and Section VI B) suitably accounting for the constraints of constant total surface area, Eq. (19) , and enclosed volume, Eq. (20) . The obtained MEP goes through a saddle point φ αc (x) for the free energy, determining the transition barriers for the forward and backward process, respectively. Figure 2 shows the computed MEP for membranes with zero spontaneous curvature, m = 0. Since the phase-field φ reaches its limiting values ±1 with an accuracy of about 3% already at a distance of ±3 from the φ = 0 membrane midsurface, we assume that l pf = 6 represents the thickness of the diffuse interface. We have shown in Section II that the phase-field description recovers the Canham-Helfrich model in the limit of small λ ∝ pf /D ve . Our numerical experiments, Appendix B, point out that this asymptotic behavior is already achieved when pf /D ve >> ( me /D ve ) max = 1/40, the latter being the maximum thickness-to-curvature radius ratio for which the Canham-Helfrich model is accepted [15]. Since the relative distance between approaching membrane segments is relevant during the topological transition, it is crucial that the diffuse interface width matches the bilayer thickness. This requirement fixes the scale of our system. Setting l pf = l me = 5 nm, the configurations shown in Figure 2a correspond to vesicles with D ve ≈ 206 nm, thus within the range of validity of the asymptotic Canham-Helfrich model. Increasing/decreasing α corresponds to moving along the path in the direction of the fusion/fission (forward/backward) process, respectively. Proceeding forward, the two vesicles come closer to each other without deforming, get in touch, and merge together forming a narrow neck that expands until the final dumbbellshaped configuration is reached. As explained in Section I, the equilibrium states of a vesicle are determined by its reduced volume and reduced spontaneous curvature, which, in the present case, are v = 1/ √ 2 and m 0 = 0, respectively, where 1/ √ 2 is the only reduced volume compatible with a vesicle obtained from the fusion of two spheres of the same radius. As shown in [66] , with these parameters, it is possible to reach two axisymmetric configurations with the topology of a sphere, namely one oblate-discocyte shape and one prolate-dumbbell shape. The latter has the lowest energy and, in the present case, is the equilibrium state assigned to the string as the final configuration, φ α=1 (x). Figure 2b shows the free energy profile along the MEP. The free energy of the final configuration (prolate) is E[φ α=1 ]/(8πk) ≈ 1.12, which is larger than the initial energy E[φ α=0 ]/(8πk) = 1 of the two spheres. Both val-ues are in excellent agreement with the data reported in [66] . One may notice that the two-spheres configuration possesses a sequence of neutral equilibrium states, corresponding to rigid translations during which the two vesicles approach/separate from each other (configurations i from 1 to 11, as also depicted in Fig. 2a) . The saddle point consists of two spheres connected by a small narrow neck and is located between configurations i = 14 and i = 15, with the latter having the highest energy of the two, E[φ α=αc ]/(8πk) ≈ 1.45. It should be noticed that such a configuration possesses the bending energy of two spheres together with the Gaussian energy and the topology of a single sphere. Hence, the forward and backward computed free energy barriers are ∆E † 0→1 /(8πk) ≈ 0.45 and ∆E † 1→0 /(8πk) ≈ 0.33, respectively. Considering a bending rigidity k of order 20 k B T [10] , it turns out that both fusion and fission processes require further agents in order to happen, in addition to elasticity and thermal fluctuations. These agents are typically protein systems. Still in Fig. 2b , it is possible to observe a substantial asymmetry between fusion and fission, with a much steeper energy increase required to reach the transition state in the fusion process. The main plots in Fig. 2c provide the bending and Gaussian contributions to the free energy along the MEP. Apparently, the forward barrier ∆E † 0→1 is almost entirely due to the Gaussian energy jump associated with the topological change. On the other hand, the backward barrier ∆E † 1→0 builds up continuously with the progressive deformation of the prolate shape to form the narrow neck preceding the actual fission. The inset shows the evolution of the area and enclosed volume along the MEP, confirming that the constraints are perfectly satisfied at each string image. Figure 3 focuses on the region of the MEP where the most relevant events associated with the topological transition take place, images i = 11, . . . , 40. The contour plots in the lower half panels of Fig. 3 provide the structure of the phase-field as a function of radius r * and axial coordinate z * , with φ smoothly joining the inner region φ = 1 to the outer region φ = −1 through the layer of dimensionless thickness * pf = 6. As explained in Section VI C, each image of the string can be rendered a state of equilibrium by introducing a force field f = −δE/δφ∇φ that counterbalances the membrane elastic reaction. Considering the forward transition, 0 → 1, such force field from α = 0 to α = α c can be interpreted as the external force needed to drive the transition under quasi-static conditions, thus spending the minimal work W 0→1 = ∆E † 0→1 . Once the critical state is overcome, the system can be left to evolve spontaneously until it reaches the final equilibrium state α = 1. Symmetric considerations hold for the backward transition 1 → 0. The dimensionless vector fields f * α (x) are depicted as arrows in each panel of Fig. 3 , where, for the sake of better readability, they are plotted only on the φ = 0 isoline. The contour plots on the upper part of each panel display the component of the force normal to φ-isolines. It should be noticed that the scale of the arrows changes from panel to panel, at least for the upper frames, i = 11, . . . , 14. For the forward process, the latter are the configurations achieved just before the critical state. In this region, the MEP is particularly steep, requiring more intense forces, which result to be strongly localized near the vesicles contact region. On the contrary, the backward process requires a more distributed force field, as shown in images i = 15, . . . , 40. The arrows reverse their direction between configurations i = 14 and i = 15, showing that in this interval the force field vanishes, confirming that the critical state occurs somewhere between these two images. As shown in Fig. 2 , the phase-field model is able to account for the Gaussian energy and its effects on the fusion and fission processes, allowing to bridge the gap across the topological transition. The initial and final (meta)stable states of the process under consideration, i.e. the two spheres and the dumbbell-shaped vesicle, are perfectly consistent with the predictions of the sharp interface description à la Canham-Helfrich [66] , both in terms of energy and shape. The two states have the same reduced volume v = 6 √ πV /A 3/2 = 1/ √ 2, which is constant throughout the process, since both the volume and the area are accurately conserved, inset of Fig. 2c . The force field f = −δE/δφ∇φ is able to induce both the fusion and fission processes with a minimal work expenditure. This external driving force, depicted in Fig. 3 , is clearly required only during the climbing phases of the energy landscape, Fig. 2b , i.e. from i = 1 (two spheres configuration) to i = 15 (saddle point) for fusion and from i = 100 (dumbbell-shaped vesicle) back to i = 15 for fission. Proceeding in the fusion direction, the energy landscape is initially flat and concerns the apposition of the two vesicles. In the biological context, the final part of this stage corresponds to the dehydration of the two facing membranes. This effect is clearly absent in our model, although it could be included by introducing a suitable effective repulsive potential. After this flat region, the energy undergoes a sharp jump, which is related to the variation of the Gaussian energy across the topological transition, |∆E G |/(8πk) = 0.5. This jump is associated with the change of the Euler characteristic χ, with χ = 4 for the two spheres and χ = 2 for the prolate-dumbbell shape, with the latter homeomorphic to a single sphere. Therefore, the force field required to complete the fusion process is very intense and localized, Fig. 3 , and essentially due to the Gaussian energy contribution. This leads to an asymmetry between the fusion and fission processes. Indeed, in the former case, the activation energy is entirely due to the topological transition, while in the latter case the barrier progressively builds up through the continuous membrane deformation. The forces required for overcoming the fusion topo-logical barrier are too strong to be directly exerted by the sole mechanical action of proteins. For example, setting k = 20 k B T [10] , the resulting activation energy is Consistently with the present findings, Deserno [67] suggests that fusion proteins, besides a mechanical action, may contribute to lowering the energy barrier by locally modifying the Gaussian modulus in the contact region of the approaching membranes. Indeed, the introduction of a suitable, spatially dependent Gaussian modulus is expected to reduce the stiffness associated with the GB theorem, opening alternative routes to the topological change. Our results show that this scenario is actually possible since the forces associated with the Gaussian energy are localized in the region of contact between the two spheres and, therefore, it is reasonable that a variation of k G in such a region could lower the activation energy. For example, this situation is compatible with the observation that influenza virus hemagglutinin proteins, in addition to having an apposition activity, are also able to perturb the membrane lipid bilayer by insertion of their amphipathic fusion peptide [68] . Interestingly, the present phase-field approach can be easily adapted to the instance of a topological transition with a spatially dependent Gaussian modulus, a case we leave for future work. As illustrated in Fig. 2a , the topological transition is mediated by the formation of a catenoid-like neck [69] , similarly to what has been observed in the experiments [32] . Operationally, we define the neck region as the zchunk of the fused vesicle where the local contribution to the Gaussian energy is positive, k G G j(z) ≈ k G 2πr ψ G dr > 0. The approximate equality follows by considering that the phase-field functional E neck G approaches the corresponding sharpinterface Canham-Helfrich Gaussian energy, where the neck midsurface is described by the Monge representation r = r(z) and dA = j(z)dz is the corresponding area element. The computed phase-field Gaussian energy of the neck along the MEP is shown in the top panel of Fig. 4 , blue line. Proceeding from left to right, E neck G (Z)/(8πk) sharply increases to a value close to (though smaller than) 0.5 and subsequently decreases. It is worth stressing that, from a sharp interface point of view, the total curvature, namely the integral of the Gaussian curvature, is 4π for a single sharp sphere, which corresponds to E CH G /(8πk) = −0.5 in terms of Canham-Helfrich Gaussian energy. Given two initially disjoint sharp spheres (total curvature 8π, E CH G /(8πk) = −1), the formation of a joining neck changes the topology and reduces the total curvature to that of a single sphere, 4π, E CH G /(8πk) = −0.5. There are two main reasons why the phase-field model provides a E neck G (Z)/(8πk) contri-bution to the Gaussian energy that is slightly smaller than 0.5: i) close to the transition state, the curvature of the neck generatrix in the r − z plane is comparable with the finite thickness of the bilayer, so that the sharpinterface model is inappropriate; ii) the total curvature of the neck midsurface is always larger than −4π and can reach the latter limiting value only when the tangent to the generatrix gets orthogonal to the z-axis at the two edges, see note [70] . Evidently, the neck contribution to the Gaussian energy is also the main contribution to the barrier ∆E † 0→1 /(8πk) ≈ 0.45. Proceeding to the right along the MEP, beyond the saddle point, the computed Gaussian energy of the neck progressively decreases, top panel of Fig. 4 , blue line. Since, Fig. 2c , in that region the total Gaussian energy remains overall constant, E G /(8πk) = −0.5, the energy lost by the neck is redistributed to the remaining, dome-like parts of the vesicle. Figure 4 , top panel, orange line with dots, also provides the neck Gaussian energy as a post-processing based on the sharp interface Canham-Helfrich energy, Eq. 1, computed considering the φ = 0 level set as the membrane midsurface. In such a manner, the neck Gaussian energy can be evaluated as where r(±Z) is the radius of the circular neck boundary and R n is the radius of normal curvature of the vesicle cross section with the plane containing both the surface normal and the tangent to the circles, evaluated at the neck boundary [71] . In order to get a geometrical understanding of these quantities, the bottom panel of Fig. 4 shows a few vesicle configurations, with their upper neck boundary (yellow circle) and the osculating (red) circle defined by the radius of normal curvature. The yellow and red arrays depict r(±Z) and R n , respectively. A more detailed discussion about Eq. (22) is provided in the Supplemental material [61] . The agreement between the phase-field Gaussian energy of the neck and its sharp interface counterpart is excellent as long as different membrane segments do not approach each other to a distance comparable to the bilayer thickness. Indeed, the agreement progressively deteriorates when getting closer to the saddle point, due to the increased curvature of the membrane generatrix. Actually, also in this stage, one could better and better reproduce the sharp interface energy by reducing the regularizing parameter λ. On the other hand, from a physical point of view, the thickness of the bilayer is finite, making the sharp interface model inappropriate when the saddle point is approached. Overall, these results confirm the accuracy of the proposed phasefield expression for the Gaussian curvature. As anticipated, the forces at play during fission are more distributed and less intense than for fusion. The large region they act on, Fig. 3 , is consistent with the cooperation of several protein systems, like, e.g., in clathrin mediated endocytosis, which involves clathrins polymerization and the subsequent action of the constrictase dynamin [32] . One can estimate the minimal work the protein system needs to perform to induce the topological change by comparing the free energy barrier ∆E † 1→0 with the protein work W 1→0 = f p ∆r, where f p is the order of magnitude of the protein force and ∆r = r max − r 0 is the change in vesicle radius at the neck, between the equilibrium prolate (r max ) and the saddle point configurations (r 0 ). Given the scale of the system, Section III, we find ∆r = 37.4 nm which, from the barrier height, provides f p = 0.91 k pN/k B T . Interestingly, for the values of k proper of fluid lipid membranes, we obtain protein forces in fairly good agreement with the experimental estimates reported in [30] , e.g. 20 pN for dynamin, 65 pN for ESCRT-III and 80 pN for FtsZ. For example, by assuming k = 20 k B T , we obtain a protein constriction force f p of 18.2 pN. For the same bending rigidity, Fig. 5 shows, red curve with squares, the energy needed to complete the fission process as a function of the current neck radius r n , ∆E(r n ) = E(r n ) − E(r 0 ) (note that the fission proceeds from larger to smaller neck radii, i.e. from right to left along the abscissa). The corresponding image number i along the MEP is provided on the second abscissa axis on the top of the frame. The slope of the plot, d∆E/dr n , orange line with triangles, provides the estimate of the constriction force (positive when constrictive). A plateau is apparent at d∆E/dr n 20 pN in the range of radii 16 ≤ r n ≤ 21 nm. Notably, it is known from the literature [72] that, e.g., dynamin polymerizes on tubules with radius between 10 and 30 nm, exherting forces of the order of 20 pN. In order to facilitate comparison with published data, Fig. 5 also provides in blue, with dots, f p = ∆E(r n )/(r n − r 0 ). [72] . The horizontal light orange strip depicts the value of dynamin constriction force, measured in the experiments [30, 72] . It is worth stressing that during topological transitions, when the relative distance between approaching membrane segments becomes comparable to the bilayer thickness, the scale invariance of the asymptotic Canham-Helfrich Hamiltonian (1) is broken. For such a reason, Section III, we defined the scale of our system by matching the bilayer thickness with the diffuse interface width. However, far from the topological change, scale invariance is preserved, implying that, during the phases where the Gaussian energy remains constant, the results of Figs. 2 and 3 hold for vesicles larger than those considered so far, too. Furthermore, for vesicles of smaller λ = /D ve , the saddle point energy E[φ c ] ought to approach the limiting value 1.5, namely the energy of two spheres connected by a zero radius neck. So far, we have discussed the fusion and fission processes as obtained by relaxing the system to the MEP using the Allen-Cahn dynamics, endowed with the additional constraints on surface area and enclosed volume. Possible alternatives can be considered, e.g., the Cahn-Hillard dynamics where the relaxation is driven by the divergence of a flux. Both choices lead to the same free energy barriers since the critical configurations (local minima and saddle point) are the same. However, the path connecting the saddle point to the two minima depends on the specific relaxation dynamics and the Allen-Cahn relaxation does indeed provide the minimal energy path. Let us conclude with a few general remarks. We have described the full-scale process of topology change in the fusion/fission process of two large unilamellar vesicles (LUVs) with an approach that can be extended to deal with giant unilamellar vesicles (GUVs). The model thoroughly accounts for the Gaussian energy, for which a suitable phase-field expression has been provided, and, far from the topology change, recovers the Canham-Helfrich description in the limit of small bilayer thickness. During the transition, our proposal should be interpreted as a rational way to regularize the singularity, leading to a process that smoothly matches the external solution before and after the transition. The model clearly misses the many molecular details associated with the dynamics of the lipids forming the bilayer. It could be argued, however, that the correction due to such details should be small as compared to the energy barrier associated with the full-scale evolution of the vesicle. Naively, one may argue that proteins systems could have evolved to overcome the large barrier that stabilizes the vesicle topology by following a minimal energetic pathway. Hence, we have evaluated the minimal free energy path for the transition and extracted the force field able to drive the process with minimal work expenditure. The free energy profile we find confirms the strong asymmetry between the fusion and the fission processes. For fusion, the force field required to overcome the estimated barrier is too intense to be exerted by protein systems and calls for additional effects that could locally modify the Gaussian modulus during the topological change. On the contrary, the spatial scales and the forces acting during fission are consistent with the experimental estimates for typical fission protein systems. Finally, it can be noted that the diffuse interface approach can naturally be coupled with hydrodynamics [73, 74] to include the dynamic effect of the external and internal aqueous environments of the vesicles. One may also stress that the Gaussian energy functional can find a much broader scope, e.g., as an indicator of the topological genus in the context of cluster analysis [75, 76] , or as a way to provide a barrier towards undesired/unphysical fusion processes. A compelling example concerns emul-sions where surfactant-covered droplets behave much like lipid micelles [77, 78] , suggesting that the Gaussian energy is the reason why surfactants do act as emulsifiers. The numerics relies on FFT-based spectral differentiation in cell-centered grids which provide high accuracy solutions, with special regard to the estimate of the Gaussian energy. The accuracy in evaluating the Gaussian energy, Eq. (7), is shown in Table I The energy pathways of Section III are obtained by means of the string method, which is briefly described in Section VI B. The remaining simulations reported in this paper, i.e. the one shown in Fig. 1 and those in Appendix B, are carried out using the Allen-Cahn dynamics where M is the mobility coefficient and δĒ/δφ is the functional derivative of the augmented energȳ Here, the additional terms added to the energy (4) are needed when constraining to A 0 and V 0 the vesicle area (19) and volume (20) , respectively. M 1 , M 2 are two penalty constants, whereas γ and ∆p are updated at each time step according to the augmented Lagrangian method, [79] : Therefore γ and ∆p are estimates of the Lagrange multipliers that improve at every time step. Starting from an assigned initial condition, the Allen-Cahn dynamics causes the energy to monotonically decrease in time until it reaches a critical steady-state. The dimensionless time and mobility are t * = t/τ R and M * = 8πkM τ R / 3 , respectively, with τ R a suitable time scale. With the help of the PETSc library [80] , a Crank-Nicolson time-stepping scheme is employed to integrate the Allen-Cahn gradient flow, while a semi-implicit Euler single step scheme is used to solve the more computationally demanding string dynamics. The explicit form of the functional derivative δĒ/δφ is given in Appendix C. The zero-temperature string method [54] is a technique for computing free energy barriers and transition pathways on a given energy landscape. The method proceeds by evolving in time a string, namely a curve parameterized by α ∈ [0, 1]. For each α the image of the string is a phase-field function φ α (x) representing a membrane state. Given an initial guess for the pathway connecting two local minima, the string evolves in time following the dynamics where M is a mobility coefficient, δĒ/δφ α is the functional derivative of (24) evaluated on the image φ α and (δĒ/δφ α ) ⊥ is its component normal to the string. This last quantity can be computed as (δĒ/δφ α ) ⊥ = δĒ/δφ α − δĒ/δφ α τ τ , where τ = ∂ α φ α / ∂ α φ α | ∂ α φ α 1/2 is the unit tangent to the string and · | · is the standard L 2 inner product. In this way, at steady state, the string converges to a minimal energy path [63] . In order to eliminate the trouble of projecting the functional derivative and in order to use the equal arc-length parameterization, the string dynamics can be rewritten [65] as whereλ = λ + M δĒ/δφ α τ and λ is a Lagrange multiplier for the purpose of enforcing the chosen parameterization ∂ α ∂ α φ α | ∂ α φ α 1/2 = 0. The algorithm follows the steps: 1. Evolution from t to t + ∆t of the discrete string, made up of N images φ i , with the dynamics Time integration is performed in wave number space by means of the semi-implicit Euler single step scheme. The evolved images at time t + ∆t are denoted asφ i . 2. Computation of the arc lengths corresponding to the evolved images: Thus, the evolved images have parameters α i = s i /s N . 3. Linear interpolation of the evolved images in order to compute the new images at equal arcs α i = i/N . These are the actual solutions at time t + ∆t. It is worth noticing that linear interpolation conserves vesicles volume. 4. Go back to one and iterate until convergence. Given a membrane state, it is possible to compute the external force needed to balance the elastic force arising from the energy of the membrane. For this purpose, let's consider an arbitrary and infinitesimal variation δφ of the phase-field, consistent with the area and volume constraints, if present. This variation results in a spatial displacement δx of the field lines. The displacement can be thought to occur in a virtual time interval δt, within which the field lines move with a virtual velocity u such that ∂φ/∂t = −∇φ · u (null material derivative condition). By integrating in time this last equation from t to t + δt, we are left with the first order approximation Hence, the work performed by the external force field f to deform the membrane is and one can identify the force field thanks to the arbitrariness of δx. ACKNOWLEDGMENTS Support is acknowledged from the 2020 Sapienza Large Project: Dynamics of Biological and Artificial Lipid Bilayer Membranes. Concerning computational resources we acknowledge: PRACE for awarding us access to Marconi successor at CINECA, Italy, PRACE 23rd call project Nr. 2021240074; DECI 17 SOLID project for resource Navigator based in Portugal at https://www.uc. pt/lca/ from the PRACE aisbl; CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support (ISCRA-B FH-DAS). Let's assume that where x ∈ Ω, being Ω a cylindrical domain of radius R and height L in the ordinary three-dimensional space, and d(·) the signed distance from an axisymmetric surface in Ω. This assumption leads to and, moreover, we set Using the cylindrical coordinates system, it is possible to show by a direct computation [52] that one of the two principal curvatures is Therefore, remembering that ∇ · n = −(k 1 + k 2 ) and n · ∇k i = k 2 i , with n = ∇d, Eq. (15) can be rewritten as where Supposing to have a single, connected, closed surface, after letting Ω invade R 3 , and still considering relation (17) with (A2), we obtain recovering the Gauss-Bonnet theorem (2) in the axiallysymmetric case. The last equality is justified by the fact that the Dirac delta function counts the intersections of the surface with the z-axis, which is equivalent to checking whether the surface has a hole. In literature [81] , it is well known that, in the presence of the bending energy alone, two initially close-by spheres merge together during the Allen-Cahn dynamics, Eq. (23) . If the area and volume constraints are included, at the steady state, a dumbbell shape with a reduced volume v = 1/ √ 2 is obtained. This happens because the bending energy of the two spheres is greater than that of the obtained dumbbell shape. Moreover, such a numerical experiment shows that there is no energy barrier for the process. This behavior is no longer possible if Gaussian energy is also included. Indeed, in this case, the whole Canham-Helfrich energy of two spheres is less than that of the dumbbell shape. Therefore, our first numerical validation experiment is to repeat this simulation including the new Gaussian energy term, Eq. (7). As shown in Fig. 6 , two spheres of equal radius R * = 10 at distance R * /2 from each other do not merge. This simulation has been carried out in a [0, 40]×[0, 40]×[0, 66] full 3D x * − y * − z * domain with a grid of 40 × 40 × 66 nodes (grid length interval h * = 1), * = h * = 1, 1/λ = 20 √ 2, M * = 1 and time step dt * = 0.8. In Fig. 6 , the energy monotonically decreases over time, revealing the stability of the scheme. At steady state, the final computed bending and Gaussian energies are E * B ≈ 1.941 and E * G ≈ −1.023. The same simulation has also been carried out in the r * − z * plane, exploiting the axisymmetry, with * = 2h * = 1 and two different time steps, dt * = 0.8 and dt * = 0.4, respectively, still obtaining the same behavior. Convergence has also been observed setting * = 1.5 h * = 1 and 1/λ = 40 √ 2, with the two spheres at distance R * /4 from each other. With the same parameters, let's take the dumbbell shape obtained merging the two spheres in a simulation with the sole bending energy, and let's use it as a new initial condition for the Allen-Cahn dynamics where the Gaussian energy is now included. As shown in Fig. 7 , the dumbbell shape remains substantially unchanged, showing that the configuration is still a local energy minimum and that there exists an energy barrier that prevents it from dividing into two spheres. It is worth noticing that the computed energies are in excellent agreement with the ones reported in [66] . Bending and Gaussian contributions to the energy of the final configuration are E * B ≈ 1.625 and E * G ≈ −5.095 · 10 −1 . Finally, we test a toroidal topology case. The initial condition is a torus with exact circular cross section of radius R * = 10 and v ≈ 0.6. The dynamics leads to a torus with a cross section that is no more perfectly circular, in excellent agreement with [82] , both as regards the shape and the energy. Figure 8 shows the energy evolution both with and without the Gaussian energy term. The two dynamics appear to be very similar, confirming that the Gaussian energy term plays no role as long as no topological transitions occur. These axisymmetric simulations have been carried out in a [0, 40] × [0, 40] computational domain in the r * − z * plane with a grid of 60 × 60 nodes, * = 1.5 h * = 1, 1/λ = 20 √ 2π, M * = 1 and dt * = 1. With the Gaussian term included, the final computed bending and Gaussian energies are E * B ≈ 1.831 and E * G ≈ −4.813 · 10 −2 , respectively. Noteworthy, the computed Gaussian energy is greater than that reported in Table I , since higher order corrections to the tanh -solution are present, see Section II, Eq (15) . Performing the same simulation with 1/λ = 40 √ 2π, * = 1.5 h * = 1, and the same dt * , the computed energies at the same final time are E * B ≈ 1.813 and E * G ≈ −1.096 · 10 −2 . The error with respect to the data reported in Table I decreases, since, by reducing the dimensionless thickness λ, the higher order terms become less and less important. In all simulations presented in this Appendix, vesicles area and volume are conserved with the same accuracy reported in Fig. 2 . The functional derivative of the energy (4) is The bending term is well known in literature and its explicit expression is As regards the Gaussian term, noticing that an integration by parts of (7) leads to This simplifies the computation of the functional derivative, which turns out to be In a more readable form: which further simplifies in the axisymmetric case to Furthermore, the functional derivatives of the area (19) and volume (20) are In conclusion, the functional derivative of the augmented energy (24) is with and In the limit of vanishing the field energy will be shown below to converge to the sharp interface energy of the Canham-Helfrich model, see Eq. (1) in the main text, for a membrane identified with the mid surface implicitly defined by the field isosurface φ(x) = 0. In this case, the various parameters entering the expression for the energy will be interpreted as follows: m is the spontaneous curvature of the membrane, assumed to be closed, taken to be positive if the membrane bulges towards the exterior; k is the bending rigidity; k G is the Gaussian curvature modulus (for definiteness, k G = −k will be assumed, as suggested by the experimental data, see the main text for details). Hence the equilibrium configurations of a Canham-Helfrich membrane can be obtained either by minimizing the Canham-Helfrich functional with respect to the membrane configuration or, alternatively, by minimizing the phase-field energy functional in the limit of vanishing . As explained in the main text, apart from practical advantages in the minimization process, the phase-field approach presents the main advantage over the sharp interface description of allowing for the membrane topological transitions. More technically, it is shown below that an ansatz for the phase function of the form where d(·) is the signed distance from the (assigned) membrane midsurface configuration and f is such that f (±∞) = ±1 and f (0) = 0: i) does indeed allow to find a minimizer of the field energy for small enough and given midsurface geometry; ii) leads to a limiting, sharp interface form of the field energy which coincides with the Canham-Helfrich energy, which is a functional of the membrane configuration. From the discussion, it will become clear that the parameter > 0 is related to the width of the diffuse interface and, thus, to the small, but finite, membrane thickness. We choose to define the signed distance such that n = ∇d computed on Γ is the inwardpointing unit normal to the vesicle. It is worth noticing that the signed distance function satisfies the eikonal equation, namely |n| = 1. As anticipated, setting d * (x) = d(x)/ , we also require that lim d * →±∞ φ = ±1 and φ = 0 for d = 0. Therefore, ±1 are the values for the stable phases of the inside and outside bulk and the level set φ = 0 identifies the membrane midsurface. After introducing the membrane characteristic length D ve = A/π, where A is a characteritic membrane area, the sharp-interface limit corresonds to λ = /D ve → 0. Since ∇φ(x) = f (d * (x)) n/ , where the prime denotes the derivative with respect to d * (x), a direct computation leads to where the bar indicates lengths normalized with D ve . Following [1] , in order to minimize E = E B + E G as λ → 0 the phase-field function is expanded as where R λ (d * (x)) = O(λ 2 ), obtaining where the terms O(λ) and O(λ 2 ) in Eqs. (24) and (25), respectively, follow from the property Ω η(d * (x))ξ(x) dV = O(λ), for each η(d * (x)) = O(1) such that η(±∞) = 0 sufficiently fast. Denoting with k 1 and k 2 the principal curvatures of the membrane, we have ∇ · n = −(k 1 +k 2 ) = −2M and n·∇k i = k 2 i , with the result that (∇·n) 2 +n·∇(∇·n) = 2k 1 k 2 = 2G. Now, noticing that for λ → 0 where the limits should be understood in a weak sense and δ(x) is the Dirac delta function, which indeed coincide with the Cahnam-Helfrich functional given as Eq. (1) in the main text. From the above analysis, denoting φ e (x, λ) = argmin E[φ, λ], it follows that the surface Γ e = x : lim λ→0 φ e (x, λ) = 0 is a minimizer of the sharp interface energy whenever Γ e is sufficiently smooth to allow the evaluation of the Cahnam-Helfrich functional. The above reasoning can be extended to the case where the field is subject to the constraint on the vesicle areas and enclosed volumes, with the phase-field area and volume estimators given in Eqs. (19) and (20) in the main text. With reference to Fig. 4 of the main text (reported below for convenience), proceeding to the right along the MEP, beyond the saddle point, the computed Gaussian energy of the neck progressively decreases. In order to get a geometrical understanding of this decrease, for each image along the MEP, we extract the φ = 0 isoline in the r − z plane from the phase-field and find its interpolating polynomial r = r(z). By resorting to the Gauss-Bonnet theorem, the total curvature of the neck (a surface with Euler characteristic χ = 0) is (minus) the integral of the geodesic curvature k g along the neck boundary (made up of two circles of Approximating the vesicle with its sharp-interface counterpart, the neck Gaussian energy is E neck G (Z)/(8πk) ≈ r(±Z)k g /2 = 1 − (r/R n ) 2 /2, where R n = 1/k n . From the insets of the figure, near the saddle point (configuration denoted by the triangle) R n >> r(±Z), leading to k g ≈ 1/r(±Z). Moving away from the saddle point, the two curvature radii become eventually comparable, R n ≈ r(±Z), leading to k g ≈ 0. The main plot in the figure provides the actual Gaussian energy of the neck as computed from the phase-field (blue line) compared to the sharp-interface approximation (orange line with dots). The agreement is excellent as long as different membrane segments do not approach each other to a distance smaller than the bilayer thickness. Indeed, the agreement progressively deteriorates when getting closer to the saddle point, due to the increasing curvature of the membrane generatrix. Actually, one could better and better reproduce the sharp interface energy also in this stage by reducing the regularizing parameter λ. On the other hand, from a physical point of view, the thickness of the bilayer is finite, making the sharp interface model inappropriate when the saddle point is approached. Overall, these results confirm the accuracy of the proposed phase-field expression for the Gaussian curvature. Viral membrane fusion and the transmembrane domain Viral membrane fusion Membrane and organelle dynamics during cell division Protein-lipid interplay in fusion and fission of biological membranes Lipid nanoparticles-from liposomes to mrna vaccine delivery, a landscape of research diversity and advancement Lipid nanoparticles for mrna delivery Engineering precision nanoparticles for drug delivery The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell Elastic properties of lipid bilayers: Theory and possible experiments Recent developments in the field of bending rigidity measurements on membranes Fluid lipid membranes: From differential geometry to curvature stresses Elastic torques about membrane edges: A study of pierced egg lecithin vesicles Molecular asymmetry and saddle-splay elasticity in lipid bilayers Determining the gaussian curvature modulus of lipid membranes in simulations The giant vesicle book Free-energy calculation methods for collective phenomena in membranes Budding and fission of nanovesicles induced by membrane adsorption of small solutes Membrane fusion and drug delivery with carbon nanotube porins The mechanism of vesicle fusion as revealed by molecular dynamics simulations Molecular dynamics simulations of lipid vesicle fusion in atomic detail Tensioninduced fusion of bilayer membranes and vesicles Pathway of membrane fusion with two tension-dependent energy barriers Thermodynamically reversible paths of the first fusion intermediate reveal an important role for membrane anchors of fusion proteins Lipid intermediates in membrane fusion: formation, structure, and decay of hemifusion diaphragm Point-like protrusion as a prestalk intermediate in membrane fusion pathway The gaussian curvature elastic modulus of n-monomethylated dioleoylphosphatidylethanolamine: relevance to membrane fusion and lipid phase behavior The gaussian curvature elastic energy of intermediates in membrane fusion Minimum membrane bending energies of fusion pores Membrane fission: model for intermediate structures Solveig Mareike Bartelt, Seraphine Valeska Wegner, Rumiana Dimova, and Reinhard Lipowsky Budding and fission of vesicles Endocytic sites mature by continuous bending and remodeling of the clathrin coat Lipid bilayer vesicle fusion: intermediates captured by high-speed microfluorescence spectroscopy Proceedings of the National Academy of Sciences of the United States of Time scales of membrane fusion revealed by direct imaging of vesicle fusion with high temporal resolution Energetics of stalk intermediates in membrane fusion are controlled by lipid composition Secretory and viral fusion may share mechanistic events with fusion between curved lipid bilayers Low energy cost for optimal speed and control of membrane fusion Phase-field theories for mathematical modeling of biological membranes Toward a thermodynamically consistent picture of the phase-field model of vesicles: curvature energy A phase field approach in the numerical study of the elastic bending energy for vesicle membranes Modeling the spontaneous curvature effects in static cell membrane deformations by a phase field formulation Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions Dynamic model and stationary shapes of fluid vesicles Simulating vesicle-substrate adhesion using two phase field functions A two phase field model for tracking vesicle-vesicle adhesion Modelling and simulations of multi-component lipid membranes and open membranes via diffuse interface approaches Model for curvature-driven pearling instability in membranes Rheology of red blood cells under flow in highly confined microchannels: I. effect of elasticity Collective behavior of red blood cells in confined channels The dynamics of shapes of vesicle membranes with time dependent spontaneous curvature Retrieving topological information for phase field models On gaussian curvature and membrane fission String method for the study of rare events Wetting transition on patterned surfaces: transition states and energy barriers How crystals form: A theory of nucleation pathways Nucleation and growth dynamics of vapour bubbles Water cavitation from ambient to high temperatures Theory and algorithms to compute helfrich bending forces: a review Asymptotic analysis of phase field formulations of bending elasticity models See supplemental material at [url will be inserted by publisher] for: i) the whole phase-field model asymptotic analysis; ii) a detail discussion about the gaussian energy of the neck δ(x) is the Dirac delta function and W denotes a weak limit in the sense of distributions The string method as a dynamical system String method in collective variables: Minimum free energy paths and isocommittor surfaces Simplified and improved string method for computing the minimum energy paths in barrier-crossing events Shape transformations of vesicles: Phase diagram for spontaneous-curvature and bilayer-coupling models The 2018 biomembrane curvature and remodeling roadmap Common energetic and mechanical features of membrane fusion and fission machineries Gaussian curvature directs the distribution of spontaneous curvature on bilayer membrane necks Gaussian curvature is G j = n · ∂n/∂ξ 1 × ∂n/∂ξ 2 , where j is the Jacobian and n the unit normal. Introducing the Gauss map, that associates the point x = x(ξ 1 , ξ 2 ) ∈ M to a corresponding point n(ξ 1 , ξ 2 ) on the unit sphere S, one finds M GdS = S(M) dΩ, where Ω is the (signed) solid angle and S(M) is the image of M through the map Differential geometry of curves and surfaces: revised and updated second edition Membrane curvature controls dynamin polymerization Diffuse-interface methods in fluid mechanics The sharpinterface limit of the cahn-hilliard/navier-stokes model for binary fluids Cluster analysis and mathematical programming Overview on techniques in cluster analysis The biophysics and cell biology of lipid droplets Complex fluidfluid interfaces: rheology and structure A constrained string method and its numerical analysis PETSc/TS: A Modern Scalable ODE/DAE Solver Library Colliding interfaces in old and new diffuse-interface approximations of willmore-flow Vesicles of toroidal topology with g = g(d * (x)) and h a combination of functions factorized in the form η(d * (x))ξ(x).Therefore, both g and h are O(1) in λ. Expansion (9) entails the ordering of the two partial energy functionals as power series in the small parameter λ. Therefore, in order to minimize the otherwise dominating bending energy (order 1/λ 3 ), the leading order term f 0 must satisfyi.e.,which identifies as directly proportional to the interface width. Note that, once f 0 is determined, the order λ −1 contribution to the Gaussian energy is fixed. Concerning the bending energy, we are left withwhere we have usedandh = O(1) shares the same structure of h.Now, considering thatone finds thatAn integration by parts leads toSince the contribution of the above integral is of higher order, substitution in (14) shows that minimizing the λ −1 -order bending energy term is tantamount to solving the equationTwo independent solutions are easily found. Given Eq. (16), f 1a = f 0 is clearly one of them, while the other can be sought in the form f 1b (d * ) = w(d * )f 0 (d * ). Thus, Equation (18) implieswhich, after setting w = z and still considering Eqs. (12) and (15), is reduced toThe solution is z(and thereforewhere, due to the factors √ 2 in the arguments, lim d * →±∞ f 1b (d * ) = ±∞. The two constants in the general solution of (18),are determined as follows. Since f 0 (±∞) = ±1, B = 0 by requiring that lim d * →±∞ φ = ±1, and, given that f 0 (0) = 0, A = 0 by requiring that φ = 0 on the membrane midsurface.Therefore, f 1 ≡ 0, and we are left with