The influence of higher order geometric terms on the asymmetry and dynamics of membranes

Abstract

We consider membranes as fluid deformable surfaces and allow for higher order geometric terms in the bending energy related to the Gaussian curvature squared and the mean curvature minus the spontaneous curvature to the fourth power. The evolution equations are derived and numerically solved using surface finite elements. The two higher order geometric terms have different effects. While the Gaussian curvature squared term has a tendency to stabilize tubes and enhance the evolution towards equilibrium shapes, thereby facilitating rapid shape changes, the mean curvature minus the spontaneous curvature to the fourth power destabilizes tubes and leads to qualitatively different equilibrium shapes but also enhances the evolution. This is demonstrated in axisymmetric settings and fully three-dimensional simulations. We therefore postulate that not only surface viscosity but also higher order geometric terms in the bending energy contribute to rapid shape changes which are relevant for morphological changes of cells.

---

1 Introduction

Membranes are ubiquitous and essential in biology; they compartmentalize biomaterials, separate the cell from its exterior and organelles from the cytoplasm, dynamically remodel and change conformation. Geometric properties of the membrane have been identified as key players for such processes.^1–3 As the typical thickness of a membrane is orders of magnitude smaller than its lateral extension, treating the membrane as a two-dimensional surface embedded in a three-dimensional space is a reasonable approximation. This separation of length scales allows for a mesoscopic modeling where details related to membrane molecular structure are considered in effective material parameters and geometric quantities and led to the success of the classical Canham–Helfrich model,^4,5 which builds on a bending energy with mean curvature [script letter H], Gaussian curvature [scr K, script letter K], bending rigidity parameters k_2,0 and k_2,1, a spontaneous curvature [script letter H]₀ and additional (local or global) area and volume constraints. For definitions see Section 2.1. Assuming constant values of k_2,1, the second term reduces to a topological measure and thus a constant as long as the topology does not change. We will therefore neglect this term.

Equilibrium shapes, resulting from minimizing the bending energy [scr F, script letter F]_BE, have been extensively studied, see ref. 6 for a review. For [script letter H]₀ = 0 and wide ranges of the reduced volume, which is the ratio of the volume and the volume of an equivalent sphere with the same area, they are dominated by prolate and oblate shapes.⁷ The spontaneous curvature [script letter H]₀, which accounts for the asymmetry of the membrane, is able to modify these shapes.⁷ A wealth of studies also aims to produce tubes as minimizing shapes. Tubes are ubiquitous in membranes and play crucial roles in trafficking, ion transport, and cellular motility. For idealized situations this is rather simple as [scr F, script letter F]_BE = 0 if [script letter H] = [script letter H]₀, which is achieved for an infinite tube of radius r and [script letter H]₀ = −1/r. However, this solution is not unique, a sphere of radius 2r also leads to [scr F, script letter F]_BE = 0. More realistic situations with finite volume and area require further considerations, e.g., introducing a spontaneous curvature deviator.^8–10 Also, various ideas have been proposed to consider higher order geometric terms in the bending energy [scr F, script letter F]_BE to enforce the stability of tubes.^2,11–14 For example, fourth order terms proportional to [scr K, script letter K]² seem plausible as [scr K, script letter K] = 0 for tubes.

Most of these studies only focus on equilibrium shapes, comparing the bending energy [scr F, script letter F]_BE of different configurations and addressing their stability. But the dynamic evolution of the membrane and associated shape changes are also of interest, e.g. in the context of the formation of bulges with pinch-offs, which can be associated with endocytosis and exocytosis^15–17 or shape changes associated with cell motility.^18–20 All these processes require additional phenomena, which are not considered in the Canham–Helfrich model. But there is growing interest in the role of membranes in these processes. One striking example where membranes at least participate are frequently forming and retracting filopodia.²¹ This process requires rapid shape changes in cells and according to ref. 22 the membrane and the underlying cortex act as an integrated system to globally coordinate such changes in cell shape. To facilitate these rapid morphological changes, cells maintain an excess of membrane that is organized in membrane reservoirs and is available to the cell on the order of seconds.²³ To understand such processes thus not only requires unveiling the secretes of equilibrium shapes but also to consider the flow of membranes which facilitates the rapid shape changes. To model such processes is a different story and even if various models for the cellular cortex exist^24–27 and also first attempts to couple them with membrane models to form the mentioned integrated system have been proposed,^22,28 we here refrain from such couplings and only aim first for a minimal model of the membrane alone, which facilitates rapid shape changes.

Surface viscosity has been identified as a key player²⁹ and considering membranes as fluid deformable surfaces^30–33 opened new perspectives on the description of the dynamics of membranes. Fluid deformable surfaces can be viewed as two-dimensional viscous fluids with bending elasticity. Due to this solid–fluid duality, any shape change contributes to tangential flow, and vice versa, any tangential flow on a curved surface induces shape deformations. This tight coupling between shape and flow makes curvature a natural element of the governing equations. As demonstrated by numerical studies of the equations for fluid deformable surfaces, surface hydrodynamics can significantly speed up the evolution^34,35 and can enhance bulging and furrow formation in membranes.¹⁷

Combined with higher order geometric terms in the bending energy, models for fluid deformable surfaces has the potential to enable the rapid shape changes in the release and formation of membrane reservoirs and the formation and retraction of filopodia. We computationally explore the effect of these terms on the equilibrium shapes and the dynamics to reach them and analyze their impact to facilitate rapid shape changes. Any coupling with the cortex,²² interaction with proteins that induce curvature³⁶ or to forces exerted on the membrane^37,38 are not considered. We also neglect any interaction with the surrounding bulk phases. This results from the theoretical interest in exploring the membrane properties without any additional influence and the limit of a large Saffman–Delbrück number. This number describes the relation between the viscosities of the membrane and the typically less viscous bulk fluid and if large allows the decoupling of surface and bulk flows.³⁹

