# Phases, morphologies, and transitions in a membrane model for the
endoplasmic reticulum

Jaya Kumar Alageshan1, Yashodhan Hatwalne2, Rahul Pandit1
1 Centre for Condensed Matter Theory, Department of Physics,
Indian Institute of Science, Bangalore 560012, India.
2 Raman Research Institute, C.V. Raman Avenue,
Sadashivanagar Bengaluru 560080, India

###### Abstract

We introduce a novel model, comprising self-avoiding surfaces and
incorporating edges and tubules, that is designed to characterize the
structural morphologies and transitions observed within the endoplasmic
reticulum (ER). By employing discretized models, we model smooth membranes
with triangulated surfaces, and we utilize numerical variational methods to
minimize energies associated with periodic morphologies. Our study obtains
phases, their morphologies, and their transitions and examines their
dependence on the membrane chemical potential, the line tensions, and the
positive Gaussian curvature stiffness. By starting with diverse topological
structures, we explore shape variations by using Surface Evolver, while
maintaining fixed topology. Notably, we identify the region of parameter space
where the model displays lamellae, with a lattice of helical edges connecting
the layers; this resembles structures that have been observed in the rough ER.
Furthermore, our investigation reveals an intricate phase diagram with
periodic structures, including flat lamellar sheets, sponge phases, and
configurations comprising tubules with junctions, which are akin to the
morphology of the smooth ER. An estimation of lattice parameters is achieved
through fluctuation analysis. Significantly, our model predicts a transition
between homotopically equivalent lamellae, with helical edges and
configurations featuring tubules with junctions.

##  1 Introduction

The endoplasmic reticulum (ER) surrounds the nucleus in most eukaryotic cells
and plays a crucial role in cellular structure and function. Its intricate
structure and diverse functions make it a cornerstone of cellular activity.
The basic building blocks of the ER are a network of membranous tubules and
flattened sheets called cisternae, which occur throughout the cytoplasm of
eukaryotic cells [1] and exhibit distinct curvature, lipid composition, and
protein organization [2]. The interconnected tubular structures extend from
the cisternae (see Fig. 1). The tubules are dynamic and flexible structures
that branch out and fuse to form the complex three-dimensional network of the
ER. The biophysics of the ER encompasses the study of its structure, physical
properties, and the molecular dynamics underlying its functions. In the ER,
the following two main regions have been identified:

  * •

The Rough Endoplasmic Reticulum (RER) often appears as flattened cisternae and
has a rough appearance under an electron microscope because of ribosomes on
its outer surface [3]. These ribosomes are responsible for protein synthesis
and the RER helps in the folding and also the modification of proteins within
its lumen [4]. Furthermore, the newly synthesized proteins’ structural
integrity and correct folding are ensured by chaperone proteins within the ER
[5]. The lumen of the RER is contiguously connected to the lumen of the
nuclues [1].

  * •

Away from the nucleus, the ER breaks down into a filamentous structure that
forms the Smooth Endoplasmic Reticulum (SER). The SER lacks ribosomes on its
surface, so it has a smooth appearance [3] and consists predominantly of
tubules. It is involved in various metabolic processes, such as lipid
synthesis [6] and detoxification of toxins [7]. In particular, in muscle
cells, the SER stores and releases calcium ions that are necessary for muscle
contraction [8]. Calcium ions also play a pivotal role in cell signalling.
therefore, the ER plays a crucial role in their regulation and the proper
functioning of the cell [10]. The SER produces phospholipids and steroids that
are essential for cell-membrane structure and hormone production [11].

Figure 1: A schematic diagram illustrating the smooth and rough endoplasmic
reticula (SER and RER), ribosomes, the cell nucleus, nuclear pores, and the
nuclear envelope.

The ER is dynamic and continually undergoes remodelling via the fusion and the
fission of tubules and cisternae; thus, it adapts to changes in cellular
physiology and structure. Proteins such as reticulons and REEPs are
responsible for ER shaping and dynamics [15], and they maintain the ER
morphology and influence its biophysical properties. The protein insertion,
trafficking, and signalling processes within the lipid bilayer of the ER
membranes influence its mechanical properties, fluidity, and curvature [12].
In turn, these biophysical properties influence the movement of vesicles and
molecules within the ER.

The ER membrane is subject to active fusion and fission processes [16, 17].
There is a traffic of protein-carrying vesicles in the ER-Golgi-plasma
membrane, through the secretory pathways [18]. Typically, about $1000$
vesicles leave the ER every second [19]. Despite this strong remodelling
activity, which is energetically driven by the consumption of ATP/GTP, the ER
has a well-defined stable lamellar structure [20]. The fusion and fission
processes occur predominantly at the outer interface between the smooth ER and
the cell plasma [1]; so we hypothesize that local equilibrium is valid in the
bulk of the ER. Our findings demonstrate that this hypothesis leads naturally
to the structures of the ER that have been uncovered in recent experiments
[13].

Freeze-fracture studies [13] have established that the RER comprises stacks of
flat sheets connected by screw dislocations with helical edges (see Fig.
8(a)). To avoid geometric frustration, the constituents of the RER must have
an equal number of left- and right-handed dislocations; a two-dimensional (2d)
section, normal to the axes of the dislocation yields a square lattice with
alternating handedness, as shown in Fig. 8(b). A morphological transition,
from the flat lamellar sheets to lamellae with helical edges occur as the
lattice parameter, i.e., the distance between the helical edges goes to
infinity. We note, in passing, that this transition is akin to the TGB-SmA
(Twist Grain Boundary-Smectic-A) transitions in lyotropic liquid crystals
[14].

The ER is an essential and versatile cellular structure that plays a pivotal
role in various cellular functions. Thus, understanding the biophysics of the
ER is of paramount importance as it governs the organelle’s function and
interactions with other cellular components. It is important to delve into the
physical properties of ER membranes, protein dynamics, and molecular
processes, which occur within this organelle; this is essential (a) for
gaining insights into cellular physiology and disease mechanisms, particularly
those that are related to ER stress disorders, and (b) for the development of
targeted therapies. The ER works in conjunction with other organelles, such as
the Golgi apparatus, vesicles, and mitochondria, to orchestrates a complex
array of molecular activities, significantly contributing to the overall
functionality and health of eukaryotic cells.

In this study, we concentrate on the structural aspects of the ER, and
organise the remaining part of this article as follows: In Sec. 2 we present
our model. Then, in Sec. 3, we give the details of the numerical shape
variation that we employ to find the minimal-energy configuration for each of
the ER morphologies by using Surface Evolver [25]. Finally, in Sec. 4, we
present the phase diagrams and morphologies we obtain. The concluding Sec. 5
is devoted to a discussion of the significance of our results and possible
limitations of our study.

