# Mean field elastic moduli of a three-dimensional cell-based vertex model

Kyungeun Kim  Physics Department, Syracuse University, Syracuse, NY 13244 USA
Tao Zhang  Department of Polymer Science and Engineering, Shanghai Jiao Tong
University, Shanghai, 200240 China  J. M. Schwarz  Physics Department,
Syracuse University, Syracuse, NY 13244 USA  Indian Creek Farm, Ithaca, NY
14850 USA

###### Abstract

The mechanics of a foam typically depends on the bubble geometry, topology,
and the material at hand, be it metallic or polymeric, for example. While the
foam energy functional for each bubble is typically minimization of surface
area for a given volume, biology provides us with a wealth of additional
energy functionals, should one consider biological cells as a foam-like
material. Here, we focus on a mean field approach to obtain the elastic
moduli, within linear response, for an ordered, three-dimensional vertex model
using the space-filling shape of a truncated octahedron and whose energy
functional is characterized by a restoring surface area spring and a restoring
volume spring. The tuning of the three-dimensional shape index exhibits a
rigidity transition via a compatible-incompatible transition. Specifically,
for smaller shape indices, both the target surface area and volume cannot be
achieved, while beyond some critical value of the three-dimensional shape
index, they can be, resulting in a zero-energy state. As the elastic moduli
depend on curvatures of the energy when the system, we obtain these as well.
In addition to analytically determining the location of the transition in mean
field, we find that the rigidity transition and the elastic moduli depend on
the parameterization of the cell shape with this effect being more pronounced
in three dimensions given the array of shapes that a polyhedron can take on
(as compared to a polygon). We also uncover nontrivial dependence on the
deformation protocol in which some deformations result in affine motion of the
vertices, while others result in nonaffine motion. Such dependencies on the
shape parameterization and deformation protocol give rise to a nontrivial
shape landscape and, therefore, nontrivial mechanical response even in the
absence of topology changes.

##  I Introduction

Different types of objects fill three-dimensional space. Such objects include
foams, metallic grains, and biological cells. There are indeed constraints on
the shapes of these objects to be space-filling. Example shapes include the
cube, the dodecahedron, and the truncated octahedron, all of which are convex
shapes. See Fig. 1 for a three-dimensional example as well as a two-
dimensional one. Efforts to put further constraints on the shapes that, for
example, have the least possible surface area for a given volume, otherwise
known as Kelvin cells, also exist [1]. In 1994, Weaire and Phelan discovered
an space-filling arrangement with two types of cells whose surface area was
0.3% smaller than the originally-proposed polyhedron by Kelvin [2, 1].

As all objects are deformable to greater or lesser extents, it, therefore,
makes sense to explore the shape-changing capability of such objects and ask
how does this capability affect the mechanics of the material? More
specifically, how does the shear modulus differ between a tessellation of
truncated octahedra and a tessellation of cubes? We will address this question
for a more recent energy functional, as compared to foams, that is inspired by
biological cells. More specifically, this energy functional is quadratic in
both surface area and volume, each with a respective target surface area and
target volume [3, 4, 5]. The area ”spring”, if you will, captures the
elasticity of the cell cortex lying just beneath the cell membrane that is, in
part, responsible to the structural integrity of a cell. The volume ”spring”,
on the other hand captures the volume compressibility of a cell given that it
contains water, proteins, fats, etc. For foams, on the other hand, the energy
functional scales linearly with the surface area [1, 6].

|
---|---

Figure 1: Space-filling packings. Left: A packing of hexagons in two
dimensions. Right: A three-dimensional packing built out of truncated
octahedrons.

Some of the earliest work considering a space-filling collection of cells as a
foam-like material in two dimensions, otherwise known as vertex models, was
pioneered by Honda and collaborators to quantify, for example, the cellular
response to a wound in the eye of a cat [7, 8]. Later work began to consider
an energy functional with a quadratic contribution of the surface area [9,
10]. Certainly, two-dimensional, open-source code bases has helped accelerate
this work [11, 12]. Such two-dimensional models have demonstrated rather rich
collective behavior, to name a few [10, 13, 14, 15, 16]. To be more specific,
in a ordered system and in the absence of cellular rearrangements, there
exists a rigidity transition from a rigid phase to a floppy phase as a
dimensionless version of the target surface area, otherwise known as the shape
index, is increased [10]. Such a transition signifies the breakdown of
nonlinear elasticity with an emergent, anomalous coupling between compression
and shear [13]. There also exists another rigidity transition from a rigid
phase to a fluid phase when cellular rearrangements are allowed [3]. In
addition, such models exhibit the phenomenon of micro-demixing [14], shear-
driven solidification [15] and even fracture [16]. Several aspects of the
theoretical richness have been confirmed experimentally make the vertex model
approach a compelling one [17, 14], though other complimentary theoretical
approaches, such as Voronoi models [18, 19], do exist.

There also exists much progress on three-dimensional vertex model approaches.
For instance, over ten years ago Okuda and collaborators outlined an algorithm
to simulate such models that included reconnection events, or cells moving
past each other [20]. They also explored cellular swelling to determine how an
initially flat monolayer of cells deformed into a curved monolayer of cells
[20]. This exploration was done using a different energy functional from the
one that will be implemented here [20]. More recently, Okuda and Sato have
demonstrated that polarized interfacial tension can induce motion of clusters
of cells within a model tissue [21]. Combining Turing models with three-
dimensional vertex models as well as branching in three-dimensional vertex
models has also been explored, to name a few [22, 23].

Given these initial results, one wonders how the energy functional affects the
physics of such three-dimensional vertex models. Very recent work begins to
address that very question using the restoring surface area and volume springs
to find that there exists a rigidity transition in a cellular-based vertex
model [5]. This work also identified a new boundary-bulk effect in the spatial
organization of cells in which boundary cells are much more aligned than bulk
cells [5]. Even more recent work on multiscale modeling addresses how such a
three-dimensional, cellular-based vertex model can be coupled with a minimal
mechanical model for a cell nucleus to address foundational questions in
biology [24].

While there are a multitude of directions that can be explored in this three-
dimensional, cellular-based vertex model, inspired by recent work focusing on
elastic moduli in an ordered version in two dimensions [25], here, we focus on
a mean field approach to determine the elastic moduli in three dimensions. To
do so, we will revisit the two-dimensional calculations for hexagonal cells
and expand the results to include octogonal cells. This expansion will help us
better understand our three-dimensional results for truncated octahedra. We
will determine the shear modulus, bulk modulus, and Young’s modulus as a
function of a dimensionless shape index in both two and three dimensions for
the three different models. In doing so, we will be able to formulate some
closed form estimates for the location of the rigidity transition. The
manuscript, thus, continues with describing the model and the parameterization
of the cell shape in both two and three dimensions, followed by a description
of deformation protocol in both dimensions. We then present our results in
both dimensions and conclude with a discussion of the implications of our
results.

##  II A cellular-based Vertex Model

###  II.1 The two-dimensional version

Consider a cross-section of a monolayer of cells. The $i$th cellular cross-
section is represented as as a polygon with $A_{i}$ area and $P_{i}$
perimeter. Each cell is assumed to have a target area of $A_{0}$ and a target
perimeter of $P_{0}$. Assuming all edges of a polygon are equivalent, the
elasticity of the cortical tension can, in its simplest form, be expressed as
a perimeter spring (with a target perimeter), while an area spring encodes the
bulk compressibility of a cell. One can, therefore, express the energy of a
collection of $N$ cells, or the tissue energy, $E_{T}$ as
$E_{T}=\frac{1}{2}K_{A}\sum_{i=1}^{N}(A_{i}-A_{0})^{2}+\frac{1}{2}K_{P}\sum_{i=1}^{N}(P_{i}-P_{0})^{2}$
where $K_{P}$ denotes cell perimeter stiffness/contractility and $K_{A}$, cell
area stiffness [3]. Note that the edges between the cells are shared, as
indicated in Fig. 1, representing the adhesiveness of the cell.

To build more intuition about this model, one can readily implement a mean-
field theory approach for an ordered cellular packing in which each cell is
equivalent and, therefore, contributes equally to the packing. In other words,
$E_{T}=N\,E_{S}$, where $E_{S}$ represents the energy of a single cell.
Moreover, we use a dimensionless form of cell energy
$e_{s}=\frac{E_{s}}{K_{A}A_{0}^{2}}$ to arrive at

$e_{S}=\frac{1}{2}(a-1)^{2}+\frac{1}{2}k_{r}(p-p_{0})^{2},$ (1)

with $a=\frac{A}{A_{0}},p=\frac{P}{\sqrt{A_{0}}}$, and
$p_{0}=\frac{P_{0}}{\sqrt{A_{0}}}$. In addition, $k_{r}$ stands for the
rigidity stiffness and is expressed as $\frac{K_{P}}{K_{A}A_{0}}$.
Furthermore, the single cell’s target dimensionless area is unity. We will
assign $p$ to be the dimensionless shape index in two dimensions. Now that the
energy functional for a single cell has been defined, let move to further
characterize its shape.

####  II.1.1 Hexagonal Model

We begin with a single cell as a hexagon, which is parameterized by two edge
lengths as illustrated in Fig. 2. The cell is embedded in a box of area
$L_{x}\times L_{y}$, where $L_{x}$ and $L_{y}$ are the lengths of the box in
the $x$ and $y$ directions, respectively. Note that the box is indicated as a
shaded gray in Fig. 2.

$l_{1}$$l_{1}$$l_{2}$$l_{2}$$l_{2}$$l_{2}$$\theta$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$O$
Figure 2: Parameterization of the hexagonal cell model embedded in an $L_{x}$
by $L_{y}$ box.

Moreover, the angle between the two lines $l_{1}$ and $l_{2}$ only is denoted
by $\theta$. The vector of each vertex in relation to the origin $O=\\{0,0\\}$
can be stated as follows (noting that $\cos\theta\leq 0$ as $\theta$ is an
obtuse angle):

$\displaystyle v_{1}$
$\displaystyle=\bigg{\\{}\frac{l_{1}}{2}-l_{2}\cos\theta,0\bigg{\\}},v_{2}=\bigg{\\{}\frac{l_{1}}{2},l_{2}\sin\theta\bigg{\\}}$
$\displaystyle v_{3}$
$\displaystyle=\bigg{\\{}-\frac{l_{1}}{2},l_{2}\sin\theta\bigg{\\}},v_{4}=\bigg{\\{}-\frac{l_{1}}{2}+l_{2}\cos\theta,0\bigg{\\}}$
$\displaystyle v_{5}$
$\displaystyle=\bigg{\\{}-\frac{l_{1}}{2},-l_{2}\sin\theta\bigg{\\}},v_{6}=\bigg{\\{}\frac{l_{1}}{2},-l_{2}\sin\theta\bigg{\\}}.$

