# Patterning of multicomponent elastic shells by Gaussian curvature

Curt Waltmann  Department of Materials Science and Engineering, Northwestern
University, Evanston, IL, 60208  Ahis Shrestha  Center for Computation and
Theory of Soft Materials, Northwestern University, Evanston, IL 60208
Department of Physics and Astronomy, Northwestern University, Evanston, IL,
60208  Monica Olvera de la Cruz  Department of Materials Science and
Engineering, Northwestern University, Evanston, IL, 60208  Department of
Chemistry, Northwestern University, Evanston, IL, 60208  Department of Physics
and Astronomy, Northwestern University, Evanston, IL, 60208  Center for
Computation and Theory of Soft Materials, Northwestern University, Evanston,
IL 60208

###### Abstract

Recent findings suggest that shell protein distribution and morphology of
bacterial microcompartments regulate the chemical fluxes facilitating
reactions which dictate their biological function. We explore how the
morphology and component patterning are coupled through the competition of
mean and Gaussian bending energies in multicomponent elastic shells that form
three component irregular polyhedra. We observe two softer components with
lower bending rigidities allocate on the edges and vertices while the harder
component occupies the faces. When subjected to a non-zero interfacial line
tension, the two softer components further separate and pattern into
subdomains that are mediated by the Gaussian curvature. We find that this
degree of fractionation is maximized when there is a weaker line tension and
when the ratio of bending rigidities between the two softer domains $\approx
2$. Our results reveal a patterning mechanism in multicomponent shells that
can capture the observed morphologies of bacterial microcompartments and,
moreover, can be realized in synthetic vesicles.

##  I Introduction

The hierarchical nature of molecular biology allows for incredible complexity
to arise from simple molecular building blocks. A striking example of this
hierarchy is the division of biological processes into specific subcellular
compartments [1]. Although many of these subcomponents share similar building
blocks such as lipids [2] and proteins [3], they have a vast range of
properties [4, 5, 6, 7, 8, 9, 10, 11], functions [12, 13, 14, 15], and
morphologies [16, 17, 18, 19, 20]. Thus, it is critical to understand how
heterogeneity in the composition [6, 21, 22] of these compartments can lead to
unique morphologies and overall functions. We study a membrane compartment
system composed of multiple components with different mechanical properties
and investigate how these properties lead to complex surface patterning that
affects shell morphology.

One class of subcellular compartments that are comprised of heterogeneous
constituents is the bacterial microcompartment (BMC) [23, 24]. Similar to
viruses, they have semi-permeable shells [25] that are composed with copies of
various shell proteins [26, 27, 28, 29, 30, 31]. Atomistic simulations of
these different shell protein have shown that there are unique mechanical
interactions between the different pairs of shell proteins, namely, the
bending interactions between two species of protein oligomers have distinct
bending energies and energy-minimum bending angles [32, 26]. Thus, at a large-
scale continuum level, the entire shell is expected to behave as a
multicomponent membrane, particularly in terms of component-specific elastic
parameters such as (mean) bending and Gaussian rigidity [33, 34, 35, 36, 37].
Understanding and quantifying the patterning mechanism in such multicomponent
shells could provide key insights into how the different shell proteins are
spatially arranged [38] and their biological functions. Moreover, BMCs house
enzymatic pathways that are involved in many other important cellular
functions including $\text{CO}_{2}$ fixation. Hence, they are of great
interest in the synthetic biology community due to their potential application
as nano-bioreactors. For such applications, it is important to control the
size and shape of these compartments as previous work [26] has shown that
reaction rates can be altered by changing the morphology of BMCs.

Fluid [39, 40, 41] and crystalline [42, 43, 44] membranes have been
successfully studied using continuum elasticity theory that contain relatively
few adjustable parameters [45, 46, 47]. For instance, viral capsids provide a
great illustration of how buckling transitions in single component membranes
are due to the competition of stretching and bending energies. Smaller virus
particles tend to be spherical minimizing the membrane bending energy while
larger particles buckle into icosahedra as the stretching energy scales with
the surface area of the shell [48, 49]. In contrast, multicomponent membranes
allow for a greater variety of shapes including irregular polyhedra [50, 51,
36, 10]. In particular, they form shapes with highly bent vertices and edges
when one component is much less rigid than the others. Unlike icosahedral
viruses studied in previous works [48, 49], BMCs have been observed with
irregular polyhedral shapes [24, 52] possibly due to diverse mechanical
properties that arise from heterogeneous bending interactions between the
shell protein components [32, 26]. The question of how these components are
patterned on the shell surface, which is highly relevant to functional
features such as BMC permeability [53, 25], remains open [38, 31]. In this
work, we address the effects of component dependent rigidity as well as
interfacial line tension on shell patterning and morphology.

The role of Gaussian curvature becomes crucial for multicomponent membranes
bearing domains of distinct rigidities [33, 35, 50, 54]. In a single
component, uniform membrane, the Gaussian curvature is only relevant when
considering a system that undergoes a change in topology. This can occur
during processes such as membrane fusion or fission events, for instance, cell
division or viral budding [55, 56]. This is due to the Gauss-Bonnet theorem,
by which the total Gaussian bending energy of a closed sphere amounts to a
topological constant. However, in multicomponent membranes with components
that have different bending and Gaussian rigidities, the total Gaussian energy
of the system depends on the contact lines of the multidomain boundaries and
the relative values of Gaussian rigidities in the domains [34, 35, 50, 54, 33,
57]. The numerical value of the Gaussian rigidity of a given component has
been shown to be directly proportional to its bending rigidity by the negative
of a dimensionless constant whose reported value is approximately 1 in
experiments and molecular dynamics simulation of lipid membranes [58], as well
as discretized mesh models [42, 9, 50, 36, 49, 59, 60]. As a result, a
component with a larger bending rigidity will have a more negative Gaussian
rigidity than the other components which can create a competition between the
mean and Gaussian curvatures leading to patterning of components on the
membrane. However, in previous works on multicomponent closed membranes, the
mechanical regimes and conditions necessary to address these patterns have not
been accessed [36, 50, 8, 37, 51]. Two component polyhedra [36, 50] do not
allow for multiple edge and vertex components needed for patterning, while
three component icosahedra [37] inherently can not form the complex shapes
dictated by this competition. Here, we investigate how the mean and Gaussian
curvatures are distributed among two components which combine to form the
edges and vertices of a buckled polyhedra whose faces are comprised of a
third, much more rigid component (see Fig. 1). This competition creates
regions of high or low mean and Gaussian curvature which correspond to domains
of the two softer components affecting the shell morphology. We investigate
how this is impacted by the strength of a line tension between the different
components, the ratio of bending rigidities between the components, and the
overall composition of the membrane, finding values which maximize the number
of domains. These low symmetry patterns give insight into how individual
components may be distributed on subcellular compartments with complex and
irregular shapes.

