Fluid membrane vesicles in confinement

We numerically study the morphology of fluid membrane vesicles with the prescribed volume and surface area in confinement. For spherical confinement, we observe axisymmetric invaginations that transform into ellipsoidal invaginations when the area of the vesicle is increased, followed by a transition into stomatocyte-like shapes. We provide a detailed analysis of the axisymmetric shapes and investigate the effect of the spontaneous curvature of the membrane as a possible mechanism for shape regulation. We show that the observed morphologies are stable under small geometric deformations of the confinement. The results could help us to understand the role of mechanics in the complex folding patterns of biological membranes.


Introduction
The folding of tissues and membranes plays a dominant role in the interplay between structure and function in living matter. Countless examples can be given both on the macroscopic (e.g. the gut [1] or the kidney [2]) and the microscopic level (e.g. the endoplasmic reticulum [3] or the cytoplasmic membrane [4]). There are various reasons for folding to occur such as differential growth in tissues [5][6][7] or the cooperative action of proteins in lipid membranes [8].
In particular, the inner membrane of mitochondria is of special interest, since its alteration and remodeling is closely linked to numerous neurodegenerative diseases [9]. It is found that this membrane consists of three parts: (i) the inner boundary membrane, which is in close proximity to the surrounding outer membrane, (ii) the invaginations of the inner membrane inside the interior of the mitochondrion, also called cristae, and (iii) the cristae junctions, which are narrow tubular necks connecting (i) and (ii) [10,11]. Experimental studies of the inner membrane suggest that the structure of the membrane results from the cooperation of different effects such as surface growth, lipid composition (such as the concentration of cardiolipin in the membrane [12]) and protein interactions [13,14]. To disentangle the different effects, one is urged to search for simple model systems which allow one to understand the influence of each contribution in a clearly defined way.
To reproduce cristae-like deformations of lipid bilayers experimentally, Khalifat et al [15] and Fournier et al [16] have changed the local pH gradient of giant unilamellar vesicles. Ponnuswamy et al [17] and Ghochani et al [18] proposed free energy models for given geometries of a crista junction to predict its size. However, the simplest and most intuitive way of understanding the generation of membrane invaginations is to consider the buckling of a closed growing membrane in a constrained volume. Applying this model in a finite element framework, Renken et al [19] concluded that minimizing the bending energy is sufficient to provoke invaginations without giving a detailed analysis of the generated shapes. In [20], we studied these shapes with the help of finite element simulations and provided a morphological phase diagram for cristae-like membrane invaginations. 3 In this paper, we explore the results of [20] in more detail and compare them to experiments. A set of differential equations describing axisymmetric invaginations of the membrane is solved numerically with a shooting method. The solutions coincide with our previous findings from finite element simulations and corroborate the emergence of the axisymmetric morphological phase. The length scales of the axisymmetric invaginations can be controlled by incorporating a non-zero spontaneous curvature into the model. Extending our description to ellipsoidal confinements results in the same morphological phases up to small local modifications. This shows that our findings are stable towards small geometric modifications such as those found in nature.

