Mem3DG: Modeling membrane mechanochemical dynamics in 3D using discrete differential geometry.

Biomembranes adopt varying morphologies that are vital to cellular functions. Many studies use computational modeling to understand how various mechanochemical factors contribute to membrane shape transformations. Compared with approximation-based methods (e.g., finite element method [FEM]), the class of discrete mesh models offers greater flexibility to simulate complex physics and shapes in three dimensions; its formulation produces an efficient algorithm while maintaining coordinate-free geometric descriptions. However, ambiguities in geometric definitions in the discrete context have led to a lack of consensus on which discrete mesh model is theoretically and numerically optimal; a bijective relationship between the terms contributing to both the energy and forces from the discrete and smooth geometric theories remains to be established. We address this and present an extensible framework, Mem3DG, for modeling 3D mechanochemical dynamics of membranes based on discrete differential geometry (DDG) on triangulated meshes. The formalism of DDG resolves the inconsistency and provides a unifying perspective on how to relate the smooth and discrete energy and forces. To demonstrate, Mem3DG is used to model a sequence of examples with increasing mechanochemical complexity: recovering classical shape transformations such as 1) biconcave disk, dumbbell, and unduloid; and 2) spherical bud on spherical, flat-patch membrane; investigating how the coupling of membrane mechanics with protein mobility jointly affects phase and shape transformation. As high-resolution 3D imaging of membrane ultrastructure becomes more readily available, we envision Mem3DG to be applied as an end-to-end tool to simulate realistic cell geometry under user-specified mechanochemical conditions.

INTRODUCTION

Computational modeling of lipid bilayer mechanics has long been accepted as a way to probe the biophysical aspects of membrane curvature generation. The ability of lipid bilayers and cellular membranes to bend in response to various applied forces has been studied extensively from the mathematical modeling perspective. However, the nonlinear system of equations that result from such modeling often leads to a computational bottleneck to generate predictions from simulations that can be tested against experimentally observed shapes. In this study, we develop a mesh-based model using discrete differential geometry (DDG) to reduce this bottleneck. To justify why our method is necessary and is a computational advance, we first describe the importance of membrane curvature generation in biology, the current state of the art in membrane mechanics modeling, and finally explicitly state the goals of our approach.

As one of the most essential and conserved structures of cells, cellular membranes perform many functions. First, they form compartments to separate chemical environments. Beyond the passive role of partitioning space, lipids in the membranes interact with proteins and other cellular components influencing cell signaling (e.g., by localizing molecules and acting as an entropic barrier) (1,2). Membrane morphology and topology changes are critical for trafficking cargo in and out of cells and are very carefully regulated (3–8). Central to these roles is the ability of the membrane to curve and adopt varying morphological configurations from spheres to highly-curved and variegated structures.

Advances in experimental studies of membrane-protein interactions (9–20), ultrastructural imaging (21–30), and image analysis (9–11,31–37) have revealed much about the molecular interactions that regulate membrane curvature. To investigate the mechanics behind these interactions, many theoretical and computational models in terms of membrane energetics and thermodynamics have been developed (7,38–52). These models, owing to the ease of in silico experimentation, have become an important tool for generating and testing hypotheses (53,54). These mechanics models and associated simulations have been used to provide intuition on the mechanical requirements for forming and maintaining complex cellular membrane shapes (55–63).

While the utility of this approach has been established and many models have been developed (38), many models are limited by critical assumptions or other technical challenges. For example, the ability to use geometries from membrane ultrastructural imaging experiments as a starting condition would improve model realism (64). With respect to computational complexity, the solver should be able to model deformations and topological changes in three dimensions and be compatible with both energy minimization and time integration for comparing with static and time-series experiments respectively. This is in contrast to the current assumptions of the existence of an axis of symmetry that is quite commonly made for purposes of ease of simulation (65). An additional feature for these solvers should be that their implementation is modular such that the addition of new physics or increasing model complexity should be straightforward. This includes the potential for coupling the membrane model with agent-based and other simulations to propagate other cellular components such as the cytoskeleton (66). Thus, new computational tools which are general, easy to use, and without restrictive assumptions are needed to bring modeling closer to experimental observations of membrane shapes in cells.

To emphasize the motivations behind our choice of extending and developing a new mesh-based membrane model, we provide a summary of the legacy literature in modeling membrane mechanics. The most common theoretical model of membrane bending is the Helfrich-Canham-Evans Hamiltonian (The Helfrich energy is related to the Willmore energy in the mathematics literature (67)), which describes the lipid bilayer as a 2D fluid-like structure that exhibits resistance to bending in the out-of-plane direction (39,40,68–70). It is a continuum model that describes the bending energy of the membrane as a function of its mean and Gaussian curvatures. The assumptions for the continuum are satisfied as long as the deformations are much larger in length scale compared with the individual lipid components.

Given the necessary material properties and boundary conditions, by minimizing the Helfrich energy, we can obtain the equilibrium shape of the membrane (39,70–72). While straightforward in concept, energy minimization requires the determination of the forces on the membrane, which is a challenging task (65). The forces on the membrane are given by the variation of the energy with respect to the embedded coordinate (i.e., shape) of the membrane (we call this variation the shape derivative, which is distinct from the chemical derivative that will be introduced later in the context of mechanochemical coupling). Taking the shape derivatives of the Helfrich energy produces the “shape equation,” so termed because solutions of this partial differential equation (PDE), with the prescribed boundary conditions, produce configurations at equilibrium (i.e., force balance).

Solving the shape equation is non-trivial since it is a PDE with fourth-order nonlinear terms. As a result, analytical solutions of the shape equation are known only for a few cases constrained to specific geometries and boundary conditions (42). For most systems, we must resort to the use of numerical methods. The simplest numerical schemes can be formulated by making restrictive assumptions such as considering only small deformations from a plane (e.g., Monge parametrization) or assuming that there exists an axis of symmetry such that the resulting boundary value system can be integrated (38). While these methods are suitable for idealized shapes, these assumptions are not consistent with the membrane shapes found in biology are and thus not general enough for advancing the field.

Solvers of membrane shape in 3D have also been developed and can be categorized into three groups: 1) phase field or level set methods (73–77), 2) finite element method (FEM) (78–85), and 3) discrete surface mesh models (60,86–99). These methods and others, reviewed in detail by Guckenberger et al. (100), differ in the strategy used to discretize the membrane domain and compute the relevant derivatives. We compare the aforementioned general, 3D models with our established model criteria in Table 1 and elaborate below.

TABLE 1

Comparison of common mathematical frameworks for modeling membrane mechanics with specifications to advance the mission of computational membrane mechanobiology

	Phase field/level set	FEM	Discrete mesh/Mem3DG
General 3D	✔	✔	✔
Statics + dynamics	✔	✔	✔
Membrane heterogeneity	✔	✔	✔
Incorporation of agents/particles			✔
Incorporation of stochastic dynamics (e.g., DPD or MC)			✔
Explicit surface parametrization		✔	✔
Coordinate-free evaluation			✔
Ability to support topological changes	✔		requires mesh surgery
Error analysis	✔	✔

A general framework will permit the easy transfer of inputs and results between model and experiments. Models that can be coupled with other modeling schemes representing other cellular components can help address the complexity of cell biology. Discrete mesh models have many desirable traits, with respect to these specifications, at the cost of forgoing rigorous error analysis.

Phase field and level set methods solve the shape equation by propagating a density field on an ambient volumetric mesh. The membrane shape is implicit in these models and can be found by drawing an isosurface or level set of the model. While this is ideal for modeling membrane topological changes, the implicit representation of the membrane adds complexity for interfacing with data generated using modern methods of visualizing membrane ultrastructure. The meshes output from ultrastructural studies must be converted into a density or phase field prior to input to the model. While this conversion is possible, representing the dynamic and variegated shapes of cellular membranes would require a dense volume mesh, which reduces computational tractability. The implicit surface representation also complicates the addition of new in-plane physics for end users.

FEM and discrete mesh models use an explicit surface parametrization (i.e., a mesh). Thus the meshes output from ultrastructural imaging datasets can be used in these frameworks with minor modifications (32,101). FEM relies on elementwise interpolation functions and is commonly derived from the weak formulation of boundary value problems. Comparing FEM methods with our specifications we identify a few key challenges. First, the numerical evaluation of smooth geometric measurements on arbitrary manifolds in an FEM framework requires non-intuitive tensor algebra to translate the shape equation in coordinate where it is ready to be solved. After this formulation, solving the shape equation can require the use of high-order function basis such as the C1 conforming FEM based on subdivision scheme (78,79) or isogeometric analysis (IGA) (81–83,85), which adds code complexity and run-time cost. Extending an FEM framework to incorporate new physics, topological changes, or interfaces with other models requires advanced mathematical and coding skills. This can restrict the usage to the computational math community and prevent broad usage by the biophysics community.

