Model for ring closure in ER tubular network dynamics

Tubular networks of endoplasmic reticulum (ER) are dynamic structures whose steady-state conformations are maintained by a dynamic balance between the persistent generation and vanishing of the network elements. While factors producing the ER tubules and inter-tubular junctions have been investigated, the mechanisms behind their elimination remained unknown. Here we addressed the ER ring closure, the process resulting in the tubule and junction removal through constriction of the network unit-cells into junctional knots followed by the knot remodeling into regular junctions. We considered the ring closure to be driven by the tension existing in ER membranes. We modeled, computationally, the structures of the junctional knots containing internal nanopores and analyzed their tension dependence. We predicted an effective interaction between the nanopores facilitating the knot tightening and collapse of additional network unit cells. We analyzed the process of the pore sealing through membrane fission resulting in formation of regular junctions. Considering the hemi-fission as the rate-limiting stage of the fission reaction, we evaluated the membrane tensions guarantying the spontaneous character of the pore sealing. We concluded that feasible membrane tensions explain all stages of the ER ring closure.


27
Endoplasmic reticulum (ER) is one of the most vital membrane bound organelles of eukaryotic 28 cells serving as a platform for protein synthesis, folding, modification and transport to intra-29 cellular destination (Alberts et al., 2017, chap. 12). The ER functions determine its structure 30 consisting of the double-membrane nuclear envelope, and the peripheral ER, which spans most 31 of the cytoplasm of mammalian cells (Westrate et al., 2015). The peripheral ER includes 30-32 50nm thick flat membrane sheets and a highly intricate network of membrane tubules with 25-90 33 nm diameters interconnected by three-way junctions (Schroeder et al., 2018;Terasaki, 2018;Voeltz 34 et al., 2002;Westrate et al., 2015). In spite of their complex architecture and topological features, 35 the membranes of all ER elements are interconnected into a continuous system that encloses a 36 common luminal space.

37
Several decade-long studies by the methods of diffraction-limited optical(Lee and Chen, 1988;38 Terasaki et al., 1986), electron (Palade, 1956), and, more recently, super-resolution(Nixon-Abell  (Westrate et al., 2015). The ER nanoscopic structures described so far include ER Movie). As follows from these and other simulations, polygonal unit-cells transform into 115 equilateral triangular unit cells, which then undergo collapse into junctional knots ( Fig.1 A-D). 116 Here we focus on the second stage of the ring closure ( Fig.1 E-H  (1)

124
The bending energy, , is associated with the shape of the membrane midplane that is characterized at 125 each point by the mean, , and Gaussian, , curvatures (Spivak, 1999). We use Helfrich model of 126 membrane bending elasticity (Helfrich, 1973) according to which the membrane material 127 parameters setting the value of the bending energy are the spontaneous curvature, , the bending 128 modulus, , and the modulus of Gaussian curvature, ̅ (Helfrich, 1973). The spontaneous 129 curvature, , which accounts for the intrinsic tendency of the membrane to acquire a bent shape 130 and determines the intra-membrane torque, is generated for ER tubular membranes by the 131 curvature generating proteins of the reticulon/REEP family and has values in the 0.01 -0.1nm -1 132 range, as inferred from the measurements of the ER tubule diameters (Schroeder et al., 2018;133 Terasaki, 2018;Voeltz et al., 2002). The bending modulus, , which sets the scale of the 134 bending energy, while varying in a certain range (Dimova, 2014), has a typical value of, ≈ 10 −19 Joule (Dimova, 2014). The modulus of Gaussian curvature, ̅ , whose direct measurements 136 have been, largely, prohibited by geometrical restrictions, is, typically, negative and has been 137 estimated to have values in the range between 0 and − [see e.g. (Templer et al., 1998;Terzi et 138 al., 2019)].

139
The tension energy, , is associated with the membrane area, , and the system parameter 140 setting the scale of this energy is the tension value, .

141
We define the system energy as the thermodynamic work of creating the junctional knot (or the 142 regular junction) out of an effective membrane reservoir represented by the surrounding tubular 143 ER. Based on the ER structure, we consider the reservoir to consist of cylindrical membrane 144 tubules characterized by a spontaneous curvature, , and subject to a tension, . The radius, , 145 and the corresponding mean curvature, , of the reservoir tubules are given by 147 as set by the mechanical equilibrium between the tension and the intra-membrane torque (see SI 148 A)). We consider the junctional knot to be connected to the reservoir by three tubular arms, 149 which merge with the reservoir along circular border-lines of radius, (Eq. 2), and denote the 150 distance from the border-lines to the unit-cell center by (Fig. 1I). The boundary conditions 151 imposed on membrane shape at the border-lines require a smooth transition between the 152 membranes of the tubular arms and those of the reservoir. Taking the distance, , to be 153 substantially larger than the characteristic length-scale of the system, ≫ √ (see SI A)), we guarantee a negligibly small influence of the boundary conditions on the junctional knot 155 configurations.