A minimal model
Lipid bilayer membranes are composed of two layers of polar lipids held together by their amphiphilic properties. In larger membrane structures such as cell membranes or mitochondria, the lateral size exceeds the bilayer thickness by several orders of magnitude, allowing a representation of the bilayer membrane in terms of a surface embedded in three-dimensional space. On this level of description, the mechanics of the membrane is governed by three different constituents: firstly, bending energy, which microscopically is due to the imbalance of lipids in the inner and the outer layer when the membrane is curved; secondly, surface energy from extension or compression of the membrane. Due to the free movement of the lipids in the membrane, the surface energy is free of shear and can be understood as a mechanism for the conservation of area. Finally, an osmotic pressure arises from concentration differences of solutes in the enclosed and outer volume. The energy scales of the surface and osmotic parts are considerably larger than the bending energy [21]. A useful and common simplification therefore is to treat them as constraints on the total surface areaĀ and enclosed volumeV . The classical model for describing the bending contribution was proposed independently by Canham [22], Helfrich [23] and Evans [24], and relates the bending energy E b to a second order expansion of curvatures: where H is the mean and K the Gaussian curvature, and the material constants κ and κ G , respectively, are the bending rigidity and the Gaussian rigidity, or saddle-splay modulus.
The spontaneous curvature C 0 represents an intrinsic preferred mean curvature value of the membrane due to a possible asymmetry of lipid density in the membrane layers 5 . For a given topology of genus g, the second term is a topological invariant of the value 4π(1 − g) according to the Gauss-Bonnet theorem. Since its variation is zero, the term does not appear in the equilibrium equations and is neglected in the following.
Our aim is to study the shapes of fluid vesicles inside a rigid cavity as a function of their target volumeV and surfaceĀ. The size of the cavity introduces an additional length scale to the problem, such that the setup is more concisely described in terms of the scaled volume v and area a, v =V /V 0 and a =Ā/A 0 , where V 0 and A 0 are the volume and the area of the cavity. Clearly, prescribing a above 1 means that spherically symmetric shapes no longer fit in the cavity, and the emergence of non-trivial solutions with lower symmetry is expected. For even higher values of a, the deformed vesicle will eventually touch itself, such that self-contact, i.e. the impenetrability of the membrane with itself, has to be included in the description. An analytical treatment of this problem is, however, in general not accessible. We rely in the following on two different numerical methods: for symmetric shapes without self-contact, we recast the energy equation (1) into a set of differential equations, essentially resulting in a classical 'shooting problem' (see appendix A). Asymmetric deformations or deformations with self-contact require a full numerical treatment of the problem, pursued in this work using the finite element technique proposed in [28] (see appendix B for details). In the finite element model, the entire surface is discretized into a mesh of triangles and minimizers of the discretized version of equation (1) are searched for. This technique guarantees that one minimum is found but does not exclude that others exist for the same value of area and volume. The enclosing cavity is modeled as a rigid container using a quadratic repulsive force penalizing the movement of the membrane beyond the boundary of the cavity. It is important to note that this approximate enforcement of the cavity constraint allows the membrane to slightly penetrate into the cavity. Self-contact is modeled with a similar approach employing a quadratic repulsive potential between intersecting polygons. The thickness of the bilayer is not taken into account by this contact formulation.

Previous results
The confinement of a fluid vesicle with zero spontaneous curvature in a spherical container has been studied by us only recently with the finite element method [20] 6 . A phase diagram was derived in the range of a = 1.1, . . . , 1.6 and v = 0.6, . . . , 0.9. For a fixed enclosed volume, several shape transitions were observed when the area was increased. As soon as one sets the area of the inner membrane to a value larger than the area of the container, the membrane has to fold into the interior and forms an invagination (see figure 1 for v = 0.8 and increasing values of a). The shapes for moderate values of surface growth were found to be axisymmetric within the numerical errors. For higher values the membrane has to break symmetry and forms an ellipsoidal invagination with a slit-like neck at which the membrane touches itself (see figure 1(b)). If the area is increased even further, the invagination folds up and a secondary invagination forms inside the first one which resembles a stomatocyte connected to the outer part via the neck (see figure 1(c)). The system still possesses one symmetry plane perpendicular to the slit-like neck. A simple theoretical model allowed us to explain these transitions qualitatively [20]: by neglecting the neck and assuming that the confined membrane consists of a spherical part (= the part in contact with the container) and a closed inner vesicle with reduced volume ν (= invagination) the transition lines were found to be defined approximately by the equation with 0.59 ν 0.65 for the ellipsoidal invagination with self-contact at the neck and ν 0.59 for the stomatocyte-like invagination [32].