##  2 The Model

The ER has a complex multi-scale structure, with diverse morphologies and
functions that span length scales ranging from nanometers to micrometres. The
ER membranes are composed of lipid bilayers that have a thickness $\simeq 4-6$
nanometers. Protein folding, chaperone interactions, and post-translational
modifications and other molecular processes occur within the ER lumen and
determine the double-bilayer thickness, which is typically $\simeq 30-50$
nanometers [26]. Furthermore, the ribosomes that decorate the outside of the
RER are $\simeq 20-30$ nanometers in diameter; they are responsible for
protein synthesis, which occurs at the nanometer scale within the ER
cisternae. Throughout the cytoplasm, the SER forms an extensive network, which
extends over micrometre scales in most eukaryotic cells. Vesicles within the
ER network are involved in the transport of proteins, lipids, and other
molecules, within the ER and between the ER and other cellular compartments,
such as the Golgi apparatus. The sizes of these vesicles can vary, but their
diameters lie typically in the range $25-50$ nanometers.

The ER’s diverse length scales, which are vital for its functions, allow it to
perform tasks ranging from molecular-level protein synthesis and modification
to coordinating cellular processes over larger distances. This ability to
cover multiple scales enables efficient communication and coordination within
the cell and underscores the ER’s importance in maintaining cellular
homeostasis. The ER is complex: it has many types of protein molecules, with
non-trivial properties, embedded within its membrane.

From the structural perspective, the principal building blocks of ER are the
fragments of sheets and tubules made up of double lipid bilayers [28]. We
model the double bilayer as a self-avoiding two-dimensional (2d) surface and
the tubules as one-dimensional (1d) curves (see fig. 2), in concurrence with
Ref. [37]. The edges of the double bilayer (red) and the tubules (blue) are
decorated with reticulons [43]. The quantitative calculation of the free
energy for a self-avoiding fluid membrane with tubules is complex: it involves
the use theoretical models that have been motivated by experimental
observations. We use a continuum-mechanics-based approach to study the phases,
transitions, and morphologies of such membranes under different biophysical
conditions.

Figure 2: (a) A colour-coded schematic diagram of the elements of the double
bilayer membrane in our model: ($\mathcal{A}$)\- 2D surface bulk;
($\partial\mathcal{A}$)\- surface edge; ($\mathcal{T}$)\- tubule; (orange)\-
junction of tubule and edge; (b) the corresponding equivalent 2d surface and
1d curve in our model.

In our study, we consider the Helfrich Hamiltonian [36] for a fluid membrane
together with line tension, and self-avoidance given below:

$\displaystyle\mathcal{H}=\underbrace{\gamma\oint_{\partial\mathcal{A}}ds\;\;+\;\;\;\Gamma\int_{\mathcal{T}}ds}_{\text{Line-
tension}}\;\;+\;\;\underbrace{\iint_{\mathcal{A}}\left(\frac{\kappa}{2}\;H^{2}+\kappa_{G}K\right)\sqrt{g}\;d^{2}\sigma}_{\text{Bending}}$
$\displaystyle+\;\;\underbrace{\mu\iint_{\mathcal{A}}\sqrt{g}\;d^{2}\sigma}_{\text{Surface
tension}}\;\;+\;\;\underbrace{4\;\mathcal{B}\iint_{\mathcal{A}}\left(\left(\frac{\Delta}{\delta}\right)^{12}-\left(\frac{\Delta}{\delta}\right)^{6}\right)\sqrt{g}\;d^{2}\sigma}_{\text{Self-
avoidance}}\,,$ (1)

where $\mathcal{A}$ is the area of the surface elements and
$\partial\mathcal{A}$ are the corresponding boundaries; $\mathcal{T}$’s refer
to the tubules [see Fig. 2]. We model the bending energy of the membrane by
using the Helfrich terms [36], which depend on the mean curvature $H$ and
Gaussian curvature $K$. The surface- and line-tension energies arise from the
stretching of the membrane’s constituent molecules and depend on the tension
or stress applied to the membrane. The term with coefficient $\mu$ penalizes
deviations from the membrane’s preferred surface area, $\gamma$ is the line
tension of the surface edge, and $\Gamma$ is the line tension of the tubule.
The membrane’s edges can coalesce to form tubules, as shown in Fig. 3, so a
tubule is effectively two merged edges (see Fig. 3), whence $\Gamma=2\gamma$.
We neglect the energetics of the tubule-edge junction,  orange-region in Fig.
2, because it involves higher-order gradients than those we retain in Eq. 1.
The self-avoidance interactions between membranes are captured by the Lennard-
Jones-type potential [38], with coefficient $\mathcal{B}$, where $\delta$ is
the inter-layer spacing and $\left(2^{1/6}\Delta\right)$ is the preferred
spacing.

Figure 3: Tubule formation: (a) A schematic where the edges of the double
bilayer join together to form a tubule; (b) The effective process when
modelled with a 2d surface and 1d curve preserves the Euler characteristic
$\chi$, with $V$-vertices, $E$-edges, and $F$-faces.

For a free-standing membrane, with $\delta\gg\Delta$, a variation of the
Hamiltonian (1), with respect to its shape, leads to [39, 40]

$\displaystyle 4\kappa H(H^{2}-K)-2\mu H+\kappa\nabla^{2}(2H)=0\,,$ (2)

in the bulk; and the corresponding boundary conditions are:

$\displaystyle\left.\left[2\kappa
H+\kappa_{G}\kappa_{n}\right]\right|_{\partial\mathcal{A}}$ $\displaystyle=$
$\displaystyle 0\,,$
$\displaystyle\left.\left[-2\kappa\partial_{\perp}H+\mu\kappa_{n}+\kappa_{G}\partial_{\parallel}\tau_{g}]\right]\right|_{\partial\mathcal{A}}$
$\displaystyle=$ $\displaystyle 0\,,$ $\displaystyle\left.\left[2\kappa
H^{2}+\kappa_{G}K+\mu+\gamma\kappa_{g}\right]\right|_{\partial\mathcal{A}}$
$\displaystyle=$ $\displaystyle 0\,,$ (3)

where $\kappa_{n}$ and $\kappa_{g}$ are, respectively, the normal and the
geodesic curvatures of the boundary curve, with respect to (w.r.t.) the
surface, $\tau_{g}$ is the torsion,
$\partial_{\parallel}\equiv\hat{\mathbf{t}}\cdot\nabla$ is the derivative
tangential to the boundary curve, and
$\partial_{\perp}\equiv\hat{\mathbf{n}}_{s}\cdot\nabla$ is the derivative
directed along the outward normal to the surface boundary.