156
The bending energy of shaping a unit area of the reservoir membrane into a unit area element of 157 the junctional knot (or the regular junction) is given by The total bending energy, , is obtained by integration of over the system area , which 160 includes the area of the junctional knot (or the regular junction) and that of the connecting arms, The tension energy is given by, 163 = . (5)

164
It is convenient to relate the energy of the junctional knot (or the regular junction), , to that of 165 a hypothetical reference state, 0 , which we define as consisting of three homogeneous 166 cylindrical membrane tubules of radius (Eq. 2) pulled out of the reservoir through the 167 boundary lines up to the center of the unit cell, i.e., each having the length (Fig. 1J). The

182
The configurations and energies of the junctional knots (and the regular junctions) were obtained 183 by numerical minimization of the relative energy, , using the Brakke's Surface Evolver 184 program (Brakke, 1992). The examples of the resulting configurations computed for a specific dimensions and shapes of their external contours ( Fig. 2A,B), the former differs from the latter 188 by an hour-glass pore in the middle ( Fig. 2A). This pore is a remnant of the micron large space 189 confined within the initial network unit-cell before its contraction. The equilibrium size of the constriction and, hence, resists this process. By contrast, the tension favors the pore tightening.
Therefore, the radius of the pore waist decreases from 8nm to 4nm with increasing membrane 194 tension from 0.1 mN/m to 0.5 mN/m (Fig.2C). For even larger tensions the pore radius tends to 195 decrease beneath 4 nm which is the lower boundary of our calculation validity set by the 196 membrane thickness.

197
We further analyzed the structures of junctional knots resulting from a sequential contraction of Here we analyze sealing of the intra-knot pores, which converts the junctional knots into regular 216 junctions. This is a final step of the ring closure event, which reduces the overall number of 217 junctions in the ER network. For a junctional knot containing one pore, the sealing completes the 218 removal of two junctions from the network ( Fig. 1E-H).

219
We consider the membrane fission mediating the sealing (Kozlovsky and Kozlov, 2003) to occur 220 either within the hourglass-like membrane neck forming the rim of the pore, or within one of the 221 short neck-like membrane elements forming the sides of the junctional knot and surrounding the 222 pore. We propose that the factor driving the membrane fission is the membrane stress that 223 accumulates within the junctional knot membrane under the tension-driven constriction of the 224 pore and relaxes as a result of the pore sealing.

225
Favorability of the pore sealing. We start with analyzing the conditions required for the pore 226 sealing to be energetically favorable, i.e., leading to decrease of the system energy. To this end 227 we compare the energy, 0 , of a three-way tubular junction with no pores ( Fig. 2B) with the 228 energy, 1 , of a junctional knot with one pore ( Fig. 2A).

229
The computational results for 0 , 1 and for the energy released by the pore sealing, ∆ = 230 1 − 0 , are presented in (Fig. 3A,B) in dependence on the membrane tension, . According  Based on the previous work (Kozlovsky and Kozlov, 2003), we consider the fission reaction to knot with a pore, 1 , such that ∆ = − 1 .