The perimeter of the shape is $p=2l_{1}+4l_{2}$, while the area is equal to
$-2\frac{1}{2}(2l_{2}\sin\theta)(l_{2}\cos\theta)+2{l_{1}}l_{2}\sin\theta=2l_{2}\sin\theta(l_{1}-l_{2}\cos\theta)$.
Accordingly, the dimensionless energy of a single cell (with a unit area) is

$\displaystyle e_{S}$
$\displaystyle=\frac{1}{2}(a-1)^{2}+\frac{1}{2}k_{r}(p-p_{0})^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}2l_{2}\sin\theta(l_{1}-l_{2}\cos\theta)-1\big{)}^{2}+\frac{1}{2}k_{r}(2l_{1}+4l_{2}-p_{0})^{2}.$
(2)

Note that for a regular hexagon with unit area,
$p=2l_{1}+4l_{2}=\sqrt{8\sqrt{3}}$ and
$2l_{2}\sin(\frac{2\pi}{3})(l_{1}-l_{2}\cos(\frac{2\pi}{3}))=1$, such that
$l_{1}=l_{2}=\frac{\sqrt{2}}{3^{\frac{3}{4}}}$.

The above energy depends on three parameters, $l_{1}$, $l_{2}$, and $\theta$.
This parameterization makes it clear that area and perimeter do not alone
uniquely define cell shape, as pointed out in Ref. [25]. To look for a
compatible-incompatible transition, one can minimize the energy with respect
to these three parameters to find that for $p_{0}<p_{0}^{*}(6)$, the energy is
nonzero with the cell unable to achieve both the target area and perimeter,
where $p_{0}^{*}(6)$ is the shape index for a regular hexagon with unit area.
While for $p_{0}\geq p_{0}^{*}(6)$, the cell is able to achieve its target
perimeter and area, hence, the cell’s shape is compatible with its energy.
Please see Appendix A for details. Similar analysis was done in Ref. [25] for
a slightly different shape parameterization.

####  II.1.2 Octagonal Model

In looking towards the three-dimensional model, it will be useful to study an
octogonal cell model in two dimensions, despite the fact that the octogonal
cell alone cannot tile two-dimensional space. Interestingly, however,
octogonal cells can achieve rectangular shapes.

We parameterize the two-dimensional octagonal model in the following way (see
Fig. 3).

$l_{1}$$l_{1}$$l_{1}$$l_{1}$$l_{2}$$l_{2}$$l_{2}$$l_{2}$$v_{8}$$v_{7}$$v_{6}$$v_{5}$$v_{4}$$v_{3}$$v_{2}$$v_{1}$$\theta$$\theta$$\phi=\frac{3\pi}{2}-\theta$$\phi$$\phi$$O$
Figure 3: Parameterization of octagonal cell model.

Assuming the octagon is symmetric along $x,y$ direction about the origin
$O=\\{0,0\\}$ and setting $\theta$ as the angle between the horizontal $l_{1}$
and $l_{2}$, as previously, the remaining angles will be automatically defined
from the following coordinate vectors of each vertex:

$\displaystyle v_{1}$
$\displaystyle=\bigg{\\{}\frac{l_{1}}{2}-l_{2}\cos\theta,-\frac{l_{1}}{2}\bigg{\\}},v_{2}=\bigg{\\{}\frac{l_{1}}{2}-l_{2}\cos\theta,\frac{l_{1}}{2}\bigg{\\}},$
$\displaystyle v_{3}$
$\displaystyle=\bigg{\\{}\frac{l_{1}}{2},\frac{l_{1}}{2}+l_{2}\sin\theta\bigg{\\}},v_{4}=\bigg{\\{}-\frac{l_{1}}{2},\frac{l_{1}}{2}+l_{2}\sin\theta\bigg{\\}},$
$\displaystyle v_{5}$
$\displaystyle=\bigg{\\{}-\frac{l_{1}}{2}+l_{2}\cos\theta,\frac{l_{1}}{2}\bigg{\\}},v_{6}=\bigg{\\{}-\frac{l_{1}}{2}+l_{2}\cos\theta,-\frac{l_{1}}{2}\bigg{\\}},$
$\displaystyle v_{7}$
$\displaystyle=\bigg{\\{}-\frac{l_{1}}{2},-\frac{l_{1}}{2}-l_{2}\sin\theta\bigg{\\}},v_{8}=\bigg{\\{}\frac{l_{1}}{2},-\frac{l_{1}}{2}-l_{2}\sin\theta\bigg{\\}}.$

With the above parametrization, the dimensionless single cell energy function
(of unit target area) will be

$\displaystyle e_{S}$
$\displaystyle=\frac{1}{2}(a-1)^{2}+\frac{1}{2}k_{r}(p-p_{0})^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}l_{1}(l_{1}-2l_{2}\cos\theta)+2l_{2}\sin\theta(l_{1}-l_{2}\cos\theta)-1\big{]}^{2}$
$\displaystyle+\frac{1}{2}k_{r}(4l_{1}+4l_{2}-p_{0})^{2}.$ (3)

When $\theta=\frac{3\pi}{4}$ and $a=1$, then
$l_{1}=l_{2}=\frac{1}{\sqrt{2+2\sqrt{2}}}$. We note that there is a similar
shape transition from incompatible to compatible shapes as the target shape
index is increased with $p<p_{0}^{*}(8)$.

###  II.2 The three-dimensional version

To generalize the two-dimensional vertex model discussed in the prior
subsections, we will simply transform areas $A$ to volumes $V$ and perimeters
$P$ to areas $A$ to arrive at [5]

$E_{T}=\frac{1}{2}K_{V}\sum_{i}(V_{i}-V_{0})^{2}+\frac{1}{2}K_{S}\sum_{i}(A_{i}-A_{0})^{2}.$
(4)

Here $K_{V}$ represents the cell volume stiffness and $K_{S}$ the cell surface
area stiffness. Since $L\sim V^{\frac{1}{3}}$ and $L\sim A^{\frac{1}{2}}$ with
$L$ as the length unit in the three-dimensional scenario, we can express a
dimensionless version of tissue energy with
$e_{T}=\frac{E_{T}}{K_{V}V_{0}^{2}}$ as

$e_{T}=\frac{1}{2}(v-1)^{2}+\frac{1}{2}k_{\rho}(s-s_{0})^{2},$ (5)

where $v,s,s_{0}$ are given by
$v=\frac{V}{V_{0}},s=\frac{A}{(V_{0})^{\frac{2}{3}}}$, and
$s_{0}=\frac{a_{0}}{(v_{0})^{\frac{2}{3}}}$. Finally,
$k_{\rho}=\frac{K_{S}}{K_{V}{V_{0}}^{\frac{3}{2}}}$. Note that we are now
using $s$ to denote the dimensionless surface area (as opposed to $a$) and it
is also referred to as the shape index.

To further parameterize the model, we implement the following criteria:

  1. 1.

A single polyhedron whose copies fill three-dimensional space.

  2. 2.

The polyhedron must have flat, non-curved surfaces.

  3. 3.

To ultimately compare with earlier work, each vertex in the tiling has four
edges attached to it.

  4. 4.

Consider an ordered packing.

The initial requirement fulfills the confluency requirement implemented in
standard vertex models. The second criterion is also typically implemented in
standard vertex models and helps make for easier computations. The third
requirement originates from the model with reconnections that we would like to
compare [20, 5]. As for the last requirement, measuring mechanical properties
of disordered cellular packings is more complex than for ordered cellular
packings and readily allows us to implement a mean field theory approach,
though an effective medium theory has been implemented for disordered packings
in two dimensions [26].

A truncated octahedron shape (see Fig. 1) meets the three criteria above and
so we form an ordered packing of them (fourth criterion). Furthermore, the
planar view of the truncated octahedron matches the 4-8 Lattice outlined in
Ref. [10].

Figure 4: A graphical representation of the bottom half of a truncated
octahedron through stereographic projection. Instead of curvy lines, it is
depicted with straight lines between projected points.

A stereographic projection of a bottom half of a truncated octahedron shows a
shape with one square and four hexagons, as illustrated in Fig. 4. The
projection is conformal, meaning it maintains angles and shapes of the curve.
By stretching the four vertices diagonally, we observe a lower hemisphere of a
truncated octahedron embedded in an octagon. Fig. 5 visualizes this process.

Figure 5: (a) Through stretching the four vertices diagonally, the bottom half
of truncated can be expressed as four hexagons and a square $S$ within an
octagon. $D$ is used to symbolize a diamond shape on the side connecting the
top and bottom halves. (b) Parameters to identify a hexagon tilted in space
from a planar view.

We observe that the right side of Fig. 5a is analogous to a truncated
octahedron projected onto the $xy$-plane. Observe that the hexagon is smaller
than the octagon due to the fact that it is tilted in space. As in Sec.
II.1.2, the octagon produced from a projection for a regular truncated
octahedron has two separate edge lengths ($l_{1},l_{2}$). Additionally, the
side edge length $\frac{l_{2}^{\prime}}{2}$ of the hexagon on a plane can be
calculated using vertical coordinate $z$ and parameters $l_{1},l_{2},\psi$ if
points $A,O$ are known. Therefore, we can use simple parametrization from the
octagon shown below.

$l_{1}$$l_{1}$$l_{2}$$l_{2}$$\theta$$\phi$ Figure 6: Octagonal cell model for
a truncated octahedron.

As defined in Section II.1.1, it is possible to define a hexagon using two
edge lengths and one angle. The octagonal cell model from Fig. 5b, however, is
characterized by $l_{1},\frac{l_{2}}{2}$, and $\theta$ (or
$\phi=\frac{3\pi}{2}-\theta$) which coincides with the octagonal model from
Section II.1.2. The octagonal cell model from Fig. 5b is managed by
$l_{1},\frac{l_{2}}{2}$, and $\theta$ (or $\phi=\frac{3\pi}{2}-\theta$), which
coincides with the octagonal model from Section II.1.2. The $z$ directional
coordinate is chosen to satisfy the regular octahedron at
$l_{1}=1,l_{2}=\sqrt{2},\theta=\frac{3\pi}{4}$, and $z=\sqrt{2}$. Within mean
field, we are now ready to specify the dimensionless single cell energy as

$e_{S}=\frac{1}{2}(v-1)^{2}+\frac{1}{2}k_{\rho}(s-s_{0})^{2},$ (6)

where

