A mechanical model reveals that non-axisymmetric buckling lowers the energy barrier associated with membrane neck constriction

TL;DR

This study uses a mechanical model to show that non-axisymmetric buckling significantly reduces the energy barrier for membrane neck constriction, highlighting the importance of symmetry loss and helical forces in biological membrane scission.

Contribution

The paper introduces a mechanical model demonstrating how relaxing symmetry constraints and applying helical forces lower the energy barrier for membrane neck constriction.

Findings

Non-axisymmetric buckling lowers energy barriers for constriction.

Relaxing symmetry restrictions dramatically reduces the energy required.

Helical proteins further decrease the energy barrier compared to cylindrical proteins.

Abstract

Membrane neck formation is essential for scission, which, as recent experiments on tubules have demonstrated, can be location dependent. The diversity of biological machinery that can constrict a neck such as dynamin, actin, ESCRTs and BAR proteins, and the range of forces and deflection over which they operate, suggest that the constriction process is functionally mechanical and robust to changes in biological environment. In this study, we used a mechanical model of the lipid bilayer to systematically investigate the influence of location, symmetry constraints, and helical forces on membrane neck constriction. Simulations from our model demonstrated that the energy barriers associated with constriction of a membrane neck are location-dependent. Importantly, if symmetry restrictions are relaxed, then the energy barrier for constriction is dramatically lowered and the membrane buckles at lower values of forcing parameters. Our simulations also show that constriction due to helical proteins further reduces the energy barrier for neck formation compared to cylindrical proteins. These studies establish that despite different molecular mechanisms of neck formation in cells, the mechanics of constriction naturally leads to a loss of symmetry that can lower the energy barrier to constriction.

Tables

Table 1. Table 1 : Model parameter values

Parameter	Value	Reference
1. Boundary membrane tension (${\lambda }_0$)	10-2 - 10-1 pN$\cdot$nm-1	[51, 52]
2. Bending rigidity of bare membrane ($\kappa$)	320 pN$\cdot$nm	[53]
3. Length scale for non-dimensionalization ($R_0$)	20 nm	[9]

Equations

$$W=\kappa H^2+{\kappa }_{G}K,$$

$$E={\int }_{\Omega }(\kappa H^2+{\kappa }_{G}K+\lambda )dA.$$

$$\nabla \cdot \mathbf{\sigma }+p\text{n}=\text{f},$$

$$\nabla \cdot \mathbf{\sigma }={\mathbf{\sigma }}_{;\alpha }^{\alpha }=(\sqrt{a}{)}^{-1}(\sqrt{a}{\mathbf{\sigma }}^{\alpha }{)}_{,\alpha },$$

$${\mathbf{\sigma }}^{\alpha }={\mathbf{T}}^{\alpha }+S^{\alpha }\text{n},$$

$${\text{T}}^{\alpha }=T^{\alpha \beta }{\text{a}}_{\beta },T^{\alpha \beta }={\sigma }^{\alpha \beta }+b_{\mu }^{\beta }{M}^{\mu \alpha },S^{\alpha }=-M_{;\beta }^{\alpha \beta }.$$

$${\sigma }^{\alpha \beta }=\rho (\frac{\partial F(\rho ,H,K;x^{\alpha })}{\partial a_{\alpha \beta }}+\frac{\partial F(\rho ,H,K;x^{\alpha })}{\partial a_{\beta \alpha }}),$$

$$M^{\alpha \beta }=\frac{\rho }{2}(\frac{\partial F(\rho ,H,K;x^{\alpha })}{\partial b_{\alpha \beta }}+\frac{\partial F(\rho ,H,K;x^{\alpha })}{\partial b_{\beta \alpha }}),$$

$$H=\frac{1}{2}{a}^{\alpha \beta }{b}_{\alpha \beta },K=\frac{1}{2}{\epsilon }^{\alpha \beta }{\epsilon }^{\lambda \mu }{b}_{\alpha \lambda }{b}_{\beta \mu }.$$

$$F(\rho ,H,K;x^{\alpha })=\tilde{F}(H,K;x^{\alpha })-\frac{\gamma (x^{\alpha },t)}{\rho }.$$

$${\sigma }^{\alpha \beta }=(\lambda +W)a^{\alpha \beta }-(2H{W}_H+2\kappa W_K)a^{\alpha \beta }+W_H{\tilde{b}}^{\alpha \beta },$$

$$M^{\alpha \beta }=\frac{1}{2}{W}_H{a}^{\alpha \beta }+W_K{\tilde{b}}^{\alpha \beta },$$

$$\lambda =-(\gamma +W).$$

$$p+\text{f}\cdot \text{n}=\Delta \frac{1}{2}{W}_H+(W_K{)}_{;\alpha \beta }{\tilde{b}}^{\alpha \beta }+W_H(2H^2-K)$$

$$+2H(K{W}_K-W)-2\lambda H,$$

$$N_{;\alpha }^{\beta \alpha }-S^{\alpha }{b}_{\alpha }^{\beta }=-({\gamma }_{,\alpha }+W_K{k}_{,\alpha }+W_H{H}_{,\alpha })a^{\beta \alpha }$$

$$=(\frac{\partial W}{\partial x_{|exp}^{\alpha }}+{\lambda }_{,\alpha })a^{\beta \alpha }=\text{f}\cdot {\text{a}}_s.$$

$$\text{r}(s,\theta )=r(s){\text{e}}_r(\theta )+z(s)\text{k}.$$

$${\kappa }_{\nu }={\psi }^{{}^{\prime }},{\kappa }_{\tau }=r^{-1}\sin \psi ,\tau =0.$$

$$H=\frac{1}{2}({\kappa }_{\nu }+{\kappa }_{\tau })=\frac{1}{2}({\psi }^{{}^{\prime }}+r^{-1}\sin \psi ),K={\kappa }_{\tau }{\kappa }_{\nu }=\frac{{\psi }^{{}^{\prime }}\sin \psi }{r}.$$

$$r^{\prime }$$

$$r{\psi }^{\prime }$$

$$\frac{L^{\prime }}{r}$$

$$-2\left(H-C\right)\left[H^2+{\left(H-r^{-1}\sin \psi \right)}^2\right],$$

$${\lambda }^{\prime }$$

$$a(s)=2\pi {\int }_0^{s}r(\xi )d\xi \to \frac{da}{ds}=2\pi r.$$

$$\begin{array}{r}\alpha =\frac{a}{2\pi R_0^2},x=\frac{r}{R_0},y=\frac{y}{R_0},h=H{R}_0,c=C{R}_0,l=L{R}_0,\\ {\lambda }^{\ast }=\frac{\lambda R_0^2}{{\kappa }_0},p^{\ast }=\frac{p{R}_0^3}{{\kappa }_0},f^{\ast }=\frac{f{R}_0^3}{{\kappa }_0},{\kappa }^{\ast }=\frac{\kappa }{{\kappa }_0},\end{array}$$

$$x\dot{x}$$

$$x^2\dot{\psi }$$

$$\dot{l}$$

$$-2\left(h-c\right)\left[h^2+{\left(h-x^{-1}\sin \psi \right)}^2\right],$$

$$\dot{{\lambda }^{\ast }}$$

$$\zeta ={\alpha }_{bp}\frac{\alpha -{\alpha }_{bp}}{{\alpha }_{tot}-{\alpha }_{bp}},$$

$$x_1\frac{d{x}_1}{d\alpha }$$

