Coupling anisotropic curvature and nematic order: mechanisms of membrane shape remodeling

Abstract

This study theoretically investigates how anisotropic curved membrane components (CMCs) control vesicle morphology through curvature sensing, nematic alignment, topological defects and volume constraints. By comparing arc- and saddle-shaped CMCs, we identify a rich spectrum of steady-state phases. For fully CMC-covered vesicles, arc-shaped components drive a pearling-to-cylinder transition as nematic interactions strengthen, while on partially CMC-covered vesicles the saddle-shaped CMCs stabilize necks between the convex regions of bare membrane. We map the steady-state shapes of vesicles partially covered by arc- and saddle-shaped CMCs, exposing how different vesicle shapes depend on the interplay between nematic interactions and volume constraints, revealing several novel phases. By investigating the in-plane nematic field, we find that topological defects consistently localize to high-curvature regions, revealing how intrinsic and deviatoric curvature effects cooperate in membrane remodeling. These findings establish a unified framework for understanding how proteins and lipid domains with anisotropic intrinsic curvature shape cellular structures—from organelle morphogenesis to global cell shape.

---

1. Introduction

Membranes are fundamental structures in cells that act as barriers between compartments and facilitate essential biological processes. These membranes dynamically change shape to accommodate various cellular functions, influenced by interactions with membrane proteins and lipids.^1–6 Understanding how these shape changes occur is crucial to understanding the intricate workings of a cell. An integral part of lipid membranes is curved membrane components (CMCs), aggregates of proteins or lipid rafts⁷ that can move laterally along the membrane surface, bend the membrane locally, and are sensitive to local membrane curvature.

There is ample experimental evidence to suggest that CMCs have different curvature preferences along different directions, which means they are intrinsically anisotropic. A notable example of anisotropic CMCs is the family of BAR domain proteins.^1,8 Amphiphysin 1, an N-BAR protein, concentrates in membrane nanotubes and induces tubulation.^9,10 Another example is IRSp53, a protein that contains an I-BAR domain and exhibits a strong preference for negatively curved membranes and is found on the inner leaflet of membrane nanotubes pulled from giant unilamellar vesicles.^3,11–13

The binding of BAR domains to lipid bilayers is mediated by electrostatic interactions. The concave, positively charged surface of BAR domains attracts negatively charged lipids, initiating membrane contact at the tips. As hydrophobic loops insert into the bilayer, local bending occurs due to leaflet asymmetry. Stronger electrostatic attraction at the BAR domain edges promotes lipid demixing, enhancing curvature induction.¹⁴ The electric field strength increases with curvature radius, which means that BAR proteins with more flat intrinsic curvature exert greater electrostatic pull, ultimately molding the membrane to match their shape.

Furthermore, the interplay between curvature and electrostatics at protein–protein interfaces influences binding affinity and residue composition. While charged and polar residues dominate small interfaces due to pronounced electrostatic effects, larger interfaces tend to bind more due to hydrophobic interactions. Positively charged domains (e.g., BAR proteins, antimicrobial peptides) preferentially bind to negatively charged lipids, inducing local inward bending by attracting the anionic membrane surface.¹⁵ In saddle (Gaussian) curvature, electrostatics can promote pore formation where charged residues stabilize high-curvature edges via lipid demixing and dipole rearrangements, so saddle-like anisotropic CMCs can stabilize pores.¹⁶ Ionizable groups further modulate this effect: protonation/deprotonation in response to curvature changes (e.g., lysine at pore rims) can switch between stabilizing negative or saddle curvature. In this way, electrostatics not only selects the curvature sign but also couples topological changes to chemical state.

Coupling between the non-homogeneous lateral distribution of membrane components and the local anisotropic membrane curvature has also been recently indicated in the Golgi, where some of the membrane components are concentrated on the bulbous rims of the Golgi vesicles and where the difference between the two principal membrane curvatures is very large.¹⁷ Similar phenomena have also been observed in photoreceptor discs,^18,19 endoplasmic reticulum shapes,²⁰ and flattened endovesicles of the erythrocyte membrane.²¹ These examples indicate that the coupling between the non-homogeneous lateral distribution of generally anisotropic membrane components may be a general mechanism stabilizing highly curved membrane structures and demonstrate the importance of anisotropic membrane components in driving cellular processes.^{5,7,10,22–31} Sometimes, several types of curved-like proteins may be working in a coordinated manner to induce membrane morphologies, and tubulation can be studied in the context of a mixture of positive and negative curvature proteins.³²

Numerical simulations of membranes that are not bound by constraints of axisymmetrical shapes have been a valuable tool in understanding the coupling between membrane curvature and curved membrane components (CMCs) that can drive feedback loops by curvature sensing and clustering. In a comprehensive review by Ramakrishnan et al., shape transformations of membranes with an in-plane nematic field were explored for partial and full coverage fraction.³³ Similarly, an open source tool was developed to study the effects of nematic interactions used for the analysis of real membranous systems.³⁴

However, studies of membrane deformation that include anisotropic CMCs in combination with the constraint of fixed volume remain poorly understood and to our knowledge, lacking. Here, we expand upon our previous work³⁵ and present new numerical results obtained from Monte Carlo simulations of triangulated, closed membrane shapes using anisotropic CMCs which interact nematically. We explore how the spontaneous curvature of CMCs (arc-shaped or saddle-shaped), CMC concentration, and strength of nematic ordering, together or in the absence of constraints on volume, influence the overall membrane shape. We first start with vesicles fully covered with arc- and saddle-like CMCs and find that the steady-state shapes are mostly cylindrical and globally saddle-shaped, respectively. We find a pearling-to-cylinder transition that occurs at low nematic interaction strength. When the vesicles are only half covered by non-interacting arc-shaped CMCs, with constraints of a constant volume, we find that there is an accelerated tendency toward global prolate shapes in comparison with empty vesicles, which can sustain an oblate phase up to larger relative volumes. When nematic interaction is present, we explore the phase diagram in the space of relative volume and CMC interaction strength, and characterize a boomerang-like phase that marks the transition from the oblate to the prolate phase. For saddle-like CMCs, we find that full coverage and volume constraints can lead to global saddle shapes with the formation of singular membrane protrusions. When the concentration of CMCs is decreased, the volume-binding strength phase diagram features stomatocyte, oblate-like, and protrusive phases. We find that at small saddle-like CMC concentrations, these have a tendency to aggregate in the necks between two convex parts of the membrane, which leads to the overall elongation of the vesicles. We find that simple physical interactions between anisotropic CMCs and their associated in-plane nematic fields can drive diverse membrane transformations, such as tubule formation and pore creation.³⁶

2. The model

2.1. Membrane–protein interactions

The theoretical framework is grounded in the derivation of the membrane bending energy and anisotropic inclusion interactions.^37–40 The model describes the membrane as a 2D anisotropic continuum surface covered with an arbitrary fraction of CMCs. The elastic energy of a small membrane section is lowest when its principal curvatures C₁ and C₂ align with the intrinsic curvatures of the anisotropic CMCs, C_1m and C_2m, at that location. The level of alignment can be quantified by comparing the orientations of the local membrane curvature tensor C and the intrinsic membrane curvature tensor C_m.
These two tensors—describing the directions of the membrane and CMC principal curvatures, respectively—are generally rotated relative to each other by an angle ω in the tangent plane. The mismatch tensor is defined as M = RC_mR⁻¹ − C, where R is the rotation matrix that quantifies the angular disparity:^38–40
A CMC reorients in the tangent plane by ω radians to match the membrane curvature and minimize energy, reflecting the energetic cost of deformation (see Fig. 1(b)). The approximate elastic energy, E₁, is expressed as a series expansion in the independent invariants of the mismatch tensor M. Using the trace and determinant of M as invariants yields the following expression:^38–41integrated across the entire vesicle. Substituting the bending moduli K₁ and K₂ gives:^31,38–40where D = (C₁ − C₂)/2 is the curvature deviator, while H_m = (C_1m + C_2m)/2 is the intrinsic mean curvature, and D_m = (C_1m − C_2m)/2 the intrinsic curvature deviator. The curvature deviator D can be nondimensionalized as d = DR, where R is the radius of a sphere with the same volume V and area A as the vesicle. In the case of an isotropic membrane where D_m = 0, it becomes evident that eqn (2) is equal to the Helfrich bending energy density⁴² described by E_b = k_c²(2H − C₀)² + k_GK, where H = (C₁ + C₂)/2 is the mean curvature, C₀ is the spontaneous curvature, K = C₁C₂ is the Gaussian curvature, and K₁ = k_c and K₂ = k_G, where k_c and k_G are the bending and splay moduli, respectively.^38,40,41 To compute the principal curvatures at each vertex of a triangulated membrane mesh, we adapt the method of Ramakrishnan et al.³³ with several key modifications (see Appendix B.1). While our method adopts the same CMC–CMC interaction model as Kumar et al.,⁴³ it differs fundamentally in the formulation of membrane bending energy (eqn (1)).

Fig. 1 (a) Example of an arc-shaped CMC with H_m = D_m = 0.25 (C_1m = 0.5 and C_2m = 0). (b) Example of a saddle-shaped CMC, with H_m = 0, D_m = 0.5. (c) The mismatch between the principle curvatures of a membrane protein and the local curvature of the triangulated membrane. Here, n⃑ is the orientation of the CMC and t⃑ is its perpendicular direction in the local tangent plane.

2.2. Isotropic protein–protein binding

Two neighboring CMCs may bond to each other on the membrane, thereby lowering the overall energy:Here, the type and binding strength w are vertex properties that reflect local aggregation. Throughout this paper, we refer to this interaction as isotropic bonding. Such bonding is commonly observed in surfactants and proteins, as it results from spontaneous self-assembly on the membrane driven by hydrophobic and other intermolecular effects.^44,45

2.3. Nematic protein–protein interactions