$\displaystyle v=2z(l_{1}-l_{2}\cos\theta)(l_{1}+l_{2}\sin\theta)$
$\displaystyle
s=2\bigg{[}\frac{l_{2}z^{2}\sin\theta}{\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}}-l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}+l_{2}\cos\theta^{2}\big{(}$
$\displaystyle\frac{2l_{1}l_{2}}{\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}}+\frac{l_{2}^{2}\sin\theta}{\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}}+\frac{z}{\sqrt{\cos\theta^{2}+\sin\theta^{2}}}\big{)}$
$\displaystyle+l_{2}\sin\theta^{2}\big{(}\frac{z}{\sqrt{\cos\theta^{2}+\sin\theta^{2}}}+\frac{2l_{1}l_{2}}{\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}}\big{)}$
$\displaystyle+l_{1}\bigg{(}l_{1}+2z^{2}\big{(}\frac{1}{\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}}+\frac{1}{\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}}\big{)}\bigg{)}\bigg{]}.$

Both $v$ and $s$ above have been computed using surface triangulation from
Mathematica [27]. As before, one can also find $l_{1},l_{2}$ values for unit
volume with $l_{1}=\frac{1}{2\cdot
2^{\frac{1}{6}}},l_{2}=\frac{\sqrt{2}}{2\cdot 2^{\frac{1}{6}}}$. The
dimensionless shape index $s$ is $s\simeq 5.31474$. Just in the case of the
two-dimensional hexagon, there is a compatible-incompatible transition at the
target shape index of the regular truncated octahedron. See Appendix A for the
analysis.

##  III Methods: Deformation Protocol and Determining Elastic Moduli

To determine the elastic moduli in mean field, we have implemented a similar
deformation protocol in Ref. [25], which is a basic linear deformation for an
affine transformation. In two dimensions, this can be expressed by
$F\mathbf{x_{i}}+g$ with $g=0$ and $i$ is the index of a vertex of a polygon.
To express it simply, $F\mathbf{x}$ is used to refer to all coordinate vectors
of the vertices for $n-$gons transformed by the matrix
$F=\begin{bmatrix}C_{x}&0\\\ 0&C_{y}\end{bmatrix}$
$\mathbf{x_{i}}=\begin{bmatrix}x_{i}\\\ y_{i}\end{bmatrix}$.

Figure 7: Affine deformation of the box.

To measure the elastic moduli, for a given set of parameters, each cell
initially has a unit area or volume, depending on the dimension, with
$C_{x}=C_{y}=1$. Then, a small change of $\delta=0.0001$ is made to either one
or both sides, depending on the elastic modulus under exploration. Please see
Table I for the different types of deformations. The size of box $L_{x},L_{y}$
is derived from the coordinates’ maximum and minimum.

With this deformation, the vertices move in response. Given where the vertices
are now positioned, we compute the perimeter/area and area/volume to determine
the change in the dimensionless single cell energy $e_{S}$ without any energy
minimization. As the configuration may not necessarily be an optimal one with
force-balance no longer achieved, this configuration, and its associated
energy, is dubbed the constrained state, following Ref. [25]. We also study
the energy minimized state given the hidden $\theta$ degree of freedom, for
example, which is the relaxed state. We then apply additional deformation by
increasing $\delta$. The single cell energy $e_{S}$ is then re-calculated for
both the constrained state and the relaxed state for each strain step.

We then compute the curvature of the energy function $e_{S}$ as a function of
$\delta$, which is proportional to the elastic moduli. Since the dimensionless
shear modulus $G$ is defined as
$G=\frac{1}{a(0)}\frac{\partial^{2}e_{s}}{\partial\delta^{2}}$ in two
dimensions, with $a(0)=1$, the curvature equates to the dimensionless shear
modulus. Notating $K(e)$ as a curvature of the energy function, here is how
the two states are characterized:

$\displaystyle\text{Constrained}\colon\lim_{C_{x}\to 1}$ $\displaystyle
K(e_{S}(\text{fixed vertex coordinates}))$
$\displaystyle\text{Relaxed}\colon\lim_{C_{x}\to 1}$ $\displaystyle
K(\min{[e_{S}(\text{vertex coordinates }}$ $\displaystyle\text{with a given
parameterization})]).$

In the constrained case, we simply evaluate the curvature as a function of
strain to compute the dimensionless elastic moduli for each target
$p_{0}/s_{0}$, the parameter that allows the system to toggle between
compatible and incompatible states, should the boundary exist. We also compute
the elastic moduli for the relaxed case using energy minimization. As Table 2
illustrates, in the relaxed case, some non-affine deformations may result
given the hidden degree of freedom. We chose unbounded relaxation of energy,
in contrast to Ref. [25] who measured relaxation with boundary conditions. The
number of variables will drop if we fix the boundary for a single cell, and we
can only use non-affine relaxation. When the box size is determined, our only
option is to discover a shape that meets the target
parameters($a,p_{0}/v,s_{0}$) within the box. As for the additional details on
the energy minimization procedure with Mathematica [27], we evaluated the
solutions based on these four criteria:

  * •

We assume small shifts $\delta$ will result in a continuous deformation, or
relaxation, of the structure.

  * •

Given the possibility of a transition between compatible and incompatible
states, we assume continuous energy values for the target $p_{0}$ such that
the configurations $l_{1},l_{2},\theta,\dots$ of $p_{0}$ and $p_{0}+\epsilon$
($\epsilon\simeq 0.1$) are similar.

  * •

Minimization is allowed within the convex region of the shape
($\frac{\pi}{2}\leq\theta\leq\pi$) as long as it stays within the physically
reasonable range ($l_{1},l_{2}\geq 0$).

  * •

In order to find solutions, we lowered the accuracy of minimization function
when the minimum energy value was not zero.

Since the moduli depend on curvatures of the energy, we compare two methods
for computing the curvature using a discrete method and a continuous method.
Following the parameters discussed in Sec. II.1.1, we simulate the mechanical
responses of a hexagonal cell by measuring the curvature of energy function
according to the deformation protocol shown in Table 1. We examine two
distinct methods for the constrained case to verify the shape before examining
the relaxed case.

We start by computing the discrete curvature of the energy function. To
calculate Young’s modulus, we measure $e_{s}$ by varying $C_{x}$ over
$(1-0.0005,1+0.0005)$. We then calculate the discrete curvature $K$ for every
data point. We acquire the value of $K$ for $C_{x}=1$ through linear fitting.
The second case is fully analytical. Referring to Section II.1.1, where we
have the function $e_{S}$, defined by parameters $l_{1},l_{2},\theta$, the
continuous arc curvature is calculated and displayed in Fig. 8.