$$x_1^2\frac{d{\psi }_1}{d\alpha }$$

$$\frac{d{l}_1}{d\alpha }$$

$$-2\left(h_1-c_1\right)\left[h_1^2+{\left(h_1-x_1^{-1}\sin {\psi }_1\right)}^2\right],$$

$$\frac{d{\lambda }_1^{\ast }}{d\alpha }$$

Full text

1Department of Mechanical and Aerospace Engineering, University of California San Diego, La Jolla CA 92093; 2 Department of Mechanical Engineering, University of Wisconsin-Madison, Madison, WI 53706, USA; 3 Biophysics Graduate Group, University of California, Berkeley, CA 94720, USA; 4 Department of Molecular and Cell Biology, University of California, Berkeley, CA 94720, USA; and 5 Departments of Mechanical Engineering and Mathematics, Michigan Institute for Computational Discovery & Engineering, University of Michigan, Ann Arbor, MI 48109, USA.

*∗*Corresponding Author

$|$Email: padmini.rangamani@eng.ucsd.edu $|$

Significance statement

Membrane tubule constriction is a critical step of cellular membrane trafficking processes and is thought to be mechanically regulated. Mechanical modeling techniques employing the Helfrich Hamiltonian and axisymmetric continuum frameworks have previously described energy barriers to constriction as a function of location along a membrane tubule. Recent advances in numerical modeling using spline basis functions (Isogeometric Analysis) enable us to conduct our analyses of membrane mechanics in a generalized 3D framework. Here, we implement a novel 3D Isogeometric Analysis framework and juxtapose it against an axisymmetric model to study the influence of location, symmetry constraints and helical collars on the constriction pathway. We show that an unsymmetric, “crushed soda can” neck consistently displays a lower energy barrier than a symmetric neck.

Introduction

Many cellular transport processes involving the plasma membrane including different forms of endocytosis [1, 2, 3], exocytosis [4, 5], and vesicle budding from intracellular organelles [6, 7] require mechanical deformation of the cellular membranes. The generation of membrane curvature is essential to trafficking, and the morphology of membranes has often been characterized as distinct shapes including U-and $\Omega$-shaped bud profiles [8, 9, 10] and tubulovesicular structures [11, 12]. The molecular mechanisms of these processes can be attributed to biochemical components of the protein machinery involved [13, 14]. For example, in the case of clathrin-mediated endocytosis (CME), more than 50 proteins are involved in regulating the different steps of membrane invagination such as nucleation, cargo selection, coat assembly, neck formation and scission [15, 3, 16] and contribute to the robustness and progression of endocytosis.

The formation of a membrane neck and scission are the last steps during many trafficking processes preceding vesicle formation. This neck formation is mediated by multiple biochemical mechanisms including mechanoenzymes belonging to the dynamin family [17], helix insertion due to BAR domain proteins [18] and ESCRT proteins [19, 20]. A common organizational feature of these different proteins is that they form helical assemblies at the membrane neck through oligomerization [21, 22].

Studies using reconstituted systems of lipid tubules decorated with protein assemblies have identified certain geometric and mechanical features of scission. Notably, studies of dynamin-mediated scission [23, 24, 25], the most investigated scission mechanism, have shown that the location of neck formation along a membrane tube, membrane tension, and bending rigidity play important roles in membrane tube constriction and scission [26]. Collectively, these studies support an emerging view that fundamental physical laws and geometric bounds confer a universality on membrane constriction phenomena and scission.

Crucially, neck formation occurs at a length scale of $<$ 10 nm, which is challenging to image even with high resolution electron tomography (ET) as radiation damage and low signal to noise ratio (SNR) can limit contrast [27]. Alternatively, equipped with extensive information from experiments such as those described above, mathematical and computational models can provide insight to the mechanics and energetics of membrane neck formation. Almost all of these models are rooted in the Helfrich elastic energy framework [28]. The physical principles underlying the Helfrich model are simple enough – the elastic energy of membrane deformation depends primarily on the curvatures of the membrane. Computational implementation of the governing equations resulting from this model, however, remain extremely challenging (see [29] for a detailed review). Therefore, many studies have assumed an axisymmetric configuration of the membrane for ease of computation [30, 31, 32, 33, 34, 9, 10].

In the most relevant of these studies to the present work, we and others have shown that a snap-through instability governs the first energy barrier associated with the formation of a membrane neck during CME [9, 35, 36]. An important limitation of the assumption of axisymmetry is that membrane deformation pathways associated with neck constriction that may have lower symmetries are not accessible (Figure 1B) and helicoidal protein assemblies [37, 38, 12] cannot be explicitly modeled.

In this study, we systematically investigate the energy barriers to constriction at different locations of a membrane geometry with and without symmetry restrictions (Figure 1). Importantly, we tackle the challenging problem of modeling non-axisymmetric membrane deformations with a benchmark comparison to axisymmetric modeling. We use a minimal, but fundamental, model of collar pressure-mediated tube constriction to obtain insights from a mechanical and energetic perspective. Using this model, we seek to answer the following foundational questions for the broader field of membrane deformation processes: First, how does the local pre-existing curvature along a tube influence the energy barrier associated with neck constriction? Second, how does relaxation of a priori imposed symmetry restrictions impact the energy barriers associated with constriction of the neck? And finally, how do cylindrical versus helical protein assemblies modulate this energy barrier?

To answer these questions, we have developed a computational framework for solving membrane mechanics problems on complex geometries using numerical techniques that exploit Galerkin methods, specifically Isogeometric Analysis [39]. This framework draws upon recent far-reaching advances on the use of spline basis functions in computational mechanics and brings them to the world of biological membranes, while building upon recent literature on finite element modeling of liquid shells [40]. As a result, we can now investigate membrane deformation using simulations of neck constrictions under conditions that are notably less restrictive than those adopted previously in the literature (i.e. no enforced axis of symmetry). Importantly, this allows us to probe realistic helical constriction pathways within a continuum framework, a different approach than recent efforts using coarse-grained modeling [37]. Using this framework, we applied constriction pressures at three different locations along the membrane tube (see Figure 2) – the ‘cap’(positive mean and Gaussian curvature), ‘cylindrical tube’ (positive mean and zero Gaussian curvature) , and ‘base’ (Positive- negative mean and negative Gaussian curvature). Our simulations show that the energy barriers associated with membrane neck constriction are indeed curvature-dependent, and therefore location-dependent, regardless of symmetry restrictions. Most importantly, we show that access to less symmetric shapes of membrane deformation lowers the energy barrier for scission considerably. These results suggest that loss of symmetry of the membrane neck may be an important energetic feature of successful neck formation.

Model development and simulations

Helfrich energy

The lipid bilayer is modeled as a thin elastic shell using the Helfrich energy [28] based on the assumption that the thickness of the bilayer is negligible compared to its radius of curvature [41, 23]. The Helfrich energy density is defined as

[TABLE]

where $\kappa$ is the bending rigidity, H is the mean curvature, K is the Gaussian curvature and $\kappa_{G}$ is the Gaussian rigidity. Furthermore, we assume that the membrane is incompressible (i.e the membrane area is constant ) [42]–a constraint that is implemented using a Lagrange multiplier field. Thus, while the Helfrich energy is defined entirely in terms of the geometry of the surface, the Lagrange multiplier, often interpreted as membrane tension [43, 44], is an important parameter that determines the minimum energy configuration. We ignore any fluid [45, 46] and friction [47, 48, 49] properties of the bilayer, guided by the dominance of unstable and stable equilibrium states over relaxation/rate processes. The augmented Helfrich Hamiltonian that is being minimized on the surface $\Omega$, including the Lagrange multiplier $\lambda$ is given as [44, 50, 45]