The rest of the paper is structured as follows. In Section 2 we introduce the full model and briefly mention the considered numerical approach. More details on the derivation of the model and on the numerical approach including convergence studies are provided in Appendices (Model derivation, Numerical method and Validation). Computational results are described in Section 3. They contain axisymmetric and full three-dimensional computations addressing the dynamic evolution and the equilibrium shapes. In Section 4 we draw conclusions and mention possible directions to extend the described model to enable simulations of rapid morphological changes of cells.

2 Model

2.1 Notation

The considered mesoscale model requires basic notation from differential geometry and geometric partial differential equations. Besides classical mathematical text books in these fields we refer to ref. 2 for an introduction in the context of membranes. We follow the same notation as in ref. 17 which is here repeated for convenience. We consider a time dependent smooth and oriented surface [scr S, script letter S] = [scr S, script letter S](t) without boundary, embedded in [Doublestruck R]³. The enclosed volume is denoted by Ω = Ω(t). We denote by ν the outward pointing surface normal, the surface projection is P = I − ν⊗ν, with I the identity matrix, the shape operator is [scr B, script letter B] = −∇_Pν, the mean curvature [script letter H] = tr[scr B, script letter B], and the Gaussian curvature . We consider time-dependent Euclidean-based n-tensor fields in Tⁿ[Doublestruck R]³|_{[scr S, script letter S]}. We call T⁰[Doublestruck R]³|_{[scr S, script letter S]} = T⁰[scr S, script letter S] the space of scalar fields, T¹[Doublestruck R]³|_{[scr S, script letter S]} = T[Doublestruck R]³|_{[scr S, script letter S]} the space of vector fields, and T²[Doublestruck R]³|_{[scr S, script letter S]} the space of 2-tensor fields. Important subtensor fields are tangential n-tensor fields in Tⁿ[scr S, script letter S]≤Tⁿ[Doublestruck R]³|_{[scr S, script letter S]}. Let p ∈ T⁰[scr S, script letter S] be a continuously differentiable scalar field, u ∈ T[Doublestruck R]³|_{[scr S, script letter S]} a continuously differentiable [Doublestruck R]³-vector field, and σ ∈ T²[Doublestruck R]³|_{[scr S, script letter S]} a continuously differentiable [Doublestruck R]^3×3-tensor field defined on [scr S, script letter S]. We define the different surface gradients by ∇_P p = P∇p^e, ∇_P u = P∇u^eP and ∇_C σ = ∇σ^eP, where p^e, u^e and σ^e are arbitrary smooth extensions of p, u and σ in the normal direction and ∇ is the gradient of the embedding space [Doublestruck R]³. The corresponding divergence operators for a vector field u and a tensor field σ are div_P u = tr(∇_P u) and div_C (σP) = tr∇_C (σP), where tr is the trace operator. The relations to the covariant derivatives ∇_{[scr S, script letter S]} and the covariant divergence div_{[scr S, script letter S]} on [scr S, script letter S], with Δ_{[scr S, script letter S]} = div_{[scr S, script letter S]} ·∇_{[scr S, script letter S]} the Laplace–Beltrami operator, read ∇_P p = ∇_{[scr S, script letter S]} p and div_P u = div_{[scr S, script letter S]}(Pu) − (u·ν)[script letter H], respectively.

2.2 Governing equations

The material velocity u ∈ T[Doublestruck R]³|_{[scr S, script letter S]} can be decomposed into u = u_Nν + u_T, with u_N = u·ν and u_T = Pu, the normal and the tangential part, respectively. The pressure p ∈ T⁰[scr S, script letter S] serves as Lagrange multiplier for the inextensibility constraint. The governing equations for these unknowns read

div_Pu = 0 (2)

where [∇_wu]_i = (∇_{[scr S, script letter S]}u_i,w), i = 1, 2, 3, with w = u − ∂_tX is the relative velocity and X a parameterisation of [scr S, script letter S], is the rate of deformation tensor, Re > 0 is the Reynolds number, γ ≥ 0 is a friction coefficient, λ ∈ [Doublestruck R] is a Lagrange multiplier to ensure a constant enclosed volume, and b denotes a bending force, defined bywhere [script letter H]₀ is a spontaneous curvature and k_2,0,k_4,0,k_4,2 ∈ [Doublestruck R] are bending rigidity parameters. The system of equations considers a model for fluid deformable surfaces. Such models consist of incompressible surface Navier–Stokes equations with bending forces and a constraint on the enclosed volume. The highly nonlinear model accounts for the tight interplay between surface evolution, surface curvature and surface hydrodynamics and allows for the modelling of membranes with surface viscosity. For derivations of the model (with k_4,0 = k_4,2 = 0) we refer to ref. 17, 31, 34, 40 and 41. They consider different principles and build on a nonlinear Onsager formalism,³¹ thin film limits from three-dimensional models^34,40 and a Lagrange–D’Alembert approach.^17,41 For a comparison of derivations for b = 0 and without the constraint on the enclosed volume we refer to ref. 42–44.

The mentioned solid–fluid duality of these models, which leads to a tight coupling between shape and flow, and lets any shape change contribute to tangential flow, and vice versa, can best be seen by rewriting eqn (1)–(3) as coupled equations for u_N = u·ν and u_T = Pu.^34,42,43 One out of several of the resulting coupling terms results from eqn (2), which leads to div_Pu = u_N[script letter H] + div_Pu_T = 0 and explicitly demonstrates the relation between the mean curvature and the tangential velocity.

Considering the overdamped limit, formally letting γ → ∞, leads to the classical dynamic equations

u = −∇_{[scr S, script letter S]}p − p[script letter H]ν + b − λν (5)

div_Pu = 0 (6)

for an inextensible membrane with constant volume. Using (5) in (6) provides the equation for the Lagrange multiplier for the inextensibility constraint −Δ_{[scr S, script letter S]}p + p[script letter H]² + λ = b·ν[script letter H] and corresponds to previous models, if b only contains second order geometric terms.^45–49 Further relaxing the constraint on inextensibility leads to the classical Canham–Helfrich models with area and volume constraints.⁴⁸