$\displaystyle\text{Young's modulus}\colon$ $\displaystyle
K(e_{s}(l_{1},l_{2},C_{x},1,\theta,k_{r},p_{0}))$ $\displaystyle\text{Shear
modulus}\colon$ $\displaystyle
K(e_{s}(l_{1},l_{2},C_{x},1/C_{x},\theta,k_{r},p_{0}))$
$\displaystyle\text{Bulk modulus}\colon$ $\displaystyle
K(e_{s}(l_{1},l_{2},\sqrt{C_{x}},\sqrt{C_{x}},\theta,k_{r},p_{0})).$

We varied $C_{x}$ and $p_{0}$ and held the other parameters fixed
($l_{1}=l_{2}=\frac{\sqrt{2}}{3^{3/4}},\theta=\frac{2\pi}{3},k_{r}=0.1$).
Figure 8 demonstrates that there is effectively no difference between the two
approaches, so we move forward using the discrete method.

Figure 8: Curvature for the constrained case using two different methods:
Curves A,C,E are generated using discrete curvature, while curves B,D,F are
formed using continuous curvature.

Note that when performing the computations, we typically acquire approximately
10-20 points close to $C_{x}=1$ (for instance $C_{x}=1\pm 0.0005$). In some
instances, there may be distinct solutions for $C_{x}=1+\delta$ and $C_{x}=1$.
For some cases, it is related to extreme values (such as $\theta=\pi,\pi/2$).
The higher number of samples is selected when there are two groups of
solutions from the minimizer.

We have picked $k_{r}=0.1$ for the two-dimensional model and
$K_{S}=1,K_{V}=10$ for the three-dimensional one. We then measure discrete
curvature and obtain a value at $C_{x}=1$ via linear fitting. A simple
algorithm has been put in place for the discrete curvature from a triangle.
The formula for $K(\Delta)$ is $\frac{4\text{ triangle area}}{\text{product of
$3$ edges}}$, where$\Delta$ is a triangle with neighboring points.

In addition to measuring the dimensionless shear modulus, bulk modulus, and
Young’s modulus for a given deformation protocol, Poisson’s ratio can be
calculated from
$\frac{1}{2}\frac{(C_{y}^{\prime}-C_{y})L_{y}/(C_{y}L_{y})}{(C_{x}^{\prime}-C_{x})L_{x}/(C_{x}L_{x})}=\frac{1}{2}\frac{\Delta
C_{y}/C_{y}}{\Delta C_{x}/C_{x}}$. In order to measure the deformation of one
side, we divided $\Delta C_{y}$ by 2, since $C_{y}$ applied to both sides of
the box.

Table 1 and 2 summarize the deformation protocol and control parameters.

Table 1: Deformation Protocol | Young’s | Shear | Bulk | Poisson’s ratio
---|---|---|---|---
$C_{x}$ | $1\pm\delta$ | $1\pm\delta$ | $\sqrt{1\pm\delta}$ | $1\pm\delta$
$C_{y}$ | 1 | $\frac{1}{1\pm\delta}$ | $\sqrt{1\pm\delta}$ | -
$C_{z}$ (for 3d) | - | - | $\sqrt{1\pm\delta}$ | -
Table 2: Control parameters Parameters | Purpose | Usage
---|---|---
$C_{y}$ | Relaxed (affine) | Young’s
$C_{y},\theta$ | Relaxed (non-affine) | Young’s
$l_{2},\theta$ | Relaxed (non-affine) | Shear/Bulk
$\theta$ | Relaxed (angle) | Transition

For our three-dimensional protocol, taking the truncated octahedron shape as
the reference, we set up minimal parameters to control regular polyhedron
under an affine transformation with a simple rescaling of the now three-
dimensional box. We only analyze the alteration of coordinates $\mathbf{x}$
within a box $(L_{x},L_{y},L_{z})$ defined by $F\mathbf{x}$ for
$F=\begin{bmatrix}C_{x}&0&0\\\ 0&C_{y}&0\\\ 0&0&C_{z}\end{bmatrix}$
$\mathbf{x}=\begin{bmatrix}x\\\ y\\\ z\end{bmatrix}$. The regular truncated
octahedron is linear in the $z$ direction, meaning there are just five values
($-z,-\frac{z}{2},0,\frac{z}{2},z$). The other parameters will be established
inside the octagon as indicated in Fig. 9 (b). Taking a look at Fig. 9, square
$E$ must be a rectangle with a right angle to prevent a non-convex polyhedron
when scaling in $x$ and $y$ with $C_{z}=1$. Furthermore, this polyhedron will
be symmetric along the $xy$-axis, which means two hexagons $A$ and $B$ will be
equivalent (symmetric along the $y$-axis) with $C$ and $D$ being equivalent
(symmetric along the $x$-axis). As a result, Fig. 9 (b) can be uniquely
defined by $l_{1},l_{2},l_{3},\theta$ (recall $\phi=\frac{3\pi}{2}-\theta$) up
to a translation.

To better understand this simple extension to three dimensions and demonstrate
that additional parameterization is not required, consider the following. The
truncated octahedron is made up of $24$ vertices. Taking into account affine
transformation that are rescalings, all coordinate vectors can be determined
by the manipulation of four vertices ($v_{1}\sim v_{4}$). As illustrated in
Fig. 9 (b), $8$ vertices are combined to create an octagon (black nodes). The
upper half contains eight points placed diagonally and a rectangle. Thus,
$8+8\times 2$ results in $24$. Hexagons $A$ and $B$ are the same, and we can
denote the three-dimensional coordinate of $C$ as $v_{1}=(r_{1},r_{2},r_{3})$,
$v_{4}=(r_{4},r_{5},r_{6})$, $v_{5}=(-r_{1},r_{2},r_{3})$, and
$v_{8}=(-r_{4},r_{5},r_{6})$. Using the imposec symmetries for $x,y,z$, all
other points can be described in the same way. Additionally, since the
octagon’s $8$ points are on the $xy$-plane ($z=0$), we can set
$v_{1}=(r_{1},r_{2},r_{3})$, $v_{2}=(r_{7},r_{8},0)$,
$v_{3}=(r_{9},r_{10},0)$, and $v_{4}=(r_{4},r_{5},r_{6})$. It is evident that
$r_{1}=r_{9}$ and $r_{2}=r_{8}$ due to the fact that they are determined by
$l_{1},l_{3}$. In addition, each hexagon is symmetrical in the $x,y$ plane,
thus, we get
$v_{4}=(r_{4},r_{5},r_{6})=(\frac{r_{1}+r_{7}}{2},\frac{r_{1}+r_{10}}{2},\frac{r_{3}}{2})$.
It is possible to visualize that each hexagon on the upper half has one edge
(with $v_{1}$) at $z=r_{3}$ and its opposite side (with $v_{2}$) at $z=0$.
Rewrite the four vertex points as $v_{1}=(r_{1},r_{2},r_{3})$,
$v_{2}=(r_{7},r_{2},0)$, $v_{3}=(r_{1},r_{10},0)$, and
$v_{4}=(\frac{r_{1}+r_{7}}{2},\frac{r_{1}+r_{10}}{2},\frac{r_{3}}{2})$, and to
control all vertices, only $r_{1},r_{2},r_{3},r_{7},r_{10}$ must be known,
which can be calculated from the parameters $l_{1},l_{2},l_{3},\theta,z$
relatively easily.

$A$$C$$D$$E$ (a) Truncated octahedron.

$\frac{l_{2}}{2}$$l_{3}$$l_{1}$$\theta$$\phi$$l_{2}$$L_{x}$$L_{y}$$A$$B$$C$$D$$E$$O$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$v_{6}$$v_{7}$$v_{8}$
(b) Schematic of 3d model.

Figure 9: A 3d representation of a vertex model with $l_{1},l_{2},l_{3}$ (a)
and (planar) schematic (b). The box unit indicated by ($L_{x}\times L_{y}$) is
represented by the shaded area (b).

Because working out analytical solutions in three dimensions is cumbersome, we
will not compare the discrete and analytical data for the constrained case.
The shape can easily be predicted from the two-dimensional one. We make use of
simple linear deformation, as described in the preceding section. The only
alteration is adjusting three-dimensional expansion, or
$C_{x}=C_{y}=C_{z}=\sqrt{1\pm\delta}$, to measure the bulk modulus.

##  IV Results

###  IV.1 The two-dimensional model

The moduli and Poisson’s ratio are measured for both the constrained and
relaxed case as $p_{0}$ is varied from $0$ to $4.8$.

####  IV.1.1 Hexagonal Model

Figure 10 demonstrates the mechanical behavior for both the constrained and
relaxed cases. The following plots are generated using discrete curvatures, as
specified previously.

Let us understand the shape of these curves by focusing on the constrained
case first. To measure mechanical response, we apply an affine deformation. If
we express energy function in the presence of the deformation as
$e_{S}(\delta)$, we write
$e_{S}(\delta)=\frac{1}{2}(a(\delta)-1)^{2}+\frac{1}{2}k_{r}(p(\delta)-p_{0})^{2}$
where $\delta$ implies $C_{x}\to 1+\delta,C_{y}\to 1/(1+\delta)$ for the shear
modulus, $C_{x}\to\sqrt{1+\delta},C_{y}\to\sqrt{1+\delta}$ for the bulk
modulus, and $C_{x}\to 1+\delta$ for Young’s modulus. Since $C_{x}C_{y}=1$ for
the shear modulus, $a(\delta)=1$ if $a=1$ (the initial condition). On the
other hand, $C_{x}C_{y}=1+\delta$ for the bulk modulus and the Young’s
modulus, so $a(\delta)\neq 1$. We can calculate the curvature functions for
each case as demonstrated below. Specifically,

$\displaystyle K(e_{S}(\delta))_{\text{B,Y}}$
$\displaystyle=|(a(\delta)-1)a^{\prime\prime}(\delta)+(p(\delta)-p_{0})p^{\prime\prime}(\delta)+{a^{\prime}(\delta)}^{2}$
$\displaystyle+{p^{\prime}(\delta)}^{2}|/(1+$
$\displaystyle((a(\delta)-1)a^{\prime}(\delta)+(p(\delta)-p_{0})p^{\prime}(\delta))^{2})^{3/2},$
(7) $\displaystyle K(e_{S}(\delta))_{\text{G}}$
$\displaystyle=\frac{|(p(\delta)-p_{0})p^{\prime\prime}(\delta)+{p^{\prime}(\delta)}^{2}|}{(1+((p(\delta)-p_{0})p^{\prime}(\delta))^{2})^{3/2}},$
(8)

where the primes denote derivatives with respect to $\delta.$ For the hexagon,
because $a(\delta)$ is same for both bulk and Young’s moduli, we arrive at
$(a(\delta)-1)a^{\prime\prime}(\delta)=0$ and
$a^{\prime}(\delta)=2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta$. Moreover,

$\displaystyle\lim_{\delta\to
0}(p(\delta)-p_{0})p^{\prime}(\delta)_{\text{G}}=$ $\displaystyle
2k_{r}(2l_{1}+4l_{2}-p_{0})(l_{1}+2l_{2}\cos 2\theta)),$ (9)
$\displaystyle\lim_{\delta\to
0}(p(\delta)-p_{0})p^{\prime}(\delta)_{\text{B}}=$ $\displaystyle
k_{r}(2l_{1}+4l_{2}-p_{0})(l_{1}+2l_{2}),$ (10) $\displaystyle\lim_{\delta\to
0}(p(\delta)-p_{0})p^{\prime}(\delta)_{\text{Y}}=$ $\displaystyle
2k_{r}(2l_{1}+4l_{2}-p_{0})$ $\displaystyle\cdot(l_{1}+l_{2}+l_{2}\cos
2\theta).$ (11)

Finally, we obtain

$\displaystyle\lim_{\delta\to
0}K(e_{S}(\delta))_{\text{G}}=|(2k_{r}(2l_{1}^{2}+6l_{1}l_{2}+16l_{2}^{2}-3l_{2}p_{0}$
$\displaystyle+2l_{2}(2l_{1}-4l_{2}+p_{0})\cos
2\theta+l_{2}(-2l_{1}+p_{0})\cos 4\theta))|$
$\displaystyle/(1+4k_{r}^{2}(-2l_{1}-4l_{2}+p_{0})^{2}(l_{1}+2l_{2}\cos
2\theta)^{2})^{\frac{3}{2}},$ (12) $\displaystyle\lim_{\delta\to
0}K(E(\delta))_{\text{B}}=|(4(4l_{1}^{2}l_{2}^{2}+l_{2}^{4}+k_{r}(l_{1}+2l_{2})p0$
$\displaystyle-4l_{1}l_{2}^{3}\cos\theta-4l_{1}^{2}l_{2}^{2}\cos
2\theta+4l_{1}l_{2}^{3}\cos 3\theta-l_{2}^{4}\cos 4\theta))|$
$\displaystyle/(4+(2k_{r}(l_{1}+2l_{2})(2l_{1}+4l_{2}-p_{0})$
$\displaystyle+4l_{2}(-l_{1}+l_{2}\cos\theta)\sin\theta(1-2l_{1}l_{2}\sin\theta+l_{2}^{2}\sin
2\theta))^{2})^{\frac{3}{2}},$ (13) $\displaystyle\lim_{\delta\to
0}K(E(\delta))_{\text{Y}}=|(4(4k_{r}l_{1}l_{2}\cos\theta^{4}-2l_{1}l_{2}^{3}\cos\theta\sin\theta^{2}$
$\displaystyle+l_{1}^{2}(k_{r}+l_{2}^{2}\sin\theta^{2})$
$\displaystyle+l_{2}\cos\theta^{2}(4k_{r}l_{2}+(6k_{r}l_{1}+l_{2}^{3}-k_{r}p_{0})\sin\theta^{2})))|$
$\displaystyle/(1+1/4(4k_{r}(2l_{1}+4l_{2}-p_{0})(l_{1}+l_{2}+l_{2}\cos
2\theta)$
$\displaystyle+4l_{2}(-l_{1}+l_{2}\cos\theta)\sin\theta(1-2l_{1}l_{2}\sin\theta+l_{2}^{2}\sin
2\theta))^{2})^{\frac{3}{2}}.$ (14)

Notice that effects of $C_{x}\to 1+\delta$ and
$C_{x}\to\sqrt{1+\delta},C_{y}\to\sqrt{1+\delta}$ act differently for the cell
perimeter. These differences will lead to different curves for the bulk and
Young’s moduli. For the shear modulus curvature function, if
$l_{1}=l_{2}=\frac{\sqrt{2}}{3^{3/4}}$, $\theta=\frac{2\pi}{3}$, and
$k_{r}=0.1$, then denominator becomes 1. Thus, we have

$\displaystyle\lim_{\delta\to 0}K(e_{S}(\delta))_{\text{G}}\simeq
0.2|2.79181(3.72242-p_{0})|.$ (15)

This function becomes linear and has a minimum at $p_{0}=3.72242$, as
indicated in Fig. 9A. However, $\lim_{\delta\to 0}K(E(\delta))_{\text{B}}$ and
$\lim_{\delta\to 0}K(E(\delta))_{\text{Y}}$ (Figs. (B and C)) both have curves
related to $\frac{|a+bx|}{(c+d(e-x)^{2})^{3/2}}$, i.e., they are not linear.
More precisely,