Entropic contributions arise because of the membrane’s ability to explore
various conformations via thermal fluctuations that lead to scale-dependent
renormalizations of the coefficients in Eq.( 1). Hence, we consider the free-
energy contribution associated with the fluctuations of the membrane. which
can be written as

$\displaystyle\mathcal{F}$ $\displaystyle=$
$\displaystyle-\ln{\mathcal{Z}}/\beta\,,\text{ where }\beta=1/k_{B}T\,,\text{
and }$ $\displaystyle\mathcal{Z}$ $\displaystyle=$
$\displaystyle\int\mathcal{D}[\mathcal{C}]\;\exp(-\beta\mathcal{H})\,,$ (4)

where $\mathcal{Z}$ is the functional integral over all possible
configurations of the membrane shape, with $\mathcal{D}[\mathcal{C}]$ the
measure. We assume that the free energy has the same form as the Hamiltonian
(1), with these renormalized coefficients. The inclusion of short-wavelength
fluctuations renormalizes the bending moduli at the length scale $L$ as
follows [31]:

$\displaystyle\tilde{\kappa}$ $\displaystyle=$
$\displaystyle\kappa_{0}-\frac{3}{4\pi\beta}\;\ln\left(L/a\right)\,,$ (5)
$\displaystyle\tilde{\kappa}_{G}$ $\displaystyle=$
$\displaystyle\kappa_{G0}+\frac{5}{6\pi\beta}\ln\left(L/a\right)\,,$ (6)

where $a$ is the molecular-length cut-off; $\kappa_{0}$ and $\kappa_{G0}$ are
bare quantities, and $\xi_{\kappa}=a\;\exp\left(4\pi\beta\kappa_{0}/3\right)$
and $\xi_{\kappa_{G}}=a\;\exp\left(6\pi\beta\kappa_{G0}/5\right)$ as the
corresponding thermal persistence lengths. The renormalized surface tension is

$\displaystyle\tilde{\mu}=\mu_{0}+\frac{\mu_{0}}{4\pi\beta\;\kappa_{0}}\ln\left(L/a\right)\,.$
(7)

The self-avoidance term also undergoes renormalization, with

$\displaystyle\tilde{\Delta}=a+1/\sqrt{\tilde{\kappa}|\tilde{\mu}|\beta^{2}}\,.$
(8)

Furthermore, if $\delta=2^{1/6}\delta_{0}+\Delta$, then the quadratic
expansion of the compression term is $\sim
2^{17/3}\mathcal{B}\;(\Delta/\delta_{0})^{2}$. Therefore,the renormalized
layer compressibility is [42]

$\displaystyle\tilde{\mathcal{B}}=\mathcal{B}_{0}+\frac{\pi^{2}\tilde{\Delta}}{2^{2/3}\;\kappa\;\beta^{2}\;(\tilde{\Delta}-\delta_{0})^{4}}\,.$
(9)

Explicit expressions for the renormalization of the line tension in fluid
membranes can be intricate and model-dependent, and there is no universally
accepted formula for them. So, for simplicity, we retain the bare value of
this line tension.

In Sec. 4, we consider different morphologies that we classify on the basis on
their topologies. Furthermore, to simplify our analysis, we study only
periodic structures. Therefore, to compare the stabilities of the different
phases and their morphologies, we use the free energy per unit volume. All the
energy-minimizing solutions to the total free energy in Eq. (1) are not
amenable to analytic solutions because of the geometry-related non-
linearities. It is imperative, therefore, to look for approximate numerical
solutions. In the next Section, we give a brief overview of the Surface
Evolver package that we use for such numerical simulations.

##  3 Methods: Surface Evolver

Surface Evolver (SE) is an open-source numerical package used in computational
physics and mathematics to model and simulate the behaviours of surfaces and
interfaces. Initially devised by Kenneth Brakke for mean-curvature-flow
analysis [23], this tool has gained wide usage across diverse realms,
prominently in exploring the physics of soap films, minimal surfaces, and
various interfacial phenomena.

Within SE, smooth 2D surfaces are represented as piece-wise-linear
triangulated surfaces embedded within Euclidean 3D space.

Key attributes of SE encompass:

  1. (a)

Numerical Surface Energy Minimization: Enables the discovery of minimal
surfaces and configurations, while adhering to specified constraints.

  2. (b)

Versatile Geometry Simulation: Facilitates the simulation of intricate
geometries, encompassing diverse surface types, shapes, and interactions.

  3. (c)

Dynamic Interface Simulation: Allows the study of surface evolution over time,
yielding insights into growth, deformation, and phase transitions.

  4. (d)

Parameter Control: Offers manipulation of parameters like surface tension,
boundary conditions, and constraints, thereby enabling the exploration of
surface behaviours in varying environments and physical conditions.

SE harnesses different numerical techniques [energy-minimization algorithms]
to ascertain surface-equilibrium configurations. These algorithms calculate
the most energy-efficient shapes and structures for surfaces, with specified
boundary conditions and constraints. Its applications in material science
include studies of surface energies, stability assessments, and interface
analyses. Moreover, SE serves as a valuable tool in comprehending surface
interactions within biological systems, so it can aid in the modelling of
biological membranes and cell structures.

Working with SE typically involves scripting, in its dedicated programming
language, which allows users to define surface energies, constraints, and
boundary conditions. In our specific study, we introduce an auxiliary term,
represented as the self-avoidance $(\mathcal{B})$ in Eq. 1, into the SE
calculations to account for membrane interactions. This term incorporates the
interaction between membranes by estimating the inter-layer distance, denoted
as $\delta$, within the lamellar phase. We determine this distance from
triangles along the vertex normal, as illustrated in Fig. 4(a).

(a) (b)

Figure 4: (a) Schematic diagram demonstrating the calculation of the distance
between a vertex and neighbouring layers ($\delta$), along the vertex-normal.
(b) Lennard-Jones-type potential between the surfaces.

SE minimizes the energy by numerical shape variation, which is implemented by
using a conjugate-gradient descent [53], w.r.t. the coordinate variables of
the vertices, as follows:

$\displaystyle\frac{\partial\mathbf{x}_{\alpha}}{\partial\tau}=\mathbf{\Upsilon}_{\alpha\beta}\;.\;\frac{\delta\mathcal{F}}{\delta\mathbf{x}_{\beta}}\,,$
(10)