Finally, evaluating discrete mesh-based methods, which define the system energy and/or forces using geometric primitives from a mesh, we find that they satisfy many of the requirements in Table 1. Due to the ease of use and implementation, discrete mesh models have gained in popularity and many different schemes can be found in the literature (60,86–99,102,103). These schemes differ in their approach to defining and computing geometric measurements necessary for defining the energy and forces on a discrete object. Discrete geometries have discontinuities and limited information that leads to degenerate definitions for geometric values. For example, there is no canonical definition for the normal of a vertex of a mesh as opposed to the normal of a smooth geometry (89,104,105). One challenge for selecting the suitable formulation to use is the lack of approximation error metric for which the discrete definition best matches the smooth theory. Another confounding factor is the step at which the problem is discretized. Some implementations discretize the energy of the system by constructing standalone discrete energy, which captures the behavior of the Helfrich energy (65). From this discrete energy, they take the shape derivatives to obtain an expression for the discrete force. Without careful consideration, the discrete forces derived in this manner are unstructured and there is little resemblance to expressions of force from smooth theory. A second option is to discretize the smooth force expression directly (65,100). While this preserves the geometric connection for the forces, there is no longer well-defined discrete energy. Several discrete mesh methods were benchmarked by Bian et al. (89) and Guckenberger et al. (100) who found differences in the accuracy, robustness, and ease of implementation (89,100).

In this work, we outline a discrete mesh framework for modeling membrane mechanics with the following goals in mind: 1) we do not make a priori assumptions about axes of symmetry or restrict the coordinates in any way; 2) we resolve the ambiguity in the definition of geometric measurements on the mesh and permit a direct comparison for both the energy and force expressions in smooth and discrete contexts; and 3) this framework allows for use of meshes generated from ultrastructural imaging. We begin by defining discrete energy that is analogous to the Helfrich energy. Then, using concepts from DDG, we derive discrete shape derivatives analytically and group terms to produce a discrete shape equation. We will show that our discrete shape equation has a clear correspondence between the terms of the smooth shape equations (57,67,70,71). Beyond establishing this important connection, we will show that the elegant analytical expressions for discrete variational terms from the DDG also yield improved geometric intuition and numerical accuracy (104,105).