$\displaystyle\lim_{\delta\to 0}K(E(\delta))_{\text{B}}$
$\displaystyle\simeq\frac{4|(2+0.186121p_{0})|}{(4+0.138564(3.72242-p_{0})^{2})^{3/2}},$
(16) $\displaystyle\lim_{\delta\to 0}K(E(\delta))_{\text{Y}}$
$\displaystyle\simeq 4|(0.270335+0.155101(0.248161$
$\displaystyle+0.75(0.611035-0.1p_{0})))|$
$\displaystyle/(1+0.034641(3.72242-p_{0})^{2})^{3/2}.$ (17)

The constrained case acts an upper bound on the curves. Note that we did not
compute Poisson’s ratio as both $C_{x}$ and $C_{y}$ are constrained.

Let us now comment on the less constrained cases. Three separate parameter
sets are used for the relaxed case. Changing only the box size in the $y$
direction ($C_{y}$) softens the Young’s modulus, for instance, as compared to
the constrained case. However, using two parameters ($C_{y},\theta$) allows
the system to soften even more so. In addition, relaxing both $l_{2}$ and
$\theta$ softens the shear modulus as compared to the constrained case. The
same trend observed in the bulk modulus. It is very reasonable to expect this
given the addition degree of freedom. Note that depending in the definition of
the moduli, not all parameters can be minimized. For instance, we do not
compute the shear modulus by minimizing $C_{y}$.

For all three moduli, with a two parameter minimization, a compatible-
incompatible transition is observed with the moduli vanishing for $p\geq
p_{0}^{*}(6)$. Recall that $e_{S}$ reaches a minimum at this target shape
index. The transition also leads to a drop in Poisson’s ratio. Furthermore,
relaxation also takes place with non-affine deformations as $l_{2},\theta$ are
non-affine parameters. Finally, our curves differ slightly from those in Ref.
[25] as our shape parameterization differs slightly from theirs and are only
investigating a single cell (as opposed to four cells with relevant boundary
conditions). By utilizing an unbounded single cell, our strategy allows for
non-affine/affine variables during relaxation and less restrictive geometry
constraints. With a fixed boundary, the box size ($L_{x}\times L_{y}$) is
determined by $C_{x},C_{y}$ values, allowing a cell to relax its energy using
non-affine deformation (ratio of $l_{1},l_{2}$ and a relaxation using
$\theta$). Therefore, there will be a similar mechanical response in an
incompatible state as in a constrained case. The relaxation of energy is more
achievable for a cell in a compatible state.

Figure 10: Mechanical response for hexagon.

####  IV.1.2 Octagonal Model

The next step is to compute the mechanical response of the octagonal cell
described in Sec. II.1.2.

Figure 11: Mechanical response for the octagon cell.

We observe that the mechanical behavior of the octagonal model is, not
suprisingly, rather similar to that of the hexagonal case. See Fig. 10. The
modulus of ($C_{y},\theta$) or ($l_{2},\theta$) will be lower than the
hexagonal case when relaxed, as a result of the greater impact of altering
$C_{y}$ or $l_{2}$ on the perimeter. To directly compare to the two shapes,
Fig. 12 plots the Young’s modulus for two different cell geometries. The
octagonal model has a transition point near $p_{0}^{*}(8)\approx 3.64072$,
which is denoted by the red dotted line, while the transition point for the
hexagonal model is near the shape index for a regular hexagon (with unit
area).

$0$$2$$4$$0$$1\cdot 10^{-1}$$0.2$$0.3$Target shape index $p_{0}$Young’s
modulusHexagonOctagon Figure 12: Young’s modulus for the two cellular models
for the relaxed case with ($C_{y},\theta$) as the minimized parameters. The
solid red line denotes $p_{0}^{*}(6)$, while the dashed red line denotes
$p_{0}^{*}(8)$.

####  IV.1.3 Octagonal model with a different parametrization

To better understand how the mechanics depends on the shape parameterization,
we include an extra edge length $l_{3}$ in addition to $l_{1}$ and $l_{2}$ to
mimic edge length alteration from an affine deformation for a standard
octagon. See Fig. 12. Two unique edge lengths will come from different scaling
in the $x$ and $y$ directions. Moreover, we alter the parameters of the
octagonal cell using the following criteria.

  * •

We are now employing three edge lengths ($l_{1},l_{2},l_{3}$)
instead of two ($l_{1},l_{2}$).

  * •

These lengths can be converted to the box-related parameters $w,h_{1},h_{2}$
that were inspired by Ref. [25]. This can be done so via ($L_{x},L_{y}\colon$
total length in x,y of octagon)

$\displaystyle L_{x}$ $\displaystyle+l_{1}=2l_{1}-2l_{2}\cos\theta=w,$
$\displaystyle L_{y}$ $\displaystyle=h_{1}+h_{2},$ where $\displaystyle
h_{1}=2l_{2}\sin\theta,h_{2}=l_{3}.$

  * •

Revise energy function by employing new variables from
$e_{S}(l_{1},l_{2},l_{3},\theta,k_{r},p_{0})$ to
$e_{S}(w,h_{1},h_{2},\theta,k_{r},p_{0})$.

We select $w=L_{x}+l_{1}$ instead of $w=L_{x}$ to form an unbalanced
parametrization; the main purpose for this is that the area function, when
reparametrized using $w,h$ specified by
$w=l_{1}-2l_{2}\cos\theta,h=l_{1}+2l_{2}\sin\theta$, has a singularity at
$\theta=3\pi/4$
($a=\frac{(h+w)^{2}\cot\theta+2hw\csc\theta^{2}}{2(1+\cot\theta)^{2}}$ and
$\cot(\frac{3\pi}{4})=-1$). With this new parametrization, the cell have the
capacity to shift between three distinct shapes. In other words, we have
modified the shape landscape and will now determine how the new shape
landscape potentially affects the moduli.

  * •

Octagon

  * •

Rectangle ($h_{1}=0$ or $\theta=\pi/2,\pi$)

  * •

Hexagon ($h_{2}=0$ or $\theta=\pi-\delta$ where $\delta$ is related to $w,h$
ratio).

The corresponding angle $\phi$ is equal to $\phi=3\pi/2-\theta$. Fig. 14
illustrates the mechanical response of a reparametrized octagon.

$\frac{h_{1}}{2}$$h_{2}=l_{3}$$l_{3}$$l_{1}$$l_{1}$$l_{1}$$\theta$$\phi$$l_{2}$$w$
Figure 13: Schematic of the reparameterized octogonal model using
$h_{1},h_{2},$ and $w$. Box ($L_{x}\times L_{y}$) is represented by the shaded
area. Figure 14: Shear modulus and Young’s modulus for the reparametrized
octagon model.

If we take $\theta$ as the sole variable for optimizing energy $e_{S}$, two
peaks show up as presented in Fig. 15. The double peak pattern occurs when
$(C_{y},\theta)$ is minimized in the relaxed case. What is the mechanism
behind this two-peaked structure? If we look at Fig. 16, the first case
$(C_{y},\theta)$ displays the cell shape transitioning to a rectangular shape,
which occurs near $p_{0}=4.1$. The emergence of the rectangular shape
constrains the angle $\theta$, thereby leads to an increase in the modulus and
so the minimization is effectively reduced to a one-parameter minimization.
Observe that the Young’s modulus continues to increase for Fig. 16 (b) as
$p_{0}$ is increased beyond approximately 4.1.

$0$$2$$4$$0$$0.2$$0.4$$0.6$$0.8$$1$$p_{0}$Young’s modulusRelaxed
($\theta$)Relaxed ($C_{y},\theta$) Figure 15: Young’s modulus as a function of
$p_{0}$ for the reparameterized octagonal model and for two different relaxed
cases.

(a) $p_{0}=3.3,3.7,3.9,4.3$.

(b) $p_{0}=3,3.7,4.2,4.8$.

$0$$2$$4$$1$$2$$3$$p_{0}$Angle $\theta$Relaxed ($\theta$)Relaxed
($C_{y},\theta$)

(c) Angle values of (a) and (b).

Figure 16: (a,b) Optimized shapes for relaxed conditions for target $p_{0}$
$3-4.8$. (a) For the relaxed case ($C_{y},\theta$), the cell is an octagonal
shape when $p_{0}=3.7$ and is a rectangular shape when $p_{0}=4.3$ (dark to
light) (b) In the relaxed state ($\theta$), when $p_{0}$ is 3.7, the cell has
an octagonal shape, whereas when $p_{0}$ is 4.8, the cell is hexagonal (dark
to light). (c) Angle $\theta$ values for Fig. 15 for the two different relaxed
cases.

###  IV.2 3D Vertex Model

####  IV.2.1 Regular truncated octahedron model

Now that we better understand some of the more subtle points behind how the
shape landscape can depend on the parameterization and, therefore, influence a
system’s mechanical response, we turn to the three-dimensional case. Given the
results of the prior subsection, we first work with our initial
parameterization of the truncated octahedron with $l_{1},l_{2},\theta$,$z$ and
now including $l_{3}$ such that $l_{1}=l_{3}$ when we first introduced the
three-dimensional model. Please note that our evaluation is based on a regular
truncated octahedron with unit volume. If we use an irregular one with unit
volume, the result would be different due to a higher $s_{0}$.

To better understand these curves, let us first do the same analysis for the
constrained case as it serves an as upper bound for the curves using the one-
or two-parameter minimization. If we set $l_{1}=l_{3}=\frac{1}{2\cdot
2^{1/6}}$ and $l_{2}=\frac{\sqrt{2}}{2\cdot 2^{1/6}}$, $K_{V}=10$, $K_{S}=1$,
and $\theta=\frac{3\pi}{4}$, for the constrained case we obtain

$\displaystyle\lim_{\delta\to 0}K(E(\delta))_{\text{G}}$
$\displaystyle\simeq\frac{|((3+4\sqrt{3})(3+6\sqrt{3}-2\cdot
2^{1/3}s_{0}))|}{3\cdot 2^{2/3}},$ (18) $\displaystyle\lim_{\delta\to
0}K(E(\delta))_{\text{B}}$
$\displaystyle\simeq|(576(20+(8138+4440\sqrt{3})^{1/3}))|$
$\displaystyle/(64+17577\cdot 2^{2/3}+8424\cdot 2^{2/3}\sqrt{3}$
$\displaystyle+72s_{0}(-3(37+30\sqrt{3})$
$\displaystyle+2^{1/3}(13+4\sqrt{3})s_{0}))^{3/2}$ (19)
$\displaystyle\text{(note that $z=\frac{\sqrt{2}}{2\cdot 2^{1/6}}$ is scaled
by $\sqrt{1+\delta}$)},$ $\displaystyle\lim_{\delta\to
0}K(E(\delta))_{\text{Y}}$ $\displaystyle\simeq|(12(160+131\cdot
2^{1/3}+42\cdot 2^{1/3}\sqrt{3})$ $\displaystyle-8\cdot
2^{2/3}(3+4\sqrt{3})s_{0})|$ $\displaystyle/(3(16+1953\cdot 2^{2/3}+936\cdot
2^{2/3}\sqrt{3}$ $\displaystyle+8s_{0}(-3(37+30\sqrt{3})$
$\displaystyle+2^{1/3}(13+4\sqrt{3})s_{0}))^{3/2}).$ (20)