In contrast to previous approaches for fluid deformable surfaces^{17,31,34,35,50} the bending force b also contains higher order geometric terms.

2.3 Bending forces

The bending force eqn (4) results from the bending energywith bending energy density f_BE. Formally we can derive f_BE via Taylor expansion at the spontaneous curvature [script letter H]₀ leading to
in terms of geometric orders n ≤ N ∈ [Doublestruck N] for different bending rigidity parameters k_n,α ∈ [Doublestruck R]. This expression corresponds to the generalized form of the classical Canham–Helfrich energy (n = 2) introduced in ref. 11 and also considered in ref. 2, 12 and 13 if gradient terms are neglected. It can be simplified assuming certain properties of the bending energy: as we are interested in variations of the energy, we can omit the constant contribution k_0,0. We will further not allow for topological changes, thus, we can omit the k_2,1[scr K, script letter K] term, applying Gauss–Bonnet’s theorem. Furthermore, we only allow for terms which guarantee boundedness from below. This excludes all odd geometric orders as well as k_4,1. Considering these points and only contributions up to geometric order N = 4 leads to the following bending energy density

f_BE([script letter H],[scr K, script letter K]) = k_2,0([script letter H] − [script letter H]₀)² + k_4,0([script letter H] − [script letter H]₀)⁴ + k_4,2[scr K, script letter K]²,

which has been considered in ref. 14 to study the stability of tubular shapes considering k_2,0 > 0, k_4,0 = 0 and k_4,2 ≷ 0 as parameters to stabilize cylindrical shapes (k_4,2 > 0) or to enforce pearling (k_4,2 < 0). The same form is also considered in ref. 12 arguing that a cylinder with radius R is stable if k_2,0 < 0 and with k_4,2 not determined. We here restrict the parameter space to k_2,0 > 0, k_4,0 ≥ 0 and k_4,2 ≥ 0. We would like to remark that more general parameter combinations are possible, still leading to well-posed bending energies and stable solutions. However, [scr K, script letter K] = 0 is a property of a tube and thus k_4,2[scr K, script letter K]² seems to be the most plausible higher order extension in the bending energy to stabilize tubular structures. k_4,0([script letter H] − [script letter H]₀)⁴ is considered for completeness. In any case the bending force b is derived as the negative of the variational derivative of the bending energy (8)

A detailed derivation of the bending force is done in the Appendix: Model derivation.

2.4 Numerical approach

The numerical approach extends the approach used in ref. ¹⁷ and ³⁵ which is based on surface finite elements^51,52 and builds on a Taylor–Hood element for the surface Navier–Stokes equations, higher order surface parametrizations, appropriate approximations of the geometric quantities, mesh redistribution, a semi-implicit discretization in time and an iterative approach to deal with the non-local constraint on the enclosed volume. Additional challenges emerge from the higher order geometric terms. In the Appendix: Numerical method we provide a detailed description and in Appendix: Validation convergence studies for these terms. The implementation is done in DUNE/AMDiS.^53,54 Throughout the paper we consider Re = 1.0 and γ = 0 and only vary [script letter H]₀, k_2,0, k_4,0 and k_4,2. Numerical parameters, such as mesh size h and time step τ are chosen to resolve the highest curvature values and to meet stability constraints, following ref. 35.

3 Results

3.1 Axisymmetric simulation without surface viscosity

We begin by describing the membrane using cylindrical coordinates (r, θ, z) and consider a rotational symmetric tube with periodic boundary conditions. With

the bending energy density f_BE([script letter H],[scr K, script letter K]) can be reformulated and considered as a function of r(z, t), the radial distance of the axisymmetric membrane from the cylindrical symmetry axis, where z measures the coordinate along that axis and t is time. For the evolution we consider the corresponding equations to eqn (5)–(7) but drop the constraints on inextensibility and volume. The resulting equation to solve reads

with b = b·ν. While lengthy if fully written down, the resulting model can be solved using standard approaches. We again use finite elements in space, a semi-implicit discretization in time, smoothing of the geometric properties and consider DUNE/AMDiS^53,54 for the realization. Fig. 1 shows the time evolution of a periodically perturbed tube for different values of k_4,2 and k_4,0 together with the evolution of the bending energy [scr F, script letter F]_BE over time. For k_4,2 = 0 and k_4,0 = 0 this is a well studied problem of the stability of a tube.^55–57 The parameters are chosen to remain within the stability region. The results clearly indicate a stronger damping of the perturbations for increasing values of k_4,2 and k_4,0. However, the behaviour slightly differs with respect to the damping of the perturbations and adjustment of r(z). This difference is also seen in the decrease of the bending energy.

Fig. 1 (a) Time evolution of axisymmetric simulations for different parameters for k_4,2 and k_4,0 starting from the same periodic solution at t₀. The evolution goes from dark to light converging to a tube with r(z) = 0.5. The depicted time instances are marked in (b). (b) Corresponding time evolution of the bending energy for different values of k_4,2 and k_4,0. Parameters are k_2,0 = 0.05 and [script letter H]₀ = −2.

Results of a linear stability analysis⁵⁸ for the case of k_4,2 = 0 and k_4,0 = 0 but considering surface viscosity using a Stokes approximation of eqn (1) indicate a similar stability region as in the overdamped limit and we assume this holds also for k_4,2 > 0 and k_4,0 > 0 and the full problem.

3.2 Equilibrium shape of tubular cell

The shape of a tubular cell can be approximated by a cylinder with hemispherical caps, as shown in Fig. 2. For this simple shape, we can compute the bending energy as
where r is the radius and l the length of the cylinder. For , the energy of the cylindrical part vanishes and the energy simplifies to