Benchmarking of our expressions was performed with our accompanying software implementation called Membrane Dynamics in 3D using Discrete Differential Geometry (Mem3DG). Mem3DG is written in C++, released under the Mozilla Public License version 2, and comes with accompanying documentation and tutorials which can be accessed on GitHub (https://github.com/RangamaniLabUCSD/Mem3DG). Beyond the computation of discrete energies and forces on a mesh of interest, we also include functionality for performing energy minimization and time integration. Using Mem3DG, we validate the exactness of the analytical expressions of force terms by numerically examining the convergence of the force terms as a function of system energy perturbation. To illustrate compliance with our tool specifications, we apply Mem3DG to a sequence of examples with increasing complexity. Finally, we outline the steps to incorporate additional physics such as membrane-protein interactions and surface diffusion into Mem3DG.

THEORY

The lipid bilayer is modeled as a thin elastic shell using the Helfrich-Canham-Evans Hamiltonian or spontaneous curvature model (39,69,106). The bending energy, E, of a smooth surface or 2-manifold, , can be expressed in terms of the mean H, Gaussian K, and spontaneous curvature with material parameters κ the bending and κ the saddlesplay moduli. Additional energy terms Es and E account for the tension area (λ–A) and pressure-volume (ΔP–V) relationships. The total energy of the bilayer is therefore

The preferred surface area and volume, and , combined with the spontaneous curvature, , characterize the zero-energy state for the system energy. In a nutshell, given the material properties, the system energy is fully determined by its geometric measurements such as volume, area, and curvatures.

Machinery to express these measurements have been a topic of extensive study in classical differential geometry (107,108). However, finding the minima of the governing energy, solving the stationary solution to the geometric PDE, can be mathematically and numerically difficult. While differential geometry provides succinct expressions to describe the measurements in a coordinate-free fashion, computational methods often require the introduction of a coordinate basis and subsequent manipulation of expressions using tensor algebra, which can obscure the underlying geometric intuition.

As an alternative, forgoing the need for a smooth geometry, one can treat a discrete geometry (such as a geometric mesh) as the input. This perspective where the discrete geometry is the actual geometry is that of DDG (109). By eliminating the burdens of treating the input mesh as an approximation of a smooth object, DDG capitalizes upon the piecewise nature of meshes to produce efficient and parallelizable finite difference-like formulae which are amenable to algorithmic implementation while maintaining clear geometric meaning. In the following sections, we use concepts from DDG to formulate a discrete analog to the smooth membrane shape problem. Following the derivation of the discrete theory, we describe the development of an accompanying software implementation called Mem3DG.

Notation and preliminaries

We assume the following notation conventions and provide a table of important symbols (Table 2). To aid the reader on how the elements of the mesh are used in the derivation, several fundamental geometric primitives (i.e., values on a mesh which are easily measurable; listed in Table 2A) are illustrated in Fig. 1 A–C.

TABLE 2

Glossary of commonly used symbols and conventions

A. Geometric primitives
ℳ	smooth or discrete 2-manifold
r→∈ℝ^3	embedded coordinate of ℳ
l	edge length
∠	corner angle
φ	dihedral angle
A	area of mesh cell,
	e.g., face Aijk, edge Aij and vertex Ai
n→	surface normal
B. Surface integral
∫a	integrated quantity over mesh
	cell; e.g., Aiai or Aijkaijk
∑_v_i	sum over all vertices vi of the mesh
∑_e_ij	sum over all edges eij of the mesh
∑_f_ijk	sum over all faces fijk of the mesh
∑_v_j∈N(a)	sum over the vertex vj in the neighborhood of a
∑_e_ij∈N(a)	sum over the edges eij in the neighborhood of a
∑_f_ijk∈N(a)	sum over the face fijk in the neighborhood of a
C. Tensors
x∈ℝ	scalar quantity
xindextype	sub- and super-script convention;e.g., ∫f→ib is
	the bending force for vertex vi
x→∈ℝ^3	vector quantity
x = {xi}	(n×l) indexed scalar quantity
x→=x_i	(n×3) indexed vector quantity
X˜	matrix or tensor quantity
D. Derivatives
∇_r→	shape derivative
∇ϕ	chemical derivative
∇_θ→	surface gradient
a˙	time derivative
Δs	Laplace-Beltrami operator
E. Physical variables
E	energy
f	force density
μ	chemical potential
H	mean curvature
K	Gaussian curvature
A	surface area
V	enclosed volume
.¯	preferred state;e.g., H¯ is the spontaneous curvature
ϕ ∈ [0,1]	protein density parameter
λ	membrane tension
ΔP	osmotic pressure across the membrane
κ	bending rigidity
κG	Gaussian modulus
KA	stretching modulus
KV	osmotic strength constant
c¯	molar ambient concentration
n	molar quantity of enclosed solute
η	Dirichlet energy constant
ε	adsorption energy constant
ξ	membrane drag constant
B	protein mobility constant

FIGURE 1

Overview of the DDG framework

For a Figure360 author presentation of this figure, see https://doi.org/10.1016/j.bpr.2022.100062. (A–C) Illustrations of geometric primitives in the neighborhood of (A) fan around a vertex, (B) diamond around an edge, and C) triangle on a face. (D) Discrete definition of scalar edge mean curvature, ∫ H, scalar vertex Gaussian curvature, ∫ K, and Laplace-Beltrami operator, ∫ Δ( · ). (E) Comparative derivation of Helfrich shape equation in both smooth and discrete formulation.

We note that, in discrete contexts, the notation, ∫ a, should be considered the discrete (integrated) counterpart of a pointwise measurement a in a smooth setting. The rationale and significance behind the usage of an integrated quantity in discrete contexts are elaborated in Appendix B and the DDG literature (104,105). Using this notation, discrete surface integrals are expressed as sums of integrated values over the discrete mesh components listed in Table 2B (e.g., is the discrete analog to ). It is possible to interchange between integrated, ∫ a, and pointwise, a, quantities by using the dual area (A),

For simplicity, we will not use separate notations for operators applying in smooth and discrete settings. The context can be inferred from the objects to which the operators are applied. Where it serves to improve our geometric or other intuition, smooth objects will be presented alongside discrete objects for comparison.

Obtaining a discrete energy defined by mesh primitives

Following the perspective of DDG, we restrict our input to the family of triangulated manifold meshes, (i.e., discrete 2-manifolds embedded in ) (We will use for both the smooth and discrete surfaces).

Paralleling the smooth Helfrich Hamiltonian (Eq. 1), a functional of geometric measurements of a surface, the discrete Helfrich Hamiltonian is composed of discrete analog of those measurements,

In comparison with Eq. 1, H and K are pointwise mean and Gaussian curvature measurements on vertices, is the integrated Willmore measure, and the smooth surface integral is replaced by its discrete analog (i.e. finite summation), (Table 2B).

The geometric properties of a given membrane configuration can be connected to the system’s energy through constitutive relations. In this work, we assume that the surface tension follows a linear stress-strain model (110), where is the preferred surface area of the membrane, and K is the stretching modulus of the membrane. The osmotic pressure can be defined based on the van’t Hoff formula as where i, R, T, , and n are the van’t Hoff index, ideal gas constant, temperature, ambient molar concentration, and molar amount of the enclosed solute. Substituting these constitutive relations (Eqs. 4 and 5) into the energy (Eq. 3), we get explicit expressions for E and E, where is the ratio of the concentrations of the ambient and enclosed solutions. Note that the preferred volume, , which is needed to evaluate the integral in Eq. 3, is related to to the parameters in Eq. 5 by . If the system is around the isosmotic condition (e.g.,), the leading order of the energy is given as, where K ≡ iRTn groups the phenomenological parameters. Mathematically, Eqs. 4 and 7 effectively prescribe a penalty-based method for area and volume control. An alternative approach is the use of Lagrange multipliers, which have been extensively adopted in the literature (40,57,72).

To compute the energy of a system, we must obtain values for several geometric measurements that appear in the discrete energy function (e.g., H, K, A, V). For measurements such as the volume and area, there are clear approaches for their evaluation on a triangulated mesh: summing the areas of the triangular faces and summing over the signed volume of tetrahedra (Fig. 1 E, osmotic pressure and surface tension). For other measurements such as the discrete mean and Gaussian curvatures, additional care must be taken. While in smooth contexts these curvatures have unique definitions, in discrete contexts there are multiple approaches for their calculation. For example, the mean curvature can be computed via the application of the cotangent Laplacian, the kernel of the heat equation, or fitting polynomials to a local patch (65). As mentioned earlier, there are challenges with these approaches that can limit their numerical accuracy and obscure the connection to smooth theory. Here using discrete exterior calculus and identification of geometric invariants, we produce theoretically and numerically consistent discrete expressions.

Similar to the polygonal curve introduced in Appendix B, a triangulated mesh has zero curvature on facets and ill-defined curvature on edges and vertices. Using the Steiner view, the sharp corners formed by vertices and edges are made smooth with portions of spherical and cylindrical shells, which have well-defined mean curvature (Fig. 1 D). Taking the limit as the radii of the cylinders and spheres decrease, the leading order contribution of total mean curvature is given by the Steiner formula on an edge, referred to as the edge mean curvature, where l is the length of edge e, and φ is the dihedral angle on e (i.e., the angle formed by the face normals of the neighboring triangles incident to e) (illustrated in Fig. 1 B) (104,105). While not necessary, a triangulated mesh is often realized in via vertex positions; thus it is conventional to prescribe data on vertices instead of edges. Summation of edgewise quantities over the “fan” neighborhood (Fig. 1 A) provides the recipe of converting an edgewise to a vertexwise quantity, where the prefactor, 1/2, accounts for fact that each edge is shared by two vertices.

While we have an integrated mean curvature, the discrete Helfrich Hamiltonian contains a pointwise mean curvature squared term. To define a pointwise mean curvature, the size of the domain occupied by the integrated mean curvature needs to be specified (cf., Appendix B for rationale). The area, A, referred to as the dual area of the vertex v, can be defined as one-third of the areal sum of the incident triangles (fan illustrated in Fig. 1 A) (104,105). Applying Eqs. 2 and 9 to Eq. 3, the pointwise mean curvature is thus,

Substituting Eq. 10 into the integrated Willmore measure term of Eq. 3, the integrated Willmore measure can be expressed as a function of the integrated mean and spontaneous curvature,

Discrete Gaussian curvature is given by the angle defect formula, which is a well-known quantity that preserves many properties parallel to the smooth theory (e.g., Gauss-Bonnet, turning number, among other invariants). One way to derive the angle defect formula is to compute the area of a spherical n-gon contained by a local Gauss map of the neighboring n faces around a vertex (104,105).

Eq. 12 provides the general geometric definition to obtain the energetic contributions from the Gaussian curvature terms. In this study, we consider only systems with uniform saddle-splay modulus which do not undergo topological changes. For these systems, the energy can be simplified based on the discrete Gauss-Bonnet theorem, which states that where , is the Euler characteristic of a topological invariant where |V|, |E| and |F| represent the number of vertices, edges and faces of the mesh respectively, and is the discrete geodesic curvature, which measures the deviation of the boundary curve from a straight line when the surface is locally flattened. In summary, for this work, the Gaussian curvature term is non-constant only when is not closed, and the energy solely involves the boundary elements.

A numerical comparison of the discrete scalar measurements with their smooth counterparts is shown in Fig. E.1. We note that for all geometric measurements (i.e., volume, area, and curvatures), unlike in smooth differential geometry where their numerical evaluation requires the introduction of coordinates, DDG measurements are functions of mesh primitives. By substituting these discrete geometric measurements from DDG into Eq. 6 and 3, we get a numerical recipe for computing the total system energy.

FIGURE E.1

Pointwise magnitude comparison of continuous and discrete measurements: (A) scalar mean curvature, (B) scalar Gaussian curvature, (C) (scalar) bi-Laplacian term ∇H based on the cotan formula, (D) vector mean curvature, (E) vector Gaussian curvature, and (F) (vector) bi-Laplacian term based on Schlafli vector. Note that the result of the cotangent Laplacian approach in (C) produces a scalar result while our approach using the Schlafli vector in (F) is a vector result, thus their direct comparison is not meaningful.

Force from discrete shape derivative of energy

We can obtain the force by taking the negative shape derivative of the energy. In continuous settings, the differentiation is an infinite-dimensional problem that requires the use of the calculus of variations and differential geometry to find analytical expressions (39,43,70,71) (Fig. 1 E, smooth). Deriving the forces from the discrete energy (Eq. 3) is a much simpler task.

Discrete forces can be obtained by taking partial derivatives of mesh primitives with respect to the vertex embedded coordinates, (Fig. 1 E, discrete). Regarding notation, despite the computational differences, the differential operations in both the discrete and smooth contexts are called (discrete) shape derivatives and denoted as due to the common geometric meaning. We note that the computation of discrete shape derivatives for membrane modeling has been described previously in the literature (87,89). Also that there are many overlapping definitions for discrete curvature, energy, and variations thereof in the graphics literature (111–113). Our work extends upon the prior art that evaluates derivatives algebraically, by introducing simplifications based upon the grouping of terms and identification of geometric objects. These simplifications have important implications for improving the geometric understanding as well as the run-time and numerical performance of an implementation.

At the high level, our goal is to express the forces on each vertex, given a set of physics, using geometric primitives and properties defined on specific mesh elements. By grouping terms, we find that the vertexwise forces arising from the different physics can be expressed as weights that are functions of the input parameters and system configuration, multiplied by basic geometric vectors. We will show that these terms have an exact correspondence to terms in the smooth shape equation (Fig. 1 E). We remark that, in some sense, the force expressions are reminiscent of finite-difference equations, which approximate differentials as a linear combination of values at discrete points. This may have implications for the suitability of modeling smooth biological surfaces with discrete meshes.

Force from osmotic pressure

For the smooth geometry, the shape derivative of the enclosed volume yields the outward-pointing surface normal with its size equal to the local area element (114). For a discrete mesh, the shape derivative of the volume is given by the face normal on triangular faces with its local area element equaling to the face area, which is referred to as the integrated face normal, (Fig. 1 E, osmotic pressure) (89,99,104,105), where (ijk) denotes the symmetry under index permutation (e.g., a means a = a). Similar to edge values, the force normal can be converted to vertex normal, where analogous to Eq. 9, the prefactor 1/3 accounts for fact that each face is shared by three vertices. The discrete vertex forces from the derivative of the pressure-volume work, , is then given by scaling it with the uniform osmotic pressure,

Forces from surface tension

Next, considering the shape derivative of the surface energy, E, in smooth contexts, the derivative of the total surface area also points at the surface normal, with its magnitude measuring the size (dA) and the curvature (2H) of the local patch (Fig. 1 E, surface tension) (114). In a discrete case, we can directly compute the derivative of total area on each vertex by summing the area gradient of incident triangles with respect to the vertex position; the sum is therefore referred to as (twice of) the integrated mean curvature vector on vertices, where we define , and is the mean curvature vector on a triangle face corner. The capillary forces from surface tension, , are given by scaling the integrated mean curvature vector by the surface tension,

Evaluating the algebraic sum of area gradients reveals the “cotangent formula” applied to the vertex positions (Fig. 1 E, surface tension). From independent derivations with unrelated frameworks (e.g., discrete exterior calculus and FEM), discretizing the smooth Laplace-Beltrami operator on a triangulated mesh produces the cotangent formula, which is called the discrete Laplace-Beltrami operator, ∫ Δ (104,105,111). Inspecting the expressions in Fig. 1 E, surface tension, we see that our discrete expression parallels smooth theory, where the mean curvature is related to the coordinates through the application of the smooth Laplace-Beltrami operator,

Forces from bending

To evaluate the shape derivative of the discrete bending energy, we consider the terms from the Gaussian and mean curvature separately. Since we do not consider nonuniform saddle-splay modulus and topological changes in this work, the total Gaussian curvature only varies if the surface has boundaries, (cf., discrete Gauss-Bonnet theorem in section “Obtaining a discrete energy defined by mesh primitives”). The shape derivative of the bending energy requires more algebra and the introduction of halfedges, (cf., Appendix C.1). Here we will focus on the key results and refer the reader to the full derivations for each term in Appendix C.2.

There are four fundamental geometric vectors on halfedges that comprise the bending force at non-boundary vertices: the mean curvature vector (see Fig. 1 B for indices), the Gaussian curvature vector, and the two Schlafli vectors, which act to smooth the profile of local dihedral angles. Note that the shape derivatives are taken with respect to different vertices (i.e., or ), such that the mean curvature , Gaussian curvature , and Schlafli vectors are asymmetric under index permutation. To account for the orientation, we refer to them as halfedge vector quantities on (Appendix C.1). A numerical comparison of the discrete geometric vector with their smooth counterparts is shown in Fig. E.1 and Fig. E.2.

FIGURE E.2

Pointwise directional comparison of continuous and discrete measurements: discrete vertex normal based on (A) volume gradient, , mean curvature vector, , and Gaussian curvature vector, , and (B) Schlafli vector,.

The bending force (Fig. 1 E, bending) can be expressed as weights, which are functions of input parameters multiplied by basic geometric measurements in scalar and vector form, which parallels the shape derivative of the smooth bending energy, where is the shape derivative in the surface normal direction.

Comparing the smooth-discrete expressions, we make a few observations:

The Schlafli vector terms, , is the discrete analog of the smooth biharmonic term, , the high-order local smoothing force. The numerical comparison of these two terms, as well as results directly obtained using cotangent formula applied on pointwise scalar mean curvature, are covered in Fig. E.2 and Fig. E.1.

Eq. 23 is the normal component of the shape derivative of the bending energy; an additional tangential component is required if surface heterogeneity exists (e.g., κ is not spatially uniform) (40,65). In contrast, the discrete shape derivative (Eq. 22) is the total derivative in , which includes both the tangential and normal components (in the smooth sense since there is no well-defined vertex normal direction in discrete geometry). Depending on the extent and symmetry of the heterogeneity, the discrete force can point in any direction in .

The coefficients in Eq. 22 show an intriguing pattern combining contributions from both v and v. From a finite-difference approximation standpoint, this results in an approximation scheme for which a rigorous error analysis has not yet been conducted.

Net force and the benefit of DDG

By summing the force terms from each physics, we obtain the net force. Through section “Obtaining a discrete energy defined by mesh primitives” and section “Force from discrete shape derivative of energy,” we identify and show a scheme where both 1) the force is analytically derived from the discrete energy, and 2) both the discrete energy and force mirror the smooth theory. The entire process of defining energy and conduction shape derivative do not involve the introduction of coordinate and the use of tensor algebra. Owing to the discontinuities and limited information contained by a discrete geometry, there are ambiguities in geometric definitions that behave otherwise in the smooth context (cf., various discrete curvature definitions for plane curve discussed in Crane and Wardetzky (115)). Intentional choices of certain definitions of basic discrete geometric measurements can reveal the connection between various definitions, preserve useful geometric invariants, and most naturally reflect the underlying physics. Here many scalar and vector definitions of geometric measurement are connected through the chain of shape derivatives (cf., Fig. A.2) (105), which justifies their role in representing either energy or forces. For example, the discrete bending energy is also commonly constructed using the mean curvature vector, in literature (90,99). As shown in section “Forces from surface tension,” the definition of the mean curvature vector is tightly correlated with the surface tension, where directionality is embedded. The inner product used in such discrete energy definition strips away the directional information. Instead, here we construct energy using the scalar mean curvature, ∫ H, because 1) the energy is inherently scalar, and 2) the discrete curvature exists on the edges. After taking the shape derivative of the energy, the mean curvature vector appears as the effective tension component and the directional information of the vector is used for representing the force. Using the directional information, an arbitrary definition of a vertex normal is avoided. When weighted homogeneously around a vertex, each fundamental geometric vector that composes the discrete force, , , , , can be used to obtain a meaningful definition of the vertex normal. The heterogeneous weighting of these vectors around the vertex represents the incorporation of functional variation in the tangential direction, which can be used to model heterogeneities in material and other properties across the membrane (40,65). Practically, the additional structure provided by the discrete force and energy expressions allows the user to inspect term-wise contributions, which can lead to additional insights and analysis. Since the terms of the discrete energy and forces are defined locally by mesh primitives at vertex neighborhoods, the algorithms are efficient and straightforward to parallelize. The numerical accuracy of these expressions is benchmarked for several scalar and vector measurements on smooth and discrete surfaces shown in Fig. E.1, Fig. E.2, and later discussed in section “Practical considerations for applying .”