[TABLE]

Simulations in axisymmetric coordinates

In axisymmetry, the membrane is modeled using coordinates defined in Figure 1A. As the membrane tubule (Fig. 2E) has three distinct shape features (Fig. 2A - Cap, Tube and Base), local membrane geometries were modeled as a hemispherical cap (Fig. 2B, Case 1), cylindrical tube (Fig. 2C, Case 2) and a curved base with negative Gaussian curvature (Fig. 2D, Case 3). Cases 1 and 2 are constant mean curvature shapes and are solved as two-point boundary value problems. Case 3 is a negative Gaussian curvature shape with an inflection point in mean curvature with respect to the arc length and is solved as a three point boundary value problem. Case 4 includes local geometric variations in both mean and Gaussian curvature and is solved as a three point boundary value problem. The third point in these cases is an additional interface point enforced at the location of constriction, such that it satisfies continuity requirements [54]. The resulting system of equations is solved using the partial differential equation solution routines in Matlab, specifically bvp4c [54, 55]. Importantly, these equations are solved using both force control (compute membrane shape for a certain applied force) and displacement control (compute applied force for a certain membrane shape). These two approaches can lead to the same equilibrium membrane shape. In the presence of membrane bending instabilities, displacement control can access regimes of the response curve that force control cannot reach. However, this requires a precise prescription of the kinematic path. In order that a system be free to find the lowest energy pathways through a region of instability in its energy landscape, it is important that, while the force and displacement vary in a coupled manner, neither quantity be fully prescribed [56, 9]. Parameters for the bending rigidity, membrane tension and non-dimensionalization length R0 are specified in Table 1. Details of the numerical methods are provided in the Supplementary Online Material (SOM).

3D numerical model development and validation

The membrane deformation problems considered in this paper can be modeled using classical thin shell theories of mechanics. However, given the geometric complexity and the associated boundary conditions, analytical solutions are inaccessible. Instead, we obtain three dimensional numerical solutions to the membrane deformation problems using the framework of Isogeometric Analysis (IGA) [39]. An IGA method-based membrane mechanics framework has been developed for this work, and is build on top of the PetIGA [57] open source library. In an IGA approach, the membrane geometry is discretized using a spline mesh and the governing equations (Fig. 1B, see thin shell formulation in the Supplementary Information) are converted to a nonlinear system of equations. This nonlinear system of equations is then solved to obtain the deformed membrane shape, and the related force and energy metrics. Of importance to our central result is that this framework naturally admits both symmetric and asymmetric deformation modes driven by the underlying physics. This framework has three key assumptions. 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 [28]. This assumption allows us to neglect the effect of transverse shear deformations and consider the classical Kirchhoff–Love shell kinematics for thin shell geometries [58]. Second, numerical solutions to the membrane shape equations (Equations 14 and 15) in general coordinates are challenging because of continuity requirements in the numerical scheme. We have overcome this challenge by adopting both B-Spline basis functions, which allow high-order continuity, and the numerical framework of Isogeometric Analysis [39]. 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 [40], thus restoring stability to the numerical model. A companion manuscript (in preparation by the authors) describes the details of the mathematical methods and numerical formulation, and establishes the validity of the computational framework by modeling a range of problems in membrane mechanics. Here, we present a validation of the 3D computational framework by comparing the output from the simulation with a known analytical solution of the classical tube pulling problem (Figure S1A, B). In addition to demonstrating good agreement with the analytical solution, the 3D model also resolves the symmetric pathways of deformation if they are indeed the energy minimizing modes (Figure S1A). Having validated the 3D numerical scheme, we then proceeded to simulate the different cases shown in (Figure 2B-E) and compared them against axisymmetric pathways. We use three key metrics to compare the two models – (1) the radial pinching load, represented by the collar pressure that drives constriction, (2) structural stiffness of the membrane, defined as the slope of the load-displacement response, and (3) membrane bending energy. We track these metrics for different pinching radii, which are defined as the shortest distances between the membrane and the center of the necking region. For fully symmetric configurations and those with lower symmetry, this distance is the radius of the smallest circle that can be fit in the necking region.

Results

We systematically investigated the role of preexisting curvature (varying with location on the membrane) in the constriction process and the associated energy landscape using both the traditional axisymmetric calculations and the 3D computational framework. The constriction process is modeled using a collar pressure (pN/nm2) applied onto a fixed membrane height (nm). In this study, we include the effect of the height of applied pressure by reporting a force per unit length, or an effective surface pressure (pN/nm), as the product of the applied collar pressure and the fixed height. Our main results can be summarized as follows – first, the energy landscape for constriction depends on the preexisting curvature of the membrane; second, 3D modes of constriction with less than full symmetry encounter lower energy barriers when compared to pathways of higher symmetry; and finally, helical constriction modes can have the lowest energy barriers of all in 3D. We elaborate on these findings in detail below.

The energy barrier associated with constriction depends on preexisting membrane curvature

We investigated the effect of local, preexisting curvature on the energy barrier associated with tubule constriction in axisymmetry. We pulled out a membrane tube by applying an external axial force (faxial) on a small patch of the membrane to mimic a point load while maintaining a membrane tension of 0.2 pN/nm [59]. We then applied a radial collar pressure at different locations on the tube (Figure 3A) while maintaining the membrane height, a setup that can be generalized to in-vitro membrane tubules pulled by optical tweezers. In the absence of a fixed height applied as a boundary condition, the membrane deforms freely in the axial direction at negligible collar pressures (Figure S2). Results from our simulations show that pinching the tube at the cap (positive mean and Gaussian curvature) and along the cylinder (positive mean and zero Gaussian curvature) results in similar force-shape relationship (Figure 3D) and the cross section of the pinched profile remains circular by construction due to the restriction of axisymmetry. (Figure 3B,C). Surprisingly, for the same range of collar pressure applied to the base (positive-negative mean and negative Gaussian curvature), we observed the existence of a snap-through instability as the membrane constriction progresses, as shown by the red line in Figure 3D. The dotted lines (Figure 3D, base) are calculated using displacement control, i.e., compute the applied force given the membrane shape. However, given our initial conditions and the mechanism of neck formation via increasing pressure, these shapes are not accessible during constriction. As in all snap-through instabilities, this pinching instability arises from a reduced energy barrier and associated reduction in neck radius, and has been reported in other membrane physical processes as well [9, 35, 36]. Despite the existence of the snap-through instability at the base, the pressure needed for further constriction becomes unbounded as the pinching radius approaches zero. This suggests that fully symmetric membrane shapes are not favorable for constriction below a certain critical radius.

Relaxation of symmetry constraints lowers the energy barrier associated with membrane constriction