262
The computation of the configuration and the energy of the hemi-fission stalk is based on the 263 previous works (Hamm and Kozlov, 2000) and, as already mentioned, uses the elastic model of and -1 (Templer et al., 1998). Another parameter is the characteristic decay length of the tilt deformation, , which is determined by the ratio between the monolayer bending and tilt moduli 278 by, = √ , and was taken to change between 1 and 2nm (Golani et al., 2021).

279
The results for the energy of the hemi-fission intermediate, , are shown in (Fig. 3A), and for 280 the activation energy, ∆ , are presented in (Fig. 3B, 4B-C), where the energy values are scaled 281 by the bilayer bending modulus ≈ 10 −19 Joule (Dimova, 2014) , the reference system is the 282 same as used for computations of the junctional knot energy (Fig. 1J), and the tension varies 283 within a biologically feasible range.

284
The dependence of Δ on the tilt decay length, , is shown in (Fig. 4B) for = −0.5. The 285 activation energy is decreasing with increasing, , which means that the smaller the tilt modulus,

286
, the smaller the activation energy of fission.

287
The dependence of Δ on is presented in (Fig. 4C)  vanishes meaning that the fission process happens spontaneously. According to this diagram, for 294 fission to be spontaneous, the tension has to be large, the monolayer modulus of Gaussian 295 curvature has to be as negative as possible, and the tilt modulus has to be small.

296
The computed membrane configurations formed along the full pathway of the junctional knot 297 remodeling into a regular junction are presented in (Fig. 5A-C).

300
Here we considered the critical final stages of the ER ring closure, the process of constriction the ER tubules. While these ER areas must be exposed to tension, the other ER regions, which 323 are outside of the zones of the pulling force action, are expected to remain tension free.

324
According to our results, the spontaneous pore sealing corresponding to a vanishing activation

331
Such tensions could be generated by more than 10 molecular motors acting on each membrane 332 tubule, which must be, biologically, a seldom situation. This implies that special fission proteins 333 might be required to facilitate the junctional knot conversion into the regular junctions.

334
Another prediction of our results is that in ER network regions exposed to ultralow tensions, the

A. Relationship between the radius and tension of a reservoir tubule
To derive the relationship between the radius, " , the tension, , and the spontaneous curvature, % , of a tubular membrane given by (Eq. 2) of the main text, we use the equation of mechanical equilibrium of a curved membrane (see "Some aspect of membrane elasticity" (Poon and Andelman, 2006, pp. 79-96)) Since the mean curvature, , of a cylindrical membrane is homogeneous so that its gradient along the membrane surface vanishes, ∇ = 0, and the Gaussian curvature vanishes as well, = 0, the equilibrium equation adopts a simple form, where is the membrane torque related to the membrane spontaneous curvature, % , and bending modulus, , by the relationship and Δ is the trans-membrane pressure difference. Substituting (Eq.A2) into (Eq.A1) and taking into account that, Δ = 0, for ER tubules, we obtain the equation for the tubule mean curvature The relationships (Eqs. A4,A5) describe the curvature, " , and the cross-section radius, " , of the reservoir given by (Eq. 2) of the main text.