Figure 1: Multicomponent elastic shell with three distinct components denoted
as hard (pink), soft (purple), and softest(cyan), referring to their different
bending rigidities such that $\kappa_{hard}>\kappa_{soft}>\kappa_{softest}$
(top panel). These components pattern the shell based on the surface
distribution of the square of twice the mean curvature $H^{2}$ and the
Gaussian curvature $G$, which are given for the red box section of the shell
with their magnitudes indicated by the color bar (lower right panel). This
leads to polyhedra with inhomogeneous vertex curvature similar to cryo
transmission electron micrographs of BMCs on TEM grids [Reproduced from Ref.
[52] under Creative Commons Attribution 4.0 International (CC BY) license] in
the lower left panel.

This article is arranged as follows. In Sec. II, we outline our continuum
elasticity model and the numerical discretization methods we employ to
minimize the total energy. Sec. III details our analysis of the resulting
energy-minimum shell patterns under varied mechanical parameters. Our key
findings and broader outcomes are highlighted in Sec. IV. Additionally, we
provide Appendices with supplemental materials and results.

##  II Model and methods

###  II.1 Continuum elasticity theory

The elastic energy of thin shells is primarily written as the sum of
mechanical stretching and bending contributions, in which the stretching term
is given as

$E_{Stretch}=\frac{1}{2}\int_{S}\frac{Y}{1+\nu}\left(U(\vec{r})_{ij}^{2}+\frac{\nu}{1+\nu}U(\vec{r})_{kk}^{2}\right)dA\,,$
(1)

where Y is the Young’s modulus, $\nu$ is poisson’s ratio, and
$U(\vec{r})_{ij}$ is the position $\vec{r}$ dependent strain tensor [47, 36,
37]. The bending term is accounted by the Helfrich energy [45, 36, 37] given
by

$E_{Bend}=\frac{1}{2}\int_{S}\kappa(\vec{r})(H(\vec{r})-H_{0}(\vec{r}))^{2}dA+\int_{S}\kappa_{G}(\vec{r})G(\vec{r})dA\,,$
(2)

where $\kappa(\vec{r})$ is the bending rigidity, $H(\vec{r})$ is twice the
mean curvature, $H_{0}(\vec{r})$ is the intrinsic or spontaneous curvature,
$\kappa_{G}(\vec{r})$ is the Gaussian rigidity, and $G(\vec{r})$ is the
Gaussian curvature. The mean and Gaussian curvature are functions of the two
principal radii of curvature, $R_{1}(\vec{r})$ and $R_{2}(\vec{r})$, such that
$H=(\frac{1}{R_{1}}+\frac{1}{R_{2}})$ and $G=\frac{1}{R_{1}R_{2}}$. In the
special case of perfect uniform sphere for instance, the position dependence
in Eqs. (1) and (2) are lost. However, this dependence is maintained for more
general multicompnent membranes with non-uniform shapes. The competition
between the stretching and bending energies is characterized by the
dimensionless Foppl von Karman number, Fvk$=\frac{YR^{2}}{\kappa}$ where $R$
is the shell radius. At high FvK values, the stretching term is dominant and a
closed shell buckles into an icosahedron. Below the critical value
FvK$\approx$ 154 [42, 49], a smoother spherical shell remains with the bending
energy density uniformly distributed.

In this work, we investigate a system of three components (softest, soft, hard
as shown in Fig. 1) with different bending rigidities: $\kappa_{hard}$,
$\kappa_{soft}$, and $\kappa_{softest}$. We define the ratios
$\kappa_{soft}/\kappa_{softest}=\lambda_{1}\geq 1$ and
$\kappa_{hard}/\kappa_{softest}=\lambda_{2}$. We set the ratio
$\kappa_{hard}/\kappa_{soft}=50$ such that $\lambda_{2}=50\lambda_{1}$,
creating the right conditions for polyhedral buckling with the hard component
on the faces [50]. The final, softest component will have a bending rigidity
$\kappa_{softest}\leq\kappa_{soft}$ and thus will buckle as well. Determining
the exact value of the Gaussian rigidity of a given component,
$\kappa_{G,\alpha}$, is not trivial and may dependent on the geometry as well
as other intrinsic structural properties of the membrane system [42, 9, 50,
36, 49, 37, 41, 60, 58]. However, it is generally accepted that the Gaussian
rigidity is proportional to the negative of the bending rigidity, such that
$\kappa_{G,\alpha}=-c\kappa_{\alpha}$ where $c$ is a dimensionless positive
constant. Moreover, we set zero intrinsic curvature $H_{0}=0$ in the membrane
for simplicity.

We now write the mean and Gaussian bending energy as a function of
$\lambda_{1}$ and $\lambda_{2}$ by breaking up the shell into domains of the
different components, $D_{\alpha}$. We reduce Eq. (2) as the dimensionless
bending energy

$E^{*}_{Bend}=E^{*}_{Mean}-cE^{*}_{Gauss}\,,$ (3)

where $E^{*}_{Bend}=E_{Bend}/\kappa_{softest}$, and the dimensionless mean and
Gaussian energies are

$E^{*}_{Mean}=\frac{1}{2}\int_{D_{softest}}{H^{2}dA}+\frac{\lambda_{1}}{2}\int_{D_{soft}}{H^{2}dA}+\frac{\lambda_{2}}{2}\int_{D_{hard}}{H^{2}dA}\,,$
(4)

and

$E^{*}_{Gauss}=\int_{D_{softest}}{GdA}+\lambda_{1}\int_{D_{soft}}{GdA}+\lambda_{2}\int_{D_{hard}}{GdA}\,,$
(5)

respectively. Eq. 4 implies that as $\lambda_{1}$ increases, the highest mean
curvature regions should be occupied by the softest component to minimize the
bending energy in Eq. 3. However, Eq. 5 implies that the more rigid component
should be located on those regions instead. Since, these regions also tend to
be vertices of polyhedral shapes and thus there is a competition partitioning
the mean and Gaussian curvature between the soft and softest components as a
function of $\lambda_{1}$.

Given the presence of multiple, separated components at length scales where
continuum modeling is reasonable, we impose an additional energy penalty when
these components are in contact [61, 36, 37, 56, 57]. To impose a penalty on
the creation of an interface between two different components we add an
interfacial line tension term [36, 37] as follows

$E_{Interface}=\Gamma_{\alpha,\beta}\int_{\partial S}dl\,,$ (6)

where $\Gamma_{\alpha,\beta}$ quantifies this energy penalty between two
components, $\alpha$ and $\beta$. Thus, our total energy is
$E=E_{Bend}+E_{Stretch}+E_{Interface}$.

###  II.2 Discretized Monte Carlo simulations

