1 Introduction
Lipid bilayers are fundamental structural elements in biology. Vesicles are compartments enclosed by lipid bilayer membranes and they have been studied extensively over the last four decades Deuling1976a ; Seifert1997 ; Deserno2015 ; Tu2018 . Notably vesicles have been the subject of investigations that were awarded the 2013 Nobel prize in the physiology or medicine (James Rothman, Randy Schekman and Thomas Südhof). A vesicle has a relatively simple structure when compared to other entities at cellular level and they have been the subject of several detailed theoretical Canham1970 ; Helfrich1973a ; Sheetz1974 ; Evans1974 ; Svetina1989 ; Seifert1992 ; Bozic1992 ; Wiese1992 and experimental studies Sackmann1986 ; Kaes1991 . Models of vesicles are readily transferable to more complex structures, as lipid bilayer membranes enclose most cells and their nuclei as well as many viruses. Importantly, a lipid bilayer and a spectrin network constitute a red blood cell (RBC) membrane Evans1980 ; Mohandas1994 .
Lipid membranes are a few nanometers thin, but they span surfaces that are in the order of micrometers. The bending elasticity of the quasi twodimensional structures is caused by relative extension/compression of outer/inner surface, which can be expressed as a function of the membrane curvature Seifert1997 ; Deserno2015 ; Tu2018 . Canham proposed the first form of membrane bending energy, that depends quadratically on the two principal curvatures of the surface Canham1970 . This representation, along with constraints on total area and volume, constitutes the minimal model for a vesicle. Helfrich formulated a bending energy that includes a Gaussian curvature as well as a ”spontaneous” curvature Helfrich1973a that is attributed to chemical differences in two lipid layers. This is the spontaneous curvature (SC) model and by setting , it is equivalent to the minimal model. In turn, Sheetz & Singer Sheetz1974 and Evans Evans1974 accounted for the density difference of the two lipid layers. Their work led to the bilayer couple (BC) model Svetina1989 that introduces a constraint on the area difference so that , where and are the instantaneous and targeted area differences between the two lipid layers. A few years later three groups (Seifert1992 ; Bozic1992 ; Wiese1992 ) proposed that the nonlocal bending energy has a quadratic form . Here reflects reference of the area difference between the two layers. This is areadifference elasticity (ADE) model, which degenerates to the minimal, SC or BC model as a special case. The ADE model has been validated experimentally on vesicles Doebereiner1997 ; Sakashita2012 .
Despite the well established formulation of these membrane models, simulations of vesicles and redblood cells remain a challenging task. The key difficulty lies in the evaluation of the resulting forces. The bending force depends on the LaplaceBeltrami operator on the mean curvature, which further relates to of the coordinates on a surface. Hence, the evaluation of the bending force requires approximations of the fourth derivatives of the surface coordinates. A recent review showcased how computation of the bending energy and especially the force remain elusive, even for the minimal model Guckenberger2017 . These computations have similar difficulties as those met when calculating Willmore flow in differential geometry Deserno2015 ; Guckenberger2017 , which are of particular interest to the computer graphics community crane2013robust .
One approach has focused on axisymmetric shapes that reduces the calculations to solving Ordinary Differential Equations (ODEs)
Deuling1976a ; Seifert1991 ; Miao1994 ; Tu2018 . Recent works have employed either spherical harmonic expansions/quadratic approximation to represent the surface Heinrich1993 ; Khairy2008a ; Veerapaneni2011 ; Yazdani2012 or phase field Du2004 models. We note in particular level set methods for the description of the membrane surface Maitre2009 ; bergdorf2010lagrangian ; Mietke2019 thus allowing for volume based calculations of the membrane dynamics and large deformations even with topological changes.In this work, we represent the surface as a triangulated mesh, which has been used extensively to model vesicle and RBC membranes and they have been coupled to fluid solvers Noguchi2005a ; Pivkin2008 ; Fedosov2010b ; Tsubota2014 ; Li2013 ; Guckenberger2016 . However, the majority of vesicle and RBC simulations have employed only the minimal model. Exceptions are Monte Carlo simulations from groups of Seifert Jaric1995 and Wortis Lim2002 and the finite element method from Barrett et al Barrett2008a . We examine and extend four established schemes for the discretization of the general bending model on a triangulated mesh Kantor1987a ; Juelicher1996 ; Gompper1996 ; Meyer2003 . We perform a number of benchmark problems to compare their accuracy, robustness and stability.
The paper is structured as follows: we describe energy functionals of the minimal, SC, BC and ADE models in Section 2 and present an expression for the force density. In Section 3, we present the four discretization schemes (A, B, C, D) which have been previously applied to the minimal model Kantor1987a ; Juelicher1996 ; Gompper1996 ; Meyer2003 . We extend schemes B, C, D to incorporate the spontaneous curvature and nonlocal bending energy thus accessing the SC, and BC/ADE models. The details of the derivations on Sections 2 and 3 can be found in A–G. In Section 4.1, we evaluate the proposed schemes on shapes with known analytical expressions for the energy and force. In Section 4.2 we compute the phase diagrams for oblatestomatocyte, oblatediscocyte, prolatedumbbell, prolatecigar shapes described by the minimal model with the four discretization schemes, from different initial shapes such as prolate, oblate and sphere. Thereafter, we generate slices of phase diagrams of SC, BC and ADE models. We compare our results with solutions obtained by axisymmetric ODEs, spherical harmonic as well as results from the program ”Surface Evolver” and experiments. In Section 4.4 we present results from dynamic simulations. We summarize our findings in Section 5.
2 Continuous energy and force
2.1 Energy functionals
The primary model for the twodimensional membrane elasticity is based on an extension of one dimensional beam theory by Canham Canham1970 with an energy functional:
(1) 
where is the bending elastic constant, and are the two principal curvatures^{1}^{1}1We take a convention that and
are positive for a sphere. With this convention, the surface unit normal vector points inwards.
and the integral is taken over a twodimensional parametric surface embedded in threedimensional space.The Canham energy functional is scaleinvariant, with smaller vesicles corresponding to larger curvatures. This energy functional is known as the minimal model of a vesicle membrane. Its phase diagram is determined by the reduced volume between a vesicle and a sphere , with the same area and is the sphere’s radius.
A few years later, inspired by research in liquid crystal, Helfrich Helfrich1973a developed an alternative free energy functional for lipid bilayers
(2) 
where is the mean curvature and is the spontaneous curvature, reflecting the asymmetry in the chemical potential on the two sides of the lipid bilayer. is the Gaussian curvature, and is a bending elastic constant. According to GaussBonnet theorem , where is the genus of the surface. For the spherical topologies considered here, and the energy term is omitted. The remaining term in Eq. 2 constitutes the spontaneous curvature (SC) model Helfrich1973a . We note that for , Eq. (1) has the same dynamics as that of Eq. (2). The Helfrich energy with is also scaleinvariant. However, a nonzero introduces a length scale and we define .
It is important to remark that most simulations of vesicles or RBCs Li2013 ; Tsubota2014 ; Guckenberger2016 use , thus assuming no chemical differences between the two sides of the lipid bilayer. However it has been noted Seifert1997 that such an assumption lacks a physical realization. We denote the energy due to spontaneous curvature as
(3) 
The first term resembles surface tension energy encountered in models of multiphase flow Landau1987 with is a ”surface tension” constant ^{2}^{2}2Ideally, a lipid membrane is incompressible and the total area of the neutral surface is constant. However, numerically we treat the constraint on area as penalization and therefore, we also consider the energies and forces due to total area arising from SC and ADE models..
In the early 90s, three groups Seifert1992 ; Bozic1992 ; Wiese1992 proposed an additional term to the energy functional that reflects the area difference between the outer and inner leaflets of the lipid membrane. This nonlocal term is areadifference elasticity (ADE) and it is expressed as:
(4) 
Here is a bending elastic constant due to area difference and is the thickness of the bilayer. The ratio is in the order of unity Bozic1992 ; Seifert1992 ; Miao1994 ; Waugh1995 ; Lim2002 and depends on the properties of the lipid.
In summary, the Helfrich model extended by the ADE terms, is the so called ADE model Seifert1997 that can be expressed as:
(5) 
Setting in the ADE model and constrain we obtain the bilayercouple (BC) model Svetina1989
that includes the area difference as a constraint. We simplify the formulation by introducing the zero, first and second moments of the mean curvature
(6) 
We rewrite the total energy as
(7) 
Here the ”” indicates that the terms originate from either or , and we have used the expression .
2.2 Force from calculus of variation
We use the energy functional and the calculus of variations to derive (details are shown in A) the force magnitude corresponding to each moment as
(8) 
Here is the surface gradient operator and is the LaplaceBeltrami operator.
We use the above expressions to write the magnitude of force density acting along the normal as
(9) 
.
These terms can be rearranged to show the components of the force density corresponding to the energy functionals as
(10) 
3 Discretization schemes
We represent the membrane as a triangulated mesh and use the following notation (see Fig. 1):