Anisotropic CMCs exhibit nematic interactions similar to those observed in liquid crystals. The energy between neighboring CMCs will be lower if their principal curvatures are aligned, which we model using Frank's free energy density for nematic liquid crystals:⁴⁶Here, ∇ represents the covariant derivative on the curved surface, n⃑ denotes the orientation of the inclusion, and t⃑ signifies its perpendicular direction in the local tangent plane. The constants and correspond to the elastic constants governing in-plane nematic interactions. The variable r_v takes a value of 1 if a CMC is present on a vertex; otherwise, it is 0.

A discrete form of this energy is employed to facilitate implementation in the Monte-Carlo (MC) simulations. If we assume a one-constant approximation eqn (4) can be rewritten in a way that makes the implementation suitable for MC simulations (Lebwohl–Lasher model):⁴⁷

Here, is the strength of the nematic interaction. The sum is over all the nearest neighbor (i, j) vertices on the triangulated grid, promoting alignment among the neighboring orientation vectors. Here, k is the summing index in a mesh of N vertices. An even simpler form of this approximation for the in-plane orientational field is the XY model on a surface:⁴⁸

Here, w represents the strength of the direct interaction constant, and the summation runs over all protein–protein pairs. The sum of E₁ and E₂ is minimized numerically, while E₂ is given by either isotropic (eqn (3)) or nematic binding (eqn (6)). Isotropic interaction reflects the tendency of CMCs to self-aggregate, while nematic interaction also describes the tendency to align CMCs' principal curvatures. All CMCs in this work are either arc or saddle-shaped (see Fig. 1(a)).

The details of the Monte-Carlo procedure are given in Appendix A and the details of the mesh generation and finding the principal curvatures are given in in Appendix B.1.

2.4. Reduced volume

A volume constraint on the vesicle shapes can be implemented by adding an energy term:Here, the v̄ is the target reduced volume that can vary from 0 to 1. A sphere has a reduced volume of 1. A and V are the area and the volume of the vesicle, respectively. The volume modulus k_v was determined to be of the same order of magnitude as the bending and bonding terms. The sum of E₁, E₂ and E₃ is minimized numerically.

2.5. Order parameter

We want to quantify the degree of order between neighboring anisotropic CMCs. For this reason we can define the nematic order of the inclusion at vertex i bywhere θ measures the angle between directors (or principal curvature C_1m) of neighboring CMCs on the membrane. The more aligned the neighboring CMCs are, the more this value approaches unity. In our numerical approach it is calculated by a Python script: for every CMC occupied vertex, it identifies neighboring vertices and uses 8 through a dot product for the ones which are also occupied. An arithmetic average of these is then calculated and recorded in a value S_i. The average nematic order across the entire membrane is then

2.6. Mean cluster size and the gyration tensor

We can define a mean cluster size 〈N〉 for simulations where the vesicles are not fully covered by CMCs. If we index all clusters so that i has N_i vertices, 〈N〉 is calculated as:The mean cluster size is not a reliable way to differentiate between phases. While all phases contain large clusters, their structure and organization vary significantly. This measure is heavily influenced by the total number of clusters, giving too much weight to small, single-vertex clusters. As a result, it introduces too much noise to effectively separate most phases.

To clearly distinguish phases in which CMCs form large condensed clusters, we instead rely on morphological measures. The shape of the vesicle is captured through the eigenvalues of the gyration tensor, λ_i². The gyration tensor⁴⁹ is defined as the average over all the vertices, with respect to the center of mass, similar to the moment of inertia tensor for equal-mass vertices:

This can be visualized by a unique ellipsoid which has the same gyration tensorThe eigenvectors (e_i) of the gyration tensor are the directions of the semi-axes of the equivalent ellipsoid and the eigenvalues are their length squared divided by 3, ordered by their size: λ₁² ≤ λ₂² ≤ λ₃².

3. Materials and methods

The simulations were run using trisurf-ng.^45,50 Briefly, trisurf-ng simulates the vesicle as a triangulated, non-interacting mesh in 3D, where each vertex in the network can either be bare or occupied with a single CMC with an intrinsic curvature. The CMCs are mobile on the membrane between consecutive simulation steps, which is enabled through random vertex translations and random edge flips at each simulation step (for details, see Appendix A). A single CMC can therefore diffuse freely over the membrane until its energy contribution to the total energy of the system is minimized. All the parameters of the simulation are included in the so-called tape file. In the tape file, all the parameters are set prior to each simulation, including the size of the mesh (given by the integer nshell, where N = 5·nshell² + 2 and N is the number of vertices in a mesh, while the number of triangles is 2N), CMC fraction, curvature and volume constraints. The fraction of CMCs on the membrane is set at the beginning of the simulation and is kept constant ρ = N_CMC/N throughout. Each vertex is assigned one-third of the area of its adjacent triangles; therefore, the CMC vertex density is equivalent to its area density. All CMCs are initially assigned to the membrane randomly. In a typical simulation, nshell = 20 (N = 2002), mcsweeps = 10⁵, iterations = 500 (MC steps) (see Appendix A). Simulations were run in parallel, each on a separate core of a 20-core PC, with one needing usually 16 hours to finish. The resulting VTU files were viewed and prepared for figures publication in ParaView, but further analysis and graph generation were done by separate Python scripts. The bending rigidity of the simulated membranes is 20k_BT throughout. The curvature of the CMCs throughout the paper is given in units of 1/l_min, which means that the spontaneous curvature C_1m = 0.25 corresponds to the curvature of the sphere with radius 4l_min. Similarly, all energy values throughout the paper are given in units of k_BT, where k_B is the Boltzmann constant and T absolute temperature.

The simulations were run using trisurf-ng.^45,50 Briefly, trisurf-ng simulates the vesicle as a triangulated, non-interacting mesh in 3D, where each vertex in the network can either be bare or occupied with a single CMC with an intrinsic curvature. The CMCs are mobile on the membrane between consecutive simulation steps, which is enabled through random vertex translations and random edge flips at each simulation step (for details, see Appendix A). A single CMC can therefore diffuse freely over the membrane until its energy contribution to the total energy of the system is minimized. All the parameters of the simulation are included in an input file, including the size of the mesh (given by the integer nshell, where N = 5·nshell² + 2 and N is the number of vertices in a mesh, while the number of triangles is 2N), CMC fraction, curvature and volume constraints. The fraction of CMCs on the membrane is set at the beginning of the simulation and is kept constant ρ = N_CMC/N throughout. All CMCs are initially assigned randomly over the membrane. In a typical simulation, nshell = 20 (N = 2002), mcsweeps = 10⁵, iterations = 500 (MC steps) (see Appendix A). Simulations were run in parallel, each on a separate core of a 20-core PC, with each one needing usually 16 hours to finish. The resulting VTU files were viewed and prepared for publication using ParaView, but further analysis and graph generation were done by separate Python scripts. The bending rigidity of the simulated membranes was set to 20k_BT throughout. The curvature of the CMCs throughout the paper is given in units of 1/l_min, which means that the spontaneous curvature C_1m = 0.25 corresponds to the curvature of the sphere with radius 4l_min. Similarly, all energy values throughout the paper are given in units of k_BT, where k_B is the Boltzmann constant and T absolute temperature.

4. Results

4.1. Fully occupied vesicles with no volume constraints

We first investigate the steady-state shapes of vesicles completely covered with CMCs (ρ = 1) without imposing volume constraints. Although 100% CMC coverage is experimentally impossible, calculating this theoretical limit is valuable for understanding the maximum potential of shapes that are nearly or fully covered by CMCs.

4.1.1. Arc-shaped CMCs promote the formation of elongated cylindrical vesicles (H_m = D_m > 0). The binding energy-spontaneous curvature phase diagram for vesicles fully covered with arc-shaped CMCs, without volume constraints, is shown in Fig. 2. Direct nematic interaction energy w and intrinsic curvature of CMCs increase along the x and y axes, respectively. Generally, equilibrium shapes are cylindrical with the radius of the cylinders determined by the curvature of the CMCs. In the case of no interaction w = 0 and low spontaneous curvature, (bottom left corner of Fig. 2), all CMCs have random orientations, resulting in a sphere. The heat map of Fig. 2 shows the energy density of each shape. When the neighbor interaction is negligible (w = 0), the energy is predictably positive due to membrane bending (the leftmost column of Fig. 2), while an increase in w increases the nematic order and lowers the overall energy.

Fig. 2 (a) Vesicles fully covered (ρ = 1) with arc-shaped CMCs as function of the nematic interaction energy w and intrinsic mean curvature H_m. Steady-state shapes are generally all cylinders (shapes (B) and (C)). The pearling steady-state shapes in (A) arise as a consequence of neighboring CMCs assuming random orientations. Even in the absence of nematic interactions between neighboring CMCs, the membrane conforms to their spontaneous curvature. The heat map gives the total energy of each shape per vertex (E₁ + E₂ from eqn (2) and (6)), which decreases with w. (b) A heat map of the average nematic order 〈S〉 reveals that it increases with w, but is most pronounced when CMCs are less curved. The shapes characterized by junctions or necks have a low degree of nematic order (the pearling phase), as seen in the close up of Fig. 5(A). Orange lines show the direction of the principal CMC curvature C_1m.

Fig. 3 shows the convergence of the bending energy and volume for three shapes of Fig. 2, as function of MC steps. We see that convergence is achieved in all cases, but at different rates. At H_m = 0.5 and w = 0, the vesicle forms a shape composed of smaller spheres joined together to form a series of pearls connected by thin necks, the pearling phase (Fig. 2(a)). In the neck regions, S_i is low. An gradual increase in w from zero leads to fewer necks as the neighboring CMCs align to form adjoined smaller cylinders with large nematic order 〈S〉 (Fig. 4). Fig. 5 shows close-ups of the steady-state shapes in Fig. 2 with the color map indicating the local nematic order parameter S_i (eqn (8)) and the average nematic order of the shapes 〈S〉 (eqn (9)).