Axisymmetric solutions
The axisymmetric shapes that result from the finite element simulations can be checked using a numerical shooting method (see appendix A for details). This method confirms that the obtained shapes are truly axisymmetric. For (a, v) = (1.2, 0.8), for instance, the profile found with the shooting method (see figure A.1) coincides with the shape of the finite element method in figure 1(a) within the numerical errors. At the section which is in contact with the container, the mean curvature of the inner membrane is negative and constant, H = −1 (see blue solid line in figure 2(a)). To detach from the container, the membrane has to attain a smaller value H < −1 which changes its sign at the neck and becomes positive and greater than one at the invagination. For the idealized problem, delta forces at the contact point are necessary to equilibrate the system. In the finite element simulations these forces are smeared out slightly since a small part of the membrane penetrates into the container (see the red dotted line in figure 2(a)). For a fixed area a = 1.2 but slightly smaller volume v = 0.7 the shape obtained with the finite element method is still axisymmetric but instead of two, one finds three peaks in the container force (see figure 2(b)). The third peak indicates that the membrane leaves the container at the bottom as well. Consequently, we do not find any shapes with the shooting method for these values because the underlying model does not allow the membrane to penetrate into the container (see again appendix A). One should, however, note that the container forces for (a, v) = (1.2, 0.7) are rather small compared to the other forces in the system. The mean curvature of the membrane in contact deviates only slightly from −1. This changes when the volume is decreased even further. For lower v the invagination touches the rest of the membrane at the bottom of the container and both, invagination and the membrane in contact with the container, are pushed out of the cavity.

Symmetry breaking
If one fixes the volume and increases the area, a second invagination forms inside the first one as we have seen in section 3.1. If we were to increase the area further, we would expect a third invagination forming in the interior, then a fourth, a fifth and so forth. A similar structure is observed in experiments on mitochondria [33]. In these experiments the inner membrane protein mitofilin was suppressed. This protein plays an important role in the correct folding of the internal structure of the mitochondrion. When it is suppressed, onion-like structures of the inner membrane form, which resemble cuts of our numerical shapes rather closely (see figure 3). The topology of the experimental shapes is different from the spherical one of our model but the basic structure is the same. This is not surprising, since in vivo self-contacts and curvature sensing proteins are known to induce topology changes. Although not included in the present model, we can speculate that these mechanisms would lead to fusion and consequently topology changes at the neck of figure 3 where the curvature is highest and self-contact occurs. In fact in all of our shapes with high area, self-contacts arise eventually which could induce topology changes in vivo: the ellipsoidal invagination in figure 1(b) could, for instance, fuse at the slit-like neck as well and form two tubular necks that connect the invagination with the rest of the membrane.