B. Elastic energy of splay and tilt of membrane monolayer
Hemi-fission and hemi-fusion intermediates of membrane remodeling introduce deformations of inhomogeneous tilting of the hydrocarbon chains of lipid molecules Kozlov, 2000, 1998). Since the in-plane distributions of these deformation are different for the two monolayers, the elastic energy has to be computed separately for each monolayer, which can be done by using the generalized Helfrich model accounting for the tilt, splay, and saddle splay of the hydrocarbon chains as introduced in detail in the earlier works Kozlov, 2003, 2002). Here we sketch the main features of this model.
To describe the monolayer geometry, we use the dividing surface coinciding with the interface between the regions of the polar heads and hydrocarbon tails of the monolayer. The local orientation of the lipid molecules at each point of the dividing surface is characterized by the lipid director, A⃗, which is defined as a unit vector pointing from the locally averaged position of the ends of the hydrophobic lipid tails to those of the polar heads (Fig. S1). Lipid tilt, splay, and saddle splay are defined based on the variation of A⃗ along the dividing surface and the deviation of A⃗ from the unit normal vector to the dividing surface, A A⃗ .
The lipid tilt is defined as and its consequence for the monolayer internal structure can be understood as shearing of the lipid hydrocarbon chains in the direction perpendicular to the monolayer plane Kozlov, 2000, 1998;May et al., 2004) (Fig. S1).
Lipid splay, H , referred also to as the modified mean curvature, is defined as the two-dimensional divergence of the lipid director, A⃗ • A⃗, determined along the, generally, curved dividing surface. In covariant form, the lipid splay is the trace of the tensor, K L M = ∇ L M , where the sub and superscripts denote, respectively, the co-and contravariant components in the local coordinate basis of the surface. K L M is referred to as the modified curvature tensor. Therefore, Lipid saddle splay is the determinant of the modified curvature tensor, In the absence of the tilt deformation, the lipid director and the surface normal are aligned, A⃗ = A A⃗ , so that the splay, H , is reduced to the two-dimensional divergence of the surface normal, = A⃗ • A A⃗ , whose geometrical meaning is the mean curvature of the dividing surface, = 5 + -, where 5 andare the two principal curvatures of the surface. Saddle splay, N , is reduced in this case to Gaussian curvature of the surface, = 5 • -, the product of the two principle curvatures.
The elastic energy of the tilt, splay and saddle-splay deformations per unit area of the membrane plane is given by (Hamm and Kozlov, 2000), where Q is the monolayer bending modulus whose value was measured for pure lipid monolayers and varies in the range between 5 − 20 \ (Dimova, 2014;Niggemann et al., 1995;Pontes et al., 2013) ( \ ≈ 4.11 • 10 _-5 is the product of the Boltzmann constant, d , and the absolute room temperature); ̅ Q is the monolayer saddle splay modulus, which has been estimated to be negative and have an absolute value smaller but of the same order of magnitude as the bending modulus; W,Q is the monolayer tilt modulus; e,Q is the monolayer spontaneous splay, which equals the molecular intrinsic curvature of the constituting lipids. We define the ratio between the bending and saddle splay moduli as = ̅ Q Q ⁄ , which can have values in the range between -1 and 0 (Templer et al., 1998;Terzi et al., 2019). As a measure of the tilt modulus, we use = h Q W,Q ⁄ , which has a meaning of a characteristic decay length of the tilt deformation and whose estimated values vary in the range between 1nm and 2nm (Doktorova et al., 2017;Hamm and Kozlov, 1998;May et al., 2004;Nagle et al., 2015;Terzi et al., 2019). The ER membrane is composed mostly of phosphatidylcholine (PC) and phosphatidylethanolamine (PE) (Zambrano et al., 1975), both having negative intrinsic molecular curvature (Chen and Rand, 1997;Leikin et al., 1996;Szule et al., 2002), and a small amount of lysolipids, such as lysophosphatidylcholine (LPC), characterized by positive intrinsic molecular curvature (Fuller and Rand, 2001;Zambrano et al., 1975). We estimate the ER monolayers spontaneous splay by calculating the weighted average of the intrinsic molecular curvatures of the constituting lipids to be e,Q = −0.1 _5 (Kozlov and Helfrich, 1992). Since, to the best of our knowledge, there is no evidence for asymmetrical distribution of lipids components between cytosolic and luminal monolayers of ER membrane, we assume the two monolayer compositions to be identical.
The total elastic energy of the membrane is calculated by integrating the energy density (Eq. 4) over the areas of the luminal and cytosolic monolayers, where the subscripts and denote the cytosolic and luminal monolayers of the ER tubule, respectively. The energy of membrane tension, p , is given by a product between the lateral tension of each membrane monolayer, Q , and its area, which is l and o for the cytosolic and luminal area respectively, The total elastic energy is: C. Structure and energy of the hemi-fission stage of the pore sealing.
We consider the pore sealing to proceed through an intermediate stage of hemi-fission of the membrane neck constituting the rim of the pore (Fig.1E-H of the main text). The hemi-fission intermediate includes a non-bilayer structure in the middle (Fig.1F of the main text) referred to as the membrane stalk. Our goal here was to compute the elastic energy of the hemi-fission intermediate, qrs , which could be compared to the energies of a junctional knot with one pore, q5 , and of a regular junction with no pore, qt , as described in the main text.
To perform the computation of qrs , we used the results of the previous works Kozlov, 2003, 2002) in which the structure and energy of a membrane stalk were computed for various processes of membrane remodeling: abscission of an endocytic vesicle (Kozlovsky and Kozlov, 2003), hemi-fusion of flat  and curved (Martens et al., 2007) membranes. In spite of different overall configurations of these systems, the computed structure of the membrane stalk had a common feature: the middle part of the stalk could be considered as a localized sub-system referred to as the stalk-core (Fig.S1). The tilt angles in the stalk-core monolayers were equal u v in the middle of the structure, and relaxed along the mid-plane over a characteristic length of few nanometer (Kozlovsky and Kozlov, 2003), which sets the stalk-core boundary (Fig.S1).
Based on that, we approximated the hemi-fission intermediate of the pore sealing as the stalk-core embedded into the three-way tubular junction. The boundary between the stalk-core and the surrounding membrane was chosen to lie on a sphere co-centered with the stalk-core, having a radius as \ = 8 , and called below as the boundary sphere. According to the earlier works Kozlov, 2003, 2002) and verified by the present computations for the stalkcore monolayers, at this boundary the lipid tilt, practically, vanishes and the monolayer curvatures becomes small with respect to the monolayer width.
The computation of the configuration and energy of the hemi-fission intermediate included the calculation of the stalk-core energy and the energy of the surrounding membrane of the junction and minimization of the sum of the two energies with respect to the position of the interface between the stalk-core and the surrounding membrane on the boundary sphere.
By computing the stalk-core energy, we assumed that reticulons are excluded from the stalk-core because of the small dimension and non-bilayer structure of the latter so that the spontaneous splay of the stalk-core monolayers is determined solely by its constituting lipids. We further assumed that the monolayer hydrocarbon moiety is incompressible and that the monolayers dividing planes are parallel to the membrane mid-plane, meaning the lipid tail length is related to the tilt by the relationship = t √1 + -. The computation of the stalk-core energy consisted in integration of the monolayer energy density along the monolayer surfaces (Eq. B7) and minimization of the obtained overall energy with respect to the tilt distribution in each monolayer upon a requirement of the monolayer tilt angles remaining equal u v in the center of the structure (Fig.S1). The details on the parametrization of the monolayer shapes and tilt distributions we used here can be found in (Golani et al., 2021).
The computation of the energy of the surrounding membrane was performed analogously to that of the junctional knot or a regular junction, i.e. by integration of the regular density of the bilayer elastic energy ((Eq.3) of the main part) and minimization of the result with respect to the membrane shape using the Brakke's Surface Evolver (Brakke, 1992).
The computations of the stalk-core and the surrounding membrane energies were coupled through fulfillment of the boundary conditions at the interface between these two parts of the system, which included the requirement of smoothness of the membrane mid-plane profile, the continuity of the mid-plane curvature and the vanishing of the lipid tilt in the two monolayers of the stalk-core.

Figure legend
Hemifission stalk and definitions. The vertical, , axis represents the axis of rotational symmetry. Left -drawing of a cross section of a hemifission stalk, bordered with by a sphere representing the interface e between the stalk-core and the surrounding membrane. rightcomputation results of stalk core shape. Blue and red lines represent averaged lipid director of cytosolic and luminal monolayers, respectively. Definitions: A⃗ lipid director, A A⃗ Normal to dividing plane, ⃗ tilt vector, angle lipid director and ~ and angle between tangent vector to mid-plane and ~. The subscript and represent the cytosolic and luminal monolayers, respectively. Right panel represent the minimal energy configuration for the parameters: = 1.5 , = −0.5, Q = −0.1 _5 , \ = 50°, \ = 31° and \ = 8 .