We discretize the shell surface into a triangular T=192 mesh in a spherical
configuration with components placed on random mesh points according to the
given fractions of each component [36, 37, 50]. It has a radius of
approximately 11.5$r_{0}$ where $r_{0}$ is the equilibrium length between
vertices. We then minimize the discretized version of the total elastic energy
as developed by Seung and Nelson [42] with the components on the mesh points
in order to include the line tension term [36] using a Monte Carlo protocol.
We attempt two different kinds of Monte Carlo moves. One type moves a mesh
point 5% of the equilibrium mesh spacing and the second swaps the identity of
two mesh points. This allows us to minimize the morphology and patterning of
the membrane.

In the Seung-Nelson (SN) discretization, we assign harmonic springs between
each pair of nearest neighbor vertices [36, 37, 50]. When the components are
defined on the vertices, this gives us the SN discretized stretching energy as

$E^{SN}_{Stretch}=\frac{1}{2}\sum_{i,j}\epsilon_{\alpha,\beta}(|\vec{r}_{i}-\vec{r}_{j}|-r_{0})^{2}\,,$
(7)

where $\alpha$ and $\beta$ are the identities of vertices i and j respectively
and when $\alpha\neq\beta$ one must define a mixing rule for the harmonic
spring constant $\epsilon_{\alpha,\beta}$. However, in practice, we set
$\epsilon_{\alpha,\beta}=\epsilon$, as a constant. We also set $r_{0}=1$, the
unit of length. In the continuum limit this potential leads to a Young’s
modulus, $Y=\frac{2}{\sqrt{3}}\epsilon$ and Poisson ratio, $\nu=\frac{1}{3}$
[42].

The SN discretized bending energy is computed in terms of the angle between
triangles which share edge i,j where $k_{\alpha}$ is the discrete bending
rigidity, and $n_{i}$ and $n_{j}$ are the unit normal vectors of triangles
$t_{i}$ and $t_{j}$, for which there is dihedral angle, $\theta_{i,j}$. Here,
the components are defined on the vertices instead [36, 37, 50] and the
expression becomes

$E^{SN}_{Bend}=\frac{1}{2}\sum_{t_{i},t_{j}}k_{\alpha,\beta}|\hat{n_{i}}-\hat{n_{j}}|^{2}=\sum_{t_{i},t_{j}}k_{\alpha,\beta}(1-\cos(\theta_{i,j}))\,,$
(8)

where $k_{\alpha,\beta}$ is the discrete bending rigidity of an edge
connecting vertices with identities $\alpha$ and $\beta$. We choose the mixing
rule to be an average of the two discrete bending rigidities
$k_{\alpha,\beta}=\frac{k_{\alpha}+k_{\beta}}{2}$. We note that such a
parameter is not defined the continuum modeling which only defines the bending
rigidity as a property of only one component. Here and in the single component
system, the relation between the discrete bending constant, k, and the
continuum bending rigidity, $\kappa$, is assumed to be $\kappa=\sigma k$ with
$\sigma=\frac{\sqrt{3}}{2}$ following the works [42, 49, 60]. For this
discretization, we estimate that $c\approx 1.28$ by using it as a fitting
parameter for our data (see Appendix B). This is in good agreement with
previous works [49, 60] which show that c=$\frac{4}{3}$ is a valid solution
for this mesh for both spheres and cylinders. We choose the discrete bending
rigidities such that FvK${}_{hard}\approx 18$, FvK${}_{soft}\approx 900$, and
FvK${}_{softest}\approx 900\lambda_{1}$.

Following previous work of Sknepnek et al. [36], we define the discrete
interfacial line tension as a sum over the neighboring mesh points i,j

$E^{Sknepnek}_{Interface}=\frac{1}{2}\sum_{i,j}\gamma_{\alpha,\beta}(1-\delta_{\alpha,\beta})\,,$
(9)

where $\delta_{\alpha,\beta}$ is the Kronecker delta and
$\gamma_{\alpha,\beta}$ is the discrete line tension parameter with units of
energy. Thus, it only penalizes the number of interfacial contacts and not the
interface length as in the continuum model and it can not lead to domain
budding. The factor of $\frac{1}{2}$ is included to avoid double counting
these contacts. Here, we make $\gamma_{\alpha,\beta}=\gamma$ a constant such
that all interfaces are penalized equally. Throughout this work, we show the
magnitude of $\gamma$ relative to the bending energy using
$\frac{\gamma}{k_{hard}}$ where $k_{hard}$is the discrete bending constant for
the hard component.

##  III Results and discussion

Membrane morphologies are shown for different component fractions of the
softest and soft components, $F_{softest}$ and $F_{soft}$ respectively, with
and without line tension in Fig. 2 along with the corresponding energies.
These configurations are generated after an annealing protocol, which we
explore in the Appendix. We note the mechanical energy, that is the sum of
only the stretching and bending energies, drops quickly as softer components
are added to a one component shell made only of the hard component. This
accompanies a polyhedral buckling transition as shown in previous work [36,
50]. However, once a certain amount of the softer components is added
($F_{softest}+F_{soft}>0.1$) the mechanical energy seems to slowly approach a
minimum. This seems to signify that the buckling transition is ”complete”
(Fig. 2(a)). A related trend is seen in the line tension energies, which
increase as more of the softer components are added in Fig. 2(b). The line
tension energy linearly for smaller $F_{soft}+F_{softest}$ as adding the
softer components creates a line or edge that has an associated energy with
it. However, at the point where this buckling transition is completed, the
relationship between the amount of the softer components and the line tension
energy is no longer linear and instead the line tension increases much more
slowly. At this point, the bending energy is no longer dominating and the line
tension energy can be optimized as well.

Figure 2:  Typical energies and shell configurations with and without line
tension, where the rigidity ratio,
$\frac{\kappa_{soft}}{\kappa_{softest}}=\lambda_{1}=2$. (a) The mechanical
energy, $E_{Bend}+E_{Stretch}$, for elastic shells with different fractions of
the hard (pink), soft (purple), and softest (cyan) components. $F_{Hard}$, the
fraction of the hard component is given by $1-F_{softest}-F_{soft}$.
Characteristic shell patterning is shown for all combinations where
$F_{soft}=0.0,0.1,0.2$ and $F_{softest}=0.0,0.1,0.2$. (b) The line tension
energy, $E_{Interface}$ for an elastic shell with line tension. Characteristic
shell patterning is shown for all combinations where $F_{soft}=0.0,0.1,0.2$
and $F_{softest}=0.0,0.1,0.2$.