We next asked if relaxation of symmetry constraints alters the energy landscape of location-dependent constriction. To answer this question, we used our 3D model. For these simulations, we initialized the computation as a pre-formed membrane tubule to limit computational complexity (see SOM). Strikingly, we observed that once the symmetry constraints are relaxed, membrane constriction at all three locations requires a lower collar pressure by more than an order of magnitude when compared to the axisymmetric deformation (compare Figure 4D and Figure 3D). To verify this result, we enforced axisymmetry constraints in the 3D model and repeated our calculations for the “cap” (Figure 2B) and “tube” (Figure 2C). We observed that the collar pressures increased by an order of magnitude (Figures S3 and S4) when symmetry is imposed, resulting in comparable pinching profiles between the axisymmetric and 3D models. However, without the imposed axisymmetry, the collar pressures reduced significantly (Figure S5). These results allow us to conclude that absence of enforced axisymmetry alone is responsible for the significant decrease of collar pressure. We next analyzed the shapes of the membrane cross sections during 3D constrictions, which we found to be distinctly reminiscent of buckling phenomena that are observed in thin walled elastic structures [60, 61]. The pinching profiles shown in Figure 4C correspond closely to the classical result of the first buckling mode of a thin ring subjected to inward pressure on its walls. These deformation modes of buckling/pinching in thin-walled elastic rings and tubes have been known since the early twentieth century in the context of structural engineering applications [62, 63, 64, 65, 66, 67], and are also observed experimentally [68].

Finally, we observed that the base, with the preexisting negative Gaussian curvature, needed lower collar pressure to undergo constriction. This result is consistent with the observation that membranes with a negative Gaussian curvature are more amenable to constriction. Furthermore, contrary to the case of axisymmetric deformation, the collar pressure associated with increasing constriction at the base does not continue to grow with constriction (compare Figure 4D, red line and Figure 3D, red line). While the collar pressure increases sharply for initial constriction (Figure 4D, inset), it only rises gradually for the tube and cap geometries as constriction continues to increase several-fold. The initial increase in compressibility represents the natural stiffness of the membrane, after which symmetry breaking occurs leading to near spontaneous constriction. This result suggests that without the arbitrary restriction to axisymmetric deformation, near-spontaneous collapse of the neck is possible after a critical collar pressure is reached. The sudden drop in stiffness observed for the base geometry (Figure 4D) is associated with symmetry-breaking. As shown, this symmetry-breaking and loss of stiffness is not immediate, but occurs after a small amount of constriction has occurred. Thus, we predict that the energy landscape of membrane neck constriction is location dependent, but more importantly, predict that lower symmetry shapes attained by 3D constriction can significantly lower the energy barrier at the base to promote easier constriction.

Easier constriction at the base is accompanied by reduced membrane stiffness.

We further investigated the energy landscape at the base of the tubule (Figure 2D) to identify the mechanisms associated with easier constriction. For the axisymmetric pinching pathway, an increase in collar pressure results in progressive transformation of the tubule into a half-catenoid-shaped membrane. Indeed, this is the shape that is commonly seen in schematics of membrane pinching (Figure 5A). On the other hand, application of increasing collar pressure in 3D demonstrates that the membrane base is quick to break symmetry, and assumes the iconic shape of a soda can crushed by radial pinching (Figure 5B). We thus observe that the axisymmetric and 3D models invoke different constriction pathways; the axisymmetric model yields uniform pinching (Figure 5A), but the 3D model captures an asymmetric flattened tubule geometry (Figure 5B).

As before, the axisymmetric mode shows a snap-through instability (Figure 5C, Figure S6B). This instability can be attributed to a build up of negative tangential stress or tension work (Figure S6E). The sharp increase of both the bending energy (Figure S6D) and tension work (Figure S6E) at large constriction corresponds to a sharp increase in the pressure required for constriction (Figure 5C, Figure S6B). However, in 3D, a narrow constriction radius is accessible at a much lower pressure when compared with the axisymmetric mode. This result can be understood by analyzing the relationship between stiffness and constriction of the membrane in different modes of deformation (Figure 5D). The stiffness of the membrane is significantly reduced in the lower symmetry mode attained in 3D when compared with the axisymmetric mode (compare black circle lines and red diamond lines in Figure 5D). Similar comparisons for the tube (Figure S5) show a significantly reduced stiffness in lower symmetry pathways of 3D constriction when compared with axisymmetric pathways. Comparison of the bending energy in the axisymmetric and 3D modes of deformation shows that while the bending energies in both cases are similar in magnitude, the energy landscape is different (Figure 5E). All intermediate energy states along the constriction of the neck radius are accessible in the 3D pathway of deformation while in the axisymmetric pathway the energy states associated with the snap-through regime are inaccessible (dashed grey region in Figure 5E). From these analyses, we conclude that easier constriction at the base of the tube, revealed by full 3D computations, is accompanied by a reduced membrane stiffness and accessibility to all intermediate energy states along the constriction path.

Helical force collar further reduces the barrier to membrane neck constriction

The 3D model of membrane deformation allows us to probe the response of the membrane to non-symmetric force distributions such as those exerted by helical arrangements of proteins that cannot be modeled in the axisymmetric framework. Membrane scaffolding proteins such as dynamin [33, 38] and ESCRT-III [69] self-assemble into helical collars that can constrict the neck. Recent Cryo-EM maps of human dynamin-1 (dyn-1) polymer report detailed structural and molecular information on its helical geometry [12]. With GTP hydrolysis, the helical polymer actively constricts the membrane from a diameter of more than 20 nm to below 3.4 nm [12]. Disassembly of dynamin is thought to promote scission via hemifission [37, 70, 71]. Despite the well-established mechanisms of dynamin-mediated constriction, the response of the membrane to constriction and subsequent scission is not yet fully understood. This led us to investigate the role of a helical collar pressure, which is a mimic of force generated by helical protein assemblies, on membrane constriction. We explore the role of two geometric parameters of a helical collar – the pitch, defined as the distance along the axis for a complete helical turn, and the number of rings.

We first consider a single helical ring exerting a collar pressure on a cylindrical tube (Figure 6A) with different values of the pitch, non-dimensionalized by the height of the collar. Increasing helical pitch corresponds to an increased span (with same collar area) over which the load is distributed on the geometry. Qualitatively, we observe that the cross section of the neck is non-axisymmetric for different values of the pitch, confirming that the lower symmetry modes of deformation are still preferred for neck constriction with helical rings of pressure (Figure 6B, C). Quantitatively, the collar pressure associated with constriction decreases as the pitch increases (Figure 6D). Correspondingly, the stiffness also decreases for increased pitch (Figure 6E), indicating that the bending energy (Figure 6F) becomes slightly less steep. The pink and blue shaded regions in Figure 6E indicate regions of high and low stiffness (easier constriction) respectively. From these observations, we conclude that helical collars have the ability to further reduce the energy barrier to membrane constriction by a decrease in collar pressure and stiffness associated with this process, and that the pitch of the helix is an important determinant of this barrier.

To our knowledge, this is the first numerical characterization of the effect of helical squeezing forces on membrane constriction in a continuum framework. However, from a soft matter perspective, it is well-known that helical structures are known to exert squeezing forces more effectively on their support; an excellent example of this behavior can be found in the twining of plant vines [72, 73, 74] and other naturally occurring soft materials. Like cylindrical forces, helices also exert tangential and radial forces. Additionally, because of the pitch of the helix, they also exert axial forces.