Confined vesicles with spontaneous curvature
For fusion or fission to occur in cells it is not even necessary that the membrane touches itself at the beginning of the process. One classical example of vesicle formation and budding is the clathrin-mediated endocytosis in which the budding is initiated by BAR domain proteins, stabilized by the polyhedral clathrin scaffold coating the membrane [34], and scised by the recruitment of dynamin to the neck [35]. A similar budding/scission mechanism can be found in the endolysosomal pathway. However, in this pathway, ESCRT machinery operates in an opposite fashion, forming and cleaving the endosomal buds from the inside of the neck [36]. The remodeling of the structure in mitochondria depends on proteins as well [14,33].
Regardless of the exact mechanism, from a theoretical point of view one can try to understand all these processes as a curvature regulation of the membrane with the help of lipids, proteins, or a direct action of the cytoskeleton. Such an induced curvature in the membrane can be accommodated in our model via the spontaneous curvature term C 0 of the energy (1). Until this point we have considered confined membranes with zero spontaneous curvature, i.e. they prefer to be flat if there are no further constraints. By assuming a constant non-vanishing value of C 0 we can incorporate effects from membrane-deforming molecules globally. Any spontaneous curvature value different from zero changes the ground state of the membrane and sets a value that the mean curvature wants to adopt. However, due to the constraints on surface area, enclosed volume and confinement the exact shape cannot be predicted ad hoc.
Consider, for example, figure 4 where we depict the equilibrium surfaces for different spontaneous curvature values C 0 = 0, 0.1 and 1 at fixed scaled area and volume. One observes that the axisymmetry of the system is preserved for increasing spontaneous curvature even though the shape itself changes: the size of the neck shrinks gradually with increasing C 0 ; the invagination moves up and becomes more spherical. At the same time the membrane has to stay in contact with the container which fixes its mean curvature to H = −1 there.
The shapes for C 0 = 0 and C 0 = 0.1 were found with the finite element simulations and could be confirmed with the shooting method (see again appendices A and B). If C 0 is set to higher values such as C 0 = 1, the shooting method does not give a solution any more. This is understandable when looking at the corresponding simulations: similar to the example of figure 2(b) in section 3.2 the simulated membrane leaves the container at the bottom, a situation that has been ruled out in the shooting method. For even higher values of C 0 no equilibrium shape can be found at all. This indicates that a pinch-off of the neck might take place in agreement with what can happen in vivo. Even a global change of the spontaneous curvature, such as the one considered here, is thus an effective way of controlling length scales of an invagination such as the diameter of necks, a fact which has already been observed before for unconfined membranes [26,[37][38][39]. As a next step one might want to study how local changes of the spontaneous curvature influence the system. This would go beyond the scope of this paper but can in principle be achieved by incorporating the composition dynamics of a multi-phase fluid membrane in the finite element model [40].

Influence of the container geometry
Membrane invaginations in vivo such as the cristae of mitochondria usually occur in a non-spherical confinement. The shape of the outer membrane spans a wide range of geometries from small spheres to large tubular networks [41].
For a model membrane in spherical confinement with no spontaneous curvature, the observed invagination shapes can be explained by a simple reduced volume effect if the neck of the invagination is neglected (see equation (3) and [20]). One can incorporate a modified geometry of the container into the approximative model: we again consider a system where the invagination is a closed vesicle with no connection to the part of the inner membrane in contact with the container. We can then express the scaled surface a of the inner membrane as the sum of the areas of the invagination, A i , and the boundary membrane, A o , divided by the container area 7 . If we assume that the boundary part is a closed vesicle that coincides with the container, both have the same area and we obtain a = A o +A i A o . Similarly, the scaled volume is 9 given Using the definition of the reduced volume, ν = 6 √ π V A 3/2 , to combine a and v, we find the following relation: where the parameters ν i and ν o denote the reduced volume of the invagination and of the boundary membrane (i.e. the container), respectively. Note that this equation reduces to equation (3) for a spherical container since ν o = 1 in this case. For an arbitrary confinement, equation (4) predicts a shift of the transition lines in the phase diagram of the spherical system, which depends on the reduced volume of the container. One can test this prediction by applying any kind of deformation to the container in the simulations. We choose a prolate ellipsoid because, firstly, it is a comparatively simple deformation of the sphere with a tunable reduced volume and, secondly, it is a good approximation for tubular organelles. The container surface is now given by the equation where we will tune the axis parameter ξ 1, while setting the two others equal to one. Increasing the value of ξ increases the ellipticity ε = 1 − ξ −2 of the container as well. The area and volume of such a prolate spheroid are given as [42] A respectively. The reduced volume follows as This function decreases only slowly if ξ is increased. For small eccentricities, say ξ = 1.5 for example, the reduced volume is ν o = 0.96, which is still close to one. One could expect that the resulting shapes do not change a lot compared to the case of the spherical confinement. And this is indeed what happens: in figure 5, we show the results of the finite element simulations of a membrane inside a prolate container with ξ = 1.5. The scaled volume is fixed at v = 0.8 and a is set to 1.2, 1.5 and 1.6, respectively. Consistent with the argument above, we observe shapes that are similar to those inside a spherical container (see figure 1). The invaginations pop up on almost planar regions of the ellipsoid. This is not only valid for the membranes presented in figure 5 but has been observed for other values of the parameters as well. We also note that the slit-like necks of non-axisymmetric invaginations tend to align themselves along the major axis of the container. This is not totally unexpected since the reverse situation would require to pack one ellipsoid into another of comparable size with their major axes at a right angle. Such a conformation requires a large amount of surface intersection which is penalized by the contact forces and, consequently, not observed in the simulations.
For low values of a the membrane detaches from the container at the extremities of the major axis (see figure 6). The reason for this is relatively simple: the membrane has to form an invagination to fit into the container while accommodating the volume constraint v = 0.8 at the same time. This already eats up a non-negligible amount of the membrane area. The remaining part is not sufficient to fill out the container, at least for small values of a such as a = 1.1, 1.   If we increase the length ξ of the major axis, our simple theoretical model (4) predicts a shift of the transition lines (see above and figure 7). Let us consider the transition from an axisymmetric to an ellipsoidal invagination for the case of fixed v = 0.8 where ν i = 0.65. The corresponding scaled area for the spherical confinement is given by a = 1.46. For ξ = 1.5 this transition point shifts only slightly to a = 1.44 as expected. But even if we increase ξ to 3, we obtain a = 1.39, i.e. the model predicts only a small change of the transition point if the scaled enclosed volume v is not too small. In fact, equation (4) implies that a remarkable shift of the transition lines can only be found for very high ellipticities. However, for high ellipticities the contacts of the membrane with itself at the bottom of the container will become very important. These contacts will modify the resulting shapes and render our simple approximative model invalid. Further research on the subject is thus necessary to incorporate this effect as well.

Conclusions
In this paper, we studied the role of mechanics in the morphogenesis of fluid membrane invaginations. We approached the problem by simplifying the external constraints as a closed rigid confinement. In the case of a spherical container we observed axisymmetric invaginations for moderate values of surface growth. We confirmed this result using the well-known shooting method. Higher values of surface growth reduce the symmetry of the system resulting in ellipsoidal invaginations. If we increase the surface area further, stomatocyte-like invaginations, which consist of a secondary invagination inside the first one, emerge as equilibrium solutions. Random cuts through such stomatocyte-like invaginations bear similarities to electron tomography images of mitochondria lacking the membrane protein mitofilin, which is believed to be responsible for the organization and stabilization of cristae. In order to relate our work further to experimental observations, we also studied the effects of a non-zero spontaneous curvature of the membrane. We showed that already a small change of the spontaneous curvature modulates the length scales of the invagination and its connecting neck. Moreover, we found that these morphologies are stable with respect to small geometric perturbations of the confinement as was shown for ellipsoidal deformations of the spherical container.
Despite its simplicity, the model is capable of producing a basic set of invaginations which appear to be the ground states of membrane folding under external confinement. The finer ultrastructural details in vivo could be due to the local control of the spontaneous curvature and dynamical effects arising from the hydrodynamical interactions of the fluid membrane with the surrounding viscous environment. The former can be incorporated into the model by adding a Ginzburg-Landau-type interface between two kinds of lipid phases [40,43]. The latter requires the simulation of viscous fluids which can be handled either by a mixed finite element formulation [44,45] or by coupling the system to a lattice Boltzmann framework [46]. The presented model does not take into account topology changes induced during the folding process. Our approach, however, could be easily extended to study the result of such a transition once it occurred by considering a surface with higher genus. Furthermore, it is an open question if the stabilizing effect of the proteins involved in the cristae formation in vivo could be macroscopically accounted for by adding an attractive part to the membrane self-contact forces.
Following these works we use a Hamiltonian formulation and divide the shape of the membrane into two parts: a spherical cap of unit size in contact with the outer container and an upper free section which forms the invagination. The shape of the spherical cap and the container coincide and are given by the height function z(x) = ± √ 1 − x 2 (see figure A.1). To describe the shape of the free profile, we will use the angle-arc length parameterization, ψ(s), where ψ is the angle between the tangent vector and the horizontal x-axis and s is the arc length. The mean curvature in this parameterization is then given by H = − 1 2 (ψ + sin ψ x ). At the container, H = −1.
The scaled energy of the free part can be written asẼ := E/(πκ) = s s dsL with are the scaled tension and pressure, respectively. The Lagrange multiplier function λ x fixes the geometrical constraint thatẋ = cos ψ everywhere along the profile. The conjugate momenta are We switch to a Hamiltonian formulation with the scaled Hamiltoniañ from which we obtain the Hamilton equationṡ To find the shape of the free membrane section, these equations have to be solved numerically with the appropriate boundary conditions: Equation (A.9) takes into account that the membrane must not have kinks: at the contact point s = s the membrane detaches from the container and ψ equals the angle α (see again figure A.1); at the z-axis, s =s, and the profile is horizontal due to the symmetry. The contact curvature condition (A.10) results from an energy balance at the contact point [50,63], whereas the vanishing of the Hamiltonian (A.11) is due to the fact that we do not fix the total arc length s − s.
The Hamilton equations are solved using a standard shooting method [64]: for a fixedσ andP and a trial angle α the equations can be integrated with a fourth-order Runge-Kutta method. For a given α ∈ {0, π}, the values of ψ, x, p ψ and p x can be calculated at s = s with the help of the boundary conditions. The integration is stopped as soon as the vertical axis is reached. Every α which results in ψ(s) = π corresponds to the profile of a closed axisymmetric membrane in a spherical cavity.
Note that this method does not allow us to determine the shape of the membrane directly as a function of the membrane area and enclosed volume. These can only be calculated as soon as a profile is found: the area of the membrane in contact with the container is given by 2π s 0 ds sin ψ = 2π(1 − cos α) and the volume of the corresponding lower spherical cap yields π s 0 ds x 2 sin ψ = π 3 (2 − 3 cos α + cos 3 α). Integrating up the equationsȦ = 2π x anḋ V = π x 2 sin ψ from s tos gives the area and volume of the free part, respectively. Sweeping through the parameter space (σ ,P) thus allows us to find confined vesicles with variable area and volume.

Appendix B. Finite element simulations
Our aim is to find equilibrium solutions of the energy functional (1) augmented with container and self-contact constraints. One notes that the curvature integral in equation (1)  We approximate the surface by a set of E triangles with N nodes. Using the trial functions N a , the position x is then interpolated by the weighted sum x h where x a is the position of node a. For N a we employed the Loop subdivision trial functions (see, e.g., [65] for their precise form).
To make equation (1) more amenable to the finite element method, we enforce the volume and surface constraints using the penalty method, which offers improved convergence properties and easier implementation than Lagrangian multipliers. In this method, the constraints are only where µ A and µ V are constants chosen large enough to approximatively fulfill the constraint conditions. The internal energy as formulated above is not yet suitable for the finite element method: in contrast to elastic shell theory, material points in the Canham-Helfrich-Evans model can freely flow on the surface as long as the shape of the surface is not altered, corresponding to the in-plane rearrangement of lipids. For a discretization into a polygonal mesh, this degree of freedom leads to arbitrary tangential movements of mesh nodes resulting in severe mesh distortion. Two methods exists to avoid this problem: in the first, the area constraint is recast into a local form (local area constraint method). The second method enforces the area constraint globally (global area constraint method) using equation (B.3), but uses a mesh-regularization potential to stabilize local mesh quality. We employed both methods to validate our results.
In the local area constraint method, equation (B.3) is recast into the equivalent formulation in which √ a is the surface Jacobian such that A = √ a ds 1 ds 2 . While physically identical, the area constraint here is enforced locally, i.e. at the element level.
Alternatively, the global area constraint method introduces an artificial spring force, which penalizes spurious tangential motion of nodes in an iterative process, such that the equilibrium solution is not affected [40]. At the start of a simulation, we calculated the reference length l 0 = l 0 1 , . . . , l 0 M of all edges. Deviation from the reference length is penalized with a linear spring force of the form where the sum runs over all edges i of node a, l i is the current edge length and n i is the unit vector along edge i. The reference lengths l 0 are reset every 1000 iterations to the current edge lengths. This prevents the regularization forces from suppressing the fluid-like nature of the membrane and, consequently, does not affect the equilibrium solution. In addition, severely distorted elements are stabilized by an additional regularization force which favors equilaterality of the triangle. For a node a of such a triangle, the corresponding force is calculated according to where the 'equilateral reference length'l 0 i is given byl 0 i = 2 A √ 3 , i.e. the edge length of an equilateral triangle with the same area as the distorted one. For both forces, f R 1 a and f R 2 a , the constant k is set to 10.
For either area constraint method, the finite element procedure proceeds in the same way: equation (B.3) or (B.4) is expressed as a function of the nodal position coefficients x a . Taking the gradient with respect to x a then yields the corresponding nodal forces f M a (we denote by M the contribution from the membrane forces). In the so-called stress resultant formulation, f M a can be concisely written as where the stress resultant n α and the moment resultant m α measure the in-and out-of-plane forces due to bending, whereasn α and f account for the constraints. They are given by The coefficient A depends on the type of the surface area constraint chosen; it is A = √ a − √ā in the case of the local area constraint method and A = A −Ā for the global area constraint method. In the latter case, the regularization forces f R 1 a , f R 2 a as described above are added to the resulting force vector as well. Values of µ A = 10 4 and µ V = 5 × 10 4 were used for the global constraint method, whereas the local constraint method had µ A = 10 5 and µ V = 5 × 10 4 .
Contact between the membrane and the cavity is modeled by a quadratic force law applied to those nodes that are penetrating the cavity. From the current position x a of a node, we first calculate the penetration depth d a and then apply the force in the direction n normal to the container surface.
f C 1 a = k 1 d 2 a n, (B.12) where a stiffness constant k 1 = 1500 was chosen. Self-contact is implemented according to the algorithm presented in [67]. In contrast to typical contact handling methods, this algorithm does not apply contact forces based on the distance between surface segments before they penetrate, but calculates the necessary forces to disentangle intersections once they occurred: whenever parts of the surface penetrate, the contour line of the intersecting polygons is first constructed. In a second step, for each node a of a triangle involved in an intersection, the gradient G of the contour length w.r.t. the position of node a is calculated. Based on this gradient, a contact force f C 2 a = k 2 G is applied to node a. In our simulations, we chose a stiffness constant of k 2 = 9, which was sufficiently large to quickly resolve all intersections during a simulation within a few timesteps.
The resulting force f a at node a consists of the contributions from the Helfrich, constraints and contact forces, and the equilibrium is characterized by the balance of forces, i.e.
where f R a = f R 1 a + f R 2 a is only non-zero if the global area constraint method is used. In contrast to Klug and coworkers, we do not employ a nonlinear gradient descent method to find equilibrium solutions. Instead, we integrate in time the nodal forces f a until a balance of forces is reached according to Newton's equation of motion: where M is the (diagonal) mass matrix,ẍ = (ẍ 1 , . . . ,ẍ N ) is the generalized acceleration vector and C = M · diag(c, . . . , c) is the mass-dependent viscous damping matrix. We chose c in the range of c = [50,200] depending on the amount of internal energy. Numerical integration in time was performed using a standard Newmark method, see, e.g., [68]. Our simulations started from spherical meshes with N = 2562 nodes. For shapes with large surface growth, we repeatedly optimized the mesh during the simulation using the free software package ReMESH [69]. The remeshing procedure redistributes a given number of nodes on the initial surface for improved accuracy and mesh quality. We employed up to three remeshing steps within one simulation, resulting in either N = 3000 or N = 4000 mesh nodes.