When the line tension is present, we observe segregated domains of the softer
components in Fig. 2(b), while without the line tension, we observe patterns
that look like the soft component (purple) is ”reinforcing” the softest
component (cyan) which is located at the center of the edges. In Fig. 3, we
quantify the number of these softer domains
$N_{SofterDomains}=N_{SoftDomains}+N_{SoftestDomains}$, that is, the sum total
of the domains which contains both the soft and the softest parts. More
precisely, a softer domain is defined as a set of mesh points of either the
soft or softest component where, for all vertices in the set, there exists a
path between mesh points that is only connected by neighboring mesh points of
that same component. Fig. 3(a) shows that the number of softer domains is
maximized when the fraction of the soft component, $F_{soft}$ is higher then
that of the of softest component, $F_{softest}$, but when the total amount of
softer components is not too high ($F_{soft}+F_{softest}\leq 0.25$).
Representative morphologies shown in Fig. 3(b) suggest that the soft (purple)
component prefers to be located on the vertices of the polyhedra, while the
softest (cyan) component tends to be located on the edges. This creates the
right conditions for many domains to form since the soft vertices will be
separated from each other by the softest edges. Adding more of the soft
component fuses the vertex domains leading to a decrease in the number of
independent domains. The soft vertex domains also seem thicker than the cyan
edges, leading to inhomogeneous mean curvature of the vertices and edges based
on which component the edge or vertex is made of. Another way to fuse the
domains should be to increase the magnitude of the line tension parameter,
$\frac{\gamma}{k_{hard}}$. This penalizes the formation of independent
domains. Indeed, Fig. 3(c) shows a decrease in the number of domains for
almost every composition by doubling $\frac{\gamma}{k_{hard}}$. The one
composition which seems unaffected by the increased line tension is
$F_{soft}=F_{softest}=0.05$, which is addressed later in the text and in Fig.
5.

Figure 3: Analysis of domain patterning for different compositions when the
rigidity ratio, $\lambda_{1}=2$. (a) The number of softer domains for
different fractions the soft (purple) and softest (cyan) components for the
line tension parameter, $\frac{\gamma}{k_{hard}}$ = $10^{-3}$. A softer domain
is a set of mesh points of either the soft or softest component where, for all
mesh points in the set, there exists a path between mesh points that is only
connected by neighboring mesh points of that same component. Representative
snapshots for the same component fractions and line tension parameter are
shown in (b). (c) The number of softer domains for different fractions the
soft (purple) and softest (cyan) components for the line tension parameter,
$\frac{\gamma}{k_{hard}}$ = 2* $10^{-3}$. Representative snapshots for the
same line tension parameter are shown in (d). Figure 4: Analysis of the
relationship between domain formation and the distribution of the mean and
Gaussian curvatures. All values shown are with line tension paramter
$\frac{\gamma}{k_{hard}}=10^{-3}$ (a) The number of softer domains as a
function of the rigidity ratio,
$\lambda_{1}=\frac{\kappa_{soft}}{\kappa_{softest}}$ for three different
compositions shown in (c). (b) Analysis of the partitioning of the mean and
Gaussian curvatures also as a function of the rigidity ratio, $\lambda_{1}$.
Dashed lines are the fraction of the Gaussian curvature found on the soft
component, $\rho^{G}_{soft}$, as defined in eq. 10. Solid lines are the
fraction of the mean squared curvature found on the soft component,
$\rho^{M}_{soft}$, as defined in eq. 11. There is more Gaussian and mean
squared curvature on the softest component (cyan), however, there is a higher
fraction of Gaussian curvature on the soft component than mean curvature on
the soft component. (c) Snapshots of different compositions and different
$\lambda_{1}$, which correspond to (a) and (b). (d) The ratio of Gaussian and
mean curvature energy on the softest and softest components, $\beta^{G}$ and
$\beta^{M}$, as defined in eq 12 and 13, respectively.  Figure 5: Analysis of
the partitioning of the mean and Gaussian curvatures also as a function of the
rigidity ratio, $\lambda_{1}=\frac{\kappa_{soft}}{\kappa_{softest}}$. The
different colors are for different line tension parameters,
$\frac{\gamma}{k_{hard}}$. Dashed lines are the fraction of the Gaussian
curvature found on the soft component (purple), $\rho^{G}_{soft}$, as defined
in eq. 10. Solid lines are the fraction of the mean curvature found on the
soft component, $\rho^{M}_{soft}$, as defined in eq. 11. (a) Relative to the
mean curvature, the Gaussian curvature is partitioned preferentially to the
soft component ($\rho^{G}_{soft}>\rho^{M}_{soft}$) for
$F_{soft}=F_{softest}=0.1$. This only occurs in the presence of a line
tension, but doubling the line tension parameter has little effect on this
partitioning but does lead to domain fusion. For $F_{soft}=F_{Softest}=0.05$,
(c) when no line tension, there is only a slight preference for the Gaussian
curvature to be on soft domains ($\rho^{G}_{soft}>\rho^{M}_{soft}$) at low
values of $\lambda_{1}$, with mean squared curvature dominating for other
$\lambda_{1}$ ($\rho^{G}_{soft}<\rho^{M}_{soft}$). (d) When weak line tension
is added there is a preference for the Gaussian curvature to partition to the
soft component, but only when $\lambda_{1}>3$. (e) When the line tension is
doubled there is always a preference for the Gaussian curvature to be located
on the soft component, although the preference is maximized at
$\lambda_{1}=2$. The ratio of mean and Gaussian curvature energies between the
softest and soft components, $\beta^{G}$ and $\beta^{M}$, as defined in eq. 12
and 13 respectively for different values of line tension when (b)
$F_{soft}=F_{softest}=0.1$ and (f) $F_{soft}=F_{softest}=0.05$.

Snapshots in Fig. 3(d) visually show the fusion of soft and softest domains.
Although the preference for the soft component to be located on the vertices
is still noticeable, there are now more soft edges and these edges are not as
clearly defined as their softest (cyan) counterparts. We expect that the
persistent preference of the soft component to be located on the vertices is
due to the Gaussian curvature of the membrane, which is concentrated at the
vertices of the polyhedra. The relation between the Gaussian and bending
rigidity imples that the soft component has a more negative Gaussian rigidity
despite having a more positive bending rigidity. Thus, the soft component
should be preferred in regions of high Gaussian curvature such as the vertices
of the polyhedra, as we observe. Based on our definition of the Gaussian
energy in eq. 5, this relationship should also be a function of $\lambda_{1}$,
the ratio between the rigidity of the soft and softest components. To
illustrate this, we vary the rigidity ratio in the range $1\leq\lambda_{1}\leq
4$, tracking the number of domains for different compositions. We also
quantity the Gaussian and mean contributions located on the soft component
with respect to the softest component through the following dimensionless
quantities

$\rho^{G}_{soft}=\frac{\int_{D_{soft}}GdA}{\int_{D_{Soft}}GdA+\int_{D_{softest}}GdA}$
(10)

and

$\rho^{M}_{soft}=\frac{\int_{D_{soft}}H^{2}dA}{\int_{D_{Soft}}H^{2}dA+\int_{D_{softest}}H^{2}dA}$
(11)