are indices for vertices

is index for the edge connecting vertices and .

is index for the triangle connecting vertices and .

is the angle between the normal vectors of two triangles sharing edge .

and are angles associated with vertices and . is the angle associated with vertex within triangle .

are the numbers of vertices, edges and triangles.
The total energy is
(11) 
where is the area difference.
We introduce the discrete moments as ^{3}^{3}3We do not differentiate continuous and discrete symbols if they are clear from the context.
(12) 
so that the total energy becomes
(13) 
where is nondimensional.
It is apparent that from Helfrich energy and from ADE energy may compensate each other’s mechanic effects, although the origins of the two are different. Therefore, we take their combined effect represented by Lim2002 ; Khairy2008a . We assume and in the ADE model so that the resultant from a simulation of the ADE model can be nondimensionalized as . Therefore, for a given area , we have three parameters , , and . Other work have adopted and to nondimensionalize the area differences, where the denominator corresponds to the area difference of a sphere Seifert1991 ; Ziherl2005 . This approach has the advantage of eliminating the membrane thickness . To compare with different references, we keep both means of nondimensionalization and they are related as and .
We consider four discretization schemes which have been applied to the minimal model in the past. Note that scheme A defines the energy without defining and . In turn, for schemes B, C, and D, there are explicit definitions of and . In the present work we extend these three schemes to account for all energy terms in Eq. (13) and their corresponding forces.
3.1 Discrete force
ADE  momenta  