As in the two-dimensional case, the curvature as a function of $s_{0}$ for
shear modulus is a linear away from the rigidity transition location of the
regular truncated octahedron with unit volume. Moreover, the other two moduli
have the form of $\frac{|a+bx|}{(c+d(e-x)^{2})^{3/2}}$. The peaks observed in
the constrained case for the Young’s modulus and the bulk modulus are present
here given the choice of $K_{V}$ and $K_{S}$, which was motivated by the
choice in parameters in Ref. [5] and is different form the equivalent relative
stiffness in two dimensions, which is why such peaks were not prominent in the
earlier case. The mechanical response of the constrained bulk modulus exhibits
a higher peak compared to the response of Young’s modulus due to the selection
of $C_{x}C_{y}C_{z}=(\sqrt{1+\delta})^{3}$ instead of
$C_{x}C_{y}C_{z}=({(1+\delta)}^{\frac{1}{3}})^{3}$.

For the relaxed cases where there is at least a two-parameter minimization, we
also observe the same phenomena of a rigidity transition. See Fig. 18. Again,
we observe a drop/cusp in the Poisson’s ratio at the transition point.

Figure 17: Mechanical response for truncated octagon.

|
---|---
(a) Rectangular prism | (b) Rhombic dodecahedron
|
(c) Elongated dodecahedron | (d) Truncated octahedron
|
(e) Staggered hexagonal cell | (f) Trigonal trapezohedron

Figure 18: Transformable structures of the reparamterized truncated octahedron
cell model.

####  IV.2.2 3d Model with different parameterization

Similar to Sec. IV.1.3, the 3D model utilizes the same parametrization methods
by modifying the 3D coordinate vectors from the variables $l_{1},l_{2},l_{3}$
to $w,h_{1},h_{2}$ while $\theta,z$ stay the same. The new parameterization
for the truncated octahedron matches the octagonal model in the Figure 13.
With this new parameterization, additional shape transformations can be
explored such that shape landscape is slightly more complex as compared to the
prior parameterization. Fig. 18 (a), (b), and (f) relate to the lowest values
($h_{1}=0$ or $h_{2}=0$). The value of $s_{0}$ for the trigonal trapezohedron
in Table 3 was determined by using the regular trapezohedron with $n=3$, and
the $s_{0}$ for the staggered hexagonal cell was calculated from the regular
truncated octahedron with an angle of $\pi/2$.

|
---|---
(a) Mechanical response for the reparameterized truncated octahedron model. |
$3$$4$$5$$6$$7$$0$$10$$20$$s_{0}$Young’s modulus$(\theta)$$(C_{y},\theta)$ (a) Mechanical response of 3d model. |
Figure 19: (a) Mechanical response for the reparameterized truncated
octahedron model as a function of the dimensionless shape index. (b)
Mechanical response of 3d model for relaxed ($\theta$) and ($C_{y},\theta$).
The first red line is for a value of $s_{0}=5.28$, and the second one is for
$s_{0}=5.8$.

For the reparameterized version of the truncated octahedron, we find similar
behavior for the shear and bulk modulus as for the first parameterization.
However, for the reparameterized version, even when energy is optimized using
$\theta$, Young’s modulus remains high, since decreasing the energy requires
more than just changing $\theta$. If we measure Young’s modulus for $h_{2}=0$
using ($C_{y},\theta$) and minimize energy $e_{S}$, we can observe a shape
transition in the three-dimensional model from elongated dodecahedron to a
near rhombic dodecahedron and then to trigonal trapezohedron.

$3$$4$$5$$6$$7$$1$$1.5$$2$$2.5$$3$$S_{1}$$S_{2}$$S_{3}$$p_{0}$Angle $\theta$($\theta$)($C_{y},\theta$) (a) Angle values for relaxed ($\theta$) and ($C_{y},\theta$). |
---|---
|  |
---|---
(b) $S_{1}$ ($s_{0}\simeq 5.27$) | (c) $S_{2}$ ($s_{0}\simeq 5.53$)
|
(d) $S_{3}$ ($s_{0}\simeq 5.8$) |
Figure 20: (a) Mechanical response of 3d model for relaxed ($\theta$) and ($C_{y},\theta$). The first red line is for a value of $s_{0}=5.28$, and the second one is for $s_{0}=5.8$. (b), (c), and (d) $S_{1},S_{2},S_{3}$ of (a) for each target shape index $s_{0}$. Table 3: Dimensionless shape index $s_{0}$ for various shapes. |  $s_{0}$ for regular polyhedron
---|---
Rectangular prism (cube) | 6
Rhombic dodecahedron | 5.34539
Elongated dodecahedron | 5.49324
Truncated octahedron | 5.31474
Staggered hexagonal cell | 5.78808
Trigonal trapezohedron | 6

##  V Discussion

We have examined the mechanical response to infinitesimal strains of an
ordered, cellular-based vertex model in both two and three dimensions by
calculating the shear modulus, bulk modulus, Young’s modulus and Poisson’s
ratio within mean field. We find a compatible-incompatible transition in three
dimensions with the target shape index as the tuning parameter, similar to the
transition found in two dimensions [25]. The transition point location occurs
at the target shape index of the truncated octahedron (with $s_{0}\simeq
5.315$), which is a bit a lower than the rigidity transition found earlier
[5]. In this earlier three-dimensional work, cellular reconnection events are
allowed and the system is disordered. The trend in the shift of the location
of the transition resembles the two-dimensional case where the transition from
fluid to rigid state occurs at $p_{0}=3.813$ when T1 events allowed [3], while
the hexagonal model with no T1 events shows a transition at $p_{0}=3.72$ [10,
25].

Shape indeed has an effect on three-dimensional rigidity transitions as it can
either more readily or less readily allow a cell to achieve both the target
surface area and volume. The target shape index $s_{0}$ explores this
readiness. When the shape index is $s_{0}\geq 5.8$, we find the cell is more
likely to be a staggered hexagonal cell (for $h_{2}\neq 0$). See Fig. 21. As
the staggered hexagonal cell has co-planar faces, it can be thought of as a
hexagonal cylinder. We observe that two inner vertices located on the top and
bottom faces can move in either the $x$ or $y$ directions as long as they
remain inside of the hexagon. This can be considered as decreasing the number
of constraints as four independent planes for the top of the structure becomes
one. Moreover, this design makes it easier to maintain the same volume as the
cell area changes. The structure becomes floppy, provided there are a
sufficient number of relaxation parameters.

We have also found that for nontrivial shape transition pathways, particularly
in three dimensions, there can be interesting mechanical properties even
without changes in topology. For instance, the staggered hexagonal cell shape
is achieved via a particular parameterization of the truncated octahedron such
that there are no singularities. This change in the shape landscape, as
described by its shape transition pathways, influences the Young’s modulus,
again, depending on the dimensionality of the relaxation parameter space. We
also observe this effect in two dimensions with a particular parameterization
of the hexagon in which there is a re-emergence of a nonzero Young’s modulus
at higher shape indices.

As for a very simple argument for the rigidity transition, when we look at the
shape index of the tetrahedron, cube, and icosahedron, which are $s_{0}=7.21$,
$s_{0}=6$, and $s_{0}=5.148$ respectively, the larger the value of $s_{0}$,
the less faces the regular polyhedron has. From a mechanical standpoint, it is
easier to adjust and locate the simpler polyhedron in the space as it has
fewer faces. For example, the hexagonal cylinder is more advantageous than the
truncated octahedron since it requires less symmetries to fill space.
Moreover, hexagonal cylinders can readily be transformed into scutoids, which
are found in many biological cellular structures [28, 29, 30]. In comparison,
the truncated octahedron requires symmetry in the diagonal direction to fill
space, and it is more commonly seen in solid structures like [31] or regularly
packed droplets [32]. The previous section shows that, if there is an
intermediate shape able to reduce stress, the cell will transform its shape to
relax as can be seen in Fig. 20. To go beyond this simple argument in terms of
the interplay between geometry and rigidity, as is done in Refs. [33, 34, 35],
as well as higher-order rigidity [36], more work is being done.

Figure 21: Staggered hexagonal cell can be morphed into various structures
with the same volume. (a) Circled vertices can move up and down under the same
edge length. Moreover, the top and bottom can change its size while the
surface remains flat (scutoids). (b) Vertices can merge together and construct
a hexagonal cylinder.

What have we learned here that will help us understand better the rigidity
transition when cellular reconnection events are allowed? Such rules have been
studied using metallic solids, as discussed in [37] and [38], and have also
been applied to biology [39, 20]. Interestingly, Williams in 1968 [40]
demonstrated the reconnection event sequence of a truncated octahedron for
vegetable cells (or soap bubbles, metal grains) based on the cell’s own edges
and vertices. Honda et al. [39] chose a square situated on the side of the
truncated octahedron. In the ordered vertex model, the dimensionless shape
index $s$ $\simeq 5.45477$ structures with two sides, as demonstrated in Fig.
22, and $\simeq 5.58941$ for those with four sides. Assuming a reconnection
event from one side in Fig. 22, the value of $s$ is $\simeq 5.38742$ for a
transition point. Note that $s\simeq 5.38742\sim 5.4$ is supported by the 3D
Voronoi model analysis [4, 41] and the 3D vertex model analysis [5]. However,
the task of achieving complete space filling of a single shape cell is
hampered by the formation of gaps. The use of top and bottom squares can help
achieve space filling in a single vertex model.

As for another reconnection event configuration where gaps are not an issue,
consider Fig. 23, a transition is shown from horizontal face (red line) to a
vertical face instead of a horizontal to vertical line transition, which is a
T1-like event in three dimensions. Instead of the hexagon to pentagon
transition in two dimensions, we chose a hexagon to rectangular transition in
three dimensions (Refer to Fig. 23, hexagonal and rectangular faces are
connected through a red line.). We can augment the T1-like transition with
more steps, as depicted in Fig. 24 in accordance with the reconnection model
articulated in Ref. [20]. We envision cells contracting from face, line, to a
point to alter its direction. In addition, by changing the lengths of the
edges (when $h_{2}\neq 0$) in the Type I state, horizontal cells are converted
into a truncated octahedrons and vertical cells become elongated
dodecahedrons. Moreover, the same process can be done with hexagonal cylinder
(staggered hexagonal cells), implying there are numerous reconnection choices
in three dimensions [20]. Indeed, we are only beginning to understand the
nuances of rigidity transitions in three dimensions that occur in cellular-
based vertex models with and without reconnection events.