FIGURE A.2

Steiner’s formula in continuous and discrete geometry: chain of smooth and discrete shape derivatives of integrated geometric measurements (144).

Defining metrics for simulation and error quantification

For monitoring simulation progress, exactness of force calculations with respect to the discrete energy, and convergence studies of computed quantities upon mesh refinement, we introduce the following norms.

L2 norm

From a PDE perspective, the vertex forces are also called the residual of the shape equation, whose solution represents the equilibrium solution. The simulation is terminated when the residual is smaller than a user-specified threshold. The rationale for using the L2 norm is justified by perturbing the system configuration and conducting an expansion on the system energy, where we refer the inner product of the force matrix as the L2 norm of the forces. Computationally, this is equivalent to the standard Frobenius matrix L2 norm,

Using the L2 norm and Eq. 24, we can perform a numerical validation of the exactness of the discrete force calculation with respect to the discrete energy. We expect the force to approximate the energy up to second order with respect to the size of a perturbation. This validation will be further elaborated in section “Membrane dynamics with full mechanochemical feedback.”

L1 norm

A scale-invariant L1 norm is well suited to quantify the magnitude of the error on varying domain size and mesh resolution. Given a vertexwise local scalar measurement, ∫ , or a vector measurement, , and their reference values, , and , where the normalizing factor, the total surface area A, is used to obtain a pointwise estimate of the error. The L1 norm is applied in the local comparison of discrete and smooth measurements, which we further elaborate in section “Practical considerations for applying .”

SOFTWARE IMPLEMENTATION: Mem3DG

Along with the theoretical developments, we have developed an accompanying software implementation written in C++ called Mem3DG. Our goal in developing this software is to enable the easy use and application of the corresponding theory developed above to biological problems of interest.

Mem3DG is a library that contains several components to support this goal. Fig. 2 provides a synopsis of Mem3DG. The input to Mem3DG includes a triangulated mesh with its coordinate embedded in . Users can choose to use Mem3DG to construct idealized meshes (e.g., icosphere, cylinder, or flat hexagonal patch) as an input or to read in meshes from several common mesh formats. Meshes are stored and manipulated in Mem3DG using the halfedge data structure provided by Geometry Central (116). The supported input file formats are those which are readable by hapPLY and Geometry Central (116,117). Once a mesh and parameters are loaded, Mem3DG can evaluate the discrete energy and forces of the system. Mem3DG adopts a modular design that facilitates the use of different energy and force components and has utilities which help the user to specify the physics and governing parameters. Mem3DG also supports local system simulations where the input mesh has boundaries. Additional details about the supported boundary conditions are given in section “Prescribing boundary conditions with force masking.”

FIGURE 2

Overview of data flow within Mem3DG. The user provides Mem3DG with an initial condition in the form of a triangulated mesh and vertexwise state and kinematic variables (green box). The main loop (black loop) of Mem3DG evaluates the discrete energy and forces and propagates the trajectory, among other supporting steps. Modules in dashed lines are optional depending on whether the system mesh has boundaries and if external forces are specified. User-specified options and possible extensions of Mem3DG to accommodate various physics are highlighted in yellow boxes. Mem3DG automatically exits the simulation when the system converges or the maximum time step is reached.

To perform energy minimization and time integration of the system, various schemes have been implemented. These schemes are described in section “Time integration and energy minimization.” As discussed further in section “Practical considerations for applying ,” when a user wishes to use Mem3DG to represent complex biological membrane geometries, additional care regarding the quality of the mesh is necessary. Mem3DG includes algorithms for basic mesh regularization and remeshing, which can be toggled by the user to support their applications. The simulation terminates when it reaches the time limit or the system reaches equilibrium, whose criteria is determined using the energy L2 norm introduced in section “Defining metrics for simulation and error quantification.” A user can choose between several formats to output a trajectory over time or the configuration of the local minima from Mem3DG. In addition to the mesh outputs supported by Geometry Central, we have also developed a scheme for outputting mesh trajectories in NetCDF format (118). Mem3DG can read and visualize the output trajectories and mesh configurations using Geometry Central and Polyscope (116,119).

