Biomembranes undergo complex, non-axisymmetric deformations governed by Kirchhoff-Love kinematics and revealed by a three dimensional computational framework

Biomembranes play a central role in various phenomena like locomotion of cells, cell-cell interactions, packaging of nutrients, and in maintaining organelle morphology and functionality. During these processes, the membranes undergo significant morphological changes through deformation, scission, and fusion. Modeling the underlying mechanics of such morphological changes has traditionally relied on reduced order axisymmetric representations of membrane geometry and deformation. Axisymmetric representations, while robust and extensively deployed, suffer from their inability to model symmetry breaking deformations and structural bifurcations. To address this limitation, a 3D computational mechanics framework for high fidelity modeling of biomembrane deformation is presented. The proposed framework brings together Kirchhoff-Love thin-shell kinematics, Helfrich-energy based mechanics, and state-of-the-art numerical techniques for modeling deformation of surface geometries. Lipid bilayers are represented as spline-based surfaces immersed in a 3D space; this enables modeling of a wide spectrum of membrane geometries, boundary conditions, and deformations that are physically admissible in a 3D space. The mathematical basis of the framework and its numerical machinery are presented, and their utility is demonstrated by modeling 3 classical, yet non-trivial, membrane problems: formation of tubular shapes and their lateral constriction, Piezo1-induced membrane footprint generation and gating response, and the budding of membranes by protein coats during endocytosis. For each problem, the full 3D membrane deformation is captured, potential symmetry-breaking deformation paths identified, and various case studies of boundary and load conditions are presented. Using the endocytic vesicle budding as a case study, we also present a"phase diagram"for its symmetric and broken-symmetry states.