Given our observation that the base of the tube is energetically favorable to constriction (Figure 5), we next added a helical collar pressure to the base to investigate the effect of coupling the three key design elements – negative Gaussian curvature, access to non-axisymmetric modes, and helical collar pressure (Figure 7A) – on membrane neck constriction. We found that this combination also results in reduced collar pressure (Figure 7D) and energy barrier (Figure 7F) with increasing constriction when compared to a ring of collar pressure (Figure 7D,F for zero pitch). However, a ring of lower pitch was advantageous in terms of a lower stiffness; the stiffness for Pitch 4 is greater than that associated with Pitch 0 (Figure 7E). This is possibly the result of recruiting the stiffer cylindrical portion of the membrane tube with greater pitch.

Comparing the effect of helical pinching on the tube (Figure 6) versus the base (Figure 7), we arrive at the following conclusions. For a single helical ring, the relationship between helix geometry and the membrane geometry is non-trivial. It appears while both the tube and base geometry show a clear relationship between helical pitch and collar pressure (Figure 6D, Figure 7D), the same is not true for the stiffness (Figure 6E, Figure 7E). This suggests a complex interaction between the geometries of the membrane and the helical ring, possibly due to a mechanical feedback between membrane curvature and the constricting action of the helicases.

Effects of multiple helical collars are also location-dependent

Since a helical ring of collar pressure offers a lower energy barrier to constriction, we next asked if an increase in the number of helical rings can further promote scission. This question is motivated by observations that an increase in the number of dynamin rings is a predicted response to delayed scission and higher membrane tension [75, 24, 76, 77]. To answer this question, we simulated 3 rings of helical collar pressure reminiscent of dynamin rings assembled on a membrane tubule [76, 38, 12]. More than 3 rings of a dynamin collar are unlikely to exist in vivo due to disassembly of the dynamin oligomer [78]. Collar pressure due to multiple helical rings not only leads to membrane constriction as expected (Figure 8A) but also appears to stabilize the membrane tube against a sideways wobble that is observed with one ring (compare Figure 6B with Figure 8B). We also observed that while the values of the collar pressure are of the same order of magnitude for 3 rings as with 1 ring (compare Figure 6D and Figure 8C), the stiffness profile is different (Figure 8D). During the early part of the constriction, the stiffness values go from high to low, quickly leading to a region of easier constriction (Figure 8D, blue shaded region). However as constriction increases, the stiffness increases again over a finite range of constriction and at a lower value of constriction (Figure 8D, pink shaded region), unlike the very narrow range of stiffening in Figure 6E. This may be due to the increase in the surface area that is constricted by three rings as compared to the surface area constricted by one ring, increasing the structural resistance to constriction. Thus, multiple rings assist neck formation on a tubule during an initial constriction region (Figure 8D, 4 - 15 nm of constriction), after which disassembly and possibly additional proteins are required.

Interestingly, the presence of 3 rings at the base resulted in an increase in both the collar pressure (Figure 8G) and stiffness (Figure 8H), such that the region of easier constriction (Figure 8H, blue shaded region) is much smaller than for a single helical ring (Figure 7E, blue shaded region). The membrane then transitions into a region of high stiffness at a smaller value of constriction (Figure 8H, pink shaded region) due to the larger surface area of the three rings that recruits more of the cylindrical tube to resist constriction.

Furthermore, multiple helical rings achieve easier constriction at shorter constriction distances for a cylindrical geometry (Figure 8D, blue region) and at larger constriction distances for the base geometry (Figure 8H, blue region). However, they resist further constriction at narrow radii independently of the pre-existing curvature (Figure 8D and Figure 8H, pink region) . Given these observations, it is possible that helical polymers might preferentially undergo conformational rearrangements such as a change in pitch or number of rings based on feedback with the underlying membrane curvature so as to achieve a lower energy barrier to constriction. Such structural rearrangements in dynamin have also been reported in experiments [12].

Discussion

Membrane constriction and subsequent scission are universal to membrane remodeling processes in vitro and in vivo. While the molecular machineries may differ across systems, these deformation processes likely share the same common physical principles. In this study, using computational modeling, we show that there are three key design elements that play important roles in promoting membrane constriction – (1) location i.e. preexisting curvature of the membrane being constricted, (2) access to lower-symmetry modes of deformation, and (3) access to helical loading.

From a mechanical standpoint, membrane constriction can be interpreted as a deformation mechanism driven by a radial collar pressure applied by the scission proteins in the vicinity of the necking region. For axisymmetric constriction, the pinching pressure needed to cause membrane constriction increases with the narrowing of the neck radius. This monotonic growth of the radial pressure results in a high energy barrier for pinching. 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 deformations (e.g. wrinkling of thin membranes and graphene sheets [79], surface tension induced buckling of liquid-lined elastic tubes [80], snap-through of elastic columns [81], barrelling modes of thin cylinders [60, 61], etc.). Notably, they have lower symmetry than the fully axisymmetric deformations. If such modes exist, and are accessible in cell membranes, their being triggered would naturally lead to a reduction in the energy barrier to constriction and scission. Building on this conventional understanding of buckling analysis of thin-walled structures, we predict the existence of lower energy modes of constriction in membrane tubules. The conclusions from our simulations provide insight to a number of recent experimental studies and suggest new experimental design as discussed below.

Dynamin and dynamin-related proteins (DRPs) have been shown to be essential for scission events during mitochondrial division [82] and during clathrin-mediated endocytosis via mechanical feedback with actin in both yeast [83, 84] and mammalian cells [85]. In dynamin-mediated fission during endocytosis [86, 87], dynamin preferentially interacts with curved membranes [88, 89], indicating a curvature dependence. Morlot et al. [23] showed that the local energy barrier to constriction is lower at the edge of the dynamin helix (large curvature) in optical tweezer experiments of dynamin-mediated fission. More recent experiments and models suggest that fission can also occur in the middle of the dynamin-coated region [90, 37]. While our results cannot confirm where fission will occur, we predict two important effects - (1) constriction is indeed curvature-dependent and (2) the membrane shape at the center of a given helical pitch is highly curved in 3D. These predictions are consistent with observations [91, 23, 90]. For example, Dar et al. [90] showed that dynamin1 polymers cause membrane constriction with high probability when the tubule radius approaches 16 nm or less, consistent with predictions from our model (see Figure 2F of [90] and compare against Figure 6E, J).

A central conclusion from this study is that a crushed soda can shape of the neck is energetically favorable for constriction over radially symmetric pinching. This prediction suggests that mechanisms such as those proposed in Figure 6 of Dar et al. can be revised to include lower degrees of symmetry (compare Figure 6 of [90] with Figure 6B). With advances in 3D imaging methods such as electron tomography, it should be possible to examine the cross-sections of necks during the progression of constriction by different molecular machines and quantify the relationship between membrane tubule symmetry and the particular protein assembly. Furthermore, determining the curvature-dependent rate constants for these proteins binding to the membrane will be important to quantify the relationship between the shape of the buckled membrane and the disassembly of monomers from polymerizing helical filaments such as dynamin [88, 92]. We predict that this feedback between membrane curvature and kinetics of helix assembly-disassembly is particularly important for the membrane curvatures where our simulations determine that it is energetically expensive for multiple rings to achieve the progression to scission.

Our results also apply to cases where dynamin is not involved in the scission process. In the absence of dynamin, BAR domain proteins and actin are thought to work closely in the formation of long tubular necks [93]. Indeed, in dynamin and clathrin-independent endocytosis, actin is the primary driver of scission of tubular invaginations via a constriction force [94]. These observations suggest that while scission may be less efficient, it is still functional in the absence of dynamin. Our results show that cylindrical collars, such as those enforced by actin, can promote constriction in the absence of dynamin.