For rapid prototyping and enumeration of simulation conditions, we have also developed a Python API called PyMem3DG. The functionality in C++ is exposed in Python using bindings from pybind11 (120). Illustrative examples of using both Mem3DG and PyMem3DG are provided in the online tutorials. For the experiments discussed in this work, all of the simulations were performed using PyMem3DG and the accompanying code and initial configurations are on GitHub: https://github.com/RangamaniLabUCSD/Mem3DG.

Defining properties of a membrane reservoir for systems with open boundaries

To facilitate correspondence with wet experiments and to support the reduction of computational cost, it is possible to construct systems using meshes with open boundaries in Mem3DG. For example, when modeling the formation of a small endocytic bud from a large cell, the deformation is small compared with the broader system. If we assume that the bulk of the cell is invariant with respect to bud formation, the computational burden can be reduced by modeling only the local deformation; we can assume that the modeled patch is attached to an implicit membrane reservoir. To define this coupled system, the constant area (A) and volume (V) of the reservoir must also be provided. The total area and volume of the broader system is given by A = Apatch + A, and V = Vpatch + Vr, where Apatch and Vpatch are area and “enclosed volume” of the mesh patch respectively. In our models, we enforce that all elements of a boundary loop are on the same plane; this way Vpatch can be unambiguously defined as the enclosed volume when each boundary loop is closed by a planar sheet. The capability to model systems attached to a reservoir reduces the modeled degrees of freedom while enabling intuitive physics to simplify the process of mimicking experimental conditions using Mem3DG.

Prescribing boundary conditions with force masking

Mem3DG supports modeling membranes with and without boundaries: a sphere (with no boundaries), a disk (with one boundary), and an open cylinder (with two boundaries). For systems without boundaries, the discrete forces conserve angular and translational momentum of system (as was also noted by Bian et al. (89)). Because the (discrete) potential energy is invariant under rigid body motions (i.e., the energy of the membrane is given only by the geometry), and the discrete forces are analytically derived from the energy, the discrete forces will not contribute to any rigid body motions since these components do not change system energy. In other words, the forces that lead to rigid body motion are, by construction, orthogonal to the shape derivative of the potential energy. To study systems with boundaries, Mem3DG currently supports three types of boundary conditions:

Roller, where the movement of boundary vertices is restricted along a given direction or plane.

Pinned, where the position of boundary vertices are pinned while the curvature is allowed to vary.

Fixed, where both the position and the boundary curvature are fixed for vertices on the boundary.

The different boundary conditions are achieved by masking the elements of the force matrix corresponding to the boundary vertices and their neighborhood. For example, to apply roller boundary conditions, we mask the Z-component of the force on the boundary vertices, therefore constraining their movement to the X-Y plane; pinned boundary conditions mask all force components for the boundary vertices to fix their position; fixed boundary conditions mask all force components for the outermost three layers to fix both their position and curvature.

Time integration and energy minimization

In this work, we use the forward Euler algorithm to integrate the system dynamics and the nonlinear conjugate gradient method to solve for equilibrium conditions. Both solvers are complemented by a backtracking line search algorithm, which satisfies Wolfe conditions to support adaptive time-stepping and robust minimization (121).

The forward Euler scheme was chosen as the simplest dynamical propagator; physically it represents over-damped conditions where the environment of the membrane is too viscous for the system to carry any inertia. Mathematically, the physics is described by, where ξ is the drag coefficient. From an optimization perspective, forward Euler is equivalent to the gradient descent method for minimizing an objective function, which is the discrete energy in our case.

A second propagator is the nonlinear conjugate gradient method for locally minimizing the discrete energy to yield the equilibrium shape of the membrane. Since the system is nonlinear, we periodically perform forward Euler (gradient descent) steps after several conjugate gradient steps. This approach of iterating between conjugate gradient and gradient descent steps is commonplace in the literature for solving nonlinear systems (121).

We note that other time integrators and energy minimizers are also compatible with Mem3DG. Included in the software are reference implementations of velocity Verlet integration (for symplectic time integration), and limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS, a quasi-Newton method to the equilibrium shape for large-scale problems where fast computation is needed). We do not discuss these additional solvers in this work.

Practical considerations for applying Mem3DG to biological problems

As we have noted above, in the DDG perspective, the mesh is the geometry and thus the formulation of the discrete forces and energies is exact. There are therefore very few restrictions on the resolution and quality of the input mesh. However, in biophysics, we often consider biological membranes as smooth systems. We expect that many users of Mem3DG may wish to approximate a smooth system using our discrete model. In doing so, they make an implicit assumption that such an approximation is reasonable. Although the relationships between geometric objects and their shapes are preserved between the smooth and discrete contexts, our ability to approximate a smooth problem with a discrete mesh is not guaranteed. Similar to finite differences and FEM, additional constraints on mesh quality and resolution must be imposed. To verify and understand the limitations of the assumption that the discrete mesh is the geometry and includes all of the geometric information, we numerically test the convergence of the discrete quantities under variation of resolution on an oblate spheroid mesh. The additional details regarding these numerical experiments are presented in Appendix D.

Setting the characteristic length scale of the finest mesh to be h = 1, as the mesh coarsens (i.e., mesh size increases) h increases. Fig. 3 shows the scaling relationship of the deviation in magnitude between the smooth and discrete quantities. Fig. 3 A shows the convergence property of global measurements that determines the energy (Eqs. 1 and 3), including the total area, A, enclosed volume, V, and total Gaussian curvature and mean curvature (squared), , and , respectively. Except for the total Gaussian curvature being an exact topological invariant, all integrated quantities exhibit approximately second-order convergence rate.

FIGURE 3

Comparison of discrete quantities with their smooth counterparts on spheroid shape. (A) Convergence plot of global quantities, including total area, volume, mean curvature (squared), and Gaussian curvature; and (B) Convergence plot of L1 norm of scalar and vector local quantities, including the mean curvature, Gaussian curvature, and the biharmonic term.

We acknowledge that convergence of global measurements does not imply that local measurements will also converge. To validate the convergence of local measurements, which determines the convergence of local forces on the membrane (e.g., Eqs. 15, 17, and 22), we utilize the L1 norm (Eq. 26) to evaluate the deviation of local quantities from their smooth counterparts. Fig. 3 B shows the local convergence plot. Similar to their global counterparts, local scalar mean and Gaussian curvature, ∫ H, and ∫ K, converge at . Fig. 3 B also shows the convergence of vector quantities, which not only contribute to the magnitude of the force but also set the direction of the force. The test shows that most vector quantities converge slightly slower than their scalar counterparts. Two terms exhibit poor convergence, the Schlafli vector term in Eq. 22, H ∫ S, and a scalar counterpart, ∫ ΔH. The latter term corresponds to the direct application of the cotangent Laplacian (Eq. 17) to the pointwise scalar mean curvature field; this approach is not used in our framework but is common in the literature (65). Both non-convergent expressions are discrete representations of the biharmonic term, ΔH, which have been noted to be sensitive to noises of vertex coordinates in the prior literature (100). Recall that the biharmonic term is the fourth-order derivative of the embedded coordinates. Although the traditional approximation theories suggest that higher-order derivatives often exhibit slower rates of convergence (122), to the best of our knowledge, there is not yet a rigorous study that connects DDG with an approximation theory. Nevertheless, we anticipate that similar principles hold. Two spatial plots comparing local measurements between smooth and discrete contexts are provided in the appendix (Fig. E.1 and Fig. E.2); each test is conducted using the finest mesh size (h = 1). Based on this numerical validation, we conclude that the computation of energy converges with a second-order rate (Fig. 3 A). While most components of the forces converge, the biharmonic term remains a limiting factor.

One other practical consideration for our models is that the Helfrich Hamiltonian, matching the in-plane fluidity of biological membranes, has no resistance to shearing. Without additional constraints, the mesh is susceptible to shearing motions, which can deteriorate mesh quality in some conditions (83). This can degrade the implicit assumption that the discrete mesh represents a smooth geometry. To ensure that such an approximation can remain valid throughout a trajectory, we have incorporated algorithmic solutions to adaptively maintain an isotropically well-resolved discrete geometry. This is achieved by two operations: 1) mesh regularization using local force constraints, which are common in FEM (78,82,83,85) (Appendix E.2); and 2) mesh mutations such as decimating, flipping, and collapsing edges. Beyond regularization, these local force constraints can also be used to model underlying physics of a problem of interest. For example, similar restoring forces between vertices (Eq. E.1) have been adopted to model actin-spectrin cortex in red blood cell (95). Mesh mutations are also a common practice to cope with deterioration and a means to perturb system configuration in other mesh simulations that use a Monte Carlo integration (60,89–95). The algorithms for mesh regularization and mutation are further described in Appendix E.