These is computed using the per node definitions of the Gaussian and mean
curvatures found in eq. 16 and 17. In Fig. 4(a), we observe a maximum in the
number of softer domains for rigidity ratios of 2 or 3, depending on the
specific compositions. In Fig. 4(b), we show that both $\rho^{G}_{soft}$ and
$\rho^{M}_{soft}$ are always less than .5, meaning that more of both the mean
and Gaussian contributions are on the softest component. However, there is
always a higher Gaussian contribution than mean contribution on the soft
component ($\rho^{G}_{soft}>\rho^{M}_{soft}$). This creates a maximum in the
number of domains when the rigidity ratio is large enough for there to be a
difference in Gaussian rigidity of the two components but not so large that
the mean curvature term dominates completely. This is illustrated by the
snapshots at different $\lambda_{1}$ values in Fig. 4(c). As $\lambda_{1}$
increases the softest component creates edges that are skinnier relative to
the soft component allowing the softest component to cover more of the
surface, increasing its fraction of the Gaussian and mean curvatures as shown
in Fig. 4(b). At first, this leads to more segregated domains as it mostly
occupies the edges, leaving segregated vertices occupied by the soft
component. However, it eventually takes up more of the vertices as well,
leading to a decrease in the average number of domains.

One can approximate the relationship between $\rho^{G}_{soft}$ and
$\lambda_{1}$ goes as $1/(1+\beta^{G}\lambda_{1})$, where $\beta^{G}$ can be a
fitting parameter that describes how fast $\rho^{G}_{soft}$ decays. It can be
shown through Eq. (12) that in this approxtimation, we find

$\beta^{G}=\frac{\int_{D_{softest}}GdA}{\lambda_{1}\int_{D_{Soft}}GdA}\,,$
(12)

which means that here $\beta_{G}$ is the ratio of Gaussian energy in the
softest component to that of the soft component. Similary, for the mean
curvature using Eq. (13), we get

$\beta^{M}=\frac{\int_{D_{softest}}H^{2}dA}{\lambda_{1}\int_{D_{Soft}}H^{2}dA}\,,$
(13)

Interestingly, $\beta^{G}\approx 1$ for all values of $\lambda_{1}$ and for
different compositions as shown in Fig. 4(d), which means the soft and softest
components have the same Gaussian curvature energy. This is not the case for
the mean curvature energy, where $\beta^{M}$ has a stronger dependence on
$\lambda_{1}$ and seems to plateau around 1.6 for all compositions. This means
that the softest component is taking on more of the mean curvature energy. The
difference in $\beta^{G}$ and $\beta^{M}$ quantifies the energies involved in
the morphological ”splitting” of the mean and Gaussian curvatures through the
low mean curvature vertices of the soft component.

Finally, we look at the effect of line tension on the splitting of the mean
and guassian curvatures. In Fig. 5(a), we plot $\rho^{G}_{soft}$ and
$\rho^{M}_{soft}$ as a function of the rigidity ratio at
$\frac{\gamma}{k_{hard}}=0,10^{-3},2*10^{-3}$ and $F_{soft}=F_{softest}=0.1$.
We observe only a slight difference in the curves for
$\frac{\gamma}{k_{hard}}=10^{-3}$ and $2*10^{-3}$ with a slight increase in
the fraction of both mean squared and Gaussian curvature on the soft component
at higher values of $\lambda_{1}$. This suggests that the magnitude of the
line tension plays only small role in determining how the mean and Gaussian
curvatures are distributed between the two components, despite the fact that
the lower value of line tension lead to significantly more domains in Fig. 3.
The values of $\beta^{G}$ and $\beta^{M}$ also remain nearly unchanged and can
be seen in Figure 5(b). Thus, this effect was only due to domain fusion and
not a change in curvature distribution or energy. For shells without line
tension, patterning on the vertices is lost as there is more mean squared
curvature than Gaussian curvature on the soft component and both $\beta^{G}$
and $\beta^{M}$ become much larger with $\beta^{G}>\beta^{M}$ (see Fig. 5(f)).
This is the case when there is plenty of excess softer component to buckle the
shell as shown in Fig. 2. In contrast, when all of the softer components are
needed to buckle the shell, as in Fig. 5(b-d) where
$F_{soft}=F_{softest}=0.05$, the line tension has a large impact on the
partitioning of mean and Gaussian curvature. Even without line tension in Fig.
5(c), there is a slight preference for the Gaussian curvature on the soft
phase when $1<\lambda_{1}<2$. When the weak line tension is applied, this
preference is much stronger and extends to $\lambda_{1}>3$ (Fig. 5(d)). When
the line tension is doubled there is always a preference for the Gaussian
curvature to be located on the soft component, although the preference is
maximized at $\lambda_{1}=2$ and decreases afterwards as shown in Fig. 5(e).
Here the line tension and rigidity ratio seem to be competing against each
other, with line tension favoring the formation of soft domains located on the
vertices.

##  IV Conclusion and Outlook

In this article, we show that Gaussian curvature plays a crucial role in the
morphology and patterning of multicomponent elastic membranes. We observe, for
shells composed of three components with distinct bending rigidities, energy-
minimum configurations where the two softer components distributed on the
edges and vertices while the harder component forms the faces. Among the two
softer components that occupy the edges and vertices, a preference of the
component with higher bending rigidity is predominantly found at the vertices
leading to many subdomains, although with low symmetry. These vertices were
much flatter than their counterparts occupied by the component with lower
bending rigidity leading to polyhedra with inhomogeneous vertices as has been
seen in BMCs with many components. These effects required that the two softer
components on the edges and vertices have different bending rigidities but not
so different that one component will occupy the vast majority of the bent
edges and vertices.

We saw a peak in the number of domains at a rigidity ratio of $\approx 2$, in
contrast with the larger ratios necessary for the polyhedral buckling
transition. Since this ratio is relatively low, the necessary differences in
bending properties may result from only subtle changes to chemical properties
of the molecules that make the different shell proteins. This is consistent
with the BMC system where the two most common types of shell proteins found
are structurally similar and have almost identical amino acid sequences [29].
These effects also required a line tension to help segregate the components,
but not one so strong as to fuse the domains back together, again suggesting
relatively similar components can create these morphological patterns.

Besides such distinctions in their elastic properties, the shell protein
components of BMCs are also known to have distinct electrostatic properties,
namely, they carry different ratios of charged amino acids. Thus, the
resulting morphology and component patterning observed here may in turn lead
to patterned surface charge densities on the shell as seen in the icosahedral
shells that have been crystallized. Moreover, such BMC systems can facilitate
chemical reaction that produce fluxes of ionic chemical species. The combined
effect of such surface features provided with their asymmetric distribution on
the shell can induce self-phoretic movement [62, 63]. Both the elastic and
electrostatic properties arise due to small compositional variations in the
protein components, and thus, it may also be possible to regulate these subtle
chemical perturbations in synthetic vesicles that can be optimized for
specialized functionalities.

###### Acknowledgements.

CW was funded by the Northwestern University Nicholson Fellowship. AS was
supported by the Department of Energy (DOE), Office of Basic Energy Sciences
under Contract DE-FG02-08ER46539. MOdlC thanks the computational support of
the Sherman Fairchild Foundation. CW, AS and MOdlC thank the support of the
Center for Computation and Theory of Soft Materials.