From a structural mechanics standpoint, the differences in the membrane responses to helical versus cylindrical collars can also be understood by drawing analogies again with the buckling of thin cylindrical tubes. The distributed radial pressure in a helical collar creates an lateral torque that induces a bending moment on the tube. Under these conditions, the cylindrical tubes are now susceptible to both radial collapse (through pinching) and buckling under bending moment, which can cause accelerated pinching. The soda can shape also locally reduces the area moment of inertia and this can induce buckling through a process called buckling by ovalization [95]. In a completely different setting, the helical structures of twining vines are also known to exert squeezing forces on the support rods, suggesting that helical structures as force generating mechanisms are quite common in nature at different scales [73]. The helical collar mechanism opens up a wider parameter space (helical pitch, collar height, lateral bending, squeeze induced by helical twist, etc.) to optimize for achieving effective pinching. Such analogies with common engineering principles and with biological materials can help build our intuition on membrane-protein interactions; however, we note that the results presented in this work are specific to elastic, incompressible membranes only.

Based on the insights derived from our simulations, future work should include further complexities such as the influence of the structure of the helical polymer, the compositional heterogeneity of cellular membranes and the effect of contact constraints between the protein and tubule that can permit potential sliding of the protein on the tubule during the scission process. While recent molecular dynamics (MD) simulations of dynamin-mediated fission also reveal non axisymmetric pathways of constriction via the formation of transient pores [37], better connections between continuum descriptions of the lipid bilayer and membrane-protein interactions at the mesoscale need to be developed to close this gap. This is an ongoing research effort in our group.

Acknowledgements

We would like to thank David Drubin, Jasmine Nirody, and Morgan Chabanon for their feedback on the study. P.R. would like to acknowledge the Office of Naval Research N00014-17-1-2628. S.R. would like to acknowledge the Wisconsin Alumni Research Foundation (WARF) and the Grainger Institute for Engineering at UW-Madison for funding support, and thank Prof. Xiaoping Qian at UW-Madison for his advise on geometric modeling. M.A. would like to acknowledge the Arnold O. Beckman Postdoctoral Fellowship. KG would like to acknowledge NSF DMREF grant #1729166.

Supplementary Online Material

Assumptions

• •

The lipid bilayer is a modeled as a thin elastic shell. We use the Helfrich energy [28] based on the assumption that the thickness of the bilayer is negligible compared to its radius of curvature. This allows us to neglect shear deformations and consider classical Kirchoff-Love shell kinematics for thin shell geometries. Furthermore, we assume that the membrane is areally incompressible since the maximum elastic stretch is only 4 $\%$ [42]. This incompressibility constraint is numerically enforced using a Lagrange multiplier field. Additionally, we ignore any fluid [45] and friction [47] properties of the bilayer. Thus, the membrane is in mechanical equilibrium at all times.

• •

The lack of resistance to shear deformation modes in the Helfrich energy formulation can result in rigid body (zero energy) modes of deformation. To circumvent this limitation, in the 3D numerical simulations, we add shear stabilization terms to the classical Helfrich energy functional [40]. These stabilization terms are of a smaller magnitude relative to the traditional bending energy terms, and restore stability to the numerical model without significantly effecting the kinematics of bending.

• •

The membrane tubule is modeled both as an axisymmetric and 3D lipid bilayer. A pinching force is applied at different locations - ‘cap’, ‘tube’, and ‘base’ (Fig. 2). Since the tether is pulled from a membrane reservoir that can buffer changes in membrane tension [96, 51], we assume that elastic properties like membrane tension and bending rigidity are constant.

• •

Since we do not consider the fluid properties of the membrane, we cannot consider scission explicitly. We assume that the large stresses at the neck can lead to the formation of a hemi-fission intermediate [35, 32].

• •

The interaction of the constriction proteins and the membrane tubule can be numerically thought of as a contact model where the proteins apply a contact force of constriction on the tubule. Here, we do not consider a contact model but rather apply a follower load type collar pressure in the constriction region.

Thin shell formulations: Axisymmetric and 3D models

Equilibrium equations for the axisymmetric model

First, we write the force balance on the membrane as

[TABLE]

where $\bm{\sigma}$ is the stress tensor, p is the pressure difference between the inside and outside of the volume bounded by the membrane, and f is any externally applied force per unit area on the membrane. In our simulations, we assume that the tubule has equilibrated the pressure difference, and thus set p to 0. f includes both the axial and pinching forces applied on the membrane. By introducing the covariant derivative as $()_{;\alpha}$, the surface divergence in Eq. 3 can be rewritten as [44]

[TABLE]

where $a$ is the determinant of the first fundamental form metric $a_{\alpha\beta}$. The surface stresses in Eq. 3 can be split into normal and tangential component given by

[TABLE]

where

[TABLE]

The two tensors $\sigma^{\alpha\beta}$ and $M^{\alpha\beta}$ can be expressed by the derivative of $F$, the energy per unit mass, with respect to the coefficients of the first and second fundamental forms, $a_{\alpha\beta}$, $b_{\alpha\beta}$, respectively [45, 44]

[TABLE]

where $\rho$ is the surface mass density. $H$ and $K$ are mean and Gaussian curvatures given by

[TABLE]

Here $(a^{\alpha\beta})=(a_{\alpha\beta})$ is the dual metric and $\varepsilon^{\alpha\beta}$ is the permutation tensor defined by $\varepsilon^{12}=-\varepsilon^{21}=\frac{1}{\sqrt{a}},\varepsilon^{11}=\varepsilon^{22}=0$.

A reasonable assumption to make is that the membrane tubule has a fixed area. We introduce an area incompressibility ($J=1$) constraint using a general form of free energy density per unit mass given as

[TABLE]

Here $\gamma(x^{\alpha},t)$ is a Lagrange multiplier field required to impose invariance of $\rho$ on the whole of the surface (see [44] for full derivation). Substituting $W=\rho\tilde{F}$ into Eq. 10 we get

[TABLE]

where

[TABLE]

Combining Eqs. 12, 6, and 5 into Eq. 3 gives the equations in normal and tangential equations as

[TABLE]

and

[TABLE]

Here $\Delta(\cdot)$ is the surface Laplacian and $()_{|exp}$ denotes the explicit derivative respect to coordinate $\theta^{\alpha}$.

Axisymmetric model

Using the axisymmetric parametrization

[TABLE]

we define $\psi$ as the angle made by the tangent with respect to the horizontal (see Fig. 1). This gives $r^{\prime}(s)=\cos(\psi)$, $z^{\prime}(s)=\sin(\psi)$, which satisfies the identity $(r^{\prime})^{2}+(z^{\prime})^{2}=1$. Using this, we define the normal to the surface as $\textbf{n}=-\sin\psi\textbf{e}_{r}(\theta)+\cos\psi\textbf{k}$, the tangent to the surface in the direction of increasing arc length as $\textbf{a}_{s}=\cos\psi\textbf{e}_{r}(\theta)+\sin\psi\textbf{k}$, and unit vector $\tau$ = $\textbf{e}_{\theta}$ tangent to the boundary $\partial\omega$ in the direction of the surface of revolution. For more details, we refer the reader to [44, 9, 43].

The expressions for tangential $(\kappa_{\nu})$, transverse $(\kappa_{\tau})$ and twist $(\tau)$ curvatures are simplified as

[TABLE]