Fig. 3 Convergence of total energy and reduced volume (eqn (2) and (6)) the evolution of inset shapes labeled (A), (B) and (C) (in Fig. 2) are shown in red and green, respectively. Convergence is achieved within 500 MC steps, which was taken as a default number of MC steps for each simulation.

Fig. 4 (a) For vesicles that are fully covered with arc-shaped CMCs, an increase in nematic strength w from 0 to 0.8 results in an evolution of the steady-state shapes from pearls to cylinders. The average nematic order 〈S〉 (eqn (9)) is given below each shape. (b) The pearls-to-cylinders transition is accompanied by a steady decrease of not only bending, but also total energy. The parameters are the same as in panel (A) of Fig. 2(a).

Fig. 5 A closer look at CMC alignment and the local nematic order (eqn (8)), shown as a heat map for the shapes labeled (A), (B) and (C) shown in Fig. 2. The average nematic order 〈S〉 (eqn (9)) is given below each shape.

The bottom row of Fig. 2(a) reveals that for spontaneous curvature H_m = D_m = 0.25 and increasing w, the steady-state shapes change from a sphere to a prolate state of different radii, in a nonmonotonous manner. This process is analyzed in greater detail in Fig. 6. Each prolate shape is approximately an ellipsoid that is characterized by two axes along its symmetrical planes. As the area is kept constant, an increase in w from zero results in an elongation of the axis in one direction and a narrowing in another (Fig. 6(a)–(c)). As w increases from zero, the CMCs begin to align with each other, transforming the initial spherical shape into an elongated prolate phase (Fig. 6(c)), resulting in a decrease of bending energy and reduced volume (Fig. 6(d) and (e)). The in-plane rotation accompanying this change is reflected in the decrease in the average angle between neighboring CMCs (Fig. 6(f)). 〈S〉 monotonically increases with increasing w (Fig. 6(g)). Since the CMCs do not perfectly align perpendicular to the axis of the tube, the tube has a radius that is smaller than the spontaneous radius of curvature of the CMC. As w increases, the CMCs align more perfectly with respect to each other, and more perpendicular to the tube axis, which increases in radius until it equals the radius of the CMC for large w.

Fig. 6 (a) Steady-state equilibrium shapes that are fully covered (ρ = 1, H_m = D_m = 0.25) by arc-shaped CMCs result in prolate phases for w > 0. For a prolate, the major and minor axes can be determined. At w = 0, the neighboring CMCs do not interact and have no effect on the orientations of their neighbors, resulting in a spherical shape with both major and minor axes having the same length (b) and (c). Increasing w forces the neighboring CMCs to align with each other, which leads to elongation of the vesicles (b) and (c). The bending energy per vertex has a minimum at w = 1, before starting to increase due to stronger bonding (d). The elongation of the vesicles is most pronounced when w = 1.5 and v = 0.56 (e). The average angle between neighboring vertices decreases as function of w (f), while the inclusion average order increases with w (g). When w is larger than 1.5, the vesicles become less elongated and tend towards the spherical shape for large values of w (c).

4.1.2. Saddle-like CMCs promote global saddle shapes (H_m = 0, D_m > 0). Saddle-like CMCs are inherently frustrated on the convex surface of the vesicle. For these CMCs, H_m = 0 and D_m > 0. The shape of the saddle of the CMCs is determined by the curvature of the hyperbolic paraboloid which is in turn defined by the deviatoric component D_m; a flat saddle corresponds to D_m = 0, while curved saddles have non-zero values D_m. Independently varying D_m and w leads to the phase diagram shown in Fig. 7. For non-interacting CMCs, w = 0, vesicles remain spherical, as in the case of arc-shaped CMCs, due to their random orientations. By contrast to arc-shaped CMCs, there is no spontaneous budding even as the spontaneous curvature increases. As w increases, the vesicles deform from the spherical shape. In the case of D_m = 0.75, the vesicle develops globally flattened regions of the membrane with increasing w (see bottom row of Fig. 7). In the case of most curved saddles (D_m = 1.2), the vesicle features a bow-tie phase. (see top row of Fig. 7). The heat map shows that the energy density decreases with the increase of w, consistent with observations for arc-shaped CMCs (Fig. 2). Fig. 8 shows the local nematic order for the steady-state shapes in Fig. 7(b), with saddle-like CMCs and an increasing interaction w. Larger 〈S〉 leads to more localized nematic defects, shown in blue.

Fig. 7 (a) The curvature-binding strength phase diagram for steady-state shapes with saddle CMCs (ρ = 1). The curvature of the hyperbolic paraboloid is given on the y-axis by D_m (H_m = 0) and w on the x-axis. When saddle CMCs are more flat (D_m = 0.75), the vesicles mirror this by exhibiting regions of flat membrane (top row), while more curved saddle CMCs (D_m = 1.2) result in global shapes that mimic the inclusion's curvatures (bottom row). The energy per vertex is shown by the corresponding heatmap. (b) The average nematic order heatmap (eqn (9)) for the curvature-binding strength phase diagram shows a gradual increase with increasing w.

Fig. 8 The heat map of the local nematic order (eqn (9)) for the steady-state shapes of vesicles covered by saddle-shaped CMCs, shown in Fig. 7(b). The nematic defects characterized by minimal S_i are shown in blue. When there is no interaction w = 0, the steady-state is a sphere (A) and a bowtie-like morphology for w = 1 (B). The inset shows the nematic defect field, which is not fully developed. (C) The nematic defect is most notable when w = 3, at the tip of the elongated membrane. Its topological charge is +1, while there are two additional charges +1/2 on either side of the wider part (inset).

4.2. Partly covered vesicles with volume constraints

We next investigate the steady-state shapes of vesicles that are only partially covered with CMCs, while also imposing a fixed volume constraint. The CMC and bare regions of the membrane are shown in blue and gray, respectively.

4.2.1. Phase diagram for vesicles half covered with arc-shaped CMCs (H_m = D_m = 0.25). First we explore the phase diagram of steady-state equilibrium shapes for ρ = 0.5, H_m = D_m = 0.25 in the v–w plane (Fig. 9). We identify several phases: oblate, prolate, and capped, with the boomerang and dumb-bell phases as subsets of the prolate phase, and the mixed phase as a subset of the oblate phase (see Fig. 9(a)). The phases can be roughly discerned by the eigenvalues λ_i² of the gyration tensor (eqn (11)), shown in Fig. 9(b)–(d). The first eigenvalue λ₁² measures the thinness and is low for oblate shapes. The second eigenvalue λ₂² is large for the oblate shapes (and roughly equal to the largest eigenvalue λ₂² ≈ λ₃²), but is minimized for elongated prolate shapes (and roughly equal to the largest eigenvalue λ₂² ≈ λ₁²). The third eigenvalue λ₃² is largest for shapes that are elongated along one principal axis and characterizes the prolate phase. Since the separation between boomerang, dumb-bell and prolate phase is more qualitative than quantitative, as they are all marked by an elongated axis, the heatmap for λ₃² shows highest values for all three. Fig. 9(e) shows the average cluster size 〈N〉 (eqn (10)). The oblate phase is characterized by CMC clustering on the rim to form discs (as in the case of isotropic inclusions⁵¹). An increase in v results in the rim of CMCs not closing up entirely, while still keeping their approximate oblate shape. Prolate shapes are limited to values of w ≤ 1 and contain two sub-phases: the boomerang and the dumb-bell. The oblate phase transitions into the mixed phase approximately above v = 0.65, and the prolate below w = 1.5. The mixed phase is characterized by two effects: stronger CMC clustering due to high w and reduced flatness due to large v. The average nematic order 〈S〉 (eqn (9)) is shown in Fig. 9(f). Predictably, 〈S〉 increases with increasing CMC alignment and w. The capped phase is located at high v and w. Here, the clusters form separate patches, sometimes with a pronounced “cap” at the tips of shape protrusions. The patch separation results in a lower 〈N〉, but better nematic order of CMCs with higher 〈S〉, both discernible on the heatmaps in Fig. 9(e) and (f). Below w = 1 the average cluster size is ≤100 (Fig. 9e) and the CMC distribution is dominated by mixing entropy and smaller 〈N〉. A comparison of shapes with the same interaction w but no volume constraint is shown on the right of Fig. 9. At w = 0 and 1, the shape is roughly spherical, then changes to the mixed and capped phases at w = 2 and w = 3, respectively.

Fig. 9 Reduced volume-binding strength plane for arc-shaped CMCs. The CMC and bare regions of the membrane are shown in blue and gray, respectively, while the orange lines show the principal curvature direction C_1m. (a) Phase diagram as function of v and w for ρ = 0.5 CMC concentration for H_m = D_m = 0.25. The different phases are indicated by their names and a typical snapshot of the equilibrium shape is shown. The transition lines between the phases were drawn according to the measures shown in the bottom panels. The oblate phase is separated by the small eigenvalue λ₁², as explained in the main text (b). The prolate phase is roughly indicated by a small intermediate eigenvalue λ₂² (c). The prolate phase is indicated with largest eigenvalue λ₃², but the boomerang and dumb-bell phases' transition is difficult to determine from λ_i² (d). The average cluster size heatmap shows an approximate separation between the oblate and the capped phases (e). When v is increased above 0.7, the CMC patch on the rim of the oblate phase segregates into distinct CMC patches. This also results in the capped phase having the largest average nematic order 〈S〉 (f). The column right of panel (a) shows that a decrease in v is slight with no volume constraints and increasing w, resulting in a transition that goes straight from spherical to capped phases.