Introduction
Membrane curvature is ubiquitous in biology [1]. Indeed, the bending of cell membranes is a central aspect of function for cells and organelles in many cellular processes such as cell migration [2], cell membrane repair [3], membrane trafficking [4] and cytokinesis [5], as well as the maintenance of distinctive membrane shapes within internal organelles like the endoplasmic reticulum [6,7] and the Golgi complex [8]. Some important curved structures include tubules, sheets, vesicles and cisternae [9]. A number of mechanisms have been identified to influence membrane bending, including geometric confinement by protein or lipid components of the membrane (intrinsic factors) [10,11] and peripheral proteins and the cytoskeleton (extrinsic factors) [12,13]. These mechanisms are often coupled and are spatio-temporally regulated by biochemical signaling cascades, leading to the mechanochemical coupling of signaling and membrane deformations. Lipid bilayer models that assume an in-plane fluid-like behaviour and an out-of-plane solid-like behaviour have provided notable insight to investigations of such curvature generation mechanisms.
Despite the wealth of information provided by theoretical membrane mechanics models, an important restriction in several of these studies is the assumption of various degrees of symmetry for the membrane geometry and its deformation. Indeed, the computation of membrane bending phenomena is significantly simplified with the axisymmetric assumption, but as we have shown recently [21], this may come at the cost of generality and precision in identifying the underlying physics, as lower-energy, low-symmetry kinematic modes and even entire mechanisms may be overlooked. With growing interest in curvature-mediated biophysical phenomena and in 3D imaging and reconstruction methods [25,26], there is a need for general purpose computational tools to enable fully three dimensional numerical simulations.
The continuum mechanical treatment of solids considers deformation as a mapping of the geometry (3D volume, 2D surface, or 1D curve) from its reference, undeformed configuration to a deformed current configuration under the influence of internal or external loads, of which the latter also may appear as boundary conditions. In limited cases, the geometry, loads and boundary conditions result in a mathematical problem of deformation of a k-manifold immersed in a n-dimensional space (R n ). A 3-manifold is a volume, 2-manifold is a surface and 1-manifold is a curve. For k = n, modeling solid deformation is relatively straightforward and can be accomplished in the framework of Euclidean geometry using a rectilinear coordinate basis. However, deformation of shell-like surface geometries, as is the case with biological membranes, involves tracking the underlying kinematics and evolution of geometric configurations of a 2-manifold embedded in a 3-dimensional space [27]. Such a geometric embedding demands a non-Euclidean framework with a curvilinear coordinate basis. While the mathematical treatment of such a framework is well-developed (beginning with the celebrated work on differential geometry by Riemann in the 19 th century [28]), its application to three-dimensional modeling of biomembranes, which entails solving nonlinear partial differential equations in a curvilinear coordinate basis is relatively recent. Beginning with finite element models of Mindlin-Reissner plates [29][30][31][32] and Kirchhoff-Love shells [29,[33][34][35], initial efforts focused on developing numerical models in a rectilinear coordinate basis with approximated geometries and kinematics. However, the advent of spline-based geometric representations of surfaces and the more recent development of Isogeometric Analysis (IGA) 2/24 techniques [36] allow for an exact representation of surface geometries and the use of a curvilinear coordinate basis.
Such treatments are now gaining traction in modeling structural applications [37][38][39][40] and also in the context of biological materials [41][42][43][44]. We build upon these developments, especially from Sauer et al. [42], by adopting spline-based representations of surface geometries, treatments of membrane kinematics using a curvilinear basis, and the framework of IGA to develop a comprehensive computational modeling framework for studying complex deformations in biological membranes.
In this work, we present a three-dimensional, Helfrich-energy based, Kirchhoff-Love thin-shell computational framework for modeling the deformation of biological membranes in the regime of fully nonlinear kinematics and precise geometric representations. With this treatment, we are able to model membrane deformations, resolve geometric bifurcations, and explore post-bifurcation responses. The main ingredients of this framework are the classical Kirchhoff-Love membrane kinematics [27], a variational formulation of the governing equations underlying Helfrich-energy based mechanics on surface manifolds [42,44,45], and the numerical framework of IGA for solving the underlying partial differential equations. IGA methods form a numerical framework for finding approximate solutions to general partial differential equations [36], are a generalization of the classical Finite Element Method [46][47][48], and posses good numerical approximation and stability properties [49]. Crucially for accurate modeling of membrane biophysics, since IGA uses spline basis functions to represent the geometry and its deformation, it admits the continuity of slopes that is a characteristic of membranes in all states except for those of actual scission.
As a result, we can now investigate simulations of membrane deformation under conditions that are notably more general (having fewer restrictive kinematic assumptions) than those considered previously in the literature [20,[50][51][52][53][54].
The computational framework is implemented as an open-source software library and provided as a resource to the biophysics community [55].
To demonstrate the scope of the computational framework, we simulate three classical and non trivial membrane deformation phenomena (see Figure 1): (a) formation of tubular shapes and their lateral constriction, (b) Piezo1-induced membrane footprint generation and gating response, and (c) the budding of membranes due to the spontaneous curvature of the protein coats during endocytosis. For each case, three dimensional membrane deformation is tracked, symmetry-breaking deformation pathways identified, and a few case studies of boundary conditions and loading are presented to exhibit the fidelity and modeling potential of the proposed methodology. We also present a phase diagram of symmetric and broken-symmetry states of membrane budding during endocytosis.
In the following sections, we present an outline of the mathematical framework and the model development, followed by a presentation of the three boundary values problems considered, their modeling results and biophysical implications. Finally, a discussion of the framework and its utility is presented in the conclusion.

Materials and Methods
The mathematical framework consists of surface geometry parametrization, Kirchhoff-Love membrane kinematics, Helfrich-energy based mechanics of lipid bilayers and surface partial differential equations governing mechanical deformation. Key ingredients of this framework are described in the mathematical framework subsection below, while the detailed mathematical derivations are provided in the SI. Using the IGA apparatus, the mathematical treatment is then cast into a numerical formulation that allows for solving the governing equations to obtain the spatial evolution of 3/24

Mathematical framework
The mathematical treatment introduced here follows from Sauer et al. [42]. Only the important results are summarized in this section, and the detailed derivations are presented in the SI.

Surface parametrization and kinematics
Consider a lipid bilayer represented as a surface (2-manifold) embedded in a 3D volume, as shown in Fig. 3. Let the reference (undeformed) configuration and the current (deformed) configuration of the surface geometry be denoted by Ω 0 and Ω, respectively. The configurations Ω 0 and Ω are parametrized by the coordinates ξ 1 and ξ 2 that map a flat 2D domain to the surface coordinates X and x: The (covariant) tangent vectors in the reference and current configuration are given by: 4/24 Here, u x , u y and u z are the displacement components, t is the surface traction and t y its component along the y-axis,ñ is the normal to the boundary curve, φ is the boundary slope, γ is the surface tension applied on the membrane boundary, and h 0 is the instantaneous curvature. Blue and orange colors identify the outer and inner rims, respectively.

5/24
In the expressions that follow, except when indicated otherwise, uppercase letters are associated with the reference configuration and lowercase letters are associated with the current configuration.
Using the tangent vectors we define the surface normals as follows: From the triad consisting of the tangent vectors and the normal that form the local curvilinear coordinate basis, we can obtain expressions for the metric tensor, The second order derivatives of the surface coordinates X and x are given by: and from them we obtain the components of the curvature tensor, We are now able to define the primary kinematic metrics of interest: the mean and Gaussian curvature. The mean curvature and Gaussian curvature are frame invariant measures of a surface geometry, and hence are natural choice for representing the kinematics of the surface as it deforms. Using the components of the curvature tensor, we can obtain expressions for the mean curvature, and the Gaussian curvature,

Biophysics of membrane deformation
With a focus on representing the correct deformation, a biomembrane is often modeled as a thin elastic shell governed by the classical Helfrich formulation [14,56,57] of membrane bending energy. In this treatment, the primary kinematic variables are the curvatures capturing the bending of the membrane, and the elastic energy density of the membrane is given by: where k B and k G are the bending modulus and the Gaussian curvature modulus of the membrane, and h 0 represents the instantaneous curvature induced in the membrane.
Furthermore, we assume that the membrane is area preserving (i.e the membrane area is constant) [58] -a constraint that is implemented using a Lagrange multiplier field. Enforcing the area-preserving condition results in the following expression for the elastic energy density: where λ is the point value of the Lagrange multiplier field, and J is the surface Jacobian (ratio of area in the current configuration to the area in the reference configuration). Thus, while the Helfrich energy is defined entirely in terms of 6/24 the geometry of the surface, the Lagrange multiplier, often interpreted as the membrane tension [45,59]; i.e. a non-geometric quantity, also plays a role in determining the minimum energy configuration. We ignore any fluid [60,61] and friction [62][63][64] properties of the bilayer, guided by the dominance of unstable and stable equilibrium states over relaxation or rate processes. The augmented Helfrich Hamiltonian whose extremum is sought over the membrane surface, including the Lagrange multiplier field λ is given as: where Ω is the domain of integration over the membrane surface.

Governing equations
The governing equation for quasi-static mechanical equilibrium in 3D simulations is obtained as the Euler-Lagrange condition of the Helfrich energy functional following standard variational arguments, and is given by: where ∂Ω is the membrane boundary on which surface tractions and displacement boundary conditions can be applied, as shown in Figure 3. Furthermore, δa ij and δb ij are variations of the components of the metric tensor and the curvature tensor, respectively, Here, σ ij are the components of the stress tensor, M ij are components of the moment tensor (in the current configuration), p is the pressure applied on the membrane surface (in the case of the tube constriction boundary value problem), and t is the surface traction.
For a hyperelastic material model, we can express the stress and moment components in terms of the strain energy density as [44]: For the Helfrich type strain energy density, these take the form:

Computational implementation
The governing equations of membrane deformation, given as a system of nonlinear elliptic partial differential equations, are solved using the framework of Isogeometric Analysis [36]. An IGA method-based membrane mechanics framework has been developed for this work, and is built on top of the PetIGA [65] open source library. In an IGA approach, the membrane geometry is discretized using a spline mesh and the governing equations are converted to a nonlinear system of equations, which is then solved to obtain the deformed membrane shape. Of importance to our central result is that this framework naturally admits both symmetric and asymmetric deformation modes driven by the underlying physics.
We make three key important remarks about this framework. First, a fundamental conjecture of the Helfrich model is that the characteristic length scales of the problem are much larger than the thickness of the bilayer [14]. This assumption allows us to neglect the effect of transverse shear deformations and consider the classical Kirchhoff-Love shell kinematics for thin shell geometries [27]. Second, numerical solutions to the membrane shape equations in general coordinates are challenging because of continuity requirements in the numerical scheme. Specifically, the deformation must maintain first-order; i.e., C 1 , spatial continuity of the membrane surface in all except extreme states, such as the configuration of actual scission. We overcome this challenge by adopting both spline basis functions, which allow high-order continuity, and the numerical framework of Isogeometric Analysis [36]. Finally, an inherent limitation of the Helfrich energy formulation in three dimensional simulations is the lack of resistance to shear deformation modes. The zero energy modes corresponding to shear deformation are eliminated in this framework by adding shear stabilization terms of smaller magnitude relative to the traditional bending terms in the Helfrich energy [42], thus restoring stability to the numerical model.

Results
We demonstrate the simulation framework using three classical membrane deformation problems: formation of tubular shapes and their lateral constriction, Piezo1-induced membrane footprint generation and gating response, and the budding of membranes by protein coats during endocytosis. Through these examples, we also demonstrate the emergence of increasingly complex membrane deformations that are beyond the scope of traditional axisymmetric formulations. These problems are described in detail below.

Formation of tubular shapes and their lateral constriction
Many cell organelles and cytoplasmic projections are shaped as vesicles, tubes, or elongated membrane structures.
Some examples of such shapes are the filopodia protrusions, inner mitochondrial region, endoplasmic reticulum, the Golgi complex, etc. (see Figure 1). These tubular structures play an important role in the locomotion of cells, production and folding of proteins, and in the formation of vesicles for transporting proteins and lipids among others. A typical mechanism for producing these tubular shapes involves motor proteins that attach to the cell membrane and pull it along the filaments of the cytoskeleton [66,67]. Further, as is the case with the fission of endocytic vesicles, the tubular or vesicular structures also undergo constriction by scission proteins like dynamin [68][69][70]. This constriction mediates a membrane pinch-off mechanism that leads to the formation of vesicles. From a biophysical standpoint, it is important to gain a quantitative understanding of the interaction between the proteins and the membranes by determining the deformation mechanisms, forces exerted by proteins, and kinematic constraints.
A classic benchmark problem in the understanding of elongated biomembrane structures is the analytical model of the formation and interaction of membrane tubes proposed by Derényi et al [17]. Some key results of this model are the prediction of the magnitude of protein-membrane interaction forces and tubule radius, and their dependence on the membrane bending modulus (κ B ) and surface tension (γ). The protein pulling force, t y , and the tubule radius, r, are related to the bending modulus and surface tension of the membrane as follows: t y ∝ √ κ B γ and r ∝ κ B /γ. In addition to these analytical estimates, numerical solutions to the problem of membrane tube pulling, albeit with axisymmetric constraints on deformation, are available in the literature [71,72] and in our earlier work [21]. We take advantage of the analytical estimates proposed by Derényi et al., the numerical solutions available from axisymmetric models [21], and validate the computational framework proposed in this work by comparing the load-displacement response of membrane tube pulling from these three approaches.
The boundary value problem solved, along with the spatial discretization (mesh), boundary conditions on the displacement (u) and the membrane boundary slope (φ) are shown in Figure 2 We further consider the effect of lateral constriction pressure on the tubular geometry and demonstrate the non-axisymmetric pinching deformation profile that is predicted by the computational framework. For this boundary value problem, we consider a tubular geometry (shown in Figure 1(a), under tube constriction) and apply an axisymmetric constriction pressure that would be applied by a spiral collar protein like dynamin [57,73,74]. As can be expected, an axisymmetric model would predict an axisymmetric pinching profile in the vicinity of the constriction pressure [21]. However, the fully 3D model considered in this computational framework is not limited to axisymmetric solutions, and is thus able to predict non-axisymmetric states when they are the energy minimizing solutions to the governing equations of membrane deformation. The progression of the non-axisymmetric solution with increasing 9/24 constriction pressure is shown in Figure 5. This shape of the membrane has significant implication on the force and energy barrier of protein induced pinching of membranes, as has been studied in detail in our recent work demonstrating how non-axisymmetric buckling lowers the energy barrier associated with membrane neck constriction [21]. In that study, we used an earlier version of the computational framework proposed here to study the influence of location, symmetry constraints, and helical forces on membrane neck constriction in a lipid bilayer.
Simulations from our model demonstrated that the energy barriers associated with constriction of a membrane neck are location-dependent, and if symmetry restrictions are relaxed, the energy barrier for constriction is dramatically lowered and the membrane buckles at lower values of the constriction pressure. These studies helped establish that even though there exist different molecular mechanisms of neck formation in cells, the mechanics of constriction of a cylindrical membrane tubule naturally leads to a loss of symmetry that can lower the energy barrier to constriction. This loss of symmetry may be a common mechanism for different scission processes and clearly demonstrates the need for fully three dimensional computational framework to give predictive insights into membrane deformation.

Piezo1-induced membrane footprint and gating response
We next investigate how mechanosensitive channels can deform the membrane. Mechanosensitive ion channels on the cell membrane play an important role in the mechanosensory transduction processes of the cell. These ion channels are sensitive to the forces acting on the cell membrane and respond to these forces by undergoing conformational changes.
These changes result in the opening and closing of pores in the cell membrane and thereby regulate the flow of ions and solutes entering and egressing the cell. Examples of such mechanosensitive ion channels include Piezo1, MscL and TREK-2 [75]. In the case of Piezo1, a gated ion channel made up of three protein subunits that induce a dome-shaped structure on the cell membrane, the gating mechanism is triggered by the membrane surface tension. The membrane deformation induced by the surface tension acts as a mechanical signal that activates the protein subunits and causes them to undergo a conformational change that results in pore opening and transport of ions and solutes [76][77][78].
While the exact mechanism of mechanosensory transduction effected by the Piezo1 ion channel is still an open question, the extent of the deformed shape induced by the Peizo1 dome (referred to as the membrane footprint) is understood to significantly influence the sensitivity of the gating response of the channel [79]. As observed by Haselwandter and MacKinnon [79], an extended membrane footprint amplifies the sensitivity of Piezo1 subunits to respond to changes in the membrane surface tension. At the same time, increasing membrane tension significantly reduces the membrane footprint and thereby renders the Piezo1 subunits less sensitive to detect membrane mechanical signals.
In this analysis, we model the effect of surface tension on the area of the membrane footprint induced by the Piezo1 dome. The schematic for this boundary value problem is shown in Fig 2(b), and the simulation results demonstrating the effect of surface tension on the membrane footprint are presented in Figure 6. The plots show the 3D displacement profiles and their 2D projections under the boundary conditions enforced by the Piezo1 dome. A Piezo dome effect on the membrane is modeled by rotating the membrane (slope boundary condition) at the inner rim of the annular geometry to a value of φ = 70 degrees that is chosen so as to simulate the effect of a nearly hemispherical dome (which would correspond to φ = 90 degrees). This slope boundary condition assumes that the Piezo1 protein complex is a rigid dome that enforces a rotation on the surrounding membrane to ensure slope continuity between the hemispherical dome and the connected membrane. As can be seen from the subfigures Figure 6  and MacKinnon [79] that uses the classical reduced order Monge and arc-length axisymmetric parametrization methods to model the Piezo1-induced membrane deformation. Note that the deformation profile at the inner rim is, in general, non-axisymmetric, an effect that increases with membrane tension, γ. This illustrates the power of the 3D computational framework, which while it encompasses axisymmetric deformation, also admits non-axisymmetric modes. With access to the larger space, deformation profiles that are attainable at lower energies, and since the elasticity problem results in a (local-) minimum energy configuration, are indeed attained. Thus, while the 3D model reproduces the trends predicted by the reduced order models, its true power is in identifying more complex deformation patterns that are not accessible to the reduced order axisymmetric models.

Budding of membranes by protein coats during endocytosis
Budding of membranes by protein coats is a critical process in clathrin-mediated endocytosis (CME) that transports substance from the extracellular matrix to the cell interior. Several key features, including protein-induced spontaneous curvature, membrane properties, membrane tension, and force from actin polymerization, have been identified to govern the bud formation in CME [19,80,81]. The ability to simulate the morphological progression of bud formation 12/24 in the 3D setting under different combinations of these identified features is crucial for understanding the mechanical progression of CME. To further demonstrate the predictive capability of our simulation framework, we investigate the relationships between coat area, coat curvature, and degrees of symmetry during the budding of a vesicle as part of the endocytosis phenomena. The schematic of the simulation setup is given in Fig. 2(c). We simulated bud formation under two different conditions in our framework, similar to the setup explored in [19]. In the first case (i), the coated region has a fixed spontaneous curvature h 0 = 0.02 nm −1 with progressively increasing area of the coat. In the second case (ii), the coated region has a fixed radius of 80 nm with progressively increasing spontaneous curvature. In both cases, the uncoated membrane has a radius of 400 nm and has a surface tension with γ = 0.002 pN/nm. The bending modulus (k B ) for all cases is taken to be 320 pN·nm. The slope boundary condition with φ = 0 is enforced at both the inner and outer rims of the membrane through the penalty method to ensure the continuous differentiability with the flat membrane reservoir. As illustrated in Fig. 2(c), Dirichlet boundary conditions are enforced to eliminate rigid body motions. The hyperbolic tangent function proposed in [19,82] (see SI for illustration of this function) is used to ensure sharp but smooth transitions at the boundaries of the coated region and the uncoated membrane.
The simulation results from both setups are reported in Figs Next, we conducted six simulations for each case to construct the membrane morphology evolution phase diagram.
In each simulation of case (i), a different value of h 0 is assigned to the coated region. For a given value of h 0 we progressively increase the area of the coated region at an identical increment for all the simulations to allow the bud to form. In each simulation of case (ii), the radius of the coated region was set to a different value and then h 0 of the coated region was progressively increased at an identical increment for all the simulations to allow the bud to form. The membrane morphology evolution phase diagram for both simulation setups appears in Fig. 9 with arrows indicating progressively increasing quantities, where different patterns appear in the asymmetric region. To detect the symmetry breaking in each simulation, first, we uniformly sample 20 discrete values of u y between its minimum and maximum at each increment of the coat area for case (i) or the coat h 0 for case (ii). Next, the range of the curvature h, [h min , h max ], at every discrete u y is computed. For those heights with h min > 0, the relative change of h, denoted as ∆h = 2(h max − h min )/(abs(h max ) + abs(h min )), is computed. At each incremental step, we thus have multiple values of ∆h. Then the median value of ∆h, denoted as ∆h med , is computed for that step. Now, for each simulation, we have an array of ∆h med , whose length is equal to the total number of incremental steps of either the coat area or the coat h 0 .
Our results show that symmetry-breaking usually can be detected when ∆h med is at its minimum over increments of coat area or h 0 , pointing to a close to uniform value of ∆h for that increment. 1 . The associated values of h 0 and the 1 We remark that this is a loose criterion for detecting symmetry breaking, as not all the symmetry-breaking events occur precisely at the loading step where ∆hmed is minimal. However, this criterion generally provides a consistent indication of symmetry-breaking compared with our visual observation. Detailed discussion on the chosen symmetry-breaking criteria is provided in the SI.
14/24 area of the coated region at the specific incremental step where ∆h med is at its minimum are used to construct the phase diagram in Fig. 9. The dots with an empty surrounding square in Fig. 9 indicate the cases where the proposed symmetry breaking criterion does not hold. In the asymmetry region, not all the asymmetric patterns could be captured by the numerical simulations due to loss of numerical convergence during the solution iterations. This is due to the ill-conditioning of the system Jacobian matrix, that in turn is caused by the severe geometric and material non-linearity, potentially including bifurcation points, in the vicinity of asymmetric deformation modes. We use standard unconstrained optimization methods like arc-length and trust-region to achieve convergence of the solution iterations, whenever possible, and report successfully captured asymmetric patterns. The placement of the reported patterns are shown in Fig. 9. Here, the same markers are chosen to denote similar shapes. The fact that our computational framework could capture the symmetry breaking behavior, and even the pattern changes from 2-fold to 3-fold/4-fold/5-fold, demonstrates the advantages of the proposed 3D model over a reduced order axisymmetric model [19,20,80].

Conclusion
Biomembranes play a central role in various cell-scale and organelle-scale phenomena like locomotion of cells [2], packaging and trafficking of nutrients and signalling constituents [4], maintaining organelle morphology and functionality [6][7][8], etc. In almost all these processes, these surfaces are known to undergo significant deformation through bending; and the evolution of the out-of-plane bending deformation is a key mechanism of morphological evolution, besides in-plane fluidity. Thus, many analytical and numerical approaches exist in the literature to model bending and curvature generation, especially for solving the governing equations resulting from the Helfrich-Canham [14] characterization of membrane elasticity.
While these widely used analytical and numerical approaches (e.g. Monge parametrization, arclength parametrization and asymptotic methods) yield solutions to a wide range of boundary value problems of membrane bending, they are intrinsically limiting in capturing the complete envelope of membrane deformations due to the underlying axisymmetric restrictions on the kinematics and boundary conditions. Since the study of biomembrane deformation draws heavily from the well established models of elastic shells [14,27], it is only natural to look for the validity of axisymmetric approximations and for the existence of non-axisymmetric solutions in the deformation of elastic shell geometries. Interestingly, many elastic structures have inherent modes of instability that result in enhanced deformation or even collapse in response to loading and are associated with lower energy barriers. Such modes are ubiquitous in thin elastic shells and manifest as folding, wrinkling, creasing, and buckling (e.g. wrinkling of thin membranes and graphene sheets [83], surface tension induced buckling of liquid-lined elastic tubes [84], snap-through of elastic columns [85], barrelling modes of thin cylinders [86,87], etc.). Notably, they have lower symmetry than the fully axisymmetric deformations. If such modes exist, and are accessible in biomembranes, then they would naturally lead to a reduction in the load and energy barriers to membrane deformation, and may result in heretofore numerically unexplored deformation profiles and membrane morphologies. Accessing these lower symmetry modes and predicting the complex, three dimensional deformation profiles in biomembranes provided the primary motivation for developing the computational framework presented in this work.
We note that the first application of this framework was in our recent study demonstrating how non-axisymmetric buckling lowers the energy barrier associated with membrane neck constriction in biomembranes [21]. In that study, we 15/24 understanding of membrane deformation and the undesired effects of axisymmetry restrictions observed in that study, we have further developed the framework and expanded its scope to modeling other important membrane deformation processes.
Accordingly, in this work, we model three classical biomembrane problems: formation of tubular shapes and their lateral constriction, Piezo1-induced membrane footprint generation, and budding of membranes during endocytosis.
For each of these problems, we are able to validate against results and observation available in the literature for the simpler deformation modes, and also predict the more complex, less symmetric deformation profiles that are not accessible by the traditional analytical methods and axisymmetric numerical methods. Moreover, for the problem of endocytic vesicle budding, we also map a phase diagram classifying the symmetric and less-symmetric states.
The computational framework is implemented as an open-source software library and provided as a resource to the biophysics community. It is expected that this framework will serve as a platform for exploring complex deformation mechanisms (including geometric bifurcations and post-bifurcation responses) in biomembranes, and result in an improved understanding of the mechanics underlying various biomembrane phenomena. Future extensions envisioned are support for in-plane fluidity [60], surface diffusion (to model protein transport on the membrane) , and a contact model (to model membrane-membrane interactions). In addition, the inability of the current framework to apply non-uniform Dirichlet boundary conditions and constraints on displacement degrees of freedom inside the domain (i.e. at non-interpolatory knots of the spline surface) are significant limitations and will be addressed in future developments.