Further, we calculate the mean curvature ($H$) and Gaussian curvature ($K$) as

[TABLE]

We introduce a term $L=\frac{1}{2\kappa}r(W_{H})^{\prime}$ in order to write a system of first order differential equations governing the problem [9] as ,

[TABLE]

Eq. 19 is a function of the arc length (s). This can be rewritten in terms of membrane area (a) using

[TABLE]

We choose to non-dimensionalize our system of equations using a length scale $R_{0}$ and bending rigidity scale $\kappa_{0}$ as

[TABLE]

where R0 is the radius of the flat patch of membrane in simulations of a membrane tubule (Fig. 2E), the radius of the hemisphere for simulations of the tubule cap (Fig. 2B), and the radius of the tube for simulations of the tube (Fig. 2C) and the base (Fig. 2D). $\kappa_{0}$ is the bending rigidity of the bare membrane.

Rewriting Eq. 19 using the dimensionless variables in Eq. 21, we get [9]

[TABLE]

We enforce a third boundary point for constriction simulations at the base (Fig. 2D) and the whole tube (Fig. 2E) by introducing an independent variable [54]

[TABLE]

where $\alpha_{bp}$ is the non-dimensional area of the first ‘phase’ and $\alpha_{tot}$ is the total non-dimensional area of the membrane. $\alpha$ is the variable defining the non-dimensional area along the first ‘phase’ and $\zeta$ is the variable defining non-dimensional area along the second ‘phase’. Like $\alpha$ in the first interval, $\zeta$ ranges from o to $\alpha_{bp}$ in the second interval. Thus we can redefine our system of equations (Eq. 22) for 2 phases as

[TABLE]

where 1 and 2 are the two phases represented by non-dimensional areas $\alpha$ and $\zeta$, both of which are defined in the interval [0 $\alpha_{bp}$]. Applied forces are modeled as a smooth hyperbolic tangent function. For example, the axial force at the tip of the tubule and radial force at the interface are modeled as

[TABLE]

where $f_{axial}$ and $f_{pinch}$ are the distribution of axial and radial forces per unit area along the non-dimensional membrane area $\alpha$, $f_{ax}$ and $f_{p}$ are the magnitudes for force per unit area, and $\alpha_{\text{axial}}$ and $\alpha_{\text{radial}}$ are the corresponding non-dimensional areas of applied force in the axial and radial direction respectively. g is a constant that ensures a sharp but smooth transition. In our simulations, we use g = 20 [9].

Boundary conditions

Eq 24 can be solved given a set of boundary conditions. All the axisymmetric simulations in this study were performed using the MATLAB bvp4c toolbox [54]. A useful feature of this toolbox is the estimation of unknown parameters by providing additional boundary conditions [54]. The MATLAB subroutines used in this work are available on Github [97].

• •

Whole tubule

The tubule pinching simulations in Fig. 2E were performed using the following set of boundary conditions

[TABLE]

$\mathbf{x_{p}}$ and $\mathbf{y_{p}}$ are additional constraints for the radius at the interface and height of the tubule respectively. These additional constraints are used to estimate the axial force and pinching force required to obtain a solution to the system of equations in Eq. 24. $\mathbf{\lambda_{0}}$ is the boundary membrane tension. We note here that this now represents a system of 12 equations and 2 unknown parameters with 14 boundary conditions. The 2 unknown parameters are the axial and radial force. 6 of these boundary conditions are continuity conditions for every parameter at the interface (Eq. 27).

• •

Base

The half catenoid-like base pinching simulations in Fig. 2D were performed using the following set of boundary conditions

[TABLE]

$\mathbf{x_{p}}$ is the additional constraint for the radius at the interface. $\mathbf{x_{0}}$ is the radius at y = 0. The additional constraint is used to estimate the pinching force required to obtain a solution to the system of equations in Eq. 24. $\mathbf{\lambda_{0}}$ is the boundary membrane tension. We note here that this now represents a system of 12 equations and 1 unknown parameter with 13 boundary conditions.

• •

Tube

The cylinder/tube pinching simulations in Fig. 2C were performed by solving the system of equations in Eq. 22 for a single phase of membrane. The boundary conditions used were

[TABLE]

where $\alpha$ is the non-dimensional area of the tube, $\mathbf{x_{p}}$ is the additional constraint required for estimating the pinching force. This represents a system of 6 equations and 1 unknown parameter with 7 boundary conditions.

• •

Cap

The hemisphere/cap pinching simulations in figure 2B were performed by solving the system of equations in Eq. 22 for a single phase of membrane. The boundary conditions used were

[TABLE]

where $\mathbf{x_{p}}$ is the additional constraint required for estimating the pinching force. This represents a system of 6 equations and 1 unknown parameter with 7 boundary conditions.

3D model

Thin shell formulation

Considering the classical Helfrich formulation for biological membranes, the strain energy density of a membrane in the current configuration is given by:

[TABLE]

where $k_{B}$ and $k_{G}$ are the bending modulus and the Gaussian modulus of the membrane, $H$ is the mean curvature, $\kappa$ is the Gaussian curvature and $H_{0}$ represents the instantaneous curvature induced in the membrane.

To enforce area-incompressibility, we consider the following Lagrange multiplier formulation:

[TABLE]

where q is the point value of the Lagrange multiplier field, and J is the surface stretch (ratio of area in the current configuration to the area in the reference configuration).

The governing equation for quasi-static mechanical equilibrium in 3D simulations is obtained by minimizing the Helfrich energy functional following standard variational arguments, and is given by [40]:

[TABLE]