4.2.2. Arc-shaped CMCs (H_m = D_m = 0.25) accelerate the oblate–prolate shape transition in comparison to bare-membrane vesicles. We plot the v–E diagram of vesicles at steady-state (Fig. 10(a)), that contain no CMCs, in agreement with existing literature.⁵² Starting from v = 0.3 and increasing v the steady-state shapes follow a familiar pattern of the stomatocyte, oblate and prolate phases (Fig. 10(a)). The phases are discerned from the eigenvalues of the gyration tensor λ_i². There is overlap where the oblate and prolate phases coexist, roughly in the range v = 0.6–0.8. However, in this region, E is lower for prolates, which extend to v = 1. We next examine how the steady-state vesicle shapes change when they are half covered with arc-shaped CMCs (ρ = 0.5), and how the v–E phase diagram compares to the bare membrane vesicles. The results are shown in Fig. 10(b). We now find that the oblate phase has been pushed to exist only below v ≈ 0.5, while the prolate phase is now extending over a larger range of reduced volume (compared to Fig. 10(a)).

Fig. 10 (a) The phase diagram in the space v–E for vesicles with no CMCs shows a familiar sequence of morphologies with increasing v; stomatocytes, oblates and prolates. An interesting intermediate shape transition (ellipsoid) is observed between oblate and prolate shapes, marked in green, which was reported in literature.⁵³ (b) The phase diagram in the v–E space for arc-shaped CMCs (ρ = 0.5, H_m = D_m = 0.25, w = 0). In comparison to empty vesicles, there is no stomatocyte phase, but an intermediate boomerang phase (shown in green) which is a subset of a prolate phase. The transition to prolates is accelerated and happens at around v = 0.5, possibly due to entropic mixing of the CMCs, which leads to their migration away from the rim of the oblate vesicles. The CMC and bare regions of the membrane are shown in blue and gray, respectively, while the orange lines show the principal curvature direction C_1m.

The smaller regime of oblate stability is driven by the lower mixing entropy of the CMCs on the oblate compared to the prolate shape. The oblate phase has two flat sides where the arc-shaped CMCs face a bending energy penalty, which is reduced when the CMCs accumulate on the curved rim. This is not the case for prolates, which have more uniform curvature. We can observe the effect of mixing entropy by plotting the distribution of CMCs as function of distance from the center of mass (COM) of the vesicle (Fig. 11). When w is increased, the CMC distribution is more pronounced on the rim of the oblate shapes, lowering their entropy (see SI 8.3.4).

Fig. 11 In the case of vesicles which are half-covered by arc-shaped CMCs, the oblate–prolate phase transition includes a rearrangement of the CMC distribution measured from the center of mass of the vesicle (COM). Oblate vesicles have most of their CMCs located on the rim (a), while in prolate vesicles they are more evenly distributed (b). The CMC and bare regions of the membrane are shown in blue and gray, respectively, while the orange lines show the principal curvature direction C_1m.

4.2.3. Phase diagram for vesicles half covered with saddle-like CMCs (H_m = 0, D_m = 0.98). We plot the v–w phase space for the steady-state shapes of vesicles that are half-covered by saddle-like CMCs with H_m = 0 and D_m = 0.98 (Fig. 12). The CMC and bare regions of the membrane are shown in red and gray, respectively. Saddle CMCs with negative Gaussian curvature are inherently frustrated on convex membranes. The phase diagram features a high degree of metastability, resulting in approximate phases that often coexist for the same v and w. This is also evident in the heatmaps of λ_i², which cannot reliably distinguish between phases (not shown). The poor reproducibility of these transitions warrants further investigation but lies beyond the scope of this paper. The simulations reveal four phases with ambiguous transition boundaries: the invaginated, oblate, neck, mixed, and prolate phases.

Fig. 12 Reduced volume-binding strength plane for saddle-like CMCs. (a) Phase diagram as function of v and w with ρ = 0.5 concentration for H_m = 0 and D_m = 0.98. Invaginated steady-state shapes are common under v = 0.5 with the majority of CMCs found on the invaginated rim. The oblate shapes are limited to low w and regions approximately between v = 0.5–0.6. An elongated neck-like formation is typical at w = 3, but the elongation becomes less pronounced at higher v. In comparison, steady-state shapes with no volume constraints are shown right of panel (a). The average cluster size heatmap shows largest clusters at w = 3, approximately where neck formation takes place (b). Average nematic order gradually increases with increasing w (c). The transitions between the phases are less defined than in the case of arc-shaped CMCs, with oblate and invaginated phases often having very similar equilibrium energies. The CMC and bare regions of the membrane are shown in red and gray, respectively.

The invaginated phase is found at v < 0.4, with most CMCs localized on the saddle-curved rim of the invagination. The oblate phase appears for w < 1 and v = 0.5–0.6. The prolate phase is restricted to v = 0.8 and w < 0.5. The neck phase exhibits slender cylindrical protrusions with nematically ordered CMCs, flanked by convex, CMC-free regions. Notably, this phase persists even without volume constraints at w = 2–3 (see Fig. 12(a), right column). Heatmaps of the mean cluster size 〈N〉 and mean nematic order parameter 〈S〉 (Fig. 12(b) and (c)) show a monotonic increase with w, where the neck phase corresponds to elevated 〈N〉.

At a lower CMC density (ρ = 0.16), steady-state shapes display CMC aggregation in neck regions between convex membrane segments, resembling a pearling phase (Fig. 13(a)). For ρ = 0.5 and w = 3, saddle CMCs exhibit strong nematic alignment (Fig. 13(b)).

Fig. 13 (a) A small concentration of anisotropic saddle CMCs (shown in red) (H_m = 0, D_m = 0.98, ρ = 0.16, w = 3, v = 0.54) clusters together to form necks between empty convex membrane regions (grey). Orange lines show the direction of the principal curvature C_1m of the CMCs. (b) The nematic order heatmap and orientation of saddle-like CMCs for the case of H_m = 0, D_m = 0.98, ρ = 0.5, w = 3, v = 0.61. The heat map in both panels shows local nematic order.

4.2.4. Analysis of topological defects. A topological defect (TD) on a surface arises when the order parameter S_i of the inclusions cannot be smoothly defined everywhere, leading to singular points or lines where the order parameter abruptly changes. TDs are characterized by their discrete topological charge, an additive and conserved quantity.⁵⁴ TDs with like charges repel each other, whereas those with opposite charges attract. Given these similarities to electrostatic interactions, some studies^55–57 suggest an electrostatic analogy to describe the interplay between Gaussian curvature and TD configurations. TDs with positive (negative) value of topological charge are referred to as defects (antidefects).^55,58 On 2D surfaces, the topological charge corresponds to the winding number of the nematic director field.^54,59 According to the Gauss-Bonnet and Poincaré-Hopf theorems,^60,61 TDs must exist in all non-toroidal topologies. These theorems dictate that the sum of all topological charges (the total winding number) equals 2 for 2D surfaces of spherical topology. In nematic ordering, the topological charge can be a multiple of half an integer.

It is well-established that Gaussian curvature strongly influences the location and number of TDs.^55–57 Regions with positive (negative) Gaussian curvature attract TDs of positive (negative) topological charge. While TDs are energetically costly—leading systems to avoid them, often through the annihilation of defect-antidefect pairs into locally defect-free structures⁵⁵—their presence is frequently unavoidable due to the system's topology.

The occurrence of TDs is most apparent when CMCs cover the entire membrane surface and can be analyzed by studying S_i and variations of C_1m orientation (Fig. 2–8). In the case of vesicles fully covered with arc-shaped CMCs, the TDs of positive charge are normally found on the cap of the cylindrical shapes of the prolate phase where Gaussian curvature is positive. An example of four defects with charge +1/2 can be seen in Fig. 14(a), which is a close-up of the TD found on the cap of the prolate phase in Fig. 6(c). Alternatively, when saddle-like CMCs fully cover the membrane, a single defect of charge +1 is found on the tip of the protrusion and two defects with charge +1/2 on either side of the base (Fig. 14(b)). They are positioned on opposite sides to maximize separation, consistent with the repulsion between like-charge defects. In both these cases, the topological charge sums up to 2, in line with the Gauss-Bonnet and Poincaré-Hopf theorems.

Fig. 14 (a) For a fully covered equilibrium steady-state shapes with arc-shaped CMCs (H_m = D_m = 0.25, w = 5), the four defects of charge +1/2 are found on the cylindrical caps, where the Gaussian curvature is positive. (b) In the case of a membrane that is fully covered by saddle-like CMCs (D_m = 0.98, H_m = 0, w = 3), two +1/2 defects are found on either side of the bulbous base (the other defect lies adjacently on the other side and is not shown) and a +1 defect on the tip of the protrusion. The defects are characterized by a locally low nematic order S_i, marked by blue regions.

The electrostatic analogy for TDs holds best when the intrinsic (direct interaction) term dominates.^56,57 However, extrinsic^62,63 or deviatoric^38,64 terms may also play a role, particularly on surface patches where principal curvatures differ. These terms can impose an external ordering field, favoring CMC orientations that best conform to the surface geometry. Consequently, intrinsic and deviatoric/extrinsic terms may introduce competing tendencies.⁶³

We compare the result shown in Fig. 14(b) to simulations of axisymmetric studies of closed membrane shapes with saddle-like CMCs and find good agreement with the equilibrium steady-state shape and with the position of defects (Fig. 15).

Fig. 15 Equilibrium orientational ordering profile on a fixed axisymmetric 2D surface that is topologically equivalent to the shape presented in Fig. 14(b). The nematic ordering amplitude λ is denoted by the color code, while the local orientation of molecules is presented by the lines in the (ϕ, s)-plane, where ϕ is the azimuthal angle of the axisymmetric surface and s the arc length of the profile curve characterizing the axisymmetric surface. The equilibrium nematic ordering amplitude is denoted by λ₀ and the total length of the profile curve by L_s. Topological defects are marked with capital letters. The topological charge of TD A, B and C is +1, +1/2 and +1/2, respectively. The surface is completely covered by saddle-like CMCs with H_m = 0 and D_m = 0.98. The details of the modeling used in this calculation are described in detail in ref. 58. The only difference in the modeling is that we used the deviatoric term for 2D inclusions instead of the term for 1D inclusions, which was used in ref. 58.