where $\mathbf{x}_{\alpha}$ is the 3d position vector of the vertex labeled by
$\alpha$, $\mathbf{\Upsilon}_{\alpha\beta}$ represents the vertex-dependent
mobility matrix, and $\tau$ is the iteration scale. Since we are modeling a
fluid membrane, the redundant evolution, leading to re-parameterization (in-
plane movement) of the vertices in the bulk, is eliminated by redefining the
mobility matrix as

$\displaystyle\mathbf{\Upsilon}_{\alpha\beta}\rightarrow\tilde{\mathbf{\Upsilon}}_{\alpha\beta}\;:=\;\mathbf{n}_{\alpha}\;.\;\mathbf{\Upsilon}_{\alpha\beta}\,,$
(11)

where $\mathbf{n}_{\alpha}$ is the normal at the vertex labeled $\alpha$. The
above process enforces only the normal motion of surfaces. Furthermore, we set
all the vertex mobilities to be equal, i.e.,
$\tilde{\mathbf{\Upsilon}}_{\alpha\beta}=\Upsilon\;\delta_{\alpha\beta}\;\vec{\mathbf{n}}_{\alpha}$,
where $\delta_{\alpha\beta}$ is a Kronecker delta. The scale factor $\Upsilon$
is estimated by using Newton’s method; and it is proportional to the inverse
of the Hessian [56]. We also note that the Hessian of the edge vertices
dominates over that of the bulk, leading to slow the convergence of edge
modes. To adjust for the offset, after every iteration of the conjugate
gradient descent we perform the variation of only the edge vertices, while
keeping all the vertices in the bulk fixed.

We explore various morphologies detailed in Sec. 4. These distinct shapes
possess unique topologies, which can be classified by their homology groups
[22]. Our study is initiated with triangular meshes that represent fundamental
mesh structures corresponding to each morphology. The process of energy
minimization employs iterative triangle and edge refinements, alongside vertex
evolution guided by Eq. (10). This iterative refinement seeks stable surface
configurations, ensuring that vertex-coordinate changes remain within $5\%$ of
the unit cell size. Triangle refinements enhance surface resolution by
subdividing large triangles into smaller ones, effectively capturing intricate
surface details. We achieve triangle refinement and reshaping by using Pachner
moves [24], depicted in Fig. 5. These transformations preserve the
triangulation’s topological properties.

Figure 5: A schematic diagram of Pachner moves that we use to refine and re-
triangulate the mesh: (a) refinement; (b) triangle weeding; (c) edge weeding;
(d) edge flip.

The systems that we consider are periodic in x, y, and z directions, with
orthorhombic or cubic unit cells; we implement this in SE by using the
SYMMETRY option [52]. We obtain the energy per unit volume to compare the
stabilities of different phases of the system.

##  4 Results

Fluid membranes, which avoid self-intersections and incorporate edges, possess
the capacity to self-organize into diverse structural configurations, to adapt
to distinct environmental conditions, and thus yield various phases [see,
e.g., Refs. [46] and [57]]. Our model introduces, in addition, the possibility
of tubule formation; this enriches the spectrum of possible phases in the
system.

In our subsequent explorations of membrane stability, we consider cases in
which $\kappa>0$ significantly surpasses $k_{B}T$. We focus on
$\kappa_{G}=\kappa$; although the stability of fluid membranes allow for both
positive and negative values of $\kappa_{G}$, we exxplore only the region
$\kappa_{G}>0$.

We investigate the stability of diverse morphologies, by using our continuum
model (1) with renormalized coefficients (as discussed above). We utilize
Surface Evolver (SE) to capture the mean shapes. We focus on the region
$\mathcal{B}\sim\beta$; and $\delta_{0}$, which is equal to the double-layer
width, induces a repulsive potential between surfaces. This enables the
fluctuations to counterbalance attractive forces and allows surfaces to exist
independently without condensation.

Our exploration of the stabilities of different morphologies begins with the
initialization of minimal triangular meshes, with either orthorhombic or cubic
symmetry, that are designed to yield the desired topology. These meshes evolve
via SE to minimize the energy [Eq. (1)]. As the vertex-coordinate change falls
below a $5\%$ threshold relative to the unit-cell size, we iteratively refine
and re-triangulate the mesh by using the Pachner moves described above. We
continue this until further refinement fails to induce deformations [given the
stipulated threshold].

We measure energies in units of $\beta$ and normalize lengths by
$\xi_{\kappa}$. Figure 6 illustrates the phase diagram that we have obtained
for model (1) with the representative value $\kappa=5\beta$ in the two-
dimensional parameter space $[\exp{-(\beta\gamma)},\beta\mu]$; this shows a
variety of phases with non-trivial structures, that minimize the energy in
different regions of this parameter space. We proceed to present an in-depth
analysis of these phases and the phase diagram; and we elucidate the interplay
between $\mu$ and $\gamma$ that stabilises different phases.

Figure 6: The phase diagram for model (1) with the representative value
$\kappa=5\beta$ in the two-dimensional parameter space
$[\exp{-(\beta\gamma)},\beta\mu]$; this shows a variety of phases with non-
trivial structures, that minimize the energy in different regions of this
parameter space. considered.

  * ✮

Lamellar phases: The system exhibits a strong inclination to form lamellae
when $\mu$ has large negative values; this favours the and filling of the
domain with layered surfaces (or lamellae). In describing the membrane stack
within this phase, a persistence length ($\xi_{\kappa}$) signifies the average
distance over which neighbouring membranes collide because of thermal
fluctuations [50]. These collisions result in a steric repulsion between the
membranes and, consequently, renormalize the inter-layer spacing to
$\Delta\approx\xi_{l}\sim 1/\sqrt{\kappa|\mu|\beta^{2}}$.

    * ❏

Flat lamellae (FL): When the energy ($\gamma$), associated with edges, becomes
excessively high, the system actively avoids edges and the lamellae adopt a
flat mean profile (see Fig. 7) to minimize the bending energy. In terms of
topology, the phase comprises disconnected sheets, displaying a long-range
orientational order, of the surface normal, and quasi-long-range positional
order along this normal direction.

If $\lambda_{\perp}$ and $\lambda_{\parallel}$ are the unit cell dimensions in
the directions along the layer’s normal and lateral directions, then the free-
energy and the energy density reduce to:

$\displaystyle\mathcal{F}_{FL}$ $\displaystyle=$
$\displaystyle\tilde{\mu}\>\lambda_{\parallel}^{2}+4\>\tilde{\mathcal{B}}\left(\left(\frac{\tilde{\Delta}}{\lambda_{\perp}}\right)^{12}-\left(\frac{\tilde{\Delta}}{\lambda_{\perp}}\right)^{6}\right)\lambda_{\parallel}^{2}\,;$
$\displaystyle f_{FL}$ $\displaystyle=$
$\displaystyle\mathcal{F}_{FL}/(\lambda_{\perp}\lambda_{\parallel}^{2})\,;$
(12)

$\lambda_{\parallel}=L=\xi_{\kappa}$; and we minimize $f_{FL}$ w.r.t.
$\lambda_{\perp}$.

(a) (b)

Figure 7: Structure of flat lamellae: (a) an orthorhombic unit cell; (b) the
repetition of unit cells that form the final configuration of a flat lamellar
phase.

    * ❏

Lamellae with helical edges (LH): In the lamellar phase, as the $\gamma$ is
reduced to a large negative value, the lamellae prefer to break and form
sheets with edges. But simple circular holes or edges are unstable, so the
system develops topologically non-trivial helical edges, exhibiting a
morphology that resembles a staircase (see Fig. 8(a)). This morphology
consists of repeating, alternating layers, or terraces, connected via the
helical edge. The helices are chiral; and opposite chirality helical edges can
condense to form a square lattice as we show in Fig. 8(b).

We use the continuum elasticity model (1) to obtain the free energy. The
lattice spacing, $\alpha$, of the helical edges is stabilized by the
fluctuations; and it is determined self-consistently. Energetically, smaller
values of $\alpha$ would lead to a higher density of edges and hence a lower
energy. In our model, the lattice-spacing limit is set by the validity of the
lamellar structure, which is $\alpha\geq\xi_{\kappa}$, and we saturate this
limit. Furthermore, by approximating the surface to be a helicoid, which is a
minimal surface, we satisfy the force balance in the bulk [Eq. (2)]. The
boundary conditions in Eq. (3) relate the radius of the helical boundary to
the interlayer spacing via [54]

$\displaystyle-\frac{\tilde{\kappa}_{G}\;(\tilde{\xi}_{l}/2\pi)^{2}}{(r^{2}+(\tilde{\xi}_{l}/2\pi)^{2})^{2}}+\frac{\tilde{\gamma}\;r}{(r^{2}+(\tilde{\xi}_{l}/2\pi)^{2})}+\tilde{\mu}=0\,,$
(13)

where $r$ is the core radius (see Fig. 8(a)), and
$\xi_{l}=1/\sqrt{\kappa|\mu|\beta^{2}}$. We begin with a coarse triangulated
surface, as in Fig. 9(a), with $\lambda_{\parallel}=2\alpha$; we then let SE
evolve the vertices with multiple mesh refinements to reach the stable
structure, shown in Fig. 9(c).

(a) (b)

Figure 8: Schematic diagrams: (a) (right panel) the Helicoid core; (left
panel) the helicoid, with the core removed and showing the corresponding
helical edge; (b) a unit cell of a pair of left(L)- and right(R)-handed
helicoids arranged on a square lattice [left panel: 3d view; right panel: top
view].

(a) (b) (c)

Figure 9: Schematic diagrams showing the structure of lamellae with helical
edges: (a) An orthorhombic unit cell with a triangulated surface that has two
left- and two right-handed helical edges, as in Fig. 8(b). This surface serves
as the input initial condition for the SE. (b) The repetition of unit cells
with orthorhombic symmetry. (c) The final configuration after energy
minimization via SE.

Studies such as those of Refs. [13] and [39] associate the LH structures with
the rough ER, but they do not take into account the self-avoidance of the
membranes. Instead, to enforce a preferred inter-layer spacing in the case of
the staircase-like structure, they use a spontaneous geodesic-curvature term
of the form $\int_{\partial\mathcal{A}}(\kappa_{g}-\bar{\kappa}_{g})^{2}\;dl$,
where $\kappa_{g}$ and $\bar{\kappa}_{g}$ are the geodesic curvatures of the
boundary and the corresponding spontaneous value induced by the protein that
occupies the helical boundaries (and acts as springs with a preferred pitch).
But for a lattice of helicoids with lattice spacing greater than
$\zeta_{\kappa}$, thermal fluctuations render the structure unstable without
self-avoidance. Furthermore, the terms in $\mathcal{F}_{geo}$ is of higher
order compared to $\int K\;d^{2}\sigma$, as a consequence of the Gauss-Bonnet
theorem, which relates the geodesic curvature and Gaussian curvature through
the relation,

$\displaystyle\int\int_{\mathcal{A}}\mathcal{K}\;\sqrt{g}\;d^{2}\sigma=2\pi(1-g)+\oint_{\partial\mathcal{A}}k_{g}ds\,,$
(14)

where $\mathcal{K}$ is the Gaussian curvature, $g$ is the genus of the
surface, and $k_{g}$ is the geodesic curvature. Now, if we use the $\oint
k_{g}^{2}ds$ term, we must include, for consistency, higher-order terms in the
Gaussian curvature term in the bulk. Furthermore, the inclusion of the
quadradic geodesic-curvature term stabilizes stacks of membranes with holes of
radius $\sim 1/\kappa_{g}$; these were not reported in Ref. [13].

  * ✮

Sponge phase: If $\kappa_{G}>0$, and we use a high value of $\gamma$ and a
small value of $|\mu|$, the interlayer spacing decreases notably.
Consequently, lamellae develop the propensity to merge, giving rise to
intricately interconnected membrane networks that create a porous,
labyrinthine structure called a sponge phase. In contrast to the flat lamellar
phase, the sponge phase likes to bend and maximize negative Gaussian
curvature. This phase adopts shapes akin to minimal surfaces with zero mean
curvature without edges. The resulting sponge phase is characterized by a
disordered and porous structure devoid of long-range order.

We examine periodic sponge structures with cubic symmetry; these are often
referred to as a plumber’s nightmare because of its complex morphology [46,
57] (see Fig. 11(c)). To find the optimal plumber’s-nightmare structure, we
start with a triangular mesh, as shown in Fig. 11(a), and let the SE evolve
the surface, iteratively refining the structure until the convergence criteria
are satisfied. These structures exhibit three-dimensional crystalline order
and are such that their area $\sim\lambda^{2}$, with a factor which we
determine by using the plot in Fig. 10, that we obtain by using SE. Since it
is a minimal surface, $H=0$, and the integral of the Gaussian curvature per
unit cell is independent of $\lambda$ [it is a function only of the genus].
Furthermore, for the plumber’s-nightmare structure, the entropic contribution
is constant [46]. Hence, the total free energy density for the
plumber’s-nightmare structure for a unit cell of size $\lambda(=L)$ is