A Kantor1987a  
B Juelicher1996  
C Gompper1996  
D Meyer2003 
For schemes A, B, and C, the force at a vertex is
(14) 
We may convert between force and force density associated with each vertex as
(15) 
where is the area associated with vertex . The former is used by scheme A, B, and C, to convert force to force density for comparison with calculus of variation in Section 2.2. For scheme D, the force density is computed directly and converted to force to use in the energy minimization process. A summary of the discretization schemes is shown in Table 1 and we elaborate each scheme in the following.
3.2 Scheme A
This scheme was invented by Kantor & Nelson Kantor1987a to simulate planar polymeric network in Monge form. This model was later broadly adopted to model generally curved membranes, in particular for red blood cells Discher1998 ; Li2005 ; Pivkin2008 ; Fedosov2010b ; Vliegenthart2011 . The energy depends on each angle formed by the normal vectors of two triangles, which share the same edge (Fig. 1(a)). The total energy is
(16) 
where is the bending modulus. We split the local energy evenly onto two vertices and . Alternatively, splitting the energy evenly onto the four vertices , has marginal differences. There is no difference on the force between the two ways of splitting. is the spontaneous angle, that is, for , . If , by Taylor expansion . Therefore, the approximation of the energy reads Krueger2012
(17) 
A few notes are in order; for a triangulated cylinder surface, can be derived Seung1988 , but it is not true for other shapes Gompper1996 , as demonstrated later. For a mesh of uniform equilateral triangles on a sphere, the is related to the via a resolution parameter (e.g, the edge length of the triangle). However, and have no correspondence in general. Furthermore, there is no definition of or in this scheme and therefore, it is ambiguous to extend it to include ADE.
We will consider scheme A in its original form with (and its approximation) and compare it with the minimal model. The expression for the discrete force is given in B and the only primitive derivative to compute is .
3.3 Scheme B
This scheme was proposed by Jülicher Juelicher1994a ; Juelicher1996 to calculate bending energy of a vesicle with topological genus more than zero. To motivate the scheme we consider two triangles parallel to the original two neighboring triangles as sketched on Fig. 2. Each parallel triangle is separated from its counterpart with distance in the normal direction. The shared edge of the two triangles and the two neighboring edges of the imaginary triangles form a fraction of a cylinder, as indicated on Fig. 2. We define the first moment of mean curvature on the fraction of the cylinder as
(18) 
The two principal curvatures of the cylinder are and . The surface area for the fraction of the cylinder is , where is the length of the edge and as indicated on Fig 1(a). We note that does not depend on and has a proper limit as . On the contrary, if we had defined local mean curvature or second moment of mean curvature directly instead, they would both behave singular as 1/h. The total first moment of mean curvature reads
(19) 
For the local first moment associated with a vertex instead of edge we have
(20) 
where the prefactor is to account for the total summation. For each vertex , we have relation . This leads to a definition of as
(21) 
where summation runs over all neighboring edges around vertex . The local area associated with vertex is based on the barycentric centers of the surrounding triangles as sketched on Fig. 3(a) and it reads
(22) 
where is the area of a neighboring triangle . This means that each triangle split its area evenly onto its three vertices. The summation runs over all neighboring triangles around vertex .
3.4 Scheme C
This scheme was proposed by Gompper & Kroll Gompper1996 to simulate a triangular network of fluid vesicles and later applied to model general membranes Noguchi2005a . Its discrete definition of LaplaceBeltrami operator reads Itzykson1986 ; Pinkall1993 ; Meyer2003 :
(24) 
where is an arbitrary function. The Voronoi area associated to vertex indicated on Fig. 3(b) and is defined as
(25) 
If we consider the discrete counterpart of an expression for mean curvature known from differential geometry, that is,
(26) 
and take function as the surface coordinate in Eq. (24), we have an explicit discrete definition of . When only the Helfrich energy/force with is considered, Gompper1996 ; Guckenberger2016 ; Hoore2018 , as is along direction. However, for the general energy/force in the SC and BC/ADE models considered in this work, we need an explicit discrete definition of . We consider a discrete Thuerrner1998 ; Jin2005a ; Guckenberger2016 . where the unit normal vector for vertex is the sum of neighboring normal vectors weighted by incident angles Thuerrner1998 ; Jin2005a .
(27) 
where is the angle at vertex of triangle , is the unit normal vector of triangle .
3.5 Scheme D
This scheme was derived from an effort of developing a ”unified and consistent set of flexible tools to approximate important geometric attributes, including normal vectors and curvatures on arbitrary triangle meshes” Meyer2003 . It has been adopted by several groups to study the bending mechanics of the minimal model for vesicles and RBCs Boedec2011 ; Tsubota2014 ; Guckenberger2016 . The definition of the discrete LaplaceBeltrami operator reads as Meyer2003 ,
(28) 
which has almost identical expression as the one from scheme C, except that the area for vertex is calculated as a mixture of two approaches Meyer2003 (see Fig. 3(c)).If the triangle is nonobtuse then local areas , , and on three vertices take the contribution from Voronoi tessellation of , If the triangle is obtuse then the tessellation relies on the middle point of each edge. Therefore, the triangle contributes to , and to and each Meyer2003 .
The key difference of scheme D from the other three schemes is in terms of the force calculation. The force density of scheme D relies on the variational expression in Eq. (9) or (10). Once the discrete mean curvature is calculated at each vertex, we apply Eq. (28) second time on discrete values of to get . Furthermore, only scheme D needs a discrete definition of Gaussian curvature at vertex Polthier1998
(29) 
which employs all incident angles around vertex and its area definition. Therefore, each discrete force density in Eq. (8) is readily available and thereafter the total force density can also be obtained.
4 Results
We present the results of the comparative study for all four schemes and the proposed extensions. We distinguish applications on prescribed shapes and on dynamically equilibrated shapes.
4.1 Prescribed configurations: sphere and biconcave oblate
We first examine the performance of the four schemes on the calculation of the Helfrich energy for a sphere and a biconcave oblate described by an empirical function Evans1972 using and . The sphere is approximated with an icosahedron and higher resolutions are obtained by applying Loop’s subdivision scheme Loop1987 for triangulated meshes with , , , . The same resolution is obtained for a biconcave oblate by transforming the surface coordinates of the sphere to the empirical function of the biconcave shape. We present the results in Fig. 4 verifying the total energy converges (albeit to different values !) with increasing the number of triangles. For the sphere, all four schemes converge with increasing the number of triangles (Fig. 3(a)). The firstorder approximation of scheme A deviates significantly from scheme A when when the angle between two neighboring triangles is not small. Schemes B, C and D converge to a correct value .
For the biconcaveoblate configuration we compute the reference solution analytically (see E). We find that all four schemes converge to the reference solution with increasing number of triangles ( Fig. 3(b)). We note that the results from the firstorder approximation for scheme A deviate significantly from those of the complete scheme A for . Moreover, with , scheme A and its simplified version again converge to a higher value than the one given by the reference solution. The value of the total energy for the prescribed shapes does not depend on the details of the triangulated mesh, as long as the mesh is regular.
Subsequently, we compute total energy with , and using scheme B, C and D ( Fig. 5). All three schemes converge with increasing the number of triangles to the reference values. As shown in Fig. 5, for low resolution , scheme B deviates much less than schemes C and D from the reference values.
Furthermore, for the biconcaveoblate configuration we plot energy density versus the radial distance from the axis of symmetry ( Fig. 6). Results from scheme B, C and D with coincide with the reference line except a few discrepancies, which disappear with high resolution . Scheme A and its simplified version show no difference, but deviate from the reference values. Increasing the resolution from to does not improve the results of scheme A or those of its simplified version. The results indicate that scheme A is not an effective discretization of the Helfrich energy,consistent with the findings of two recent works Tsubota2014 ; Guckenberger2016 .
Finally, we examine the energy density for the Helfrich energy with , as computed by schemes B, C and D. We present the results in Fig. 7
. All three schemes capture accurately the profile of the reference line. Some outliers showing for
also disappear for higher resolution . Note that the ADE energy is nonlocal and has only one global value, which was already presented in Fig. 5.4.2 Numerical results on equilibrium shapes
Here we consider the equilibrium shapes of a closed membrane. The energy minimization process is initialized from a prescribed shapes such as a sphere, prolate and oblate ellipsoids, or from the equilibrium shape (e.g. stomatocyte) of an optimization using different parameters.
The minimization process relies on the calculation of the forces corresponding to the associated energy. The energy and the force are both calculated on the triangulated mesh. We set and Lim2002 ; Khairy2008a . We penalize the constraints on global area with (F) and volume with (G) Du2004 ; Li2005 ; Pivkin2008 ; Fedosov2010b
. We also add an inplane viscous damping force between any two neighboring vertices connected by an edge in order to dampen the kinetic energy. To explore the energy landscape efficiently, we also add a stochastic force of white noise between any two neighboring vertices connected by an edge. The pair of viscous and stochastic forces between vertices resemble the pairwise thermostat
Espanol1995 , which conserves momenta and has a proper thermal equilibrium. The time integration is performed by the explicit velocity Verlet method. As the minimization evolves, the configuration of the triangulated mesh may deteriorate, e.g., elongate in on direction or generate obtuse angles. We remedy this mesh distortion by introducing two regularization schemes: by introducing a constraint on the local area with as penalization (F) Li2005 ; Pivkin2008 ; Fedosov2010b or by triangle equiangulation Brakke1992 . We note that as we perform multiple runs from different initial shapes we may reach different local minima of the energy landscape. Thereafter, we select the smallest among all available local minima as the global minimum, that is, the ground state.4.2.1 The minimal model
We first consider the minimal model described by Helfrich energy with Deuling1976a ; Seifert1991 . As the equilibrium shapes are axisymmetric, the computations are reduced to solving twodimensional EulerLagrangian ODEs. We aim to reproduce the known phasediagram with configurations and energies, in particular for small reduced volumes Peng2010 . Furthermore, as the minimal model is widely used to simulate vesicles and RBCs we examine the performance of all four schemes in this situation.
We present the phase diagram of the reduced volume and the normalized energy of the vesicles as generated by scheme B in Fig. 8. Results from schemes A, C and D are quite similar to those obtained from scheme B, so in the following we only emphasize their discrepancies. We adopt the reference lines from the work of Seifert et al Seifert1991 and denote our results with symbols. Each symbol corresponds to the energy minima as obtained from our minimization procedure.
The reference lines correspond to three types of configurations. With decreasing in the prolate branch, the shape changes from sphere to prolate, dumbbell, and long caped cylinder. With decreasing in the oblate branch the shape changes from a sphere to a famous biconcaveoblate shape at around . For the oblate branch bifurcates to form two stomatocyte subbranches. One subbranch has constant energy for while the other has energy which depends on volume.
The simulation results denoted by symbols in Fig. 8, are obtained by three initial shapes: sphere, prolate and oblate ellipsoids. The volume for the initial prolate and oblate ellipsoids are the same as that of the target . However, the prolate and oblate branches generated in this work are different from the reference work Seifert1991 . The reference restricts the prolate branch to only prolatelike shapes and the oblate branch to oblatelike shapes. Here, we do not impose such constraint on the shape during the minimization procedure. More specific, all the squared symbols are from initial shape of prolates and they all stay on the prolate branch when reaching the local energy minimum. Similarly, all the circle symbols are from initial shapes of oblates, with the corresponding , that stay on the oblate branch for and jump onto the prolate branch for . This results corroborates an earlier finding Jaric1995 where the oblate branch is only locally stable for . For the minimization from a spherical initial shapes (up triangle symbols), the energy minima are found on the prolate branch for , on the oblate branch for and scatter between prolate and oblate branches for .
For the constant energy subbranch of the stomatocyte, we can not explore the partial () energy landscape directly. Therefore, we take two final shapes of stomatocyte, which are obtained by minimization with initial shapes of oblates. One final shape is from and the other is from . We use these two shapes as initial shapes to run minimization procedure with targeted reduced volume ranging from to . The results from are denoted with solid down triangles while from are solid diamonds.
We demonstrate the convergence of the minimization process for the calculated energies, by examining the results in , as this range contains the regime of oblatelike biconcave shapes () that is relevant for the simulations of RBCs. We present results from schemes B, C and D with two resolutions and in Fig. 9. It is apparent that scheme B is superior, as it already captures the reference values with , whereas scheme C and D reproduce the same results with . Moreover, the energies at delimiting points, that is, for and for are reproduced quite well by all three schemes. The minimal energies computed by scheme A are removed, as they are completely out of range. This should not be surprising, given its performances on two static configurations in section 4.1. Nevertheless, we present the equilibrium shapes by scheme A along with the other three schemes.
A  