Figure 22: The 3d reconnection model (Type O) developed by Honda et al. [39]
was implemented on the truncated octahedrons, represented by the yellow cells
($s\simeq 5.45477$). In order to fill space, it is necessary to replace the
gray-colored cells with other shapes.

|
---|---
(a) Cells before T1 swap | (b) interface contraction
|
(c) Cells after T1 swap |

Figure 23: Constructing truncated octahedrons that have a corresponding
T1-like transition. (a) Cell-cell interface indicated by a red line that is
contracted to a point (b), and subsequently stretched in the perpendicular
direction (c).

Figure 24: Modified reconnection model in an ordered truncated octahedron
model.

###### Acknowledgements.

JMS acknowledges financial support from NSF-PHY-PoLS-2014192. The authors also
acknowledge insightful discussion with Shabeeb Ameen, Danilo Liarte, James
Sethna, and Stephen Thornton.

##  VI Appendix A: Finding minima of the energy functional $E$

###  VI.1 Two dimensions

Suppose energy $E$ is a function of $x$ with $a(x)$ and $p(x)$, then we have

$\displaystyle E(x)$
$\displaystyle=\frac{1}{2}(a(x)-1)^{2}+\frac{1}{2}k_{r}(p(x)-p_{0})^{2},$
(A.1) $\displaystyle\frac{\partial E}{\partial x}$
$\displaystyle=\frac{\partial a(x)}{\partial x}(a(x)-1)+k_{r}\frac{\partial
p(x)}{\partial x}(p(x)-p_{0}).$ (A.2)

If $\frac{\partial a(x)}{\partial x}\neq 0$, then the function has to satisfy
$a(x)=1$ to make a $E(x)=0$. We can apply same strategy for $\frac{\partial
p(x)}{\partial x}\neq 0$ and find necessary conditions to achieve $E(x)=0$.
For hexagon, we have

$\displaystyle p$ $\displaystyle=2l_{1}+4l_{2},$ (A.3) $\displaystyle a$
$\displaystyle=2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta.$ (A.4)

Notice that $\frac{\partial p}{\partial\theta}=0$ and nonzero for other cases
($\frac{\partial p}{\partial l_{1}},\frac{\partial p}{\partial
l_{2}},\frac{\partial a}{\partial\theta},\frac{\partial a}{\partial
l_{1}},\frac{\partial a}{\partial l_{2}}$). We can write more explicitly

$\displaystyle E$
$\displaystyle=\frac{1}{2}(2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta-1)^{2}+\frac{1}{2}k_{r}(2l_{1}+4l_{2}-p_{0})^{2}$
$\displaystyle=\frac{1}{2}(a-1)^{2}+\frac{1}{2}k_{r}(p-p_{0})^{2},$ (A.5)
$\displaystyle\frac{\partial E}{\partial\theta}$
$\displaystyle=2l_{2}(l_{2}\cos 2\theta-
l_{1}\cos\theta)(1+2l_{2}(l_{2}\cos\theta-l_{1})\sin\theta)$
$\displaystyle=2l_{2}(l_{2}\cos 2\theta-l_{1}\cos\theta)(1-a),$ (A.6)
$\displaystyle\frac{\partial E}{\partial l_{1}}$
$\displaystyle=2k_{r}(2l_{1}+4l_{2}-p_{0})$
$\displaystyle+2l_{2}\sin\theta(-1+2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta)$
$\displaystyle=2k_{r}(p-p_{0})+2l_{2}\sin\theta(a-1),$ (A.7)
$\displaystyle\frac{\partial E}{\partial l_{2}}$
$\displaystyle=4k_{r}(2l_{1}+4l_{2}-p_{0})$
$\displaystyle+(-2l_{2}\cos\theta\sin\theta+2(l_{1}-l_{2}\cos\theta)\sin\theta)$
$\displaystyle\cdot(-1+2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta)$
$\displaystyle=4k_{r}(p-p_{0})$
$\displaystyle+(-2l_{2}\cos\theta\sin\theta+2(l_{1}-l_{2}\cos\theta)\sin\theta)(a-1).$
(A.8)

Therefore, we need $a=1$ for $\frac{\partial E}{\partial\theta}=0$ and
$a=1,p=p_{0}$ for $\frac{\partial E}{\partial l_{1}}=0,\frac{\partial
E}{\partial l_{2}}=0$.

###  VI.2 Three dimensions

For the truncated octahedron, we have
$E(x)=\frac{1}{2}K_{V}(v(x)-1)^{2}+\frac{1}{2}K_{S}(s(x)-s_{0})^{2}$ and

$\displaystyle s$
$\displaystyle=2\bigg{(}l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}-l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)},$
(A.9) $\displaystyle v$
$\displaystyle=2z(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta).$ (A.10)

Both $s$ and $v$ satisfy $\frac{\partial s,v}{\partial l_{1}},\frac{\partial
s,v}{\partial l_{2}},\frac{\partial s,v}{\partial l_{3}},\frac{\partial
s,v}{\partial\theta},\frac{\partial s,v}{\partial z}\neq 0$, we need $v=1$ and
$s=s_{0}$ for $E=0$. For $K_{S}=1$, we can also write explicitly as

$\displaystyle\frac{\partial E}{\partial\theta}=$
$\displaystyle\frac{1}{2}\bigg{(}2K_{V}\big{(}-1+2z(l_{1}-l_{2}\cos\theta)$
$\displaystyle(l_{3}+l_{2}\sin\theta)\big{)}\big{(}2l_{2}z\cos\theta(l_{1}-l_{2}\cos\theta)+$
$\displaystyle 2l_{2}z\sin\theta(l_{3}+l_{2}\sin\theta)\big{)}$
$\displaystyle+4\bigg{(}l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle-(2l_{2}^{2}l_{3}\cos\theta\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle-(l_{2}^{3}\cos\theta\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(2l_{1}l_{2}^{2}\cos\theta\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}^{3}\cos\theta^{2}\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}$
$\displaystyle\bigg{(}-s_{0}+2(l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$ $\displaystyle-
l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\bigg{)}\bigg{)}$ (A.11)
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}2K_{V}\big{(}-1+v\big{)}\big{(}2l_{2}z\cos\theta(l_{1}-l_{2}\cos\theta)$
$\displaystyle+2l_{2}z\sin\theta(l_{3}+l_{2}\sin\theta)\big{)}$
$\displaystyle+4\bigg{(}l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle-(2l_{2}^{2}l_{3}\cos\theta\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle-(l_{2}^{3}\cos\theta\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(2l_{1}l_{2}^{2}\cos\theta\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}^{3}\cos\theta^{2}\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}\bigg{(}-s_{0}+s\bigg{)}\bigg{)},$
(A.12) $\displaystyle\frac{\partial E}{\partial z}=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta)$
$\displaystyle\big{(}-1+2z(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta)\big{)}$
$\displaystyle+4\bigg{(}l_{2}+(2l_{3}z)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(l_{2}z\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(2l_{1}z)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}z\cos\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}$
$\displaystyle\big{(}-s_{0}+2(l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-
l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\big{)}\bigg{)}$ (A.13)
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta)\big{(}-1+v\big{)}+$
$\displaystyle 4\bigg{(}l_{2}+(2l_{3}z)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(l_{2}z\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}+(2l_{1}z)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}z\cos\theta)/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}\big{(}-s_{0}+s\big{)}\bigg{)},$
(A.14) $\displaystyle\frac{\partial E}{\partial l_{1}}=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}z(l_{3}+l_{2}\sin\theta)\big{(}-1+2z(l_{1}-l_{2}\cos\theta)$
$\displaystyle(l_{3}+l_{2}\sin\theta)\big{)}+4(l_{3}+2\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})$
$\displaystyle\bigg{(}-s_{0}+2(l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-l2\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\bigg{)}\bigg{)}$
(A.15) $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}z(l_{3}+l_{2}\sin\theta)\big{(}-1+v\big{)}$
$\displaystyle+4(l_{3}+2\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\big{(}-s_{0}+s\big{)}\bigg{)},$
(A.16) $\displaystyle\frac{\partial E}{\partial l_{2}}=$
$\displaystyle\frac{1}{2}\bigg{(}2K_{V}(2z(l_{1}-l_{2}\cos\theta)\sin\theta-2z\cos\theta$
$\displaystyle(l_{3}+l_{2}\sin\theta))\big{(}-1+2z(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta)\big{)}$
$\displaystyle+4\bigg{(}z+(2l_{2}l_{3}\cos\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(l_{2}^{2}\cos\theta^{2}\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}\sin\theta$
$\displaystyle+(2l_{1}l_{2}\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}^{2}\cos\theta\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}$
$\displaystyle\big{(}-s_{0}+2(l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-
l_{2}\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\big{)}\bigg{)}$ (A.17)
$\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}2K_{V}(2z(l_{1}-l_{2}\cos\theta)\sin\theta-2z\cos\theta(l_{3}+l_{2}\sin\theta))$
$\displaystyle\big{(}-1+v\big{)}+4\bigg{(}z+(2l_{2}l_{3}\cos\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+(l_{2}^{2}\cos\theta^{2}\sin\theta)/\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}\sin\theta$
$\displaystyle+(2l_{1}l_{2}\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-(l_{2}^{2}\cos\theta\sin\theta^{2})/\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}\bigg{)}\big{(}-s_{0}+s\big{)}\bigg{)},$
(A.18) $\displaystyle\frac{\partial E}{\partial l_{3}}=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}z(l_{1}-l_{2}\cos\theta)\big{(}-1+2z(l_{1}-l_{2}\cos\theta)$
$\displaystyle(l_{3}+l_{2}\sin\theta)\big{)}+4(l_{1}+2\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}})$
$\displaystyle\big{(}-s0+2(l_{1}l_{3}+l_{2}z+2l_{3}\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}\sin\theta\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}}+2l_{1}\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}}$
$\displaystyle-l2\cos\theta\sqrt{z^{2}+l_{2}^{2}\sin\theta^{2}})\big{)}\bigg{)}$
(A.19) $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}4K_{V}z(l_{1}-l_{2}\cos\theta)\big{(}-1+v\big{)}$
$\displaystyle+4(l_{1}+2\sqrt{z^{2}+l_{2}^{2}\cos\theta^{2}})\big{(}-s0+s\big{)}\bigg{)}.$
(A.20)

##  VII Appendix B: Energy functional with scaling factors

###  VII.1 Two dimensions