Fig. 2 (a) As initial surface we consider a cylinder with hemispherical caps. (b) The geometry is determined by the length of the cylinder l and the radius r of the cylinder and hemispherical caps.

Note that this energy is independent of the length of the cylinder. We consider this shape as the initial configuration and the energy as a quantity for comparison.

Fig. 3 shows the final configurations achieved by solving the full problem, eqn (1)–(3), for three different values l. In all cases, the shape deviates from the initial configuration. For k_4,2 = 0 and k_4,0 = 0 the resulting shapes correspond to the equilibrium prolate shapes in ref. 7. However, for k_4,2 > 0 and k_4,0 > 0, these shapes deviate. For k_4,2 > 0 the bending energy is reduced by increasing l and the Hausdorff distance d_H([scr S, script letter S],[scr S, script letter S]₀) to the idealized initial shape of a cylinder with hemispherical caps [scr S, script letter S]₀ is reduced by increasing k_4,2, see Fig. 4. The Hausdorff distance is a measure of the distance between sets of points, in our case the discretization points of the equilibrium shape and the initial shape. For a mathematical definition and implementation issues see Appendix: Hausdorff distance. The observed trends are far from being general. This can also not be expected due to the highly nonlinear coupling of the geometric terms and the considered constraints on volume and inextensibility. However, they confirm the intuition of the potential impact of k_4,2[scr K, script letter K]² in the bending energy on the emerging equilibrium shapes. The effect of k_4,0 > 0 differs. The increased tendency to enforce [script letter H] = [script letter H]₀ leads to undulations with a wavelength related to −2/[script letter H]₀. How well this can be achieved strongly depends on the length l. Therefore also the Hausdorff distance d_H([scr S, script letter S],[scr S, script letter S]₀) to the idealized initial shape of a cylinder with hemispherical caps [scr S, script letter S]₀ does not decrease but leads to a non-monotonic behaviour if considered as a function of l. The dependency of l on the bending energy [scr F, script letter F]_BE is negligible for the considered parameters. To summarize, within the stability region of the classical Canham–Helfrich model, higher order geometric terms related to k_4,2[scr K, script letter K]² seem to stabilize tubular shapes and terms related to k_4,0([script letter H] − [script letter H]₀)⁴ seem to destabilize these shapes.

Fig. 3 (a) Equilibrium shapes for k_2,0 = 0.125, [script letter H]₀ = −2 and k_4,2 = 0. As initial condition we consider a cylinder of radius r = 0.5 with hemispherical caps and vary the length l of the cylindrical part. From left to right the lengths are l = 2, 3 and 4. The surfaces are colored by mean curvature [script letter H] (upper row) and Gaussian curvature [scr K, script letter K] (lower row). (b) Equatorial cuts through the equilibrium shapes for different values of k_4,2 and l, keeping k_2,0 = 0.125 and [script letter H]₀ = −2. (c) Equatorial cuts through the equilibrium shapes for different values of k_4,0 and l, keeping k_2,0 = 0.125 and [script letter H]₀ = −2.

Fig. 4 (a) and (b) Bending energy of the final shape for fixed values of k_4,2 (a) and k_4,0 (b) for different lengths l. (c) and (d) Hausdorff distance between the equilibrium shape and the initial shape for different lengths l = 2, 3 and 4 as a function of k_4,2 (c) and k_4,0 (d).

3.3 Dynamic evolution

In order to further explore the influence of the higher order geometric terms related to k_4,2[scr K, script letter K]² and k_4,0([script letter H] − [script letter H]₀)⁴ on the dynamics, we consider an initial surface [scr S, script letter S]₀ as a perturbed unit sphere

with spherical harmonics Y^m_l and Legendre polynomials P^m_l and the case l = 5, m = 3 and r₀ = 0.5. The velocity field is initialized with u₀ = 0. The initial surface has non-zero mean and Gaussian curvature and is out of equilibrium. It serves as a prototypical example to study the rapid shape changes featuring moderate geometric properties of a cell with membrane reservoirs. The resulting bending force induces shape deformations in the normal direction. However, the curvature terms also induce tangential flows, which also contribute to shape deformations. This coupling between tangential flow and shape deformations is well understood and shown to enhance the evolution towards equilibrium shapes.³⁵ Here, we explore the evolution for different values of k_4,2 and k_4,0, which are k_4,2 = 0, 0.06 and 0.12 and k_4,0 = 0, 0.00625 and 0.0125, respectively. Fig. 5 and 6 show snapshots of the evolutions. The color coding corresponds to shape deformations (red – movement outwards, blue – movement inwards) and the arrows indicate the tangential velocity, where the length scales with the magnitude. Furthermore, the energy contributions are shown over time. Here, the bending energy is split into the different contributions [scr F, script letter F]_BE = [scr F, script letter F]_2,0 + [scr F, script letter F]_4,0 + [scr F, script letter F]_4,2 and the total energy [scr F, script letter F] = [scr F, script letter F]_BE + [scr F, script letter F]_K is the sum of the bending energy and the kinetic energy . All evolutions converge to an equilibrium shape, which for the cases k_4,2 = 0 and k_4,0 = 0 correspond to the associated Seifert shape,⁷ which here corresponds to an oblate. For k_4,2 > 0 the equilibrium shapes only slightly differ. The difference in orientation might result from the different dynamics. However, due to the decoupling from the surrounding bulk phases and the considered parameter γ = 0, force-free rigid body rotations are also possible.⁵⁹ While the equilibrium shapes are similar, significant changes can be observed in the dynamics. The close coupling between the bending energy [scr F, script letter F]_BE and the kinetic energy [scr F, script letter F]_K can be observed and related to significant shape changes. But their appearance differs. The plateau in [scr F, script letter F]_BE = [scr F, script letter F]_2,0 for k_4,2 = 0 between t ≈ 2 and t ≈ 6 is reduced for k_4,2 = 0.06 and 0.12 and already ends at t ≈ 4.5 and t ≈ 3, respectively. It should be noted that the absolute values of the energies cannot be directly compared as k_4,2 varies by definition. The situation changes for k_4,0 > 0. The equilibrium shapes are prolate shapes. Besides this difference in the long time behavior, the short time evolution also changes. We observe similar behavior as for k_4,2 > 0 with an enhanced influence of the kinetic energy and a faster evolution towards intermediate shapes. The final convergence towards the equilibrium shape at late times probably results form similar values for the local minima on the oblate and prolate branch as also known for the classical Canham–Helfrich model.⁶ Again, the absolute values of the energies cannot be directly compared as k_4,0 varies in the definition.