B  
C  
D 
We select four equilibrium configurations in the prolatelike/dumbbell/cigar regime with generated by the four schemes with in Table 2. We find very small differences in the shapes obtained by the four schemes.
A  

B  
C  
D 
However, in the biconcaveoblate regime for ( Table 3), shapes obtained by scheme A are different from those of the other three schemes, consistent with recent work Hoore2018 . With an oblatelike biconcave shape as initial shapes, we reproduce the same sequence of shapes by scheme A as those shown in the supplementary movie of Hoore2018 . However, the run time of Hoore2018 is too short to observe the equilibrium configurations. Moreover, the final (local) equilibrium shape by scheme A is also dependent on the initial shapes. Each shape by scheme A in Table 3 is selected as the smallest energy among steady states of three minimizations, each of which has one initial shapes among sphere, prolate and oblate. The equilibrium shapes by scheme A are always prolatelike/dumbbells/cigars, which are completely different from the analytical results in the reference study Seifert1991 . Schemes B, C, and D all produce the biconcave oblates with no distinct differences among them.
B:  

B:  
B:  
C:  
C: 
At smaller reduced volume , scheme A does not generate stomatocyte, but always ends up with vesiculations, that is, two sphere like shapes connected via thin tubes (not shown here). Schemes B and C generate equilibrium shapes consistent with what is expected by the reference study as shown in Table 4 for two and three different resolutions. Scheme B is stable even at low resolution of , although the equilibrium shapes are not as accurate. Especially for and , scheme B with does not lead to the proper stomatocyte. Instead the two spheres are connected externally via a line segment, which can be seen from two orthogonal views (the first two entries in the first row of Table 4). These two shapes are caused by the underresolved neck connecting the two spheres at this low resolution. All equilibrium shapes generated by scheme B with are accurately axisymmetric stomatocyte and show virtually almost no difference from the results of . These results are comparable to the twodimensional contours reported in the reference Seifert1991 . The results of scheme C are shown in the last two rows of Table 4 for and . Especially for the most challenging case of considered here, scheme C leads to unphysical shapes. Even for , scheme C cannot separate the inner sphere from the outer sphere with a reasonable distance, due to poor resolution at the neck. The results of scheme D are very similar to those of scheme C as the minimization procedure evolves for each target . However, once scheme D reaches quasisteady state with a proper stomatocyte shape, the overall configuration oscillates and the shape becomes unstable. We tried to employ regularization and smoothing to stabilize scheme D, but we were unsuccessful. We speculate that since scheme D does not enforce the conservation of momentum or energy this may cause this instability. To the best of our knowledge, there are no reports from previous work using triangulated mesh to simulate such low volumes. The exact reason why scheme D is unstable remain obscure.
4.2.2 Spontaneous curvature model
We compute the minimal energies determined by the SC model and the corresponding equilibrium configurations. In particular, we consider two spontaneous curvature values and so that we can compare with results of Seifert et al Seifert1991 . The possible shapes generated by the SC model are still axisymmetric. However, the SC model enables quite a few more configurations that are numerically challenging. These include pearlike shapes which break the updown symmetry, budding transition with a narrow neck connecting the neighboring compartments of the same size, budding transitions with the neighboring compartments of different sizes, vesiculation where connection between the neighboring compartments of the same size is lost and vesiculation where connection between the neighboring compartments of different sizes is lost. There are also discontinuous transitions for , but the minimization procedure of exploring the energy landscape is quite similar to the case of and we do not elaborate on this finding here.
4.2.3
B  