$\displaystyle e_{S(\text{hex})}$ $\displaystyle=\frac{1}{2}(C_{x}C_{y}\cdot
2l_{2}(l_{1}-l_{2}\cos\theta)\sin\theta-1)^{2}$
$\displaystyle+\frac{1}{2}k_{r}(C_{x}\cdot 2l_{1}$
$\displaystyle+\sqrt{{C_{x}}^{2}\cos{\theta}^{2}+{C_{y}}^{2}\sin{\theta}^{2}}\cdot
4l_{2}-p_{0})^{2}$ (A.21) $\displaystyle e_{S(\text{oct})}$
$\displaystyle=\frac{1}{2}\bigg{(}C_{x}C_{y}\cdot\big{(}l_{1}(l_{1}-2l_{2}\cos\theta)+2l_{1}l_{2}\sin\theta$
$\displaystyle-2{l_{2}}^{2}\cos\theta\sin\theta\big{)}-1\bigg{)}^{2}+\frac{1}{2}k_{r}((C_{x}+C_{y})\cdot
2l_{1}$
$\displaystyle+\sqrt{{C_{x}}^{2}\cos{\theta}^{2}+{C_{y}}^{2}\sin{\theta}^{2}}\cdot
4l_{2}-p_{0})^{2}$ (A.22)

###  VII.2 Three dimensions

Recall we have $e_{T}=\frac{1}{2}(v-1)^{2}+\frac{1}{2}k_{\rho}(s-s_{0})^{2}$,
and where $s,v$ are

$\displaystyle s$ $\displaystyle=2\bigg{(}C_{x}C_{y}\cdot
l_{1}l_{3}+2C_{y}l_{3}\sqrt{{(C_{z}\cdot z)}^{2}+{C_{x}}^{2}\cdot
l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+C_{y}l_{2}\sin\theta\sqrt{{(C_{z}\cdot z)}^{2}+{C_{x}}^{2}\cdot
l_{2}^{2}\cos\theta^{2}}$
$\displaystyle+l_{2}C_{z}z\sqrt{{C_{x}}^{2}\cos\theta^{2}+{C_{y}}^{2}\sin\theta^{2}}$
$\displaystyle+2C_{x}l_{1}\sqrt{{(C_{z}\cdot z)}^{2}+{C_{y}}^{2}\cdot
l_{2}^{2}\sin\theta^{2}}$ $\displaystyle-
C_{x}l_{2}\cos\theta\sqrt{{(C_{z}\cdot z)}^{2}+{C_{y}}^{2}\cdot
l_{2}^{2}\sin\theta^{2}}\bigg{)},$ (A.23) $\displaystyle v$
$\displaystyle=C_{x}C_{y}C_{z}\cdot
2z(l_{1}-l_{2}\cos\theta)(l_{3}+l_{2}\sin\theta).$ (A.24)

## References

  * Kelvin [1887] L. Kelvin, Phil. Mag. 24, 503 (1887).
  * Weaire and Phelan [1994] D. Weaire and R. Phelan, Philosophical Magazine Letters 69, 107 (1994).
  * Bi  _et al._ [2015] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, Nature Physics 11, 1074 (2015).
  * Merkel and Manning [2018] M. Merkel and M. L. Manning, New Journal of Physics 20, 022002 (2018).
  * Zhang and Schwarz [2022] T. Zhang and J. M. Schwarz, Phys. Rev. Res. 4, 043148 (2022).
  * Cantat  _et al._ [2013] I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer, and A. Saint-Jalmes,  _Foams: structure and dynamics_ (OUP Oxford, 2013).
  * Honda  _et al._ [1982] H. Honda, Y. Ogita, S. Higuchi, and K. Kani, Journal of Morphology 174, 25 (1982).
  * Honda and Nagai [2022] H. Honda and T. Nagai,  _Mathematical Models of Cell-Based Morphogenesis_ (Springer, 2022).
  * Farhadifar  _et al._ [2007] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, Current Biology 17, 2095 (2007).
  * Staple  _et al._ [2010] D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, The European Physical Journal E 33, 117 (2010).
  * Mirams  _et al._ [2013] G. R. Mirams, C. J. Arthurs, M. O. Bernabeu, R. Bordas, J. Cooper, A. Corrias, Y. Davit, S.-J. Dunn, A. G. Fletcher, D. G. Harvey, _et al._ , PLoS Computational Biology 9, e1002970 (2013).
  * Sussman [2017] D. M. Sussman, Computer Physics Communications 219, 400 (2017).
  * Moshe  _et al._ [2018] M. Moshe, M. J. Bowick, and M. C. Marchetti, Physical Review Letters 120, 268105 (2018).
  * Sahu  _et al._ [2020] P. Sahu, D. M. Sussman, M. Rübsam, A. F. Mertz, V. Horsley, E. R. Dufresne, C. M. Niessen, M. C. Marchetti, M. L. Manning, and J. M. Schwarz, Soft Matter 16, 3325 (2020).
  * Huang  _et al._ [2022] J. Huang, J. O. Cochran, S. M. Fielding, M. C. Marchetti, and D. Bi, Physical Review Letters 128, 178001 (2022).
  * Chen  _et al._ [2022] Y. Chen, Q. Gao, J. Li, F. Mao, R. Tang, and H. Jiang, Physical Review Letters 128, 018101 (2022).
  * Park  _et al._ [2015] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, J. Notbohm, R. Steward Jr, S. Burger, S. H. Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore, E. Israel, D. A. Weitz, D. J. Tschumperlin, E. Henske, S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen, and J. J. Fredberg, Nature Materials 14, 1040 (2015).
  * Bi  _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Phys. Rev. X 6, 021011 (2016).
  * Barton  _et al._ [2017] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek, PLoS Computational Biology 13, e1005569 (2017).
  * Okuda  _et al._ [2012] S. Okuda, Y. Inoue, M. Eiraku, Y. Sasai, and T. Adachi, Biomechanics and Modeling in Mechanobiology 12 (2012).
  * Okuda and Sato [2022] S. Okuda and K. Sato, Biophysical Journal 121, 1856 (2022).
  * Okuda  _et al._ [2018] S. Okuda, T. Miura, Y. Inoue, T. Adachi, and M. Eiraku, Scientific Reports 8, 2386 (2018).
  * Rozman  _et al._ [2020] J. Rozman, M. Krajnc, and P. Ziherl, Nature Communications 11, 3805 (2020).
  * Zhang  _et al._ [2023] T. Zhang, S. Gupta, M. A. Lancaster, and J. M. Schwarz, bioRxiv , 2023 (2023).
  * Staddon  _et al._ [2023] M. F. Staddon, A. Hernandez, M. J. Bowick, M. Moshe, and M. C. Marchetti, Soft Matter 19, 3080 (2023).
  * Damavandi  _et al._ [2022] O. K. Damavandi, M. L. Manning, and J. M. Schwarz, Europhysics Letters 138, 27001 (2022).
  * Inc. [2022] W. R. Inc., Mathematica, Version 13.2 (2022), champaign, IL, 2022.
  * Alt  _et al._ [2017] S. Alt, P. Ganguly, and G. Salbreux, Philosophical Transactions of the Royal Society B: Biological Sciences 372, 20150520 (2017).
  * Gómez-Gálvez  _et al._ [2018] P. Gómez-Gálvez, P. Vicente-Munuera, A. Tagua, C. Forja, A. M. Castro, M. Letrán, A. Valencia-Expósito, C. Grima, M. Bermúdez-Gallardo, O. Serrano-Pérez-Higueras, F. Cavodeassi, S. Sotillos, M. D. Martín-Bermudo, A. Márquez, J. Buceta, and L. M. Escudero, Nature Communications 9, 2960 (2018).
  * Lemke and Nelson [2021] S. B. Lemke and C. M. Nelson, Current Biology 31, R1098 (2021).
  * Elechiguerra  _et al._ [2006] J. Elechiguerra, J. Reyes-Gasga, and M. Yacaman, J. Mater. Chem. 16, 3906 (2006).
  * Alcinesio  _et al._ [2020] A. Alcinesio, O. J. Meacock, R. G. Allan, C. Monico, V. Restrepo Schild, I. Cazimoglu, M. T. Cornall, R. Krishna Kumar, and H. Bayley, Nature Communications 11 (2020).
  * Cauchy [1813] A. L. Cauchy, J. Ecole Polytechnique 9, 87 (1813).
  * Alexandrov [2005] A. D. Alexandrov,  _Convex polyhedra_ , Vol. 109 (Springer, 2005).
  * Gandikota  _et al._ [2022] M. Gandikota, A. Parker, and J. Schwarz, Physical Review E 106, 055003 (2022).
  * Yan and Bi [2019] L. Yan and D. Bi, Physical Review X 9, 011029 (2019).
  * Smith [1953] C. Smith, Acta Metallurgica 1, 295 (1953).
  * Ranganathan and Lord [2008] S. Ranganathan and E. Lord, Philosophical Magazine Letters 88, 703 (2008).
  * Honda  _et al._ [2004] H. Honda, M. Tanemura, and T. Nagai, Journal of Theoretical Biology 226, 439 (2004).
  * Williams [1968] R. E. Williams, Nature , 276 (1968).
  * Nogucci [2019] H. Nogucci, Biophysics and Physicobiology 16, 9 (2019).
  * Lenne and Trivedi [2022] P.-F. Lenne and V. Trivedi, Nature Communications 13, 2041 (2022).
  * Angelini  _et al._ [2011] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg, and D. A. Weitz, Proceedings of the National Academy of Sciences 108, 4714 (2011).
  * Schötz  _et al._ [2013] E.-M. Schötz, M. Lanio, J. A. Talbot, and M. L. Manning, Journal of The Royal Society Interface 10, 20130726 (2013).
  * Friedl and Gilmour [2009] P. Friedl and D. Gilmour, Nature Reviews Molecular Cell Biology 10, 445 (2009).
  * Sadati  _et al._ [2013] M. Sadati, N. Taheri Qazvini, R. Krishnan, C. Y. Park, and J. J. Fredberg, Differentiation 86, 121 (2013).
  * Kim  _et al._ [2021] S. Kim, M. Pochitaloff, G. A. Stooke-Vaughan, and O. Campás, Nature Physics 17, 859 (2021).
  * Ingber [1997] D. E. Ingber, Annual Review of Physiology 59, 575 (1997).
  * Viens and Brodland [2007] D. Viens and G. W. Brodland, Journal of Biomechanical Engineering 129, 651 (2007).
  * Nnetu  _et al._ [2012] K. D. Nnetu, M. Knorr, J. Käs, and M. Zink, New Journal of Physics 14, 115012 (2012).
  * Mongera  _et al._ [2018] A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shelton, D. A. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio, J. Giammona, and O. Campás, Nature 561, 401 (2018).
  * Oswald  _et al._ [2017] L. Oswald, S. Grosser, D. M. Smith, and J. A. Käs, Journal of Physics D: Applied Physics 50, 483001 (2017).
  * Das  _et al._ [2021] A. Das, S. Sastry, and D. Bi, Phys. Rev. X 11, 041037 (2021).