RESULTS AND DISCUSSION

To further validate the method and to provide a sense of how Mem3DG can be used and extended to solve more complex physics, we apply Mem3DG to a sequence of examples with increasing complexity. First, we model well-studied systems with homogeneous membrane conditions. We show that Mem3DG is capable of reproducing the classical solutions without imposing the axisymmetric constraint commonly adopted by other solvers. The later examples set a blueprint for extending and modifying Mem3DG for particular systems of interest. We introduce new energy and corresponding force terms to expand the formulation for complex systems of interest. We emphasize that the goal of these examples is to illustrate the generality of the theory and software and to outline specific steps for future extensions; we do not perform rigorous experimental comparisons, nor do we extract scientific insights. Additional care must be taken to mimic specific biological experiments for model validity, which is left for future work.

Each of the following sections considers a different class of membrane biophysics problem of increasing complexity in the coupling of the in-plane protein density parameter, ϕ ∈ [0, 1]. To mimic the various influences protein-lipid interactions can have on the bilayer, the protein density can be set to influence membrane properties such as the spontaneous curvature, , and bending rigidity, κ(ϕ). More complex phenomena such as the production of in-plane interfacial forces from membrane-protein phase separation (55,59,123) can also be modeled. In our final proof of concept, we extend Mem3DG to support full mechanochemical dynamics, where proteins can mobilize in and out of plane through adsorption and lateral diffusion, based on its coupling with membrane material properties and shape transformation. These scenarios highlight the relative ease of extending Mem3DG with additional physics and the potential utility to simulate realistic experimental scenarios. Note that, for all of the examples, unless otherwise specified, the bending rigidity of membrane, κ, is assumed to be the rigidity of a bare membrane, κ = 8.22 × 10−5 μm · nN. Despite the superior performance of the nonlinear conjugate gradient method in finding an energy minimizing configuration, to maintain both static and dynamic interpretability, we perform all simulations using a forward Euler integrator unless otherwise noted. All simulations presented in this work were conducted on a standard workstation with Intel Xeon processors. Although the numerical algorithms are amenable to parallelization, Mem3DG is currently a single-threaded program. Using a single core, the simulations here complete in minutes and up to 2 hours.

Modeling spherical and cylindrical membranes with homogeneous physical properties

We begin our examples by using Mem3DG to find the equilibrium shapes of membranes with homogeneous protein density, ϕ. We ask, given an initial membrane configuration with uniform bending modulus and spontaneous curvature, what are the minimizers of the system energy? The answers are the classical equilibrium solutions to the shape equation obtained analytically (42), and numerically using many methods with different assumptions (39,124). We will show solutions obtained using Mem3DG with topologies of sphere and tube (Fig. 4). These geometries were selected not only because of their potential for comparison with the legacy literature but also because they are reminiscent of biological membranous structures such as red blood cell (97,98,125,126), cell-cell tunneling and tethering (127–129), and neuron beading (130,131), among other biological processes.

FIGURE 4

Recover typical equilibrium shapes of membranes with homogeneous material properties. (A and B) Equilibrium solutions under different osmolarity and spontaneous curvature conditions, with initial condition of (A) Oblate spheroid and (B) Prolate spheroid. We vary the osmolarity by adjusting the concentration of the ambient solution, , while holding the enclosed amount of solute, n, constant. (C) Equilibrium solutions of a tubular membrane structure under variations in osmolarity and surface tension.

Starting with closed topological spheres, Fig. 4 A and B shows the typical equilibrium shapes under osmotic stress while the surface area is conserved. The preferred area of the vesicle, μm2, represents a sphere of radius 1 μm. This constraint is achieved by prescribing a large stretching modulus, K, such that the areal strain, , is less than 1%. The strength constant of osmotic pressure, K is set to be 0.1 μm · nN. Initializing the simulations from an oblate spheroid, as the osmolarity increases (e.g., the normalized ambient solution, ), we recover the well-known biconcave red blood cell shape (97,98,106,124) (Fig. 4 A). The vesicle adopts a more convex configuration as we increase the spontaneous curvature, indicating an overall increase in its mean curvature with the concomitant decrease of areas with negative mean curvature (the dimple regions). In contrast, starting from a prolate spheroid, as the spontaneous curvature increases, the vesicle adopts a dumbbell configuration as the energetically preferred state (Fig. 4 B). The size of the beads on the dumbbell is governed by the osmolarity, . These trends with respect to the variations of the spontaneous curvature and osmolarity are consistent with the analytical and numerical results in the broader literature (42,89). Qualitatively the predicted geometries of closed vesicles with homogeneous spontaneous curvature match the predictions of a detailed benchmark of mesh-based methods performed by Bian et al. (89).

We also modeled the shapes of membranes starting from an open cylinder configuration under different osmotic and surface tension conditions (Fig. 4 C). This problem is related to a well-studied phenomenon called the Plateau-Rayleigh instability (132,133). The Plateau-Rayleigh instability describes how surface tension breaks up a falling stream of fluid into liquid droplets. Compared with a liquid stream, a lipid membrane provides additional resistance against instability due to its rigidity. Zhongcan and Helfrich (134) obtain stability regimes as a function of membrane bending rigidity and spontaneous curvature using the spectral stability analysis (134). Although osmotic pressure is often reported as an important cause of morphological instability (131,135–137), the effect of osmotic pressure is difficult to isolate in wet experiments because the change to osmolarity affects the surface tension, which is a key driver of the instability. In our simulations, the two effects are decoupled, making the investigation of individual contributions to the morphology possible. All shapes in Fig. 4 C evolve from the initial tubular mesh with radius of 1 μm and axial length of 19.9 μm, under a constant spontaneous curvature of 1 μm−1. These simulations are set up as local models (cf., section “Defining properties of a membrane reservoir for systems with open boundaries”) where the explicit mesh is assumed to be coupled to a membrane reservoir. Additional geometric information defining the membrane reservoir and boundary conditions are required to initialize the local model. The tubular structure is considered to be a cylinder that connects two otherwise detached domains (e.g., membrane reservoirs), which combined have a total reservoir volume, V = 4.19 μm−3. The strength of osmotic pressure, K, is set to be 0.01 μm · nN. To isolate the effect of osmotic pressure and surface tension on the morphology, we prescribe a specific surface tension that we assume to be invariant with respect to changes to the surface area. On the two boundary loops of the mesh, we apply roller boundary conditions, which restrict the movement of boundary vertices in the axial direction. The length of the tube is thus constrained to be 19.9 μm, while the radius of the tube including the boundaries is free to shrink or expand.

As the osmolarity increases from the reference condition ( μm−3) (Video S1), we obtain near constant-mean-curvature surfaces such as unduloid pearl structure at μm−3 (Video S2), and cylindrical tube at μm−3, which follow the trends from both analytical (42,138) and experimental observations (19,131,135). As we increase the surface tension from the reference condition (λ = 1 × 10−7 nN · μm−1) to a tension-dominated regime (λ = 1 × 10−4 nN · μm−1), we obtain the beads-on-a-string structure that minimizes the surface-to-volume ratio (Video S3). The formation of beads-on-a-string is an interesting configuration that has been identified in biological membranes and other systems (130,131). Note that our simulations revealed a symmetric metastable state where two large beads form at either end (Appendix A), connected by a thin tube, prior to adopting the asymmetric conformation shown in Fig. 4 C. We believe that discretization artifacts such as mesh mutations act as a perturbation to break the symmetry of the metastable intermediate and transition the membrane to a single bead configuration (see Fig. A.1).

FIGURE A.1

A symmetric metastable state with two beads instead of a single larger bead is observed, prior to collapsing to the solution shown in Fig. 4 C, high tension.

These examples with uniform spontaneous curvature profile prove the ability of Mem3DG to reproduce the expected classical solutions for spherical and tubular membrane geometries. Note that no axisymmetric constraint is imposed in these simulations. Mem3DG solves the system in full three dimensions and the symmetrical configurations are due to the problem physics. The ability to adapt to changing and complex curvatures of the membrane using discrete mesh is achieved using mesh mutation and other manipulations within solver steps. For example, the pinched neck regions of the tubes are automatically decimated with finer triangles than other regions of the mesh. For a global closed membrane simulation such as in Fig. 4 A, B, we do not remove any rigid body motions from the system; since the forces from DDG are exact and we used the forward Euler integrator, no artificial rigid body motions are introduced throughout the simulation. These examples show that that the derivation of the discrete energy and forces along with the software implementation are accurate and proceed to test Mem3DG with more complex examples.

Modeling endocytic budding mechanisms