$\displaystyle f_{SP}\approx
A_{1}\tilde{\mu}/\lambda+A_{2}-A_{3}\tilde{\kappa}_{G}/\lambda\,,$ (15)

where $A_{1}$ is the slope of the linear fit in Fig. 10(a), $A_{2}$ is the
constant entropic contribution, and $A_{3}$ is the slope from Fig. 10(b).

(a) (b)

Figure 10: Scaling relations for the plumber’s-nightmare phase obtained by
using SE: (a) Plot of the area vs the unit-cell length ($\lambda$); (b) Plot
of the total Gaussian curvature vs $\lambda$.

Note that, for $\kappa_{G}>0$ and $\mu<0$, $f$ decreases with $\lambda$, so we
have a UV runaway, and the stability of the structure is dictated by the
short-distance cut-off; in our numerical study, the double bi-layer thickness
provides this cut-off, i.e. $\lambda\sim a$.

(a) (b) (c)

Figure 11: Structure of the plumber’s nightmare: (a) Triangulated surface in a
cubic unit cell that serves as the input to SE; (b) the repetition of unit
cells with cubic symmetry; (c) final configuration after energy minimization
via SE.

  * ✮

Droplet phase: For positive $\mu$ and positive $\gamma$, the free energy can
be reduced by the formation of spherical droplets whose radii are determined
by fluctuations. The characteristic length scale, associated with $\mu>0$,
determines the droplet size $\xi_{\mu}=1/\sqrt{4\pi\mu\beta}$; and the shape
fluctuations determine the spacing between the droplets, which is of the order
of $\xi_{\kappa}$. In general, for $\kappa_{G}>0$, the droplets proliferate
[high genus topology]. For large $\mu$, the system prefers to avoid membranes
by breaking up into dilute non-spherical droplets. We neglect the self-
avoidance term because of the dilute-droplet assumption. To perform a self-
consistent comparison between different morphologies, while keeping the
analysis simple, we use a mono-disperse distribution of spheres with
$\xi_{\mu}$ as the cut-off radius, and the lattice spacing of the order of
$\lambda=\xi_{\kappa}$. Then, the corresponding free-energy density is [32]

$\displaystyle f_{D}\approx
4\pi\>\tilde{\mu}\>\xi^{2}_{\mu}/\lambda^{2}+4\pi\>(\tilde{\kappa}/2+\tilde{\kappa}_{G})/\lambda^{2}-\ln(\xi_{\kappa}/\xi_{\mu})/(\beta\lambda^{2})\,.$
(16)

  * ✮

Tubule phases: When $\mu$ is positive and we have a negative value of
$\gamma$, tubules dominate the morphology, which has string-like undulations
and translational entropy of the junctions [33].

    * ✧

Tubule nematic: For large positive values of $\mu$, the system avoids forming
surfaces and consists only of tubules. For the simple model that we consider,
with only line tension for the tubules, they can align to give rise to nematic
order, which optimizes packing. We consider a hexagonal packing of tubules and
estimate the free-energy density.

    * ✧

Tubules with Y-junction: In the case of negative $\mu$ close to zero and
negative values of $\gamma$, the tubules have string-like undulations and form
3-way Y-junctions, with a flat triangular patch of surface, as shown in Fig.
12(a). The Y-junctions act as particles and give rise to translational entropy
as in the Tlusty-Safran transition [33]. The break-up of the membrane into
3-way junctions is achieved through the transition in Fig. 3. In our study, we
consider the cubic arrangement of Y-junctions shown in Fig. 12.

(a) (b)

Figure 12: Schematic diagram showing the structure of tubules with
Y-junctions: (a) Triangulated surface in a cubic unit cell that serves as our
input to SE. The tubules and surface edges are highlighted as 3d red tubes,
but in SE they are 1d curves; (b) final configuration after energy
minimization via SE, such that the unit cell is repeated with cubic symmetry.

(a) (b)

Figure 13: A schematic diagram showing different possible topologies with
4-way junctions, with a Burger’s loop highlighted in black: (a) fragments of
flat membranes that connect to form a closed loop; (b) saddle membranes
connected by tubules.

    * ✧

Tubules with saddle junctions: If $\kappa_{G}>0$, for large negative $\gamma$
and small positive $\mu$, the edges in the staircase can grow and coalesce to
form a network of tubules connected by small patches of sheets, to form
saddle-junctions. The patches form pieces of saddle minimal surfaces that have
negative Gaussian curvature and zero mean curvature (see Fig. 14(c)). Here the
shape and size of the patch are determined by the competition between the
positive surface tension, $\mu$, which tries to minimize the area, and the
negative contribution from the Gaussian curvature. To embed a negative
Gaussian curvature surface in Euclidean three-dimensional space, the integral
of the Gaussian curvature scales sub-linearly in the surface area [55] and
hence leads to a finite size. We consider a lattice of such saddles and
estimate their shape and size numerically.

(a) (b) (c)

Figure 14: Schematic diagram showing the structures of saddle junctions: (a)
Triangulated surface in a unit cell that serves as our input to SE. The
tubules and surface edges are highlighted as 3d red tubes, but in SE they are
1d curves; (b) final configuration after energy minimization using SE, such
that the unit cell is repeated with orthorhombic symmetry; (c) the saddle
shape of the junction with net negative Gaussian curvature minimizes the
energy for $\kappa_{G}>0$.

The Saddle junctions undergo a symmetry-breaking transition into the
Y-junction phase, through the tubule formation process described in Fig. 3,
wherein a membrane breaks up into pieces of triangles as shown in Fig. 15.

(a) (b)

Figure 15: Symmetry breaking transition from a saddle junction (a) to a
Y-junction in (b).

The phase diagram in Fig. 6 summarizes the analysis presented above. As
observed in Ref. [13], the lamellar phase with helical edges corresponds to
the rough ER and the tubule phases correspond to the smooth ER. Through our
model, we highlight that there exists a structural transition between the
lamellae with helical edge, saddle-junctions, and Y-junctions such that the
interior of the double bilayer, the lumen, maintains contiguous connectivity.

##  5 Conclusions

Our study demonstrates that the diverse structural phases witnessed within the
endoplasmic reticulum (ER) can be effectively replicated by using our model
(1) that encompasses interacting fluid membranes featuring edges. By employing
energy-minimization techniques and variational methods, we derive the stable
shapes characterizing these distinct structures.