where $\partial\Omega$ is the membrane surface and $\Gamma$ is the membrane boundary on which surface tractions can be applied, as shown in Figure 1B. $\delta a_{ij}$ and $\delta b_{ij}$ are variations of the components of the metric tensor and the curvature tensor, respectively. $\sigma^{ij}$ are the components of the stress tensor, $M^{ij}$ are components of the moment tensor, $\bm{p}$ is the pressure applied along a collar on the membrane surface (to cause constriction), and $\bm{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 [98]:

[TABLE]

For the Helfrich type strain energy density, these take the form:

[TABLE]

Computational implementation

We solve the governing equation given by Eq.33 using a Isogeometric Analysis (IGA) based numerical framework for solving problems of membrane mechanics developed as part of this work. A companion manuscript (in preparation by the authors) describes the details of the mathematical methods and the numerical formulation. The computational implementation, along with the source code for solving the boundary value problems listed below, is available as a public code repository on GitHub [99].

List of 3D simulations

For each of the 3D simulations, we solve the governing equation given by Eq.33 using a force control or displacement control approach, with the relevant displacement, angle and traction boundary conditions. The displacement boundary conditions are applied on the components of the displacement vector, $u$, that is defined as the change in position of a point on the membrane between its current and reference configuration (\mbox{\boldmathu}(\xi_{1},\xi_{2})=\textbf{x}(\xi_{1},\xi_{2})-\textbf{X}(\xi_{1},\xi_{2})). The angle boundary conditions, where needed, are enforced through the weak formulation using a penalty approach and result in the normal vector (n) at the boundaries to align along the preferred direction. In all the simulations, we have two boundaries, and these are identified as the inner boundary ($\Gamma_{\text{I}}$) and the outer boundary ($\Gamma_{\text{O}}$) as indicated in the schematic in Figure 1B. The specific numerical simulations in this work using the 3D model are described below.

• •

Tube pulling

The tube pulling simulation shown in Figure S1(A) considers a reference circular plate geometry with an outer radius of 20 nm, and an inner radius of 0.2 nm. The boundary value problem is solved as a force control problem with a traction on the inner boundary ($\Gamma_{\text{I}}$). The displacement and traction boundary conditions are as follows:

[TABLE]

See Movie M12 in the supplementary information for the evolution of the membrane deformation.

• •

Whole tubule

The whole tubule simulations shown in Figure 4 consider pinching at three different locations, identified as the cap, tube and base locations (Figure 2A). For the tubule geometry, shown in Figure 2A, the tubule radius is 20 nm and height is 100 nm. The inner boundary ($\Gamma_{\text{I}}$) at the top of the tubule has a radius of 0.2 nm and the outer boundary ($\Gamma_{\text{O}}$) at the base of the tubule has a radius of 40 nm. The boundary value problem is solved as a force control problem with pressure applied on a collar ($\Omega_{\text{collar}}$) located at the cap, tube or base location. In addition, the displacement boundary conditions are as follows:

[TABLE]

See Movies M1-M3 in the supplementary information for the evolution of the constriction process for the cap, tube and base locations.

• •

Base

The constriction simulation shown in Figure 5 considers pinching at the base location. The tube geometry considered has a radius of 20 nm, and a height of 80 nm. The tube boundary on the top is identified as the inner boundary ($\Gamma_{\text{I}}$) and the tube boundary at the bottom is identified as the outer boundary ($\Gamma_{\text{O}}$). The boundary value problem is solved as a force control problem with pressure applied on a collar ($\Omega_{\text{collar}}$) located at the base location. In addition, the displacement boundary conditions are as follows:

[TABLE]

See Movie M3 in the supplementary information for the evolution of the constriction process for the base location.

• •

Cap

The constriction simulation shown in Figure S3 considers pinching at the cap location. The cap geometry is a hemisphere with a radius of 20 nm. The cap boundary on the top, with a small radius of 1 nm, is identified as the inner boundary ($\Gamma_{\text{I}}$) and the cap boundary at the bottom is identified as the outer boundary ($\Gamma_{\text{O}}$). This boundary value problem is solved as a displacement control problem, as the force control problem is numerically unstable due to the rigid body modes induced under the displacement boundary conditions considered. As this problem is solved as a displacement control problem, this enforces axisymmetry of the pinching profile. The displacement boundary conditions are as follows:

[TABLE]

See Movie M13 in the supplementary information for the evolution of the constriction process for the cap location.

• •

Tube

The constriction simulation shown in Figure S4 considers pinching at the tube location. The tube geometry is a cylinder with a radius of 20 nm. The tube boundary on the top is identified as the inner boundary ($\Gamma_{\text{I}}$) and the boundary at the bottom is identified as the outer boundary ($\Gamma_{\text{O}}$). This boundary value problem is solved as a displacement control problem. Like in the case of the cap simulation, as this problem is solved as a displacement control problem, this enforces axisymmetry of the pinching profile. The displacement boundary conditions are as follows:

[TABLE]

See Movie M14 in the supplementary information for the evolution of the constriction process with displacement control for the tube location.

We also solve a force control equivalent of this problem, and this is shown in Figure S5. This case is discussed below in the simulation of the helical force collar at the tube location with a zero helical pitch. See Movie M15 in the supplementary information for the evolution of the constriction process with force control for the tube location.

• •

Helical force collar at the tube location

The constriction simulation shown in Figure 6 considers pinching at the tube location due to helical collar. In Figure 6A we consider a single helical ring, and in Figure 6G we consider three helical rings. The tube geometry considered for the single helical ring case has a radius of 20 nm, and a height of 40 nm. The tube geometry considered for the three helical rings case has a radius of 20 nm, and a height of 200 nm. For both cases, the tube boundary on the top is identified as the inner boundary ($\Gamma_{\text{I}}$) and the tube boundary at the bottom is identified as the outer boundary ($\Gamma_{\text{O}}$). The boundary value problem is solved as a force control problem with pressure applied on a helical collar ($\Omega_{\text{collar}}$) located at the tube location. In addition, the displacement boundary conditions are as follows:

[TABLE]

See Movies M4-M6 in the supplementary information for the evolution of the constriction process due to a helical force collar at the tube location with a non-dimensional pitch of zero, two and four, respectively, and movie M7 for the corresponding evolution of the constriction process due to a force collar with three helical rings.

• •

Helical force collar at the base location

The constriction simulation shown in Figure 7 considers pinching at the base location due to helical collar. In Figure 7A we consider a single helical ring, and in Figure 7G we consider three helical rings. The tube geometry considered for the single helical ring case has a radius of 20 nm, and a height of 40 nm. The tube geometry considered for the three helical rings case has a radius of 20 nm, and a height of 200 nm. For both cases, the tube boundary on the top is identified as the inner boundary ($\Gamma_{\text{I}}$) and the tube boundary at the bottom is identified as the outer boundary ($\Gamma_{\text{O}}$). The boundary value problem is solved as a force control problem with pressure applied on a helical collar ($\Omega_{\text{collar}}$) located at the base location. In addition, the displacement boundary conditions are as follows:

[TABLE]

See Movies M8-M10 in the supplementary information for the evolution of the constriction process due to a helical force collar at the base location with a non-dimensional pitch of zero, two and four, respectively, and movie M11 for the corresponding evolution of the constriction process due to a force collar with three helical rings.

Analytical solution for tube pulling simulation

The equilibrium values of R0 and f0 for a membrane tube are defined as

[TABLE]

where $\kappa$ is the bending rigidity, $\sigma$ is the membrane tension. For bending rigidity of 20 pN$\cdot$nm and membrane tension 0.1 pN/nm, we get

[TABLE]

which is the equilibrium value of force obtained in Fig. S1.

List of movies

• •

Movie M1: Evolution of the constriction process for the cap location for the whole tubule geometry.

• •

Movie M2: Evolution of the constriction process for the tube location for the whole tubule geometry.

• •

Movie M3: Evolution of the constriction process for the base location for the whole tubule geometry.

• •

Movie M4: Evolution of the constriction process due to a single helical force collar at the tube location with a non-dimensional pitch of zero.

• •

Movie M5: Evolution of the constriction process due to a single helical force collar at the tube location with a non-dimensional pitch of two.

• •

Movie M6: Evolution of the constriction process due to a single helical force collar at the tube location with a non-dimensional pitch of four.

• •

Movie M7: Evolution of the constriction process due to a force collar with three helical rings at the tube location.

• •

Movie M8: Evolution of the constriction process due to a single helical force collar at the base location with a non-dimensional pitch of zero.

• •

Movie M9: Evolution of the constriction process due to a single helical force collar at the base location with a non-dimensional pitch of two.

• •

Movie M10: Evolution of the constriction process due to a single helical force collar at the base location with a non-dimensional pitch of four.

• •

Movie M11: Evolution of the constriction process due to a force collar with three helical rings at the base location.

• •

Movie M12: Evolution of membrane deformation for pulling of a tubule from a flat membrane.

• •

Movie M13: Evolution of a axisymmetric constriction profile due to a displacement control approach of constriction at the cap location.

• •

Movie M14: Evolution of a axisymmetric constriction profile due to a displacement control approach of constriction at the tube location.

• •

Movie M15: Evolution of a non-axisymmetric constriction profile due to a force control approach of constriction at the tube location.