Fig. 5 Evolution of a perturbed sphere considering parameters k_2,0 = 0.125, [script letter H]⁰ = 0 and different values of k_4,2. (a)–(c) Time instances for the evolution for k_4,2 = 0, 0.06 and 0.12, respectively. The color coding indicates movement in the normal direction (red – movement outwards, blue – movement inwards) and the arrows indicate the tangential velocity. The time evolves from left to right; the time instances except for the last one are equal and are depicted in (d). (d) Evolution of the different energy contributions: [scr F, script letter F]^2,0 and [scr F, script letter F]^4,2 are the energies linked to the corresponding bending terms, [scr F, script letter F]^K is the kinetic energy. The bending energy is [scr F, script letter F]^BE = [scr F, script letter F]^2,0 + [scr F, script letter F]^4,2 and the total energy [scr F, script letter F] = [scr F, script letter F]^BE + [scr F, script letter F]^K. Shown is the evolution for k_4,2 = 0, 0.06 and 0.12 (from left to right). Corresponding videos to the evolution in (a)–(c) are provided in the SI† using a LIC filter for visualization of the tangential flow.

Fig. 6 Evolution of a perturbed sphere considering parameters k_2,0 = 0.125, [script letter H]⁰ = 0 and different values of k_4,0. (a)–(c) Time instances for the evolution for k_4,0 = 0, 0.00625 and 0.0125, respectively. The color coding indicates movement in the normal direction (red – movement outwards, blue – movement inwards) and the arrows indicate the tangential velocity. The time evolves from left to right; the time instances except for the last one are equal and are depicted in (d). (d) Evolution of the different energy contributions: [scr F, script letter F]^2,0 and [scr F, script letter F]^4,0 are the energies linked to the corresponding bending terms, [scr F, script letter F]^K is the kinetic energy. The bending energy is [scr F, script letter F]^BE = [scr F, script letter F]^2,0 + [scr F, script letter F]^4,0 and the total energy [scr F, script letter F] = [scr F, script letter F]^BE + [scr F, script letter F]^K. Shown is the evolution for k_4,0 = 0, 0.00625 and 0.0125 (from left to right). Corresponding videos to the evolution in (a)–(c) are provided in the SI† using a LIC filter for visualization of the tangential flow.

Qualitatively the higher order geometric terms further enhance the evolution and lead to alternative pathways to dissipate energy. This is most pronounced in the inlets highlighting the initial evolution; the drastic decrease of [scr F, script letter F]^4,2 and [scr F, script letter F]^4,0 is associated with large fluctuation of [scr F, script letter F]^K. This behavior increases with increasing values of k_4,2 and k_4,0. Furthermore, while [scr F, script letter F]^4,2 is roughly one order of magnitude smaller than [scr F, script letter F]^2,0, the influence on the dynamics is dramatic. This is less pronounced for [scr F, script letter F]^4,0, which is at the same order as [scr F, script letter F]^2,0. We expect these fast shape changes to be even enhanced for less regular real geometries of membrane reservoirs, due to the appearance of larger curvature gradients.

4 Conclusions

Motivated by rapid shape changes of cells, where an excess of membrane that is organized in membrane reservoirs is made available to the cell on the order of seconds,²³ we formulated a minimal mesoscopic membrane model which helps to facilitate this behavior. This, on the one side, includes an extension of the classical Canham–Helfrich model towards higher order geometric terms, and on the other side, includes the explicit treatment of the fluid properties of the membrane by considering membranes as fluid deformable surfaces. The first aspect adds to the classical bending energy , terms proportional to the Gaussian curvature squared, , and the mean curvature minus the spontaneous curvature to the fourth power, . The effects of these higher order terms are different. While the one related to the Gaussian curvature squared not only helps to damp perturbations of tubes, it also has a tendency to stabilize them, the term related to the mean curvature minus the spontaneous curvature to the fourth power has a tendency to destabilize tubes. Both effects have been considered in idealized rotational symmetric and full three-dimensional situations by analyzing the evolution and the emerging equilibrium shapes. If combined with the second aspect, which takes the surface viscosity of the membrane into account and combines the bending in the normal direction with the properties of an inextensible surface fluid, as a fluid deformable surface, the dynamic drastically changes. Already for the classical bending energy an enhanced evolution towards the equilibrium shape has been observed if the effects of surface viscosity are taken into account.³⁵ With the higher order geometric terms this is further enhanced. The considered numerical experiments for the relaxation of a perturbed sphere showed alternative pathways to dissipate energy and strong tangential flows inducing fast shape changes. However, besides enabling rapid shape changes, the evolution and also the emerging equilibrium shapes qualitatively differ for the considered higher order terms. While the model with the Gaussian curvature squared term also converges to similar oblate-like shapes as the corresponding Seifert shapes⁶ for the classical Canham–Helfrich model, the model with the mean curvature minus the spontaneous curvature to the fourth power term converges to prolate-like shapes.