The resultant phase diagram, derived through numerical optimization, unveils
transitions from lamellar configurations to helicoidal stacks, closely
resembling the experimentally observed structures within the rough endoplasmic
reticulum (RER). Furthermore, our findings highlight the existence of
homotopy-preserving transitions, illustrating transformations between
helicoidal stacks, saddle junctions, and Y-junctions.

###  5.1 Discussion

The intricate nature of endoplasmic reticulum (ER) structures, coupled with
the diverse array of factors influencing its morphology poses a significant
challenge in constructing a comprehensive phase diagram that shows the varied
ER morphologies under diverse conditions. In this context, several critical
points require attention:

  1. (a)

The diverse morphological variations across different cell types and
physiological contexts.

  2. (b)

Understanding the interconversion and remodelling dynamics between tubules and
sheets.

  3. (c)

Investigating the influence of cellular activity on ER morphology.

  4. (d)

Unraveling the intricate interplay among lipid bilayers, proteins, and
cellular signalling pathways.

  5. (e)

Deciphering the role of proteins in sculpting and determining ER morphology.

Advancements in super-resolution microscopy, live-cell imaging techniques, and
computational modelling continuously contribute valuable insights into ER
morphology. Recent research endeavours have made strides in elucidating phase
transitions and morphological alterations within the ER. However, the
construction of a detailed phase diagram that comprehensively integrates these
findings remains a work in progress.

Limitations of our model: The ER system is, in general, active and far from
equilibrium. While our work identifies key influential conditions and factors
impacting ER morphology, crafting an exact phase diagram that integrates
dynamic activity and out-of-equilibrium components remain an ongoing
challenge. The inclusion of active non-equilibrium phenomena can be
accomplished by using coloured-noise-based membrane fluctuations, as
discussed, e.g., in Ref. [21]; such studies have not yet yielded the
structures we have obtained.

##  6 Acknowledgements

YH and JKA express gratitude to V.A. Raghunathan for valuable inputs and
insightful discussions. JKA thanks Chaitra Hegde for engaging in relevant
discussions.

We acknowledge the Department of Science and Technology, New Delhi, for their
support through the DSTO1359 grant. RP acknowledges the support received from
the Science and Engineering Research Board (SERB) and the National
Supercomputing Mission (NSM). We also extend our gratitude to SERC (IISc) for
providing computational resources.

We acknowledge the use of AI tools like Grammarly and ChatGPT for sentence
rephrasing and grammar checks. Subsequently, the material underwent meticulous
proofreading to ensure precision and rectify any errors. Our process includes
thorough reviews and edits to ensure accuracy, relevance, and coherence in the
finalized text.