In our example (Fig. 14(b) and 15), the deviatoric term (eqn (2)) exerts its strongest influence along the tubular protrusion, where the difference in principal curvatures (i.e., the deviatoric curvature) is greatest. This term aligns the CMCs parallel to the protrusion, minimizing frustration. In contrast, the less deviatoric regions—particularly the nearly spherical base in Fig. 15—remain unaffected by this term, as the surface is isotropic (equal principal curvatures). The deviatoric term also has no effect at the highly curved, isotropic tip of the protrusion.

Since the TDs in Fig. 14(b) and 15 appear in regions of high Gaussian curvature, we conclude that intrinsic and deviatoric/extrinsic terms do not conflict in this case. The deviatoric term's influence is evident along the tubular protrusion, where it enforces CMC alignment. In contrast, considering only the intrinsic term could yield configurations with CMCs tilted relative to the protrusion.

5. Discussion

Our results demonstrate that anisotropic curved membrane components (CMCs) significantly influence vesicle morphology through curvature sensing and nematic alignment. By focusing on the two simplest anisotropic CMCs – the arc-shaped and saddle-shaped – we find a wider variety of steady-state phases compared to isotropic CMCs, where only convex budding phases were found across all concentrations of CMCs.^65,66 By mapping the steady-state shapes of fully covered vesicles for both arc- and saddle-shaped CMCs, we found that in the former case, there exists a competition between nematic ordering and curvature; when the former is absent, equilibrium steady-state shapes are a series of connected pearls (Fig. 2). Nematic binding tends to break up the pearling phase and results in smooth prolate phases, with the radius of the cylinder determined by the intrinsic curvature of the CMCs. With this, we have confirmed pearling using our non-axisymmetric numerical scheme and connected it to previous studies originating from analytic and theoretical considerations.^38,67 Additionally, these findings align with experimental observations of BAR domain proteins (e.g., amphiphysin and IRSp53) stabilizing tubules or negatively curved membranes.^1,3,68 The pearling-to-cylinder transition at low nematic strength mirrors the budding and tubulation processes seen in cellular membranes, suggesting that weak interactions may suffice for initial curvature sensing, while stronger alignment drives large-scale shape changes.

We next investigated how nematic ordering and volume constraints govern phase behavior. The interplay between nematic interaction strength and reduced volume reveals distinct phases (Fig. 9 and 12). A membrane half-covered with arc-shaped CMCs can exhibit both prolate and oblate phases, depending on the binding strength between CMCs; the prolate and oblate phases occur at low and high nematic interaction strengths, respectively. Tubular structures stabilized by arc-shaped CMCs are well documented in numerous numerical studies,^22,69–71 and are also supported by experiments with giant unilamellar vesicles,⁷² curvature-stabilizing proteins in the endoplasmic reticulum (ER),⁷³ and in the daughter vesicles of the erythrocyte membrane.⁷⁴ We found that arc-shaped CMCs, in the absence of neighbor interactions, accelerate the oblate-to-prolate transition compared to bare membranes where this transition is driven only by volume constraints (Fig. 10), highlighting how CMC entropy penalizes flat vesicle regions. This entropy-driven effect is consistent with the clustering of CMCs on curved rims (Fig. 11), resembling lipid raft aggregation in Golgi cisternae¹⁷ or pore stabilization.^75,76 In contrast, saddle-like CMCs exhibit metastability and neck phases (Fig. 12), similar to tubulation phenomena.³⁶ The persistence of protrusions even without volume constraints (Fig. 12(a)) underscores the role of saddle-like anisotropic CMCs in stabilizing connecting necks with negative Gaussian curvature.

At concentrations below 10% CMC and high interaction strength, we found that saddle-shaped CMCs preferentially accumulate in the elongated necks separating convex empty membrane regions. This observation agrees with previous studies^77,78 suggesting that vesicle necks contain high concentrations of anisotropic membrane components with non-zero deviatoric curvature. The accumulation of these components in membrane necks correlates with strong in-plane nematic ordering, occurring where the difference between principal membrane curvatures (and consequently, Gaussian curvature) is maximized. This was also observed in simulations and in vitro experiments of Le Roux et al.⁶⁸ The authors of the work modeled the membrane-bound N-BAR proteins as being elliptical with a concave isotropic curvature (invaginated) in contrast to our anisotropic, saddle-shaped CMCs.

A similar result was obtained by simulations of Noguchi et al.,⁷⁹ where the authors showed that an assembly of curved rod-shaped proteins assemble in the necks and may play an important role in the formation of neck-like structures during cell division and membrane budding. The difference between rod- and arc-shaped CMCs is in the fact that the former have one and the latter have two intrinsic curvatures, respectively.

Additionally, such neck-like phases are observed also when the deviatoric term is increased in the absence of direct interactions. This leads us to conclude that the pearling phase is a ubiquitous process which always requires some CMCs to develop, but the ramifications of our simulations could provide insight into the structure of pearling or neck-like phases simply by observing the width of the neck and/or the pearls themselves. Conversely, we found that in the case of arc-shaped CMCs the average mean curvature H_m controls the radius of the pearls and that nematic order in the thin necks connecting the bare parts of the membrane is zero.

We investigated topological defects (TDs) for both anisotropic CMC types. Our results show that TDs consistently localize to regions of high Gaussian curvature (Fig. 14 and 15). The +1/2 defects on cylindrical caps and +1 defects at protrusion tips (Fig. 14) match theoretical predictions,^55,56 with their total topological charge satisfying the Poincaré-Hopf theorem. Along tubular protrusions, the deviatoric curvature term dominates, enforcing CMC alignment parallel to the tube axis, while intrinsic interactions primarily determine defect positioning. This dual behavior suggests that CMCs may utilize both curvature sensing and generation mechanisms in biological systems, as observed in ER tubules²⁰ and photoreceptor discs.¹⁹

Although area and volume constraints in combination with a nematic in-plane field have been studied before,^80,81 our systematic exploration of the phase diagram reveals the transition between distinct steady-state phases. Simulations of membranes with rod-shaped membrane-bound components by Bonazzi et al.⁸² revealed the formation of prolate phases every time the CMC concentration exceeded 40%, irrespective of their intrinsic curvature. We find that contrary to this observation, oblate steady-state shapes exist even at 50% CMC concentrations and up to 0.65 reduced volume, but only if the CMCs interact nematically. The fraction of CMC vertex coverage is directly comparable to its area coverage. This is because, in a mesh, each vertex accounts for one-third of the area of every triangle it touches, making the results comparable.

When considering nematic ordering on curved surfaces, extrinsic curvature effects play a significant role in aligning conventional (rigid rid-like) nematics along directions of least principal curvature.⁶² In our system, however, the molecules (CMCs) possess intrinsic curvature, necessitating a generalized approach. For 1D curved molecules, prior work⁵⁸ demonstrated that the extrinsic (deviatoric) bending energy drives alignment toward optimal surface fit rather than minimal curvature. Extending this to 2D CMCs, our coupling term (eqn (1)) encodes a similar interplay between extrinsic curvature and molecular shape, but with the alignment now governed by the compatibility between the intrinsic curvature of CMCs and the local surface geometry. This distinction highlights how our model inherently incorporates extrinsic effects while accounting for the additional complexity of the shape-anisotropy of the molecules.

Our findings are in line with the conclusions of ref. 83, who emphasized the importance of incorporating both shape and curvature anisotropy, as well as interaction potentials, in understanding protein sorting behavior. Our study demonstrates that curvature anisotropy and interaction strength enhance sorting efficiency, while shape anisotropy can counteract it. In our model, the in-plane nematic field was also researched in the context of non-convex CMCs. Moreover, in-plane coupling to membrane curvature can drive the emergence of complex structures such as tubes and discs, reminiscent of those induced by curvature-sensing proteins.⁸⁴ Additionally, our findings support and extend prior work showing that anisotropic curvature-inducing proteins, modeled as in-plane nematics, can drive membrane remodeling and aggregation via curvature sensing alone, without explicit self-interactions.⁸⁵ Specifically, the variation of the Gaussian modulus has previously been reported to affect shape changes even in the absence of direct CMC interactions, and its increase can facilitate neck formation even in axis-symmetrical 2D systems.⁸⁶ Our article highlights the necessity of incorporating anisotropic spontaneous curvature into membrane models to accurately capture vesiculation phenomena, as also demonstrated elsewhere.⁸⁷

Our simulations establish a framework for understanding how anisotropic CMCs (including BAR proteins and lipid rafts) shape cellular membranes. The boomerang phase (Fig. 9) and neck-driven elongation (Fig. 13) potentially model Vibrio cholerae cells or budding of yeast cells, respectively.⁸⁸ However, the metastable behavior of saddle-like CMC phases (Fig. 12) requires further investigation, particularly to distinguish between invaginated and elongated phases.

6. Conclusion

In this work, we systematically investigated how anisotropic curved membrane components (CMCs) govern vesicle morphology through curvature sensing and nematic alignment. By focusing on arc- and saddle-shaped CMCs, we uncovered a far richer spectrum of membrane shapes than previously observed with isotropic CMCs, where only budding was reported. Our simulations revealed that nematic ordering plays a crucial role in shaping vesicles: for arc-shaped CMCs, weak alignment allows pearling, while stronger alignment stabilizes smooth cylindrical phases. These findings confirm long-standing theoretical predictions and align with experimental studies of BAR domain proteins and tubule-forming systems.

We mapped the morphological transitions as a function of nematic interaction strength and reduced volume, demonstrating how these two parameters control the emergence of prolate, oblate, tubular, and metastable necked morphologies. Notably, arc-shaped CMCs can drive shape changes even in the absence of direct interactions, with entropy playing a key role in destabilizing flat membrane regions and facilitating neck formation. Saddle-shaped CMCs, on the other hand, induce and stabilize negatively curved necks, even at low concentrations and without volume constraints, reinforcing their role in neck formation observed in biological systems.

