Protein aggregation and membrane bending govern nanotube morphology

Abstract Membrane nanotubes have been identified as a dynamic way for cells to connect and interact over long distances. These long and thin membranous protrusions have been shown to serve as conduits for the transport of proteins, organelles, viral particles, and mRNA. Nanotubes typically appear as thin and cylindrical tubes, but they can also have a beads-on-a-string architecture. Here, we study the role of membrane mechanics in governing the architecture of these tubes and show that the formation of a bead-like structure along the nanotube can result from a local heterogeneity in the membrane due to protein aggregation or due to membrane composition. We present numerical results that predict how membrane properties, protein density, and tension compete to create a phase space that governs the morphology of a nanotube. We also find that the stability of multiple beads is regulated by an energy barrier that prevents two beads from fusing. These results suggest that the membrane-protein interaction, protein aggregation, and membrane tension closely regulate the tube radius, number of beads, and the bead morphology.


Introduction
Cell-cell communication is fundamental for material and information exchange in multicellular organisms and has been studied for many decades. Various mechanisms for short range communication have been reported including plasmodesmata in plant cells [1,2], gap junctions in animal cells [3], and synapses in nervous systems [4][5][6]. Membrane nanotubes, also referred as tunneling nanotubes have been identified as intercellular structures that connect cells over long distances, i.e., hundreds of micrometers (see Fig.1) [7,8], providing a highway for the transfer of small molecules, organelles, cargoes, mitochondria [9], and mRNA [10]. Membrane nanotubes have been observed in a wide variety of cell types such as epithelial cells, fibroblasts, neurons, and immune cells [8,[11][12][13]. Artificial nanotubes have also been produced from lipid vesicles [8,14]. Although questions remain regarding the mechanisms governing nanotube formation, elongation, and fusion of the plasma membranes between the two connected cells, it is clear that the role played by membrane tension and heterogeneities along the nanotube remains poorly understood.
Membrane nanotubes are typically long and thin cylindrical protrusions with sub-micron diameters and lengths on the order of several hundred microns [7]. In contrast to other types of cellular projections, such as filopodia or endocytic invaginations, which are attached to the substrate and are often function specific [15][16][17], nanotubes are suspended in the medium and participate in a wide variety of cellular functions including embryonic development [18], calcium-mediated signal transduction [19], spreading of pathogens [13], spinal cord injuries [20], and even in • We assume that the total energy of the system includes the membrane bending energy and a contribution from membrane-protein aggregation in a dilute regime. Thus, the membrane is modeled using an augmented version of Helfrich energy for elastic manifolds including the membrane-protein interaction contributions [42,[53][54][55][56]. The aggregation energy of the proteins, which is proportional to the gradient of the protein density, is modeled in our system by a hyperbolic tangent function. A detailed discussion of the full energy can be found in [56]. We do not explicitly model the diffusion of these proteins and the dynamics of the domains in this study; rather we assume that diffusion has run its course and we are left with a protein heterogeneity on the surface. We vary this heterogeneity along the nanotube. A schematic depicting membrane-protein interactions that could lead to the formation of beads along a nanotube. Proteins (shown in red) can aggregate along the membrane to induce local curvature and create a tension heterogeneity, which results in the formation of beads of various sizes and shapes. We assume the dispersed proteins to be cone-shaped such that their meridian makes an angle ' with the normal vector to the surface. (E) The coordinate system to define a surface by tangent basis a 1 , a 2 and normal vector n.
• We do not consider the role of forces applied by actin-mediated nanotube formation, so that we can focus only on membrane nanotube deformation due to membrane-protein interactions [57,58].
• For simplicity in the numerical simulations, we assume that the membrane in the region of interest is rotationally symmetric (Fig. 1C).

Membrane energy and equations of motion
We model the membrane with two contributions to the strain energy -one from membrane-protein interactions and the other from membrane bending. The entropic interaction of transmembrane proteins with the lipid bilayer is written as a constitutive function of the protein density (number per unit area). While the exact form of this energy is yet to be experimentally verified, based on thermodynamic arguments, a quadratic dependence of the energy on the local protein density has been proposed as [42,53,56,59,60], where W is the energy per unit area, and ↵ 2 is an indicator of the membrane-protein interaction strength [56]. can depend explicitly on the surface coordinates ✓ ⇠ ( ⇠ 2 {1,2}) to allow for local heterogeneity. Also, these proteins are assumed to be transmembrane, conical insertions such that the meridian of each protein is at an angle ' with the normal vector to the membrane surface (n) (See Fig. 1D). In the dilute regime, the local induced-spontaneous curvature due to membrane-protein interaction can be modeled as a linear function of the surface protein density as [43,56]  Bending modulus of rigid domain [66][67][68] 320 -9600 pN · nm Protein density [69] 0 -2.5⇥10 4 nm 2 ↵ Protein aggregation strength [43] 3.57 ⇥ 10 3 nm 1.5 · pN 0.5 µ' Constant length scale [43] 200 nm A membrane Total area of the membrane 8⇡ µm 2 p Transmembrane pressure 0 pN/nm 2 where µ is a length scale that represents the lipid-protein specific entropic interactions. Other forms and relationships between spontaneous curvature and protein density are explored in [42]. The energy associated with membrane bending due to the spontaneous curvature is given by the Helfrich Hamiltonian, modified to include heterogeneous membrane properties as where H is the local mean curvature, and K is the local Gaussian curvature.  and  G are bending and Gaussian moduli respectively, which in the general case can vary along the surface coordinate ✓ ⇠ [39,61,62]. Hence, the total energy of the membrane including both bending and protein contributions is given by Protein aggregation (4) Assuming that the membrane is at mechanical equilibrium at all times, we minimized the total energy using a stress balance approach [48], and solved the governing equations numerically. A complete derivation of all equations is provided in the SOM and the notation and values of parameters used in the model are summarized in Tables S1, S2, and S3.

Numerical implementation
In the axisymmetric coordinates, the equations of motion (Eq. S14 and Eq. S15) simplify to a system of first-order differential equations (Eq. (S26)) with six prescribed boundary conditions (Eq. (S27)). In order to solve these coupled equations, we used the commercially available finite element solver COMSOL MULTIPHYSICS R 5.3. In this work, we assume that the total area of the membrane is conserved and to focus on the net effect of membrane tension, we set the transmembrane pressure at zero (p = 0).

Results
Heterogeneity in membrane properties can lead to formation of bead-like structures along a nanotube For membrane nanotubes, it has been shown that the composition of the lipid bilayer is a critical factor in determining the shape and radius of nanotubes [70][71][72]. A change in lipid composition also corresponds to a change in membrane stiffness [70]. To understand how a heterogeneity in membrane properties due to lipid composition affects the shape of a nanotube, we conducted simulations on cylindrical nanotubes with a significant aspect ratio of radius 35 nm and length 100 µm, and set the boundary tension to 0 = 0.064 pN/nm. The effect of boundary tension on the initial nanotube radius and length is shown in Figs. S1 and S2. On this nanotube, we prescribed a region of higher bending rigidity representing a change in the composition but not a change in the asymmetry between the leaflets. This defines a domain at the center of the nanotube that has a higher bending moduli with respect to the rest of the surface, with a rigidity ratio  ratio =  rigid / lipid [67,72] (Fig. 2A). To have a smooth continuous transition, we implemented the difference in bending moduli using a hyperbolic tangent function Eq. (S28), such that the area of the rigid domain is fixed (A rigid = 1.6⇡ µm 2 ), and the bending rigidity ratio ( ratio ) is varied between 1-30 [67,72,73]. With increasing  ratio , the membrane bends significantly in the area of the rigid segment and an ellipsoidal bead-shaped structure forms along the nanotube (Fig. 2B). For  ratio < 10, the radius of the bead increases almost linearly. However, for  ratio > 10 the size of the bead remains almost constant (Fig. 2B, inset). The local variation of membrane tension corresponding to the membrane bending in the beaded domain and the applied incompressibility constraint is shown in Fig. S3. For completeness, we also used a hyperbolic tangent function to parameterize a variation in the Gaussian modulus (  G ) between -20 to 20. We found that the effect of varying Gaussian modulus alone leads to a very insignificant membrane deformation (see Fig S4).  Beads formation along a nanotube due to membrane-protein interactions Motivated by our numerical observation that a heterogeneity along the membrane can lead to bead-like structures, we next explored if localized surface protein aggregation could induce a similar deformation along the nanotube. To explore how protein domains can affect the tubular shape of a nanotube, we applied a distribution of the protein density -prescribed as a hyperbolic tangent function Eq. (S28) -along the surface of the elongated nanotube in 5 Fig. 1C with 0 = 0.064 pN/nm and  ratio = 1. Here, we assumed that the number of proteins per area increases ( 0 = 0 to 0 = 1.25 ⇥ 10 4 nm 2 ), while the area covered by the proteins at the center of the nanotube remains constant (A protein = 1.6⇡ µm 2 ) (Fig. 3A). As shown in Fig. 3B, the membrane bends outward in the area of the protein aggregation in response to the protein-induced spontaneous curvature. As the number of proteins in the defined area is increased, the membrane deformation also increases such that it resembles a bead-shaped structure forming along the nanotube (Fig. 3B, inset). The local membrane tension distribution along the nanotubes and the bead-shaped structures is shown in Fig. S5.

Morphological landscape of bead-shaped structures along a nanotube
Thus far, simulation results demonstrate that two unrelated mechanisms, heterogeneity in membrane rigidity and protein-induced curvature, independently lead to formation of bead-like structures along the membrane nanotube. In order to explore how these two mechanisms might interact and modulate the shape of a nanotube, we conducted simulations where the heterogeneous domain has effects both from protein-induced spontaneous curvature and from increased bending rigidity. Initially, we repeated the simulations in Fig. 3 but this time assuming that  ratio = 11. Interestingly, we found that the interaction of these two mechanisms leads to the formation of beads with different shapes. Based on the magnitude of the protein density, three different oblate spheroid shapes were obtained; (i) an ellipsoidal bead at 0 = 0, (ii) a flat cylindrical bead at 0 = 4 ⇥ 10 5 nm 2 , and (iii) a large unduloid-shaped bead at 0 = 5.5 ⇥ 10 5 nm 2 (Fig. 4A). Thus, coupling the two modes of spatial heterogeneity along a membrane nanotube not only increases the radius of the bead (see Fig. S6) but also changes the shape of the bead. This shape transition suggests that the energy landscape of the membrane, which is now modulated by a heterogeneity in  and a heterogeneity in , can lead to a range of shapes. To further identify over what range of protein density and  ratio these shape transitions occur, we performed the same simulation over a range of protein densities, ( 0 = 0 2.5 ⇥ 10 4 nm 2 ), as well as over a range of  ratio = 1 11, encompassing soft protein domains to very stiff clusters. This variation allowed us to construct a phase diagram to identify the regions of different bead morphologies (Fig. 4B). The red region represents the formation of ellipsoidal bead-shaped structures, the blue region denotes the cylindrical beads, and the green region indicates the unduloid-shaped beads configuration.
These changes of the bead morphology can be understood better by comparing the two "induced" length scales by the rigid domain (L k = 1/2 p  rigid / ) [74] and the protein aggregation (L = 1/µ' ) [37,46]. The "induced" length scale by the membrane-protein interaction 1/µ' is the preferred radius of curvature of the protein-enriched domain and increases with decreasing (with = 0, both spontaneous curvature and the effect of protein aggregation are eliminated). The phase diagram in Fig. 4B shows that these two length scales act in tandem to regulate the bead size and shape. In order to compare the length scales, we consider a case with a fixed value of L k (e.g.  ratio = 2), and vary L by increasing the protein density (Fig. 4B). First, we see that the ellipsoidal beads form for low protein density, i.e. when L L k (Fig. 4B). Second, the cylindrical beads are energetically favorable for average values of protein density when the two length scales become comparable (L ⇠ L k ). Finally, at very large values of protein density, when the "induced" length scale by the rigid domain becomes dominant ( L k L ), the unduloid-shaped bead forms along the nanotube (Fig. 4B).
Regardless of which mechanism dominates, the edge tension 0 implicitly governs the length scale of the membrane; therefore the natural question that arises is how does this tension govern the shape and length scale transitions of the beads? We explore these questions by conducting two sets of simulations. First, for a homogeneous membrane ( ratio = 1) we varied the edge membrane tension ( 0 ) and protein density ( 0 ) in a range ( 0 = 0.004 0.064 pN/nm and 0 = 0 2.5 ⇥ 10 4 nm 2 ) (Fig. 4C). We observed that high edge tension shifted the transition of ellipsoidal to cylindrical and unduloid-shaped beads to large protein densities. Second, we fixed the protein density ( 0 = 6 ⇥ 10 5 nm 2 ) and varied the edge tension ( 0 ) and the rigidity ratio ( ratio ) between 0.004-0.064 pN/nm and 1-11 respectively (Fig. 4D). In this case, all three possible shapes of beads are accessible with varying the bending rigidity ratio. In fact, by increasing the edge membrane tension either at a constant protein density or a fixed rigidity ratio, we decreased the value of the the induced length scale by the rigid domain ( p  rigid / ) and moved from the cylindrical and unduloid-shaped beads to the ellipsoidal bead-shaped region of the phase space.

Interaction between multiple beads along a nanotube
Often, multiple stable beads are observed along a membrane nanotube, suggesting that multiple domains of heterogeneity exist along the nanotube [24,27]. Using bioorthogonal click chemistry to label sialic acids on surface cell glycans in COS7 cells [27], a series of well defined and temporally stable bead structures along membrane nanotubes has been previously observed ( Fig. 1A and B). These observations lead to the following question: how does the stability profile of these beaded strings depend on the different length scales associated with beaded nanotubes?
In order to answer this question, we conducted simulations for the formation of two beads along the membrane by prescribing two areas of heterogeneity for three cases -(i) varying membrane rigidity alone in each domain in the absence of protein aggregation, (ii) varying protein density for uniform rigidity, and (iii) varying protein density for domains with different rigidity.
We first found when there are two regions of proteins far from each other with an end-to-end area given by A separation = 1.6⇡ µm 2 (Fig. 5A, top), two independent beads form ( Fig. 5A bottom). The size and the shape of the beads are independent of the number of domains as long as the domains are far away from each other ( Fig. S7 and Fig. S8). Having established that the regions of heterogeneity are independent when they are far from each other, we next asked under what conditions might these beads interact with one another? In other words, what length scales govern the stability features of multiple beads knowing there is a certain relaxation length between the beads and cylinder? To answer this question, we repeated the simulation of two beads by varying the rigidity ratio ( ratio ) and end-to-end area (A separation ) between 1-11 and 0 1.6⇡ µm 2 respectively. Based on the results, we constructed a phase diagram separating the two possible morphologies; (i) two distinct beads represented by color blue, and (ii) one single bead denoted by color red (Fig. 5B). We found that when the two domains are bought together by decreasing the area between the heterogeneity, there is a smooth transition from two beads connected by a string to a single bead (Fig. 5B, C). This suggests that at close distances, the energy minimum of the nanotube is attained for a single large bead rather than for two beads connected by a string. We also observe the smooth transition in number of beads is accompanied by the shape transition from an unduloid-shaped bead to a large ellipsoidal bead (Fig. 5C).
However, when we varied the area separation between the beads for a range of protein densities, we found a sharp transition from two beads to one bead accompanied by a snap-through instability (Fig. 5C,D,E). For small protein densities, when the area between the two beads is decreased, the nanotube appears to have one large bead, while the transition from two beads to one bead is continuous indicating that there is no energy barrier to go from one state to another (Fig. 5D, E, dashed, black line). However, as the protein density increases, we found the emergence of a snap-through instability for small separation areas (Fig. 5E) corresponding to bead shapes that show a distinct transition from two ellipsoids to a flower-shaped bead to a large ellipsoidal bead (Fig. 5E, green dashed line, F). This means that the landscape between two beads and one bead at high protein densities is governed by a large energy barrier. And finally, when we repeated these calculations for a rigid protein domain such as  ratio = 4 ( Fig. 5 G,H,I), we found that the area and protein density still govern the energy and stability landscape but the transition point -where the snap-through instability occurs -is shifted towards lower protein densities, with no change in bead shapes (Fig. 5G, H, I). Thus, protein-mediated beads along nanotubes potentially tend to be stable because any fusion event of two beads may often be accompanied by energy barriers and nontrivial shape transitions. This energy barrier corresponds to the length scale of domain separation being one order of magnitude larger than the length scale of the heterogeneity (see Fig. S10).

Discussion
Tunneling nanotubes, membranous projections between cells, have been identified as cellular transport highways [7][8][9]. Much of the biophysics associated with these dynamic structures are only now beginning to be explored but it is becoming increasingly clear that the cellular membrane and membrane-protein interactions play a critical role in maintaining these cellular highways [24,27,70,75,76]. In this study, we used systems biophysics to explore the energy landscape of nanotube formation and how heterogeneity in the membrane either due to material properties or protein aggregation affects these structures. Our results can be summarized as -membrane heterogeneity due to either membrane rigidity (Fig. 2) or protein-induced spontaneous curvature can result in the formation of bead-like structures along a nanotube (Fig. 3). Additionally, the interaction between these two modes of heterogeneity can lead to the formation of beads with distinct shapes while the transitions between these shapes from ellipsoidal to cylindrical to unduloid-shaped beads are consequences of competing length scales in the system (Fig. 4). Finally, the length scale competition induced by these effects coupled with the length scale between two regions of heterogeneity also governs the stability landscape of multiple beads (Fig .5). Although a similar interplay between adhesion and bending energy has been explored in shaping the plasma membrane and the shape of integrin adhesions, [77,78], the shaping of membrane nanotubes resulting from the interplay between membrane heterogeneity due to a spatially varying membrane rigidity and protein-induced spontaneous curvature had not been previously been explored.
Our simulations allow us to make the following experimentally testable predictions. First, tension at the edge of the nanotube not only governs the nanotube radius but also its response to heterogeneity. Therefore, manipulating cell tension and evaluating how it affects the morphology of the nanotube will provide information on how tensional homeostasis of cells affects the membrane nanotubes. Additionally, we found that there is an energy barrier that governs the stability landscape of the transition from two beads to one bead. This energy barrier, governed by a snap-through instability suggests that the fusion of two beads depends on the membrane composition and its material properties such that under high protein density or high rigidity conditions, there is a large crossover energy for fusion. These predictions can be extend to multiple beads as well.
Although our model has provided some fundamental insights into the role of membrane composition in nanotube morphology, we have made simplifying assumptions that may need to be revisited as more experimental evidence is gathered regarding these structures. Additionally, the dynamics of membrane-protein diffusion and what could be the underlying mechanisms that govern protein aggregation along the nanotube as suggeste addressed by our model and are not yet fully understood. While there is evidence for strong curvature-mediated feedback for protein aggregation [56], it is possible that proteins in the lumen of the nanotube coupled with biochemical signaling can lead to the formation of protein microdomains. For exmaple, it is known that phase separation between two main components of the membrane -clusters of sphingolipids and cholesterol molecules -can result in formation of lipid rafts [67,79]. Second, the role of the cytoskeleton (actin and microtubule filaments) and motor protein transport along the nanotube is known to be important [7,8,80] but a correlation with the beading morphology is yet to be established. And finally, the role of active transport versus flow of cytosolic components in governing the stability of these nanotubes is not captured by our model and remains an active focus of our research and modeling efforts. Despite our simplifying assumptions, our model provides insight into the morphological complexity of membrane nanotubes. The area separation between two beads versus the bending rigidity ratio phase diagram, 0 = 0. There are two possible shapes, (i) two separated beads denoted with color blue, and (ii) one single bead marked by color red. The transition from two to one bead is smooth everywhere in this parameter space. (C) The smooth membrane shape transition from two beads to one bead with decreasing area separation for  ratio = 11 and membrane profiles for points marked in panel (B). (D) Membrane profiles for the black dashed lines in panel (E) show the smooth evolution of membrane shape from two beads to one bead at 0 = 0.7 ⇥ 10 4 nm 2 for  ratio = 1 and increasing A separation . (E) Phase diagram for the area separation between two beads versus the protein density at  ratio = 1. The colors represent same morphologies as panel (B). Smooth transition from two beads to one bead at low protein density and a snap-through instability in the transition from two beads to one bead for large protein densities are observed.

Equations of motion
In this section, we present a brief derivation of the generalized equations of motion for biological membranes under the assumptions stated in the main text. We then restrict the equations to axisymmetric coordinates and solve them in a prescribed membrane area domain. Details of all derivations can be found in [1][2][3].
The local force balance in the absence of inertia for any material point on the membrane is where ( ;⇠ ) is the surface divergence, ⇠ is the stress vector, p is the pressure difference between the inside and outside of the volume bounded by the membrane, and n is the unit vector normal to the membrane. The surface divergence in Eq. (S1) can be rewritten as [1] where a is the determinant of the first fundamental form metric a ⇠⌘ . The stress vector in Eq. (S1) can be decomposed into normal and tangential components as For elastic surfaces for which the energy density per unit mass F (a ⇠⌘ , b ⇠⌘ ) responds to the first and second fundamental forms, the normal and tangential components of the stress vector in Eq. (S3) can be expressed as [1,4] Here, ⇢ is the surface mass density, and H and K are the mean and Gaussian curvatures given by Here (a ⇠⌘ ) = (a ⇠⌘ ) is the dual metric and " ⇠⌘ is the permutation tensor defined by " 12 = " 21 = 1/ p a, " 11 = " 22 = 0.
In the case of area incompressibility (J = 1), the general form of free energy density per unit mass can be rewritten as where (x ⇠ , t) is a Lagrange multiplier field required to impose invariance of ⇢ on the whole of the surface (see [1] for full derivation). Substituting W = ⇢F into Eq. (S7) we get whereb ⇠⌘ is the co-factor of the curvature tensor, and = ( + W ) can be interpreted as the membrane tension [5]. Combining Eqs. (S9), (S4), and (S3) with Eq. (S1) gives the equation of motion in normal and tangential directions as and Tangential: N ⌘⇠ Here (·) is the surface Laplacian, and () |exp denotes the explicit derivative respect to coordinate ✓ ⇠ .

Helfrich energy and mechanical equilibrium
For the local energy density of a lipid bilayer membrane, we use an augmented version of the Helfrich energy accounting for an energy associated with the local protein density ( (✓ ⇠ )) [3, 6-8] where W is the local energy density, C is the induced spontaneous curvature due to lipid-protein interactions,  is the bending modulus, and  G is the Gaussian modulus. In Eq. (S12), ↵ is a positive constant representing derivative of the chemical potential with respect to the protein density at constant pressure and temperature [8]. It should be noted that Eq. (S12) is different from the standard Helfrich energy by a factor of 2. We take this net effect into consideration by choosing the value of the bending modulus to be twice that of the standard value of bending modulus typically used for lipid bilayers [9]. For low protein densities (dilute regime), the induced-spontaneous curvature C( ) in Eq. (S12) can be expressed in term of protein density as [6] C( ) = (µ') , where ' is the angle between cone-shaped proteins meridian and the normal vector to the surface (n), and µ is a length-scale correction representing the lipid-protein specific moietic interaction (see Fig. 1) [10]. Substituting the modified version of Helfrich energy function (Eq. (S12)) and Eq. (S13) into the first functional normal variation of total energy (Eq. (S10)) gives the so-called "shape equation," [11] [ (H (µ')  (S14) 3 A consequence of spatial variation of membrane properties and protein density is that is not homogeneous along the membrane [5,12,13]. Substituting Eq. (S12) into the balance of forces tangential to the membrane Eq. (S11) gives the spatial variation of membrane tension, where r is the partial derivative with respect to the coordinate (✓ ⇠ ).

Equations of motion in axisymmetric coordinates
We define a surface of revolution (Fig. 1C) by where s is the arc length along the curve, R(s) is the radius from the axis of rotation, Z(s) is the height from the base plane and (e r , e ✓ , k) form the basis coordinate. Defining as the angle made by the tangent with respect to the vertical gives satisfying the identity is the partial derivative with respect to the arc length. Using this, we can define the normal and tangent vectors to the surface as n = sin e r (✓) + cos k, a s = cos e r (✓) + sin k. (S18) Using this parameterization, we can now write the tangential ( ⌫ ) and transverse ( ⌧ ) curvatures as The mean curvature (H) and Gaussian curvature (K) are obtained by summation and multiplication of the tangential and transverse curvatures Finally, we define L = 1 2 R(W H ) 0 such that we obtain a system of first-order differential equations with six unknowns R, Z, , H, L, and [11,13], In order to solve this system of equations (Eq. (S21)), we need to provide six boundary conditions. We consider an axisymmetric cylindrical membrane with (1) fixed radius at both boundaries, (2) = ⇡/2 at both ends to ensure the edges remain vertical, (3) fixed height at one end and prescribed tension ( 0 ) at the other end (see Fig. 1C). These boundary conditions can be summarized as follow, (S22) One advantage of an asymmetric coordinate system is that the manifold area (A) can be expressed in term of arc length [13], Thus, Eq. (S23) allows us to convert the system of equations in Eq. (S21) that are in term of arc length to a system in terms of area. So, we can prescribe the total area of the membrane instead of arc length where (˙) denotes the derivative with respect to the area. The boundary conditions Eq. (S22) remain unchanged, and the upper limit of the domain changes to A max instead of s max . To non-dimensionalize the system of equations (Eq. (S24)), we use two parameters, the radius of the initial cylindrical membrane (R 0 ), and lipid bilayer bending modulus ( lipid ). Using these constants, we can define Rewriting Eq. (S24) in terms of the dimensionless variables in Eq. (S25), we get [13] xẋ = cos , xẏ = sin , x 2˙ = 2xh sin , and the boundary conditions in Eq. (S22) simplify to where ⇣ max = A/2⇡R 2 0 is the total dimensionless membrane area .

Numerical implementation
We solved the system of first-order differential equations (Eq. (S26)) with boundary conditions Eq. (S27) using the finite element software COMSOL MULTIPHYSICS R 5.3, using the "General Form PDE" module. Here, we summarize the steps and assumptions that we used for each simulation.
• All the simulations were performed for a fixed area of the membrane (A membrane = 8⇡ µm 2 ). • The membrane patch was initialized to be a perfect cylinder ( = ⇡/2) with radius R 0 = 200 nm. • The membrane domain (⇣) was discretized equally with mesh size = 0.001.
• To have a sharp but smooth transition in heterogeneous properties, we used a hyperbolic tangent function given by where denoted the membrane property such as bending modulus (), Gaussian modulus ( G ), or protein density ( ), g is a constant and a 0 represents the area of the heterogeneity. • The applied tension at the boundary ( 0 ) in Fig. S1, the bending modulus ratio ( ratio ) in Fig. 2, and the number of proteins per unit area ( 0 ) in Fig. 3 were progressively increased such that each solution was used as an initial guess for the next step. 6 2 Tables   2.1 Table of Notation Figure S3: Membrane tension distribution along nanotubes corresponding to local bending rigidity variation shown in Fig. 2. The region of higher rigidity and bead-shaped deformations correspond to a lower membrane tension. The inset shows the reduction of the local membrane tension versus the bead radius.  [13,15]. There is a negligible membrane deformation with increasing the Gaussian moduli difference from  G = 20 to  G = 20 such that the bead radius is almost equal to the nanotube radius (r n = 35.334 nm) at  G = 0.  Figure S5: Membrane tension distribution along nanotubes corresponding to local protein aggregation shown in Fig. 3. The region of protein aggregation and bead-shaped deformations corresponds to a lower membrane tension. The inset shows the reduction of the local membrane tension versus the bead radius. Here the negative membrane tension can be interpreted as a surface pressure [5].

Cylindrical bead
Unduloid-shaped bead Figure S6: Bead radius as a function of the protein density for a rigid protein-enriched domain ( ratio = 2). The three different observed bead morphologies in Fig. 4 are separated with same colors. The slope of bead radius versus the protein density in cylindrical and unduloid-shaped beads is smaller compared to the ellipsoidal bead.  Figure S7: Formation of two beads far away from each other with no interaction. (A) Three different shapes of two beads in the presence of two rigid domains of proteins with  ratio = 11, and 0 = 0.064 pN/nm. Here, we set the protein density distribution the same as Fig. 5A. Similar to a single bead (Fig. 4A), as protein density increases, each bead shape transforms independently from an ellipsoidal bead to a cylindrical one and finally a unduloid-shaped bead. (B) The bending rigidity ratio versus the number of the proteins per unit area phase diagram for 0 = 0.064 pN/nm. As expected, the phase diagram is exactly same as Fig. 4B because the beads are far away from each other and therefore their shapes evolve completely independently of one another.   Figure S9: Transition between two beads to one bead with decreasing area separation from A separation = 4.8 µm 2 to A separation = 0. (A) The membrane radius at the middle of nanotube -where A membrane = 4⇡ µm 2 -increases monotonically with decreasing area separation between the beads at 0 = 0.7 ⇥ 10 4 nm 2 . (B) With decreasing area separation between two beads, initially the membrane radius at the middle increases and then suddenly jumps to almost 1200 nm with formation of one bead at 0 = 2.5 ⇥ 10 4 nm 2 . The radius of the middle point decreases after snap-through transition.  Figure S10: Length scale competition determines the minimal distance between beads before they overlap. In Fig.  5C we varied the the area separation between the beads as an independent parameter. By using Eg. S23, we can covert this area separation to arc length (L d ). (A) Contour plot of the Log of the L d /L k ratio for A separation varying between 0 and 1.6⇡ µm 2 as a function of the bending rigidity ratio  ratio . We find that the transition from two beads to one bead with decreasing the distance between domains occurs when the arc length distance L d is one order of magnitude smaller than the induced length scale L k . (B) Contour plot of the Log of the L d /L k ratio as a function of protein density and area separation between beads for a homogeneous membrane  ratio = 1. As expected, the green and blue regions, where Log(L d /L k )  1, overlap with the one bead region in Fig. 5D. However, the yellow region with Log(L d /L k ) > 1 corresponds to the two beads configuration in Fig. 5D. (C) Contour plot of the Log of the L d /L k ratio as a function of protein density and area separation between beads assuming that the bending rigidity along the protein domain is 4 times larger than the bare lipid membrane ( ratio = 4).