However, the focus of this paper is not to classify the increased phase space of equilibrium shapes resulting from the higher order terms, but to address additional mechanisms which facilitate rapid shape changes. While certainly more research is needed to fully explore the potential of the higher order geometric terms, e.g., with respect to the stability of tubes extending the analysis in ref. 58 the numerical studies already clearly indicate the potential for rapid shape changes. Even if only passive contributions of a homogeneous membrane are considered, and a full model for morphological changes of a cell requires the taking of inhomogeneities, active processes of the underlying cortex, adhesion between the membrane and the cortex and probably even more phenomena into account, the study contributes to identifying underlying general mechanical principles which might help to predict and control the dynamics of cells.²⁰ Surface viscosity and higher order geometric terms in the bending energy provide mechanical cues and probably support active processes to enable rapid shape changes of cells.

A full model able to address the mentioned example of frequently forming and retracting filopodia and the associated fast shape changes of membrane reservoirs will require the coupling of models for the cellular cortex^24–27 with membrane models of the considered type, which requires additional modeling and numerical efforts. One intermediate step to link this research closer to biology could be efforts to resolve the dynamics of membrane reservoirs and compare them with simulations of the proposed model.

Data availability

Data are available from Zenodo at https://doi.org/10.5281/zenodo.14503545.

Conflicts of interest

There are no conflicts to declare.

Appendices

Model derivation

We here compute the variational derivatives of the higher order geometric terms in the bending energy. We therefore write the bending energy as
with bending density components f^n,α_BE = ([script letter H] − [script letter H]₀)^n−2α[scr K, script letter K]^α. The corresponding forces b^n,α ∈ T[Doublestruck R]³|_{[scr S, script letter S]} are given by
for all W ∈ T[Doublestruck R]³|_{[scr S, script letter S]}. The result for f^2,0_BE is well known
For the remaining terms we use the deformation derivative ð_W⁶⁰ and obtain the deformation formulaMoreover,

P(ð_W[scr B, script letter B])P = ∇_{[scr S, script letter S]}(ν·∇_CW) − [scr B, script letter B]∇_CW, (14)

ð_W[script letter H] = tr(P(ð_W[scr B, script letter B])P) = div_{[scr S, script letter S]}(ν·∇_CW) − [scr B, script letter B] : ∇_CW (15)

are valid.⁶⁰ As a consequence, the deformation formula (13) and integrations by parts result in

Another application of integrations by parts w.r.t. ∇_C, yields

where

div_C(σP) = div_{[scr S, script letter S]}(PσP) − ν·σ[scr B, script letter B] + (div_{[scr S, script letter S]}(ν·σP) + [scr B, script letter B] : σ)ν (16)

holds for all σ ∈ T[Doublestruck R]³|_{[scr S, script letter S]}.^41,60 Therefore, the tangential part of b^4,0 cancels out and we obtain
With (14) and (15), the deformation derivative of [scr K, script letter K] = 2([script letter H]² − ‖[scr B, script letter B]‖²) reveals
As a consequence, the deformation formula (13), integrations by parts, and [scr B, script letter B]² = [script letter H][scr B, script letter B] − [scr K, script letter K]P, result in
Since div_{[scr S, script letter S]}([script letter H]P − [scr B, script letter B]) = 0 holds, integrations by parts yields
Using (16), the tangential part of b^4,2 cancels out and we obtain

b^4,2 = −(2div_{[scr S, script letter S]}(([script letter H]P − [scr B, script letter B])∇_{[scr S, script letter S]}[scr K, script letter K]) + [script letter H][scr K, script letter K]²)ν.

Putting everything together yields eqn (4). For the derivation of the other parts of the model we refer to.^17,41

Numerical method

We consider a surface finite element method (SFEM)^51,52 to solve the highly nonlinear set of geometric and surface partial differential equations (1)–(3), using the approaches in ref. 17 and 35.

We combine the system (1)–(3) with a mesh redistribution approach.⁶¹ These are equations for the parametrization

∂_tX·ν = u·ν (17)

[script letter H]ν = Δ_CX, (18)

which generate a tangential mesh movement to maintain the shape regularity and additionally provide an implicit representation of the mean curvature [script letter H]. We consider a discrete k-th order approximation [scr S, script letter S]^k_h of [scr S, script letter S], with h the size of the mesh elements, i.e. the longest edge of the mesh. We use the DUNECurvedGrid library⁶² and consider each geometrical quantity like the normal vector ν_h, the shape operator [scr B, script letter B]_h, the Gaussian curvature [scr K, script letter K]_h, and the inner products (·,·)_h with respect to the [scr S, script letter S]^k_h. In the following, we will drop the index k. We define the discrete function spaces for scalar functions by V_k([scr S, script letter S]_h) = {ψ ∈ C⁰([scr S, script letter S]_h)|ψ|_T ∈ [scr P, script letter P]_k(T)} and for vector fields by V_k([scr S, script letter S][h]) = [V_k([scr S, script letter S][h])]³. Within these definitions T is the mesh element and [scr P, script letter P]_k are the polynomials of order k. We consider u_h,X ∈ V₃([scr S, script letter S]_h), [script letter H]_h ∈ V₃([scr S, script letter S]_h), and p_h ∈ V₂([scr S, script letter S]_h), which leads to an isogeometric setting for the velocity and a [scr P, script letter P]₃ − [scr P, script letter P]₂ Taylor-Hood element for the surface Navier–Stokes-like equations. We discretize in time using constant time stepping with step size τ. In each time step we solve the surface Navier–Stokes-like equations and the mesh redistribution together. We define a discrete surface update variable Yⁿ = Xⁿ − Xⁿ⁻¹, which is considered as unknown instead of the surface parametrization Xⁿ. The system to solve reads:

Find (uⁿ_h,pⁿ_h,[script letter H]ⁿ_h,Yⁿ) ∈ [V₃×V₂×V₃×V₃]([scr S, script letter S]ⁿ⁻¹_h) such that:

for all (ν_h,q_h,h_h,Z_h) ∈ [V₃×V₂×V₃×V₃]([scr S, script letter S]ⁿ⁻¹_h), where .