Our investigation of topological defects (TDs) and local nematic order provides a novel perspective on the coupling between curvature, nematic alignment, and vesicle topology. We also compared our curvature-energy framework with prior models, emphasizing our focus on mean and deviatoric curvature, which better aligns with biological observations of anisotropic curvature sensing. Our phase diagrams are, to our knowledge, the first to fully chart this behavior in the space of binding strength and reduced volume.

Altogether, our findings underscore the importance of anisotropic spontaneous curvature in modeling membrane remodeling and suggest a unifying framework for understanding how proteins such as those containing BAR domains, lipid rafts, and neck-stabilizing proteins shape cell membrane vesicle topology in vivo.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Our code is readily found on github as cited in the paper, but should the situation arise, we are also keen on sending you the original vtu files or python scripts upon request to the corresponding author.

Appendices

8. Appendix

A. Monte-Carlo procedure. The membrane is represented by a set of N vertices that are linked by variable length tethers l to form a closed, dynamically triangulated, self-avoiding two-dimensional network of approximately 2N triangles and with the topology of a sphere.^89,90 The lengths of the tethers can vary between a minimal and a maximal value, l_min, and l_max, respectively. Self-avoidance of the network is ensured by choosing the appropriate values for l_max and the maximal displacement of the vertex s in a single updating step.

One Monte-Carlo sweep (MCs) consists of individual attempts to displace each of the N vertices by a random increment in the sphere with radius s, centered at the vertex, followed by R_BN attempts to flip a randomly chosen bond. We denote R_B as the bond-flip ratio, which defines how many attempts to flip a bond are made per one attempt to move a vertex in one MCs. Note that the bond-flip ratio is connected to the lateral diffusion coefficient within the membrane, i.e. to the membrane viscosity. In this work we have chosen R_B = 3, s/l_min = 0.15 and l_max/l_min = 1.7. The dynamically triangulated network acquires its lateral fluidity from a bond flip mechanism. A single bond-flip involves the four vertices of two neighboring triangles. The tether connecting the two vertices in diagonal direction is cut and reestablished between the other two, previously unconnected, vertices. The self-avoidance of the network is implemented by ensuring that no vertex can penetrate through the triangular network and that no bond can cut through another bond.^65,91

B. Anisotropic code details.
B.1. Representing the membrane as a mesh. Our solver is called Trisurf. It models the vesicle as a closed, triangulated surface: a graph with vertices i ∈ V and edges e_ij ∈ E, and an auxiliary set of triangles t_ijk ∈ T. The triangles make the approximation of the surface, but it is the vertices which are the principal dynamical entities that hold the properties of the membrane (intrinsic curvature, membrane composition, nematic director, etc.) and move in space.

From the position of the vertices, x⃑_i we can compute the normal vector to each triangle N⃑_ijk and circumcenter O⃑_ijk (center of the circle containing the position of x⃑_i, x⃑_j, x⃑_k). This allows us to divide the triangle to six parts and assign two pieces to each vertex.

For vertex i, these are the sub-triangle between the vertex position x⃑_i, the triangle center O⃑_ijk, and the middle of one of the edges and a similar sub-triangle for the other edge ik. We can denote the vector for the side between the vertex and the edge middle as half the edge length e⃑_ij = x⃑_j − x⃑_i and the side between center and the edge middle, which is half of the dual (Voronoi) edge, as [small sigma, Greek, macron]ⁱ_jk (the other half of the voronoi edge is on the neighboring triangle σⁱ_{[small script l]j}). We can compute this σ from the circumcenter and the middle of x⃑_i, x⃑_j. This allows us to divide the triangle to six parts and assign two pieces to each vertex (Fig. 16).

With this, we assign an area A(i) and a normal N(i) for each vertex i, by running over the neighbors j − 1, j, j + 1… and the adjacent triangles (i, j − 1, j), (i, j, j + 1)…where the j vertices are the counterclockwise ordered neighbors of i.

Fig. 16 Schematic of a vertex and it's neighbors. The vertex i has neighbors, [small script l], j, k, …, which are connected by edges e_ij,…, where we explicitly show e_ij. Each triangle has a circumcenter O, and the area is assigned to each of the 3 vertices, which is shown for triangle ijk (red, green, and blue section). The total area assigned to i is shown in red across the entire cap, and the dual edge σⁱ across e_ij are shown σⁱ_jk and σⁱ_{[small script l]j}. The curvature of vertex i is essentially computed on the red area, for example, the normal of the vertex N(i) is the average of the local normal for every point in the red surface.

The fluidity of the surface is achieved by bond-flips, where a bond ij and the two triangles that share it is ijk, ji[small script l] are replaced by a cross bond k[small script l] and two triangles i[small script l]k, jk[small script l], which allows vertices to change neighbors.

B.2. Anisotropic curvature on the vesicle. We use the shape operator, which is a discreet version of the principal curvatures on a mesh, but we calculate it by projecting it onto a plane d (director) and t (tangent) to get the 2 × 2 matrix C(v). By this, we have 2H given by tr(C) and det(C) gives K.

To get the anisotropic bending energy of the surface, we use the method by ref. 48 to estimate the shape operator matrix S on each vertex, which represents the shape of the surface at the point. We then calculate the mismatch tensor M = S − C_m, where C_m is the intrinsic curvature tensor whose direction is determined by the director and the other tangent vector to the local normal (t̂ = N̂ × d̂).

where H_m = (C_1m + C_2m)/2, D_m = (C_1m − C_2m)/2 are the spontaneous curvature and spontaneous deviator at the vertex, respectively, which reflects the physical characteristics of local membrane composition. This is a change from the original method by Ramakrishnan et al.,³³ where the P_v was used as a matrix projection, and could be the detail that makes the calculations more robust and efficient.

The bending energy is calculated by inserting the mismatch tensor in the Hamiltonian

where K₁ and K₂ are the bending moduli of the vertex, again reflecting physical parameters due to local composition.

To calculate the shape curvature of a vertex i based on,⁴⁸ each edge ij is assigned a shape tensor estimation

S_ij = h_ijb⃑⊗b⃑ (18)

b⃑ is the binormal at the edge N̂_ij × e⃑_ij, where e⃑_ij = x⃑_j − x⃑_i is the edge vector and N̂_ij is the normal of the edge N̂_ij = (N̂_i,j−1,j + N̂_i,j,j+1)/|⋯| which is the sum of the normal of the two triangles sharing the edge, normalized. h_ij is a factor representing the directional derivative of the area ≈∇_pAwhere Φ is the dihedral angle (angle between the two triangle sharing the edge). Luckily there is a simple triple product formula for this factorThe cross product of the edge direction and a triangle normal gives a vector on the triangle which is perpendicular to the edge, which is at an angle from the edge normal.

The full vertex-shape tensor is a sum of the edge tensor, weighted by the match of the normal of the edge to the normal of the vertex.

We then project this 3 × 3 matrix in real space x̂, ŷ, ẑ to the tangent plane of the vertex d̂, t̂

The mean curvature H at the vertex is half the trace, while the Gaussian curvature K is the determinant, which are two of the degrees of freedom in the Hamiltonian.

The mismatch matrix can be calculated

The angle between the director and the eigenvectors of the shape matrix is what results in the ω angle dependence, which is the final degree of freedom. The mismatch tensor is simply inserted into eqn (17) to compute the bending energy of the vertex. If we integrate it across the whole membrane, we get eqn (1).

Let us shortly comment the main differences in numerical methodology. While Kumar and colleagues⁴³ project the principal curvatures—and their corresponding bending rigidities—along both principal directions in the tangent plane of each membrane vertex, we use their sum and difference, namely the mean curvature and mean deviator (and correspondingly, bending and splay stiffness). When comparing the two approaches, one might notice a discrepancy in the absolute value of the energies, but this difference is not important for the method itself, as the Monte Carlo steps and the convergence of both methods rely only on the energy difference between the former and latter states, ΔE. Our methodology is distinct because we explicitly apply the theoretical results posited nearly thirty years ago in the works of Kralj-Iglic and Iglic.^37–40 To make this example more concrete, let us directly compare our method with that of Kumar and colleagues⁸⁰ for a simple case of a tubulating membrane. To investigate the steady-state shapes of membranes fully covered by nematic CMCs, the method of ref. 80 uses c_‖ and c_⊥, while our method uses H_m and D_m. For an outward tubulation of the membrane, the first method uses c_|| < 0 and c_⊥ > 0, while our method uses H_m = D_m > 0. Similarly, for inward tubulation, the first method uses c_‖ > 0 and c_⊥ < 0, and our method H_m = D_m < 0. Of course, the length scales of the CMCs relative to the membrane size would need to be adjusted and matched in both cases. Anisotropic CMCs interact with both principal curvatures of the membrane, so in any case, the degrees of freedom are two; what differs is the parametrization.

C. Supplementary information.
C.1. Pearling transition up close. We study the transition from cylinders to pearl-like steady-state shapes by populating a previously empty cylinder with arc-shaped CMCs and running the simulation (Fig. 17).

Fig. 17 (a) Populating an empty cylinder with arc-shaped CMCs (H_m = D_m = 0.5) and setting w = 0 results in pearl steady-states shapes. (b) This transition is accompanied by a drop in bending energy, average nematic order and volume, but an increase in average deviator.

C.2. Pearling transition by varying H_m. At w = 0 and changing H_m, the steady-state shape changes from a sphere to a pearling state as shown in Fig. 19.

Fig. 18 With no nematic interaction (w = 0), but increasing K₂, the steady-states of vesicles form saddle-like necks between alternate convex regions and show perfect alignment between neighboring saddle-like CMCs at low concentrations. The x-axis is the constant K₂ in eqn (1), while 〈cos(2ω)〉 is the average orientation between neighboring CMCs. Parameters here are for simulations with H_m = 0, D_m = 0.98, ρ = 0.06, w = 0.