##  Appendix A Itzykson Discretization

In the main text, we use the SN dicretization [42] of the Helfrich bending and
stretching energy, but define the components on the vertices instead of the
edges [36]. This makes it simple to define the line tension, however, the
edges can now have a mixed identity at the interface and one must define
mixing rules, which are somewhat arbitrary and do not exist in the continuum
equations. One could also use the Itzykson discretization [43, 37]. Here, the
components are located on the vertices, however the bending energies are
computed by calculating the mean and Gaussian curvatures at the vertex instead
of the angle across triangles. In Fig. 8, we show that this discretization is
does not perform well for our annealing protocol and thus we do not show
results using the Itzykson discretization in the main text. However, we do use
the definitions of mean and Gaussian curvature at a vertex to analyze our
morphologies generated by the SN discretization.

The issue of mixed edges is resolved in the Itzykson [43] discretization by
computing the local mean and Gaussian curvature of the Voronoi cell around
each vertex. The potential, summed over all vertices, takes on the form

$E_{Bend}=\sum_{i}\left(\frac{\kappa_{\alpha}}{2}(H_{i}-H_{0,\alpha})^{2}+\kappa_{G,\alpha}G_{i}\right)A_{i}\,,$
(14)

where $\kappa_{\alpha}$ and $\kappa_{G,\alpha}$ have the same value as the
continuum parameters and thus we use the same $\kappa$ to denote them as
opposed to the SN discretization, which uses k that must be converted to
$\kappa$. $H_{i}$ and $G_{i}$ are the mean and Gaussian curvature at vertex i
and $A_{i}$ is the dual-lattice vertex area defined as

$A_{i}=\frac{1}{8}\sum_{j\in
N(i)}[\cot(\phi_{i,j})+\cot(\psi_{i,j})](r_{i}-r_{j})^{2}\,,$ (15)

where the sum is over the neighbors of i, indexed by j. The angles
$\phi_{i,j}$ and $\psi_{i,j}$ are the interior angles opposite from edge i,j
which is shared by triangles, $t_{i}$ and $t_{j}$. The mean curvature is
computed as

$H_{i}=\frac{sign}{2A_{i}}\sum_{j\in
N(i)}[\cot(\phi_{i,j})+\cot(\psi_{i,j})](r_{i}-r_{j})\,,$ (16)

where $sign$ denotes whether the curvature is concave or convex, 1 or -1
respectively. The Gaussian curvature is found by summing over the angles
adjacent to vertex i in the triangles which share vertex i as follows

$G_{i}=\frac{1}{A_{i}}\left(2\pi-\sum_{t}\theta_{t}\right)\,,$ (17)

where $\theta$ is the angle at vertex i. In Figure 8, we show this discretized
form of the energy is not adequately minimized on our mesh by our Monte Carlo
protocol. However, we use this discretization to calculate the mean and
Gaussian curvature of the mesh minimized with the SN discretization in Figures
4 and 5 in the main text.

##  Appendix B Estimation of $c$ and $\sigma$

Figure 6: Using a fitting scheme where the continuum solution is $c=1.28$ as
described in the main text we estimate c for different compositions and the
values of $\frac{\gamma}{k_{hard}}$ and $\lambda_{1}$ shown in (a), (b), (c).
We also display the average value of c for the different compositions at the
bottom of (a), (b), and (c) which also show reasonable agreement with the
continuum solution.

We attempted to estimate the continuum $c$ and $\sigma$ for our three
component meshes in a way that is consistent with the continuum result for the
one component sphere. In the SN discretization of a one component sphere, the
discrete bending energy is $4\pi k\frac{\sqrt{3}}{3}$. We can set this equal
to the corresponding continuum bending energy

$4\pi k\frac{\sqrt{3}}{3}=\frac{\kappa}{2}\int H^{2}dA-c\kappa\int
GdA=4\pi\kappa(2-c)\,,$ (18)

and thus

$\kappa=\frac{\sqrt{3}}{3(2-c)}k=\sigma k\,.$ (19)

We will use this relationship to reduce c and $\kappa$ to a single fitting
parameter. Lidmar et al. [49] previously used the $\kappa=\frac{\sqrt{3}k}{2}$
for a cylinder with no Gaussian curvature to arrive at the result that
$c=\frac{4}{3}$; a solution which is correct for both geometries. We use

$E^{discrete}_{Bend}=\kappa_{softest}E^{*}_{Bend}\,,$ (20)

where $E^{discrete}_{bend}$ is the discretized bending energy of the mesh and
$E^{*}_{Benc}$ is the same as defined in eq. 3. We found that the best fit for
the data, including different compositions, line tensions, and rigidity ratios
was $c\approx 1.28$ and $\sigma\approx 0.80$. This is in fairly good agreement
with the values suggested by Lidmar et al. In Figure 6, we show c for some
specific parameters.

##  Appendix C Monte Carlo annealing

Figure 7:  Results for annealing using (a) Seung-Nelson discretization and (b)
Itzykson discretization, both with a lower maximum temperature, $T^{*}=0.01$.
Typical energies and shell configurations with and without line tension, where
the rigidity ratio, $\frac{\kappa_{soft}}{\kappa_{softest}}=\lambda_{1}=2$ and
$\frac{\gamma}{k_{hard}}=10^{-3}$. The line tension energies and
characteristic shell patterns are shown for all combinations where
$F_{soft}=0.0,0.1,0.2$ and $F_{softest}=0.0,0.1,0.2$. Figure 8: One
representative shell morphology for annealing using the Itzykson
discretization with the standard maximum temperature, $T^{*}=0.1$. The
rigidity ratio, $\frac{\kappa_{soft}}{\kappa_{softest}}=\lambda_{1}=2$,
$\frac{\gamma}{k_{hard}}=10^{-3}$, and $(F_{softest},F_{soft})=(0.2,0.1)$.
Instabilities at vertices composed of the softest (cyan) component in (a)-(c)
with (a) showing some preference of the soft component to occupy the vertices.
(d) Shows a magnified image where it is apparent that some of the hard
component is trapped and there are many convex bending angles.

Our annealing protocol involves 7 different temperature cycles. The first 5
cycles start at an annealing temperature, $T^{*}$, of $0.1$ and then is
reduced linearly in 10 steps to $T^{*}=10^{-7}$. The 6th and 7th cycles start
at $T^{*}=0.01$ and $T^{*}=0.001$ respectively before being cooled to
$T^{*}=10^{-7}$ in the same linear fashion over 10 steps. For each $T^{*}$, we
run 5000 Monte Carlo steps for a total of 350000 Monte Carlo steps. In each
Monte Carlo step, we attempt two different kinds of Monte Carlo moves. (1) For
each vertex, we attempt to move it a distance of $0.05r_{0}$ and (2) For each
vertex, we attempt to swap the identity of two random vertices. The
probability of move acceptance, $P$, is given by

$P=min\left(1,e^{-\Delta E/T^{*}}\right)\,,$ (21)