In the above formulation we used the identity (−∇_{[scr S, script letter S]}pⁿ_h − pⁿ_h[script letter H]ⁿ_hν_h,ν_h)_h = (pⁿ_h,div_Pν_h)_h. Note that the Lagrange multiplier λ is unknown, which leads to an underdetermined system. To resolve that problem, we follow the approach introduced in ref. 35. In order to fulfill the volume constraint, λ has to be chosen such that

We consider Φ(λ) = 0 as an equation in λ and apply a Newton iteration λ^j+1 = λ^j − Φ(λ^j)/Φ′(λ^j). After convergence is achieved the new surface [scr S, script letter S]ⁿ_h needs to be computed by updating the parametrization Xⁿ = Xⁿ⁻¹ + Yⁿ, lifting the solutions uⁿ_h, pⁿ_h and [script letter H]ⁿ_h to the new surface and computing the remaining geometric quantities νⁿ_h, [scr B, script letter B]ⁿ_h, ∇_{[scr S, script letter S]}[script letter H]ⁿ_h, [scr K, script letter K]ⁿ_h and ∇_{[scr S, script letter S]}[scr K, script letter K]ⁿ_h for the new surface. While this approach showed the (expected) optimal order of convergence,³⁵ with respect to computational and numerical analysis results for the underlying surface Stokes equations on stationary surfaces^63,64 for bending terms up to second order, the approach is not sufficient if higher order geometric terms are included.

We therefore introduce a smoothing step of the surface quantities [script letter H]ⁿ_h, ∇_{[scr S, script letter S]}[script letter H]ⁿ_h, [scr K, script letter K]ⁿ_h and ∇_{[scr S, script letter S]}[scr K, script letter K]ⁿ_h. For each surface quantity or its components a_h = [script letter H]ⁿ_h, [∇_{[scr S, script letter S]}[script letter H]ⁿ_h]_i, [scr K, script letter K]ⁿ_h and [∇_{[scr S, script letter S]}[scr K, script letter K]ⁿ_h]_i, i = 1, 2 we solve one time step of the diffusion equation

a_s − εΔ_{[scr S, script letter S]}a_s = a_h,

where ε > 0 is a smoothing parameter and a_s the smoothed surface quantity.

Validation

Instead of a full convergence study of the numerical approach we only test the smoothing of surface quantities. We consider a surface for which the surface quantities can be computed analytically. The surface is parametrized by X : [0,2π) × [0,π) → [Doublestruck R]²,
We compute the L²-error of the surface quantities [script letter H], [scr K, script letter K], ∇_{[scr S, script letter S]}[script letter H] and ∇_{[scr S, script letter S]}[scr K, script letter K] for different grid widths h. The smoothing parameter ε is chosen experimentally such that it shows optimal results. Fig. 7 shows that for this test case the surface quantities converge with the optimal orders and that the additional smoothing step for all quantities improves the approximation. Together with the convergence studies in ref. 17 and 35 these results provide enough confidence in the numerical approach for the full problem including the higher order geometric terms. They require an appropriate resolution of [scr K, script letter K] and ∇_{[scr S, script letter S]}[scr K, script letter K], which is achieved in O(h²) and O(h), respectively. This motivated the considered discrete function spaces.

Fig. 7 (a) Reference surface for the tests. Depicted are the mean curvature [script letter H] and Gaussian curvature [scr K, script letter K]. (b) Errors of geometric quantities with and without additional smoothing step for different grid widths h. The orders of convergence are indicated by the dashed lines and are optimal orders for the numerical implementation.

Hausdorff distance

For two sets X, Y ⊂ [Doublestruck R]³, we consider the Hausdorff distance given by
where d(x, M) = sup_m∈M‖x − m‖ and ‖·‖ is the Euclidean norm. In the implementation we use the VTK Hausdorff distance point set filter⁶⁵ with the target distance method point to cell.

Acknowledgements

We acknowledge fruitful discussions with Patrick Zager (UCSF) and Veit Krause (TUD). This work was funded by DFG within FOR3013 “Vector- and tensor-valued surface PDEs”. We further acknowledge computing resources at FZ Jülich under Grant No. MORPHO and at ZIH under Grant No. WIR.

References