Fig. 19 The curvature of the arc-shaped CMCs determines the radius of steady-state for a fully covered membrane. (a) Without nematic interaction (w = 0), but varying curvature H_m, the steady-states of vesicles gradually change to pearls. No nematic interaction between CMCs ensures their random orientations on the membrane, leading to spherical and pearling morphologies for low- and high-curvature CMCs, respectively. (b) Total energy as a function of H_m shows that bending energy initially increases, but is reduced once the number of pearls exceeds 2.

C.3. Spontaneous ordering in necks when K₂ increases. For saddle-like CMCs at w = 0 and increasing K₂ (see eqn (1)), the steady-state shapes form saddle-like necks between empty convex membrane regions (see Fig. 13). This is due to the dominance of the deviatoric term (Fig. 18).
C.4. Arc-shaped CMCs cluster on the rim to form oblates even in absence of nematic ordering. Fig. 20 shows the total energy, average nematic order, average deviator and CMC distribution during the sphere-oblate transition for ρ = 0.5 of arc-shaped CMCs (H_m = D_m = 0.25) at v = 0.4 and three values of w. Vesicle energy (eqn (2)) converges to a minimum in all three cases, but is lower at greater values of w. Average nematic order is largest for w = 2 and lowest for w = 0 when no explicit nematic ordering is present.

Fig. 20 Oblate shapes v = 0.4 at ρ = 0.5 of arc-shaped CMCs (H_m = D_m = 0.25) and three values of w (from 0 to 2). Even without nematic interaction between CMCs (c), these are not distributed homogeneously over the membrane, but tend to accumulate on the rim, as is reflected in the inclusion distribution from the COM (row (c)). This slight ordering results in a positive average nematic order 〈S〉. Slight nematic ordering of w = 1 (b) results in the majority of CMCs being located at the rim (row (b)) and nematic order 〈S〉 increasing in the process. Deviator and radial dependence for H_m = D_m = 0.25. The CMCs profile from the centre of mass (COM) is made for the last state, namely at 500 MC steps. The CMC and bare regions of the membrane are shown in red and gray, respectively, while the orange lines show their principal orientation C_1m.

The average curvature deviator for each shape is calculated and shown in the third column of Fig. 20. The membrane deviator is defined for each vertex as D = (C₁ − C₂)/2, where C₁ and C₂ are its two principal curvatures. Dimensionless deviator is given by d = RD, where R is the radius of the sphere with the same area as the vesicle. The average deviator for each shape is calculated as .⁹² We find that the deviator slightly decreases with increasing w, as shown in the third column of Fig. 20; all shapes are oblates with varying degrees of flatness, with most flat regions corresponding to w = 2, where – due to strong nematic ordering – most of the CMCs are located at the rim (see the upper figure of the fourth column of Fig. 20). Note that a sphere would have d = 0.

Acknowledgements

Slovenian Research and Innovation Agency core founding No. P2-0232, P3-0388, projects No. J2-4427, J2-4447 and J3-60063, the University of Ljubljana interdisciplinary preparative project Nanostructurome 802-12/2024-5, and the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie Staff Exchange project “FarmEVs” (grant agreement No. 101131175). This publication is partially based upon work from COST Action COST-CA22153, supported by COST (European Cooperation in Science and Technology). N. S. G. is supported by the Lee and William Abramowitz Professorial Chair of Biophysics (Weizmann Institute) with additional support from a Royal Society Wolfson Visiting Fellowship (UK), and acknowledges support by the Israel Science Foundation (Grant No. 207/22). The views and opinions expressed in this publication are solely those of the authors and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible.

References