where $\Delta E$ is the difference in total energy of the system after the
attempted Monte Carlo move. This protocol is run 4 separate times with
different random initial conditions for all sets of parameters reported in the
main text. All values reported are an average of these 4 trials. Using a
maximum $T^{*}=0.1$ is important as if a similar protocol is used but the
maximum $T^{*}=0.01$, the shells end up in much higher energy states as shown
in Fig. 7(a). This is for the same parameters as in Fig. 2(b) in the main text
and one can see the much higher line tension energies in these cases as well
as a lack of smooth trends in Fig. 7(a).

We also ran this annealing protocol using the Itzykson discretization and
achieved very similar results in Fig. 7(b). Neither Fig. 7(b) nor Fig. 7(a)
show lowest energy states do the lower maximum annealing temperature, $T^{*}$.
When we attempted to use the higher annealing temperature, $T^{*}=1$, with the
Itzykson discretization we found that the mesh became unstable as shown in
Fig. 8. This issue appears at vertices formed by the softest component, making
it impossible to quantitatively compare the Itzykson and Nelson
discretizations for our parameters. We do note that the soft component still
seems to prefer the vertices as was the case with the Nelson discretization in
the main text. This is noticeable in Fig. 8(a).

##  Appendix D Larger line tension

Figure 9:  Representative snapshots for the shells with a strong line tension
parameter. The rigidity ratio, $\lambda_{1}=2$. There are always 2 softer
domains which are now completely separated from the hard component and the
polyhedral shape is lost. Similar effects have been seen previously in a two
component system with high line tension [36].

In the main text, we show that the line tension is necessary for patterning
the polyhedron, but also that doubling the line tension will fuse domains back
together. In Fig.9, we show that as the line tension is increased by an order
of magnitude, the polyhedral morphology is lost altogether and a sphere with
segregated regions of the softer component are formed instead. A very similar
transition has been observed in two component membranes [36]. This also
supports the idea that line tension in these systems must be relatively weak
or else the polyhedral morphology will not be observed.