Our goal is to highlight the potential of Mem3DG and its associated framework for investigating mechanical phenomena relevant to cellular biology. Endocytosis is a cellular process in which cells uptake cargo from the extracellular environment; the transported material is engulfed by the cell membrane, which then buds off to form a vesicle (13). Endocytosis occurs through various mechanisms, including clathrin-mediated endocytosis (13,139). It has been shown that clathrin aggregates on the plasma membrane, helping to deform the membrane and form a spherical bud (9,13,59). However, the specific mechanisms of how membrane-clathrin interactions facilitate membrane curvature generation remain unresolved. While there is significant literature investigating the many proposed mechanisms, here we develop models to demonstrate the bud formation via spatially localized spontaneous curvature, combined with a line tension term arising from phase separations on the membrane (140).

We model endocytic budding on a circular patch with radius 1 μm (a disc with one boundary loop). We assume that the patch is a local system which is coupled to a large vesicle (section “Defining properties of a membrane reservoir for systems with open boundaries”). A heterogeneous protein density, ϕ ∈ [0, 1], is applied to mimic the distribution of clathrin and other scaffolding proteins. Shown in Fig. 5 A, the protein density is high (ϕ = 1) toward the center of a geodesic disk with radius 0.5 μm) and decreases toward the boundaries (ϕ = 0). During simulation, the geodesic distance to the center of the patch is periodically computed using the heat method (141). Vertexwise ϕ is assigned based on the stair-step profile smoothed by the hyperbolic tangent function applied to the geodesic distance. Each experiment is initialized as a flat patch and the displacement of boundary vertices is restricted using a fixed boundary condition. Since the patch is viewed as a small piece within a larger closed vesicle reservoir, we assume that the surface tension is constant.

FIGURE 5

Budding dynamics by robust mechanisms of protein scaffolding and interfacial line tension constriction. (A–C) Control group, (D–F) bending-driven scaffolding mechanism, and (G–I) Interfacial line tension assisted budding. (A) Spontaneous curvature distribution, , on initially flat patch. (D) Normal projection of the bending force at T = 15. (G) Normal projection of the line tension force at T = 7. (B, E, H) Shape evolution through time-series snapshots of the Y-Z cross sections of the membrane, corresponding to the vertical dash lines in (C), (F), (I) trajectory plots of system energy and its competing components.

A common model to account for the preferential bending owing to protein-membrane interactions is through the spontaneous curvature; we assume , where μm−1 is the spontaneous curvature imposed by the membrane protein coat. Proteins such as clathrin are known to form stiff scaffolds on the membrane. Similar to the spontaneous curvature, we can assume a linear relationship between bending rigidity and protein density, κ(ϕ) = κ + κ ϕ, where constant κ is the rigidity of the bare membrane, and κ is additional rigidity of the protein scaffold.

Shown in Fig. 5 A–C and Video S4, is the control simulation where we set the contribution to the rigidity from protein to be the same as that of the raw membrane, κ = κ. Fig. 5 A shows the initial flat configuration of the control experiment; the color bar shows the heterogeneous spontaneous curvature resulting from the prescribed protein density profile. In the control experiment, the bending force is resisted by the surface tension (Fig. 5 C) until, at the final frame in Fig. 5 B (t = 5), the membrane reaches the equilibrium configuration where the surface tension cancels with the bending force. In a second model, we assume that the scaffolding proteins are much more rigid than the bare membrane, κ = 3κ. Fig. 5 D–F and Video S5 show the bud formation due to this increased protein scaffolding effect. The greater rigidity results in an increase of initial bending energy, which outcompetes the resistance from the surface tension (Fig. 5 F). Fig. 5 E shows the shape evolution from a flat patch to a successful bud with a pinched neck. Fig. 5 D shows the signed projection of the bending force onto the vertex normal, , at T = 15. (Outward-pointing angle-weighted normal; the same applies to the interfacial line tension.) We can see an “effective line tension” driven by the heterogeneous spontaneous curvature that constricts the neck. This phenomenon is theoretically explored in detail by Alimohamadi et al. (58).

For our third model, based on the prior observations that protein phase separation on surfaces can lead to a line tension (140), we incorporate a Ginzburg-Landau interfacial energy into the system, where η, referred to as the Dirichlet energy constant, governs the strength of the energy penalty, and is the discrete surface gradient of the protein density profile. The term is similar to previous modeling efforts by Elliott and Stinner (80) and Ma and Klug (79) using FEM; because we use the protein phase separation as a prior, we exclude the double-well term, which models the thermodynamics of phase separation, and incorporate only the Dirichlet energy component that penalizes the heterogeneity of membrane composition.

Defined as the slope of the linearly interpolation of ϕ on faces of the mesh, f, the discrete surface gradient of the protein density is, where following illustration in Fig. 1 C, is the vector aligned with the halfedge , with its length of l, and (·)⊥ represents a 90° counterclockwise rotation in the plane of f. The resulting line tension force is then the shape derivative of the Dirichlet energy, , which acts to minimize the region with sharp heterogeneity. The detailed derivation of the shape derivative is elaborated in Appendix C.3, where we follow the formulaic approach by taking geometric derivatives of basic mesh primitives shown in Eq. C.13. Note that despite bearing the same name, this line tension force differs from from those resulting from line energy, which prescribes the line tension energy based on the interfacial edge length (55,142). The line tension force from Dirichlet energy is used to model the out-of-plane component resulting from either entropic or enthalpic repulsion at the interface between heterogeneous membrane protein aggregates. The Dirichlet energy is based on a 2D field variable and the line tension is only effective when phase separation of the field variable occurs, where the interfacial line, or, more precisely, thin areal band, between phases exists. With the introduction of protein evolution later in sections “Protein aggregation on the realistic mesh of a dendritic spine” and “Membrane dynamics with full mechanochemical feedback,” thickness depends on the competition between the Dirichlet energy and other competing aggregational potential.

Fig. 5 G–I and Video S6 show the trajectory where we used control bending rigidity, κ = κ, and the additional interfacial line tension, η = 5 × 10−4 μm · nN. We find that the interfacial line tension, jointly with the bending force, lowers the system energy and helps the formation of a spherical bud (Fig. 5 I and H). Fig. 5 G shows the snapshot (t = 7) with the color map representing the signed normal projection of the interfacial line tension that acts to constrict the neck. These examples of endocytic bud formation help to illustrate the utility of Mem3DG and the accompanying theoretical framework. Since physical parameters are assigned on a per-vertex basis, it is straightforward to incorporate heterogeneity such as the nonuniform bending rigidity and spontaneous curvature. In smooth theory and its derived discrete mesh models, when the membrane is heterogeneous, it is required to decompose the force separately in normal and tangential direction (40,65). In contrast, the general derivation of the discrete bending force following the formalism of DDG permits modeling membrane with heterogeneous material properties without any modification to its formulation (section “Forces from bending”). The introduction of Dirichlet energy and line tension force serves to highlight the relative ease to extend the modeled physics.

Protein aggregation on the realistic mesh of a dendritic spine

While the prior examples have focused on the mechanical response of the membrane given a bound protein distribution, we can also model the inverse problem. Given the membrane shape, how do curvature-sensing proteins diffuse in the plane of the membrane and distribute over the domain? And how does the resultant protein distribution influence the stresses of the system? To model the protein dynamics, we use three terms corresponding to protein binding, curvature sensitivity, and lateral diffusion.

To model the binding of proteins to the membrane, we assume that the energy of adsorption, ε, is constant and uniform across the surface such that the discrete adsorption energy is, where ϕ is an order parameter representing the area density of protein at each vertex. Taking the derivative with respect to ϕ, referred to as the chemical derivative, we obtain the adsorption component of the chemical potential. To account for protein curvature sensitivity, we find the chemical potential of the bending energy, where we assume that , and where κ and are constant parameters defined in Section “Modeling endocytic budding mechanisms.” The first term of Eq. 32 endows the protein with curvature-sensitive binding. The second term of Eq. 32 is the shape mismatch penalty; considering the binding of a rigid protein that induces a significant spontaneous curvature change, if the curvature of membrane is far from this new spontaneous curvature, then the shape mismatch between the membrane and proteins will prevent binding. Alternatively, if the protein is more flexible, a shape mismatch results in a small energetic penalty.

The in-plane diffusion of the protein is accounted for by the chemical derivative of the smoothing Dirichlet energy, E, where η is the same Dirichlet energy constant introduced in Eq. 28 that governs the strength of interfacial line tension,. The total chemical potential captures the bending, adsorption and diffusion components. A mobility rate constant, B, determines the time scale of the chemical response,

We investigate the influence of curvature-dependent binding to a realistic dendritic spine geometry, which was reconstructed from electron micrographs and curated using GAMer 2 (Fig. 6 A) (32). A summary of the parameters used in the simulation is shown in Table 3. The mean curvature of the spine geometry is shown Fig. 6 C. We isolate the effect of curvature-dependent binding by assuming that the shape of the spine is fixed and impose Dirichlet boundary conditions at the base on the spine to fix the protein concentration, ϕ = 0.1 (Fig. 6 A).

FIGURE 6