1. M. Simunovic, et al., When physics takes over: BAR proteins and membrane curvature, Trends Cell Biol., 2015, 25(12), 780–792.
2. C. Has and S. Lal Das, Recent developments in membrane curvature sensing and induction by proteins, Philos. Trans. R. Soc., A, 2021, 1865(10), 129971.
3. D. H. Johnson, et al., Protein–membrane interactions: sensing and generating curvature, Trends Biochem. Sci., 2024, 49(5), 401–416.
4. C. Has, P. Sivadas and S. Lal Das, Insights into membrane curvature sensing and membrane remodeling by intrinsically disordered proteins and protein regions, J. Membr. Biol., 2022, 255(2), 237–259.
5. J. Gómez-Llobregat, F. Elıás-Wolff and M. Lindén, Anisotropic membrane curvature sensing by amphipathic peptides, Biophys. J., 2016, 110(1), 197–204.
6. A. E. Hafner, J. Krausser and A. Šarić, Minimal coarse-grained models for molecular self-organisation in biology, Curr. Opin. Struct. Biol., 2019, 58, 43–52.
7. B. Schamberger, et al., Curvature in Biological Systems: Its Quantification, Emergence, and Implications across the Scales, Adv. Mater., 2023, 35(13), 2206110.
8. J. L. Gallop and H. T. McMahon, BAR domains and membrane curvature: bringing your curves to the BAR, Biochem. Soc. Symp., 2005, 72, 223–231.
9. M. Simunovic, et al., Curving cells inside and out: roles of BAR domain proteins in membrane shaping and its cellular implications, Annu. Rev. Cell Dev. Biol., 2019, 35, 111–129.
10. M. Simunovic, et al., Physical basis of some membrane shaping mechanisms, Philos. Trans. R. Soc., A, 2016, 374(2072), 20160034.
11. F.-C. Tsai, et al., Activated I-BAR IRSp53 clustering controls the formation of VASPactin–based membrane protrusions, Sci. Adv., 2022, 8(41), eabp8677.
12. M. Drab, et al., Inception mechanisms of tunneling nanotubes, Cells, 2019, 8(6), 626.
13. C. Prévost, et al., IRSp53 senses negative membrane curvature and phase separates along membrane tubules, Nat. Commun., 2015, 6(1), 1–11.
14. D. Kabaso, et al., Attachment of rod-like (BAR) proteins and membrane shape, Mini-Rev. Med. Chem., 2011, 11(4), 272–282.
15. P. J. Kundrotas and E. Alexov, Electrostatic properties of protein-protein complexes, Biophys. J., 2006, 91(5), 1724–1736.
16. M. Kandušer, et al., Effect of surfactant polyoxyethylene glycol (C12E8) on electroporation of cell line DC3F, Colloids Surf., A, 2003, 214(1–3), 205–217.
17. A. Iglič, et al., Coupling between vesicle shape and the non-homogeneous lateral distribution of membrane constituents in Golgi bodies, FEBS Lett., 2004, 574(1–3), 9–12.
18. D. Corbeil, C. A. Fargeas and W. B. Huttner, Rat prominin, like its mouse and human orthologues, is a pentaspan membrane glycoprotein, Biochem. Biophys. Res. Commun., 2001, 285(4), 939–944.
19. D. E. Chandler, et al., Intrinsic curvature properties of photosynthetic proteins in chromatophores, Biophys. J., 2008, 95(6), 2822–2836.
20. J. Guven, G. Huber and D. Marıá Valencia, Terasaki spiral ramps in the rough endoplasmic reticulum, Phys. Rev. Lett., 2014, 113(18), 188101.
21. H. Hägerstrand, et al., Endovesicle formation and membrane perturbation induced by polyoxyethyleneglycolalkylethers in human erythrocytes, Biochim. Biophys. Acta, Biomembr., 2004, 1665(1–2), 191–200.
22. V. Kralj-Iglič, et al., Stable tubular microexovesicles of the erythrocyte membrane induced by dimeric amphiphiles, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 61(4), 4230.
23. V. Kralj-Iglič, et al., Minimizing isotropic and deviatoric membrane energy – An unifying formation mechanism of different cellular membrane nanovesicle types, PLoS One, 2020, 15(12), 1–25,  DOI:10.1371/journal.pone.0244796.
24. R. Kumar Sadhu, et al., Modelling cellular spreading and emergence of motility in the presence of curved membrane proteins and active cytoskeleton forces, Eur. Phys. J. Plus, 2021, 136(5), 1–37.
25. R. Kumar Sadhu, A. Iglič and N. S. Gov., A minimal cell model for lamellipodia-based cellular dynamics and migration, J. Cell Sci., 2023, 136(14), jcs260744.
26. X. Li, et al., A nanostructure platform for live-cell manipulation of membrane curvature, Nat. Protoc., 2019, 14(6), 1772–1802.
27. J. Zimmerberg and M. M. Kozlov, How proteins produce cellular membrane curvature, Nat. Rev. Mol. Cell Biol., 2006, 7(1), 9–19.
28. D. Kabaso, et al., On the role of membrane anisotropy and BAR proteins in the stability of tubular membrane structures, J. Biomech., 2012, 45(2), 231–238.
29. F. Campelo, et al., Modeling membrane shaping by proteins: focus on EHD2 and N-BAR domains, FEBS Lett., 2010, 584(9), 1830–1839.
30. T. Idema and D. J. Kraft, Interactions between model inclusions on closed lipid bilayer membranes, Curr. Opin. Colloid Interface Sci., 2019, 40, 58–69.
31. V. Kralj-Iglič, et al., Free energy of closed membrane with anisotropic inclusions, Eur. Phys. J. B, 1999, 10, 5–8.
32. G. Kumar and A. Srivastava, Membrane remodeling due to a mixture of multiple types of curvature proteins, J. Chem. Theory Comput., 2022, 18(9), 5659–5671.
33. N. Ramakrishnan, P. B. Sunil Kumar and R. Radhakrishnan, Mesoscale computational studies of membrane bilayer remodeling by curvature-inducing proteins, Phys. Rep., 2014, 543(1), 1–60.
34. W. Pezeshkian and J. H. Ipsen, Mesoscale simulation of biomembranes with FreeDTS, Nat. Commun., 2024, 15(1), 548.
35. Y. Ravid et al., Numerical studies of triangulated vesicles with anisotropic membrane inclusions, Advances in Biomembranes and Lipid Self-Assembly. 2024, pp. 21–40.
36. M. Fošnarič, et al., The influence of anisotropic membrane inclusions on curvature elastic properties of lipid membranes, J. Chem. Inf. Model., 2005, 45(6), 1652–1661.
37. A. Iglič, T. Slivnik and V. Kralj-Iglič, Elastic properties of biological membranes influenced by attached proteins, J. Biomech., 2007, 40(11), 2492–2500.
38. A. Iglič, et al., On the role of membrane anisotropy in the beading transition of undulated tubular membrane structures, J. Phys. A: Math. Gen., 2005, 38(40), 8527.
39. V. Kralj-Iglič, et al., Deviatoric elasticity as a possible physical mechanism explaining collapse of inorganic micro and nanotubes, Phys. Lett. A, 2002, 296(2), 151–155.
40. V. Kralj-Iglič, Stability of membranous nanostructures: a possible key mechanism in cancer progression, Int. J. Nanomed., 2012, 3579–3596.
41. V. Kralj-Iglič, S. Svetina and B. Žekš, Shapes of bilayer vesicles with membrane embedded molecules, Eur. Biophys. J., 1996, 24, 311–321.
42. W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Z. Naturforsch., C, 1973, 28(11–12), 693–703.
43. G. Kumar, S. Chaithanya Duggisetty and A. Srivastava, A review of mechanicsbased mesoscopic membrane remodeling methods: capturing both the physics and the chemical diversity, J. Membr. Biol., 2022, 255(6), 757–777.
44. M. Deleu, et al., Interaction of surfactin with membranes: a computational approach, Langmuir, 2003, 19(8), 3377–3385.
45. M. Fošnarič, et al., Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces, Soft Matter, 2019, 15, 5319–5330.
46. F. C. Frank, I. Liquid crystals. On the theory of liquid crystals, Discuss. Faraday Soc., 1958, 25, 19–28.
47. G. Lasher, Monte Carlo results for a discrete-lattice model of nematic ordering, Phys. Rev. A: At., Mol., Opt. Phys., 1972, 5(3), 1350.
48. N. Ramakrishnan, P. B. Sunil Kumar and J. H. Ipsen., Monte Carlo simulations of fluid vesicles with in-plane orientational ordering, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81(4), 041922.
49. D. N. Theodorou and U. W. Suter, Shape of unperturbed linear polymers: polypropylene, Macromolecules, 1985, 18(6), 1206–1214.
50. Github cluster-trisurf. cluster-trisurf. Accessed: 2025-08-04. 2025. url: https://github.com/yoavrv/cluster-trisurf.
51. Y. Ravid, et al., Theoretical model of membrane protrusions driven by curved active proteins, Front. Mol. Biosci., 2023, 10, 1153420.
52. U. Seifert, K. Berndl and R. Lipowsky, Shape transformations of vesicles: Phase diagram for spontaneous-curvature and bilayer-coupling models, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44(2), 1182.
53. V. Kralj-Iglič, S. Svetina and B. Žekš., The existence of non-axisymmetric bilayer vesicle shapes predicted by the bilayer couple model, Eur. Biophys. J., 1993, 22, 97–103.
54. N. David Mermin, The topological theory of defects in ordered media, Rev. Mod. Phys., 1979, 51(3), 591.
55. L. Mesarec, et al., Effective topological charge cancelation mechanism, Sci. Rep., 2016, 6(1), 27117.
56. M. Bowick, D. R. Nelson and A. Travesset, Curvature-induced defect unbinding in toroidal geometries, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69(4), 041102.
57. V. Vitelli and A. M. Turner, Anomalous coupling between topological defects and curvature, Phys. Rev. Lett., 2004, 93(21), 215301.
58. L. Mesarec, et al., Coupling of nematic in-plane orientational ordering and equilibrium shapes of closed flexible nematic shells, Sci. Rep., 2023, 13(1), 10663.
59. O. D. Lavrentovich, Topological defects in dispersed words and worlds around liquid crystals, or liquid crystal drops, Liq. Cryst., 1998, 24(1), 117–126.
60. H. Poincaré, Sur les courbes définies par les équations différentielles, J. Math. Pures Appl., 1885, 1, 167–244.
61. R. D. Kamien, The geometry of soft materials: a primer, Rev. Mod. Phys., 2002, 74(4), 953.
62. G. Napoli and L. Vergori, Extrinsic curvature effects on nematic shells, Phys. Rev. Lett., 2012, 108(20), 207803.
63. R. L. Blumberg Selinger, et al., Monte Carlo studies of the XY model on two-dimensional curved surfaces, J. Phys. Chem. B, 2011, 115(48), 13989–13993.
64. V. Kralj-Iglič, et al., Quadrupolar ordering of phospholipid molecules in narrow necks of phospholipid vesicles, J. Stat. Phys., 2006, 125, 727–752.
65. M. Fošnarič, et al., Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces, Soft Matter, 2019, 15(26), 5319–5330.
66. Ž. Pandur, et al., Surfactin molecules with a cone-like structure promote the formation of membrane domains with negative spontaneous curvature and induce membrane invaginations, J. Colloid Interface Sci., 2023, 650, 1193–1200.
67. M. Denisse Rueda-Contreras, et al., On Gaussian curvature and membrane fission, Sci. Rep., 2021, 11(1), 9562.
68. A.-L. Le Roux, et al., Dynamic mechanochemical feedback between curved membranes and BAR protein self-organization, Nat. Commun., 2021, 12(1), 6550.
69. N. Bobrovska, et al., On the role of anisotropy of membrane components in formation and stabilization of tubular structures in multicomponent membranes, PLoS One, 2013, 8(9), e73941.
70. Z. Jarin, et al., Lipid-Composition-Mediated Forces Can Stabilize Tubular Assemblies of I-BAR Proteins, Biophys. J., 2021, 120, 46–54.
71. H. Noguchi, Curvature-sensing and generation of membrane proteins: a review, Soft Matter, 2025, 21, 3922–3940.
72. B. Ugarte-Uribe, et al., Drp1 polymerization stabilizes curved tubular membranes similar to those of constricted mitochondria, J. Cell Sci., 2019, 132(4), jcs208603.
73. J. Hu, W. A. Prinz and T. A. Rapoport, Weaving the web of ER tubules, Cell, 2011, 147(6), 1226–1231.
74. K. Bohinc, et al., Shape variation of bilayer membrane daughter vesicles induced by anisotropic membrane inclusions, Cell. Mol. Biol. Lett., 2006, 11, 90–101.
75. M. Fošnarič, et al., Stabilization of pores in lipid bilayers by anisotropic inclusions, J. Phys. Chem. B, 2003, 107(45), 12519–12526.
76. A. Iglič and V. Kralj-Iglič, Stabilization of Hydrophilic Pores in Charged Lipid Bilayers by Anisotropic Membrane Inclusions, Advances in planar lipid bilayers and liposomes, 2008, vol. 6, pp. 1–26.
77. A. Iglič, et al., On the role of anisotropy of membrane constituents in formation of a membrane neck during budding of a multicomponent membrane, J. Biomech., 2007, 40(3), 579–585.
78. H. Hägerstrand, et al., Curvature-dependent lateral distribution of raft markers in the human erythrocyte membrane, Mol. Membr. Biol., 2006, 23(3), 277–288.
79. H. Noguchi, Membrane tubule formation by banana-shaped proteins with or without transient network structure, Sci. Rep., 2016, 6, 20935.
80. G. Kumar, N. Ramakrishnan and A. Sain, Tubulation pattern of membrane vesicles coated with biofilaments, Phys. Rev. E, 2019, 99(2), 022414.
81. A. Behera, et al., Deformation of membrane vesicles due to chiral surface proteins, Soft Matter, 2021, 17(34), 7953–7962.
82. F. Bonazzi and T. R. Weikl, Membrane morphologies induced by arc-shaped scaffolds are determined by arc angle and coverage, Biophys. J., 2019, 116(7), 1239–1247.
83. P. Vyas, P. B. Sunil Kumar and S. Lal Das., Sorting of proteins with shape and curvature anisotropy on a lipid bilayer tube, Soft Matter, 2022, 18(8), 1653–1665.
84. N. Ramakrishnan, J. H. Ipsen and P. B. Sunil Kumar, Role of disclinations in determining the morphology of deformable fluid interfaces, Soft Matter, 2012, 8(11), 3058–3061.
85. N. Ramakrishnan, P. B. Sunil Kumar and J. H. Ipsen, Membrane-mediated aggregation of curvature-inducing nematogens and membrane tubulation, Biophys. J., 2013, 104(5), 1018–1028.
86. N. Walani, J. Torres and A. Agrawal, Anisotropic spontaneous curvatures in lipid membranes, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89(6), 062715.
87. K. Xiao, C.-X. Wu and R. Ma, Vesiculation mechanisms mediated by anisotropic proteins, Phys. Rev. Res., 2023, 5(2), 023176.
88. B. Alberts et al., Essential cell biology, Garland Science, 2015.
89. G. Gompper and D. M. Kroll, Random surface discretizations and the renormalization of the bending rigidity, J. Phys. I, 1996, 6(10), 1305–1320.
90. G. Gompper and D. M. Kroll, Triangulated-surface models of fluctuating membranes, Stat. Mech. Membr. Surf., 2004, 359–426.
91. S. Penič, et al., Bending elasticity of vesicle membranes studied by Monte Carlo simulations of vesicle thermal shape fluctuations, Soft Matter, 2015, 11(25), 5004–5009,  10.1039/C5SM00431D.
92. V. Kralj-Iglič, et al., Minimizing isotropic and deviatoric membrane energy–An unifying formation mechanism of different cellular membrane nanovesicle types, PLoS One, 2020, 15(12), e0244796.

---