1. H. T. McMahon and J. L. Gallop, Nature, 2005, 438, 590–596.
2. M. Deserno, Chem. Phys. Lipids, 2015, 185, 11–45.
3. R. Lipowsky and R. Dimova, Soft Matter, 2021, 17, 214–221.
4. W. Helfrich, Z. Naturforsch. C., 1973, 28, 693.
5. P. Canham, J. Theor. Biol., 1970, 26, 61.
6. U. Seifert, Adv. Phys., 1997, 46, 13–137.
7. U. Seifert, K. Berndl and R. Lipowsky, Phys. Rev. A:At., Mol., Opt. Phys., 1991, 44, 1182–1202.
8. N. Bobrovska, W. Góźdź, V. Kralj-Iglič and A. Iglič, PLoS One, 2013, 8, e73941.
9. N. Walani, J. Torres and A. Agrawal, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2014, 89, 062715.
10. A. Mahapatra and P. Rangamani, Soft Matter, 2023, 19, 4345–4359.
11. M. D. Mitov, C. R. Acad. Bulg. Sci., 1978, 31, 513–515.
12. J. B. Fournier and P. Galatola, Europhys. Lett., 1997, 39, 225.
13. R. Capovilla, J. Guven and J. A. Santiago, J. Phys. A: Math. Gen., 2003, 36, 6281.
14. T. Shemesh, A. Luini, V. Malhotra, K. N. Burger and M. M. Kozlov, Biophys. J., 2003, 85, 3813–3827.
15. M. Kaksonen and A. Roux, Nat. Rev. Mol. Cell Biol., 2018, 19, 313–326.
16. S. C. Al-Izzi, P. Sens and M. S. Turner, Phys. Rev. Lett., 2020, 125, 018101.
17. E. Bachini, V. Krause, I. Nitschke and A. Voigt, J. Fluid Mech., 2023, 977, A41.
18. M. Arroyo, L. Heltai, D. Millán and A. DeSimone, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 17874–17879.
19. B. Winkler, I. S. Aranson and F. Ziebert, Commun. Phys., 2019, 2, 82.
20. R. K. Sadhu, A. Iglič and N. S. Gov, J. Cell Sci., 2023, 136, jcs260744.
21. V. T. Ruhoff, N. Leijnse, A. Doostmohammadi and P. M. Bendix, Trends Cell Biol., 2024, 35, 129–140.
22. H. De Belly and O. D. Weiner, Curr. Opin. Cell Biol., 2024, 89, 102392.
23. H. De Belly, S. Yan, H. B. da Rocha, S. Ichbiah, J. P. Town, P. J. Zager, D. C. Estrada, K. Meyer, H. Turlier and C. Bustamante, et al., Cell, 2023, 186, 3049–3061.
24. S. R. Naganathan, S. Fürthauer, M. Nishikawa, F. Jülicher and S. W. Grill, eLife, 2014, 3, e04165.
25. A.-C. Reymann, F. Staniscia, A. Erzberger, G. Salbreux and S. W. Grill, eLife, 2016, 5, e17807.
26. H. B. da Rocha, J. Bleyer and H. Turlier, J. Mech. Phys. Solids, 2022, 164, 104876.
27. A. Bhatnagar, M. Nestler, P. Gross, M. Kramar, M. Leaver, A. Voigt and S. W. Grill, Curr. Biol., 2023, 33, 5096–5108.
28. T. Itoh and K. Tsujita, Curr. Opin. Cell Biol., 2023, 81, 102173.
29. H. A. Faizi, R. Dimova and P. M. Vlahovska, Biophys. J., 2022, 121, 910–918.
30. M. Arroyo and A. DeSimone, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 79, 031915.
31. A. Torres-Sánchez, D. Millán and M. Arroyo, J. Fluid Mech., 2019, 872, 218–271.
32. A. Voigt, J. Fluid Mech., 2019, 878, 1–4.
33. S. C. Al-Izzi and R. G. Morris, Semin. Cell Dev. Biol., 2021, 120, 44–52.
34. S. Reuther, I. Nitschke and A. Voigt, J. Fluid Mech., 2020, 900, R8.
35. V. Krause and A. Voigt, J. Comput. Phys., 2023, 486, 112097.
36. M. Simunovic, G. A. Voth, A. Callan-Jones and P. Bassereau, Trends Cell Biol., 2015, 25, 780–792.
37. I. Derényi, F. Jülicher and J. Prost, Phys. Rev. Lett., 2002, 88, 238101.
38. A. Roux, G. Cappello, J. Cartaud, J. Prost, B. Goud and P. Bassereau, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 5394–5399.
39. P. G. Saffman and M. Delbrück, Proc. Natl. Acad. Sci. U. S. A., 1975, 72, 3111–3113.
40. I. Nitschke, S. Reuther and A. Voigt, Phys. Rev. Fluids, 2019, 4, 044002.
41. I. Nitschke and A. Voigt, Adv. Differ. Equ., 2025, 30, 335–420.
42. S. Reuther and A. Voigt, Multiscale Model. Simul., 2015, 13, 632–643.
43. S. Reuther and A. Voigt, Multiscale Model. Simul., 2018, 16, 1448–1453.
44. P. Brandner, A. Reusken and P. Schwering, Interface. Free Bound., 2022, 24, 533–563.
45. T. Biben and C. Misbah, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 67, 031908.
46. Y. Kim and M.-C. Lai, J. Comput. Phys., 2010, 229, 4840–4853.
47. D. Salac and M. J. Miksis, J. Fluid Mech., 2012, 711, 122–146.
48. F. Haußer, W. Marth, S. Li, J. S. Lowengrub, A. Rätz and A. Voigt, Int. J. Biomath. Biostat., 2013, 2, 19–48.
49. S. Aland, S. Egerer, J. Lowengrub and A. Voigt, J. Comput. Phys., 2014, 277, 32–47.
50. V. Krause and A. Voigt, J. R. Soc. Interface, 2024, 21, 20240056.
51. G. Dziuk and C. M. Elliott, Acta Numer., 2013, 22, 289–396.
52. M. Nestler, I. Nitschke and A. Voigt, J. Comput. Phys., 2019, 389, 48–61.
53. S. Vey and A. Voigt, Comput. Visualization Sci., 2007, 10, 57–67.
54. T. Witkowski, S. Ling, S. Praetorius and A. Voigt, Adv. Comput. Math., 2015, 41, 1145–1177.
55. T. R. Powers, Rev. Mod. Phys., 2010, 82, 1607–1631.
56. R. E. Goldstein, P. Nelson, T. Powers and U. Seifert, J. Phys. II, 1996, 6, 767–796.
57. S. C. Al-Izzi, G. Rowlands, P. Sens and M. S. Turner, Phys. Rev. Lett., 2018, 120, 138102.
58. M. Janssen, S. Liese, S. C. Al-Izzi and A. Carlson, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2024, 109, 044403.
59. M. Nestler and A. Voigt, Proc. Appl. Math. Mech., 2023, 23, e202300044.
60. I. Nitschke, S. Sadik and A. Voigt, IMA Journal of Applied Mathematics, 2023, 88, 917–958.
61. J. W. Barrett, H. Garcke and R. Nurnberg, SIAM J. Sci. Comput., 2008, 31, 225–253.
62. S. Praetorius and F. Stenger, Archive of Numerical Software, 2022, 6, 1–27.
63. P. Brandner, T. Jankuhn, S. Praetorius, A. Reusken and A. Voigt, SIAM J. Sci. Comput., 2022, 44, A1807–A1832.
64. H. Hardering and S. Praetorius, IMA Journal of Numerical Analysis, 2023, 43, 1543–1585.
65. F. Commandeur, J. Velut and O. Acosta, VTK Journal, 2011.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4fd00202d

---