Protein aggregation on a realistic dendritic spine geometry. (A) Mesh of the dendritic spine and its boundary elements. (B) Trajectory of protein evolution and components of system energy. (C) Mean curvature of the geometry. The normal component of (D) the bending force at t = 0, (E) the line tension force produced by the equilibrium protein distribution, and (F) the difference in the bending force profile produced by final protein distribution as opposed to the initial distribution.

TABLE 3

Parameters used in section “Protein aggregation on the realistic mesh of a dendritic spine”

Parameters	Values
ϕ 0	0.1
κc	0 nN · μm
H¯_c	10 μm−1
B	3 nN−1 · μm−1 · s−1
η	0.01 μm · nN

Starting from a homogeneous protein distribution, ϕ0 = 0.1, Fig. 6 B and Video S7 show the evolution of the protein distribution and a trajectory of the system energy. Note that, for simplicity, we have turned off the adsorption energy term since it only shifts the basal protein-membrane interactions, which will also be set by the Dirichlet boundary condition. Mem3DG constrains the range of ϕ ∈ (0, 1) using the interior point method (121). Due to the curvature sensitivity of the protein, illustrated by the snapshots (Fig. 6 B, T = 350) representing the final protein distribution, the protein aggregates toward regions of high curvature (e.g., on the spine head).

Although the proteins can reduce the bending energy by modulating the local bending modulus and spontaneous curvature, the protein distribution at equilibrium does not cancel out the bending energy. We expect that the Dirichlet energy term, which limits ϕ to be smooth across the geometry, restricts the protein from further aggregating to the extent required to cancel out the bending energy. The components of forces on the initial and final configurations of the spine are compared in Fig. 6 D–F. The initial homogeneous protein distribution has no line tension forces and a bending force shown in Fig. 6 D. After the protein distribution reaches the steady state, line tension appears in response to membrane heterogeneity Fig. 6 E. We hypothesize that, similar to section “Modeling endocytic budding mechanisms,” the line tension constricts the neck of the spine and helps to support the cup-like structures in the spine head. While, in most regions, the distribution of proteins reduces the force, several regions experience increased stress Fig. 6 F. Note that the magnitude of the forces generated by proteins in this model is orders of magnitude smaller than the initial bending force. Thus, this example demonstrates that Mem3DG can be used on meshes imported from realistic geometries of cellular substructures.

Membrane dynamics with full mechanochemical feedback

In this section, we will demonstrate the use of Mem3DG to model the complete mechanochemical feedback of a protein-membrane system. For the following simulations, not only can proteins bind in a curvature-dependent manner but the membrane can also deform, leading to a feedback loop. We have introduced all of the force terms in previous sections except the shape derivative of the adsorption energy, which accounts for the change in the area of protein coverage (i.e., expanded coverage if ε < 0).

Revisiting the claim that all discrete forcing is exact with respect to the discrete energy, we validate by examining the convergence of the forcing terms with respect to the size of perturbation to the system configuration, ϵ (Fig. 7 A). This is based on the leading order expansion done in Eq. 24, which concludes that the forcing terms are exact if their rate of convergence is at least second order. Shown in Fig. 7 A, this is true for all forcing terms; note that, since the adsorption energy, E, is a linear function with respect to ϕ, μ can be determined to the machine precision for all perturbation sizes. A meaningful discrete-smooth comparison of all terms in the energy and forcing similar to section “Practical considerations for applying ” requires the analytical solutions of the bending force and interfacial line tension arising from the spatially heterogeneous membrane properties, which, to the best of our knowledge, are not available. In section “Modeling endocytic budding mechanisms,” we introduced a heterogeneous membrane composition as a static property. By prescribing the protein density profile, we can get hints about how to form membrane buds from a patch. Here we lift this assumption and simulate the dynamics of osmotic pressure-driven budding from a vesicle. The dynamics couples the proteinmembrane mechanochemical feedback and includes protein binding and diffusion introduced in section “Protein aggregation on the realistic mesh of a dendritic spine.” The expressions of discrete free energy and forcings do not change but we allow the membrane configuration and protein density to interact and evolve simultaneously.

FIGURE 7

Modeling budding from a vesicle driven by the full mechanochemical feedback of membrane-protein interactions. (A) Validation of the exactness of the discrete forcing with respect to the discrete energy. The terms correspond to forces from bending f, tension area f, pressure-volume f, Dirichlet f, and protein binding f. μ, μ, and μ are the chemical potential of diffusion, bending, and binding respectively. (B) The time trajectory of budding dynamics in hypertonic, isotonic, and hypotonic osmotic conditions. (C) The final snapshot of the system configuration under hypertonic, isotonic, and hypotonic conditions. (D) Similar geometries to those shown in (C) have been observed in experiments by Saleem et al. (59).

We start each simulation from a sphere with a uniform protein concentration, ϕ = ϕ0 = 0.1. We consider the evolution of the system in different osmotic conditions: hyper-, iso-, and hypotonic, , 3.95 and 4.99 μm3, respectively. Additional parameters for these simulations are given in Table 4. Fig. 7 B shows plots of the mechanical, , and chemical response, , along with the protein density, (ϕmax + ϕmin)/2, over the trajectory for each osmotic condition. Note that under hypo- and isotonic conditions, the system reaches the (quasi) steady state where both the shape and protein distribution stabilize, while, in a hypertonic solution, the system continues to experience strong mechanical force and protein mobility, which we expect to drive further morphological changes of the membrane beyond our simulation stopping point. Fig. 7 C shows the final snapshot of each simulation across the osmotic conditions with the protein density represented by the color map. In hypertonic conditions, the osmotic pressure provides sufficient perturbations to membrane morphology, which initializes the positive feedback loop between membrane curvature generation and protein aggregation. This mechanochemical feedback jointly promotes the outward bending of the membrane and results in bud formation (Fig. 7 C, hypertonic; Video S8). Under isotonic and hypotonic conditions, the osmolarity does not permit the large change in the volume required to form spherical buds with a thin neck. In hypotonic condition, the pressure-tension balance provides substantial stability to the initial spherical configuration. In comparison, in the isotonic condition, the comparable competition between the chemical and mechanical response leads to a patterned protein distribution and an undulating morphology (Fig. 7 C, hypotonic; Video S9). This example illustrates the possibility to utilize Mem3DG to model a full mechanochemical feedback between membrane and protein. Although we do not intend to replicate the exact experimental conditions and assumptions, the geometries obtained from these simulations resemble the shapes obtained by Saleem et al. (59) who investigated budding from spherical vesicles under differing osmotic conditions (Fig. 7 D) (59).

TABLE 4

Parameters used in section “Membrane dynamics with full mechanochemical feedback” for models with full mechanochemical feedback

Parameters	Values
ϕ 0	0.1
κc	8.22 × 10−5 μm · nN
H¯_c	10 μm−1
KV	0.5 nN · μm
KA	1 nN · μm−1
B	3 nN−1 · μm−1 · s−1
ξ	1 nN · s · μm−1
ε	−1 × 10−3 nN · μm
η	0.1 μm · nN
V¯	2.91, 3.95, and 4.99 μm3

SUMMARY

In this work, we introduce a new perspective for constructing a 3D membrane mechanics model on discrete meshes. The goal of our approach is to close the gap between existing discrete mesh-based models (60,86–96,99,102,103) and the smooth theory. Specifically, we seek to advance the discussion behind the choice of algorithmic approaches for computing geometric values required for the discrete energy and force (65,89,99,100). We start by writing a discrete energy, Eq. 3, mirroring the spontaneous curvature model. Then using the perspective of DDG, we show that there is a formulaic approach for deriving the corresponding discrete force terms based on fundamental geometric vectors. By identifying geometric invariants and grouping terms, the resulting discrete forces have exact correspondence to the traditional smooth theory. This helps us to facilitate the comparison between smooth and discrete contexts enabling new geometric perspectives and understanding of numerical accuracy. Moreover, the discrete force terms are functions of readily accessible geometric primitives, and the formulation is amenable for efficient implementation and extension.

We have produced a reference software implementation called Mem3DG. Using Mem3DG, we validate our theory by reproducing the solutions to the classical shape transformations of a spherical and tubular vesicle. We further demonstrate the derivation and incorporation of additional physics terms to model protein-membrane interactions. Following our formulaic approach using DDG, we derived the discrete analog of various physics, such as the interfacial line tension, surface-bulk adsorption, protein lateral diffusion, and curvature-dependent protein aggregation. To exemplify all the introduced physics, the full mechanochemical coupling between the membrane shape and protein density evolution results in protein localization, pattern formation, and budding. These examples serve to highlight the extensibility of the framework, which does not require the introduction of coordinates and advanced tensor calculus commonly needed to solve problems on arbitrary manifolds. The software implementation Mem3DG was designed to facilitate coordination between theoretical modeling and wet experiments; we hope to support the modeling of scenes with well-resolved protein-membrane interactions such as in the electron tomograms (143). We expect that as the advances in biophysical modeling and membrane ultrastructure imaging progresses, Mem3DG will emerge as a useful tool to test new hypotheses and understand cellular mechanobiology.