## References

  * [1] Park, Seong H., and Craig Blackstone. ”Further assembly required: construction and dynamics of the endoplasmic reticulum network.” EMBO reports 11.7 (2010): 515-521.
  * [2] Alberts, Bruce, et.al., Molecular Biology of the Cell: Seventh International Student Edition with Registration Card. WW Norton $\&$ Company, 2022.
  * [3] Shibata Y, Shemesh T, Prinz WA, Palazzo AF, Kozlov MM, Rapoport TA: Mechanisms determining the morphology of the peripheral ER. Cell, 143:774-788, (2010).
  * [4] Braakman I, Hebert DN. ”Protein folding in the endoplasmic reticulum”, Cold Spring Harb Perspect Biol.;5(5):a013201 (2013).
  * [5] Caramelo, J.J.; Parodi, A.J. ”A sweet code for glycoprotein folding”, FEBS Lett., 589, 3379–3387, (2015).
  * [6] Vance, J.E. Phospholipid Synthesis and Transport in Mammalian Cells. Traffic, 16, 1–18, (2015).
  * [7] Maxfield FR, Wüstner D, ”Intracellular cholesterol transport”, The Journal of Clinical Investigation. 110 (7): 891–8, (2002).
  * [8] Baumann O, Walz B, ”Endoplasmic reticulum of animal cells and its organization into structural and functional domains”, Int Rev Cytol 205: 149–214, (2001).
  * [9] oyoshima C, Nakasako M, Nomura H, Ogawa H, ”Crystal structure of the calcium pump of sarcoplasmic reticulum at 2.6 A resolution”. Nature. 405 (6787): 647–55, (2000).
  * [10] Scorrano, et.al., BAX and BAK regulation of endoplasmic reticulum Ca2+: A control point for apoptosis. Science , 300, 135–139, (2003).
  * [11] Fawcett, D. The Cell, 2nd ed.; W. B. Saunders Company: Philadelphia, PA, USA, 1981; ISBN 9780721635842.
  * [12] Voeltz GK, Prinz WA, Shibata Y, Rist JM, Rapoport TA, ”A class of membrane proteins shaping the tubular endoplasmic reticulum”, Cell 124:573–86, (2006).
  * [13] Mark Terasaki, Tom Shemesh, et.al., Stacked Endoplasmic Reticulum Sheets Are Connected by Helicoidal Membrane Motifs, In Cell, Volume 154, Issue 2, 2013, Pages 285-296.
  * [14] Kamien, Randall D., and T. C. Lubensky. ”Chiral lyotropic liquid crystals: TGB phases and helicoidal structures.” Journal de Physique II 7.1 (1997): 157-163.
  * [15] Friedman JR, Voeltz GK: The ER in 3D: a multifunctional dynamic membrane network. Trends Cell Biol, 21:709-717, (2011).
  * [16] Chen, S., Novick, P., Ferro-Novick, S. ”ER structure and function.” Current Opinion in Cell Biology 25, 428–433, (2013).
  * [17] Borgese, N., Francolini, M., Snapp, E. ”Endoplasmic reticulum architecture: structures in flux”. Current Opinion in Cell Biology 18, 358–364, (2006).
  * [18] Lodish H, Berk A, Zipursky SL, et al. Molecular Cell Biology. 4th edition. New York: W. H. Freeman; (2000).
  * [19] Felix T. Wieland, Michael L. Gleason, Tito A. Serafini, James E. Rothman, The rate of bulk flow from the endoplasmic reticulum to the cell surface, In Cell, Volume 50, Issue 2, 1987, Pages 289-300, ISSN 0092-8674, https://doi.org/10.1016/0092-8674(87)90224-8.
  * [20] Schwarz, Dianne S, and Michael D Blower. ”The endoplasmic reticulum: structure, function and response to cellular signalling”, Cellular and molecular life sciences: CMLS vol. 73,1 79-94 (2016).
  * [21] Rao, Madan and R. C., Sarasij, Active Fusion and Fission Processes on a Fluid Membrane, Phys. Rev. Lett., (87), 12, 2001.
  * [22] Nakahara, Mikio. Geometry, topology and physics. CRC press, 2018\.
  * [23] K. Brakke. ”The motion of a surface by its mean curvature”, In: Mathematical Notes, Vol. 20. Princeton University Press, Princeton. (1978)
  * [24] Pachner, Udo, ”P.L. homeomorphic manifolds are equivalent by elementary shellings”, European Journal of Combinatorics, 12 (2): 129–145, (1991).
  * [25] Brakke KA. The surface evolver. Experimental mathematics. 1992 Jan 1;1(2):141-65.
  * [26] West M, Zurek N, Hoenger A, Voeltz GK. ”A 3D analysis of yeast ER structure reveals how ER domains are organized by membrane curvature”, J. Cell Biol. 193:333–46, (2011).
  * [27] Silke Robitzsch, et.al., ”Quantitative analysis of ER morphology in plants: The endoplasmic reticulum in Arabidopsis hypocotyls”, The Plant Cell, Vol. 33, Issue 7, 2264-2281, (2021).
  * [28] Shibata Y, Voeltz GK, Rapoport TA. ”Rough sheets and smooth tubules”. Cell. 126 (3): 435–9, (2006).
  * [29] Shibata Y, Hu J, Kozlov MM, Rapoport TA, ”Mechanisms shaping the membranes of cellular organelles”. Annu Rev Cell Dev Biol 25:329–354 (2009).
  * [30] A W. Helfrich, Helical bilayer structures due to spontaneous torsion of the edges, 1986, The Journal of Chemical Physics, 1085-1087, 85, 2.
  * [31] F. David, in Statistical Mechanics of Membranes and Surfaces, edited by D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 1989).
  * [32] Safran, S.A., et.al., Phys. Rev. Lett. 57 491, (1986).
  * [33] T. Tlusty, S. A. Safran, R. Strey, Topology, Phase Instabilities, and Wetting of Microemulsion Networks, PRL, Vol. 84, No. 6, Feb 2000.
  * [34] Lubensky, T. C. and Renn, S. R., Twist-grain-boundary phases near the nematic˘smectic-A˘smectic-C point in liquid crystals, Phys. Rev. A, Vol 81, Issue 8, 1990, pp: 4392–4401.
  * [35] A. A. Abrikosov, [Sov. Phys.—JETP 5, 1174 (1957)]
  * [36] W. Helfrich, Elastic Properties of Lipid Bilayers—Theory and Possible Experiments, Z. Naturforsch. C 28 (1973) 693-703.
  * [37] Shemesh T, Klemm RW, et.al., A model for the generation and interconversion of ER morphologies. Proceedings of the National Academy of Sciences. 2014 Dec 9;111(49):E5243-51.
  * [38] Lennard-Jones. ”Cohesion”. Proceedings of the Physical Society. 43 (5): 461–482, (1931).
  * [39] Guven, Jemal, Greg Huber, and Dulce María Valencia. ”Terasaki spiral ramps in the rough endoplasmic reticulum.” Physical review letters 113.18 (2014): 188101.
  * [40] Alageshan, Jaya Kumar, Buddhapriya Chakrabarti, and Yashodhan Hatwalne. ”Equilibrium of fluid membranes endowed with orientational order.” Physical Review E 95, no. 4 (2017): 042806.
  * [41] Gelbart, William M., Avinoam Ben-Shaul, and Didier Roux, eds. Micelles, membranes, microemulsions, and monolayers. Springer Science and Business Media, 2012.
  * [42] Lianghui Gao, Leonardo Golubovic, ”Smectic phases of semiflexible manifolds: Constant-pressure ensemble” Physical Review E 66, 051918 (2002).
  * [43] Voeltz, GK; Prinz WA; Shibata Y; Rist JM; Rapoport TA. ”A class of membrane proteins shaping the tubular endoplasmic reticulum”. Cell. 124 (3): 573–586 (2006).
  * [44] Gao, Lianghui and Golubovi, Leonardo, Smectic phases of semiflexible manifolds: Constant-pressure ensemble, Phys. Rev. E, vol. 66, issue 5, pages: 051918, 14, 2002.
  * [45] Parsegian, A. and Evans, E., Proc. Nat. Acad. Sci. USA 83 7132, (1986).
  * [46] David A. Huse, Stanislas Leibler. ”Phase behaviour of an ensemble of nonintersecting random fluid films.” Journal de Physique, 49 (4), pp.605-621, (1988).
  * [47] Golubovic, Leonardo and Lubensky, T. C., Smectic elastic constants of lamellar fluid membrane phases: Crumpling effects, Phys. Rev. B, vol. 39, issue 16, pp: 12110–12133, 1989.
  * [48] Ogawa, H., Uchida, N., Numerical simulation of the twist-grain-boundary phase of chiral liquid crystals. Physical Review E, 73(6), 060701, 2006.
  * [49] Kamien, Randall D., and Tom C. Lubensky. ”Minimal surfaces, screw dislocations, and twist grain boundaries.” Physical review letters 82.14 (1999): 2892.
  * [50] Chaikin, Paul M., Tom C. Lubensky, and Thomas A. Witten. Principles of condensed matter physics. Vol. 10. Cambridge: Cambridge university press, (1995).
  * [51] David A. Huse, Stanislas Leibler. “Phase behaviour of an ensemble of nonintersecting random fluid films.” Journal de Physique, 49 (4), pp.605-621 (1988).
  * [52] Brakke, Ken A., and John M. Sullivan. ”Using symmetry features of the surface evolver to study foams.” In Visualization and mathematics, pp. 95-117. Springer Berlin Heidelberg, 1997.
  * [53] Hestenes, Magnus R., Stiefel, Eduard. ”Methods of Conjugate Gradients for Solving Linear Systems”. Journal of Research of the National Bureau of Standards. 49 (6): 409 (1952).
  * [54] Hatwalne, Yashodhan, and Murugappan Muthukumar. ”Chiral symmetry breaking in crystals of achiral polymers.” Physical review letters 105, no. 10 (2010).
  * [55] Johannes C. C. Nitsche, Lectures on Minimal Surfaces: Volume 1, Introduction, Fundamentals, Geometry and Basic Boundary Value Problems, Cambridge University Press, (1989).
  * [56] Avriel, Mordecai. “Nonlinear Programming: Analysis and Methods.” Dover Publishing (2003).
  * [57] Menon, Gautam I., Rahul Pandit, and Sriram Ramaswamy. ”Sponge Phase Transitions from a Lattice Mode.” Molecular Crystals and Liquid Crystals Science and Technology. Section A. 288.1 (1996): 93-104.