C  
D 
We first consider and present the energy diagram in Fig. 9(a) and the representative equilibrium shapes in Table 5. The references are adapted from Seifert et al. Seifert1991 . For the regime , all three schemes accurately capture the shapes of necklace (with two necks) at , dumbbells at and budding at , as shown in the first three columns of Table 5. The minimum energies in this regime computed by the three schemes are also quite close to the reference line. Especially scheme B accurately reproduces the energy line with as shown in the left part of Fig. 9(a). For the regime , we expect prolate/dumbbell shapes as suggested by the reference. All three schemes B, C and D are able to get the minimum energies accurately as shown in the right part of Fig. 9(a), and to generate the expected shapes comparable to the reference (see figure 12 of Ref. Seifert1991 ) as shown in the last two columns of Table 5. The most challenging regime resides in , since there are multiple local minima and the energy diagram bifurcates. However, all three schemes are able to produce the budding, where there is only one very thin neck connecting the two compartments of different sizes. For this very extreme regime, we do not expect accurate computation on the minimum energy. However, the energies by scheme B and D are still surprisingly close to the reference, whereas the scheme C is completely off due to the poor representation of the narrow neck with only a couple of triangles.
4.2.4
B  
C  NA  NA  
D  NA  NA  NA 
Computations for the case of are even more challenging as, besides budding transitions, vesiculations are also expected. For the regime , all three schemes are able to generate necklaces as shown in the first two columns in Table 6 and compute minimum energy accurately as shown on the left part of Fig. 9(b). There is notable discrepancy for the equilibrium shape generated by scheme C at , as it cannot sustain the axissymmetry due to its inability to resolve the narrow necks of the budding shapes. For the regime , all three schemes can generate dumbbell shapes at equilibrium, for example, at . However, only scheme B can generate vesiculation of two spheres of equal size expected from the reference Seifert1991 , as shown by one example of on Table 6. Scheme C at may generate vesiculation, but the two compartments do not have equal size and the energy differs from the reference. Perhaps the regime of is the most challenging one, since the two compartments are expected to have updown symmetry breaking for the vesiculation. We find that both schemes C and D to run stably in this regime of vesiculations. On the contrary, scheme B delivers the correct results with a satisfactory accuracy on both the energy values and configurations.
To summarize this section, scheme B with is able to generate the whole spectrum of equilibrium shapes including dumbbells, necklaces, even buddings and vesiculations where only a couple of triangles are connecting different compartments. Furthermore, the corresponding energy value for each equilibrium shape is very accurate in comparison with the reference value. Both schemes C and D with run into troubles when budding or/and vesiculations are expected. With higher resolutions , neither scheme C nor D show significant improvement (results not presented herein).
4.2.5 Bilayercouple model
We consider the phase diagram of the BC model which can be realized by taking in Eq. (5) so that the areadifference elasticity becomes a quadratic penalization. In this model, there are two nondimensional parameters and
. We consider the BC model before presenting results from the complete ADE model for a few reasons: 1. The BC model has only continuous phase transitions
Svetina1989 ; Seifert1991 . 2. The equilibrium shapes and possible resultant in BC model contains those of the ADE model Seifert1997 . 3. The phase diagram of BC has two controlling parameters, one less than that of ADE model. 4. With , we can focus on examinations the areadifference elasticity.We consider the same parameter ranges of and as in Ref. Ziherl2005 , where the program Surface Evolver was used Brakke1992 . We present our phase diagram by scheme B in Fig. 11, where almost each individual shape has onetoone correspondence to the shape on Fig. 1 of Ziherl and Svetina Ziherl2005 . As in the SC model, we obtain similar axisymmetric shapes such as stomatocytes, prolatelike/dumbbells, budding necklaces. Different from the SC model, we obtain in addition many nonaxisymmetric shapes from the constraint of the area difference between outer and inner leaflets. These shapes include elliptocytes (e.g., the fourth one at or ), rackets (e.g., the fifth one at and sixth one at ), triangle oblates with equal sides (e.g., the sixth one at ), and triangle oblates with unequal sides (e.g., the seventh one at ). All these shapes have been recently validated by experiments on selfassembled vesicles combined with image analysis Sakashita2012 .
Furthermore, we verify the corresponding energy values as shown in Fig. 12, where the reference values are adapted from Ziherl and Svetina Ziherl2005 . We observe that energies computed by scheme B with and exhibit negligible discrepancies and they both follow closely the reference values. Therefore, we omit the equilibrium configurations generated with .
We have also computed the same phase diagram by scheme C and D with and . However, scheme C and scheme D become unstable for certain values of . We speculate that the vulnerability of scheme C is due to the definition of the discrete normal vector, which is not compatible with the discrete definition of the Laplace operator. We maintain that the instability of scheme D is due to its lack of conservation properties. We consider that a more rigorous investigations on scheme C and D are beyond the scope of this work.
4.3 Areadifferenceelasticity model
Here we compute the total energy and force in the ADE model. For comparison, we consider as reference the work of Khairy et al. Khairy2008a , which employs spherical harmonics to represent the surface and calculate the total energy. Therefore, we take the same reduced volume and , which represent the parameters of a normal RBC Lim2002 ; Khairy2008a .
We present the energy versus reference area difference on Fig. 12(a), where we include results of all three schemes and each with two resolutions of and . For each considered from to , there is only small discrepancy between the schemes with a proper resolution, that is, scheme B with and and scheme C and D with . We present the resultant versus on Fig. 12(b) and observe that increases as for smaller . It stays almost at one plateau for and at another plateau for . Correspondingly, we plot again the energy versus the resultant on Fig. 12(c). We observe that the associated energies span among possible area differences for and cluster at two discrete states of area differences for . These characteristics of area differences and minimal energies reflect the possible equilibrium configurations of the vesicles, which will be presented next.
B:  
B:  
C:  
C:  
D:  
D:  NA 
We present equilibrium configurations from scheme B, C and D with and in Table 7, where values of both and are considered, as the latter corresponds to the parameter in Fig 3c of Khairy et al Khairy2008a . With resolutions of , all schemes B, C and D resemble closely the results of the reference. As increases, we observe axisymmetric stomatocytes with smaller height, as exemplified by the first two columns of Table 7. This regime corresponds to the first regime discussed above for Fig. 13, where grows as increases; As (and ) further increases, it enters to the second regime, where stays almost plateau. In this regime, we observe equilibrium shapes as biconcaveoblate and biconcaveelliptocyte, which are exemplified by the third and fourth columns of Table 7. In the last regime, stays flat and we observe prolatedumbbells as shown for example in the last column of Table 7. These configurations of agree with the work of Kairy et al Khairy2008a . However, with the anisotropy of elliptocytes with tends to disappear and the final configurations are less elongated. The results on elliptocyte with . We note that the difference between computed energies for and is very small as shown in Fig. 13. The results from the three different schemes are consistent with each other, which leads us to suspect that the number of terms in the spherical harmonics of the reference could be insufficient to resolve this subtle difference.
4.4 A dynamic simulation of of areadifferenceelasticity model
Finally, we present dynamic simulations of lipid bilayer membranes using the ADE model discretized by scheme B. We consider the same parameters as in section 4.3, that correspond to mechanics activities of a RBC. We start with a stomatoctye shape, which takes the final shape of the previous minimization procedure (see results of or on Table 7). We set the reference areadifference as or .
We report the dynamic trajectories of local energy, nonlocal energy, areadifference and corresponding shapes on Fig. 14. The inplane viscous damping force is orthogonal to the normal direction of the membrane surface and has negligible effects on the dynamic trajectories. However, a larger viscous force is needed to stabilize simulations with a higher resolution and therefore demands a smaller time step. The dynamic simulation arrives at the prolatedumbbell shape, as shown on Fig. 14 in agreement with results on Fig. 12(a) and 12(b) and on Table 7, respectively, which are obtained from independent minimization procedures. From Fig. 13(a), we observe that the local energy does not change significantly over time, whereas the nonlocal energy decreases around ten fold. This confirms that the ADE is primarily responsible for the transformation of shapes. Comparing Fig. 13(b) with 12(b), we notice that the accessible states of are more abundant and even continuous in the dynamic simulation. Accordingly, there is a rich range of shapes including triangle oblates, curvedprolatelike shapes, and prolatedumbbells, as represented on Fig. 13(b). These intermediate configurations are dynamic (nonequilibrium) and do not have a correspondence for the intermediate values on Table 7. Further comparing with the phase diagram of the BC model on Fig. 11, we also do not find corresponding equilibrium shapes for the same values of reduced volume and . The small discrepancy between the energy trajectories of and are not surprising, as they start from slightly different values.
5 Summary
We have presented a comparative study of four bending models for lipid bilayer membranes and their respective discretization on triangulated meshes. The physical models are the minimal, spontaneous curvature (SC), bilayercouple (BC) and areadifference elasticity (ADE) models and the four schemes are termed as scheme A, B, C, and D respectively.
We find that the total energies computed by scheme A on a prescribed sphere and biconcave oblate differ significantly from the analytical values. The energy computed by scheme A cannot be tuned to match the Helfrich energy for an arbitrary shape. For reduced volume , scheme A is able to generate prolatedumbbells at equilibrium, although the corresponding energies are inaccurate. However, for , scheme A does not result in biconcave oblates as equilibrium shapes and it does not sustain the shape if it is given a biconcave oblate as the initial shape. Below , scheme A has numerical artifacts of budding transitions and vesiculations. Finally, due to the lack of direct definition of mean curvature and local area, scheme A can not be directly extended to discretize the SC, BC, or ADE models.
The energies computed by scheme B, C and D on prescribed sphere and biconcave oblate match well the analytical values. However, the accuracy of the three schemes has a varying dependence on their resolution. On the phase diagram of the minimal model, for , all three schemes produce prolate dumbbells and biconcave oblates as the equilibrium configurations. They also resolve well the critical energy values at and and , with scheme B at but scheme C and D at . However, for , Scheme B produces stomatocytes as equilibrium configurations accurately with , where the inner sphere separates clearly from the outer sphere. Scheme C generates acceptable equilibrium configurations, but cannot maintain a clear physical separation between the inner and outer spheres. For the most challenging case of , scheme C is in trouble to generate a physically correct configuration. A higher resolution of does not help scheme C significantly on the overall equilibrium shape. Scheme D has stability issues in the regime of stomatocytes that may be attributed to its lack of explicit energy and momentum conservation.
We have also computed two slices ( and ) of the phase diagram for the SC model. These are challenging benchmarks, as includes the budding transitions at equilibrium, while include both buddings and vesiculations. We observe that schemes B, C and D all capture accurately the budding process, but scheme C and D have difficulties to generate vesiculations. Furthermore, with the same resolution of scheme B computes the energy more accurately than scheme C and D. Moreover, the energies computed by scheme C at budding transitions are completely off the reference values due to the poor representation of the configuration necks.
As a special case of the ADE model, we considered the BC model, where a constraint is imposed on the area difference between the two layers of lipid. We focus on the results of scheme B and reproduce the phase diagram for different resolutions. We find that results of scheme B are comparable to the reference on the equilibrium shapes and energy values, which were calculated by Surface Evovler. We have also considered the complete energy functional of the ADE model. We reproduce the sequence of stomatocytes, biconcave oblates, biconcave elliptocytes and prolate dumbbells by all three schemes, except for the extreme case of stomatocyte where scheme D runs into instability again. As the last demonstration, we consider a dynamic trajectory of the ADE model from stomatocyte to prolatedumbbell. The values of local energy, nonlocal energy, and areadifference along the dynamic trajectory may also serve as reference for other researchers.
Our results indicate that it is a challenging task to compute accurately the bending energy and in particular forces on a triangulated mesh. However, due to its generality it is still a very promising path to pursue. At the same time we note the level set formulations introduced in Maitre2009 which incorporates some of the energy models presented in this work. We believe that this is a very promising direction to account for large membrane deformations and topological changes. Future work could be directed in comparing level set and triangulated surface representations.
In summary, the relatively simple formulation of scheme B, which bypasses definitions of LaplaceBeltrami operator and unit normal vectors, and instead defines the integral of mean curvature with respect to area, exhibits excellent accuracy, robustness, stability. Moreover as the the implementation of scheme B is also almost as simple as that of scheme A we suggest that it deserves to be more broadly investigated as an alternative for vesicle and RBC simulations. While scheme B is readily applicable to the minimal model, in this work we promote and extend its potentials to further resolve SC, BC and ADE models. With ADE model as the most accurate model for lipid bilayer membrane, we envision many numerical applications on vesicles and RBCs in a dynamic context.
Appendix A Force from variational calculus of energy
In differential geometry, the Cartesian coordinates of a point on a surface may be expressed in terms of two independent parameters and as . For further derivations, we introduce some definitions DoCarmo2016 ; Tu2018
(30) 
where free index or is either or , and summation convention applies to repeated index. Index after “,” denotes the partial derivative with respect to .
is the covariant metric tensor and
is its determinant. is the contravariant metric tensor. An infinitesimal area is then(31) 
The parameterization of and are chosen in such a way the the normal direction
(32) 
points inwards of a closed surface. Therefore, is positive for a sphere.
For a small and continuous perturbation along the normal direction , the new position of the point on the parametric surface is given as
(33) 
Thereafter, the variations for and are given as
(34)  
(35) 
where is the Christoffel symbol and is the LaplaceBeltrami operator on the surface. We note a different sign convention used Tu2018 .
To have an explicit expression of , we need variational expressions for each moment. The variation of the zero moment simply reads as
(36) 
However, to calculate variations of the first and second moments we need some fundamental equalities. For any function we have
(37) 
which can be readily proven via integration by parts. We further notice
(38)  
(39)  
(40) 
Therefore,
Comments
There are no comments yet.