## References

  * Holthuis and Ungermann [2013] J. C. Holthuis and C. Ungermann, Biological Chemistry 394, 151 (2013).
  * van Meer  _et al._ [2008] G. van Meer, D. R. Voelker, and G. W. Feigenson, Nature Reviews Molecular Cell Biology 9, 112 (2008).
  * Schoch  _et al._ [2004] G. A. Schoch, J. K. Yano, M. R. Wester, K. J. Griffin, C. D. Stout, and E. F. Johnson, Journal of Biological Chemistry 279, 9497 (2004).
  * Funkhouser  _et al._ [2013] C. M. Funkhouser, R. Sknepnek, T. Shimi, A. E. Goldman, R. D. Goldman, and M. O. de la Cruz, Proceedings of the National Academy of Sciences 110, 3248 (2013).
  * Li  _et al._ [2018] S. Li, P. Roy, A. Travesset, and R. Zandi, Proc. Natl. Acad. Sci. 115, 10971 (2018).
  * Schmit  _et al._ [2021] J. D. Schmit, M. Feric, and M. Dundr, Trends Biochem. Sci. 46, 525 (2021).
  * Dahal and Olvera de la Cruz [2019] Y. R. Dahal and M. Olvera de la Cruz, Soft Matter 15, 4311 (2019).
  * Hu  _et al._ [2011] J. Hu, T. Weikl, and R. Lipowsky, Soft Matter 7, 6092 (2011).
  * Vernizzi and Olvera De La Cruz [2007] G. Vernizzi and M. Olvera De La Cruz, Proceedings of the National Academy of Sciences 104, 18382–18386 (2007).
  * Yuan and Olvera de la Cruz [2019] H. Yuan and M. Olvera de la Cruz, Phys. Rev. E 100, 012610 (2019).
  * Brunk  _et al._ [2020] N. E. Brunk, J. Kadupitiya, and V. Jadhao, Phys. Rev. Lett. 125, 248001 (2020).
  * Swift  _et al._ [2013] J. Swift, I. L. Ivanovska, A. Buxboim, T. Harada, P. D. P. Dingal, J. Pinter, J. D. Pajerowski, K. R. Spinler, J.-W. Shin, M. Tewari, _et al._ , Science 341, 1240104 (2013).
  * Crowley  _et al._ [2010] C. S. Crowley, D. Cascio, M. R. Sawaya, J. S. Kopstein, T. A. Bobik, and T. O. Yeates, Journal of Biological Chemistry 285, 37838 (2010).
  * Yildiz  _et al._ [2011] I. Yildiz, S. Shukla, and N. F. Steinmetz, Current Opinion in Biotechnology 22, 901 (2011).
  * Oh  _et al._ [2014] D. Oh, J. Qi, B. Han, G. Zhang, T. J. Carney, J. Ohmura, Y. Zhang, Y. Shao-Horn, and A. M. Belcher, Nano Letters 14, 4837–4845 (2014).
  * Zandi  _et al._ [2020] R. Zandi, B. Dragnea, A. Travesset, and R. Podgornik, Physics Reports 847, 1 (2020), on virus growth and form.
  * Greenfield  _et al._ [2009] M. A. Greenfield, L. C. Palmer, G. Vernizzi, M. O. d. l. Cruz, and S. I. Stupp, Journal of the American Chemical Society 131, 12030 (2009).
  * Leung  _et al._ [2012] C.-Y. Leung, L. C. Palmer, B. F. Qiao, S. Kewalramani, R. Sknepnek, C. J. Newcomb, M. A. Greenfield, G. Vernizzi, S. I. Stupp, M. J. Bedzyk, and M. Olvera de la Cruz, ACS Nano 6, 10901 (2012).
  * Sknepnek  _et al._ [2011] R. Sknepnek, G. Vernizzi, and M. O. de la Cruz, Phys. Rev. Lett. 106, 215504 (2011).
  * Sun  _et al._ [2023] F. Sun, N. E. Brunk, and V. Jadhao, Phys. Rev. E 107, 014502 (2023).
  * Girard [2022] M. Girard, Phys. Biol. 20, 016006 (2022).
  * Jacobs and Frenkel [2017] W. M. Jacobs and D. Frenkel, Biophys. J. 112, 683 (2017).
  * Kennedy  _et al._ [2021a] N. W. Kennedy, C. E. Mills, T. M. Nichols, C. H. Abrahamson, and D. Tullman-Ercek, Curr. Opin. Microbiol. 63, 36 (2021a).
  * Kerfeld  _et al._ [2018] C. A. Kerfeld, C. Aussignargues, J. Zarzycki, F. Cai, and M. Sutter, Nat. Rev. Microbiol. 16, 277–290 (2018).
  * Slininger Lee  _et al._ [2017] M. F. Slininger Lee, C. M. Jakobson, and D. Tullman-Ercek, ACS Synth. Biol. 6, 1880—1891 (2017).
  * Mills  _et al._ [2022] C. E. Mills, C. Waltmann, A. G. Archer, N. W. Kennedy, C. H. Abrahamson, A. D. Jackson, E. W. Roth, S. Shirman, M. C. Jewett, N. M. Mangan, M. Olvera de la Cruz, and D. Tullman-Ercek, Nat. Commun. 13, 3746 (2022).
  * Liu  _et al._ [2021] L.-N. Liu, M. Yang, Y. Sun, and J. Yang, Curr. Opin. Microbiol. 63, 133 (2021).
  * Kennedy  _et al._ [2022] N. W. Kennedy, C. E. Mills, C. H. Abrahamson, A. G. Archer, S. Shirman, M. C. Jewett, N. M. Mangan, and D. Tullman-Ercek, J. Bacteriol. , e00576 (2022).
  * Kennedy  _et al._ [2021b] N. W. Kennedy, S. P. Ikonomova, M. Slininger Lee, H. W. Raeder, and D. Tullman-Ercek, J. Mol. Biol. 433, 166721 (2021b).
  * Sutter  _et al._ [2017] M. Sutter, B. Greber, C. Aussignargues, and C. A. Kerfeld, Science 356, 1293 (2017).
  * Lassila  _et al._ [2014] J. K. Lassila, S. L. Bernstein, J. N. Kinney, S. D. Axen, and C. A. Kerfeld, J. Mol. Biol. 426, 2217 (2014).
  * Li  _et al._ [2021a] Y. Li, N. W. Kennedy, S. Li, C. E. Mills, D. Tullman-Ercek, and M. Olvera de la Cruz, ACS Cent. Sci. 7, 658 (2021a).
  * Baumgart  _et al._ [2003] T. Baumgart, S. T. Hess, and W. W. Webb, Nature 425, 821 (2003).
  * Jülicher and Lipowsky [1993] F. Jülicher and R. Lipowsky, Phys. Rev. Lett. 70, 2964 (1993).
  * Jülicher and Lipowsky [1996] F. Jülicher and R. Lipowsky, Phys. Rev. E 53, 2670 (1996).
  * Sknepnek  _et al._ [2012] R. Sknepnek, G. Vernizzi, and M. Olvera de la Cruz, Soft Matter 8, 636 (2012).
  * Li  _et al._ [2021b] S. Li, D. A. Matoz-Fernandez, and M. Olvera de la Cruz, ACS Nano 15, 14804 (2021b).
  * Tre [2022] Structure (2022).
  * Lipowsky [1991] R. Lipowsky, Nature 349, 475 (1991).
  * Seifert [1997] U. Seifert, Advances in Physics 46, 13 (1997).
  * Deserno [2015] M. Deserno, Chemistry and Physics of Lipids 185, 11 (2015).
  * Seung and Nelson [1988] H. S. Seung and D. R. Nelson, Physical Review A 38, 1005 (1988).
  * Kohyama  _et al._ [2003] T. Kohyama, D. M. Kroll, and G. Gompper, Phys. Rev. E 68, 061905 (2003).
  * J. Bowick and Travesset [2001] M. J. Bowick and A. Travesset, Physics Reports 344, 255 (2001), renormalization group theory in the new millennium.
  * Helfrich [1973] W. Helfrich, Zeitschrift für Naturforschung C 28, 693–703 (1973).
  * Winterhalter and Helfrich [1988] M. Winterhalter and W. Helfrich, The Journal of Physical Chemistry 92, 6865–6867 (1988).
  * Landau and Lifshitz [2009] L. D. Landau and E. M. Lifshitz,  _Theory of elasticity_ , 3rd ed., Course of theoretical physics (Elsevier, Butterworth-Heinemann, Amsterdam Heidelberg, 2009).
  * Zandi  _et al._ [2004] R. Zandi, D. Reguera, R. F. Bruinsma, W. M. Gelbart, and J. Rudnick, Proc. Natl. Acad. Sci. U.S.A. . 101, 15556 (2004).
  * Lidmar  _et al._ [2003] J. Lidmar, L. Mirny, and D. R. Nelson, Physical Review E 68, 051910 (2003).
  * Vernizzi  _et al._ [2011] G. Vernizzi, R. Sknepnek, and M. O. de la Cruz, Proc. Natl. Acad. Sci. U. S. A. 108, 4292 (2011).
  * Demers  _et al._ [2012] M. F. Demers, R. Sknepnek, and M. Olvera de la Cruz, Phys. Rev. E 86, 021504 (2012).
  * Kennedy  _et al._ [2020] N. W. Kennedy, J. M. Hershewe, T. M. Nichols, E. W. Roth, C. D. Wilke, C. E. Mills, M. C. Jewett, and D. Tullman-Ercek, PLoS One 15, 1 (2020).
  * Stewart  _et al._ [2020] K. L. Stewart, A. M. Stewart, and T. A. Bobik, EcoSal Plus 9 (2020).
  * Bassereau  _et al._ [2018] P. Bassereau, R. Jin, T. Baumgart, M. Deserno, R. Dimova, V. A. Frolov, P. V. Bashkirov, H. Grubmüller, R. Jahn, H. J. Risselada, L. Johannes, M. M. Kozlov, R. Lipowsky, T. J. Pucadyil, W. F. Zeno, J. C. Stachowiak, D. Stamou, A. Breuer, L. Lauritsen, C. Simon, C. Sykes, G. A. Voth, and T. R. Weikl, Journal of Physics D: Applied Physics 51, 343001 (2018).
  * Ruiz-Herrero and Hagan [2015] T. Ruiz-Herrero and M. Hagan, Biophysical Journal 108, 585 (2015).
  * Dharmavaram  _et al._ [2019] S. Dharmavaram, S. B. She, G. Lázaro, M. F. Hagan, and R. Bruinsma, PLOS Computational Biology 15, 1 (2019).
  * Shrestha  _et al._ [2021] A. Shrestha, F. Pinaud, and C. A. Haselwandter, Phys. Rev. E 104, L022401 (2021).
  * Hu  _et al._ [2012] M. Hu, J. J. Briguglio, and M. Deserno, Biophysical Journal 102, 1403 (2012).
  * Gompper and Kroll [1996] G. Gompper and D. M. Kroll, Journal de Physique I 6, 1305–1320 (1996).
  * Schmidt and Fraternali [2012] B. Schmidt and F. Fraternali, Journal of the Mechanics and Physics of Solids 60, 172 (2012).
  * Lipowsky [1992] R. Lipowsky, J. Phys. II France 2, 1825 (1992).
  * Steinkühler  _et al._ [2022] J. Steinkühler, C. H. Abrahamson, J. Agudo-Canalejo, R. Golestanian, D. Tullman-Ercek, and N. P. Kamat 10.1101/2022.05.16.492142 (2022).
  * Shrestha and de la Cruz [2023] A. Shrestha and M. O. de la Cruz,  (2023), arXiv:2305.01102 [cond-mat, physics:physics].