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

The mechanics of a foam depends on bubble shape, bubble network topology, and the material at hand, be it metallic or polymeric, for example. While the shapes of bubbles are the consequence of minimizing surface area for a given bubble volume in a space-filling packing, if one were to consider biological tissue as a foam-like material, the zoology of observed shapes of cells perhaps motivates different energetic contributions. Building on earlier two-dimensional results, here, we focus on a mean field approach to obtain the elastic moduli for an ordered, three-dimensional vertex model. We use the space-filling shape of a truncated octahedron and an energy functional containing 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. 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. This parameterization effect is more pronounced in three dimensions than in two dimensions given the zoology of shapes that a polyhedron can take on (as compared to a polygon). We also uncover nontrivial dependence of the elastic moduli 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.

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 compatibleincompatible 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 spacefilling arrangement with two types of cells whose surface area was 0.3% smaller than the originally-proposed polyhedron by Kelvin [1,2].
As all objects are deformable to greater or lesser extents, it, therefore, makes sense to explore the shapechanging 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].
FIG. 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], sheardriven solidification [15] and even fracture [16].Several aspects of the theoretical richness have been confirmed experimentally make the vertex model approach a compelling one [14,17], 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.

A. The two-dimensional version
Consider a cross-section of a monolayer of cells.The ith 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 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 = Es K A A 2 0 to arrive at with a = A A0 , p = P √ A0 , and p 0 = P0 √ A0 .In addition, k r stands for the rigidity stiffness and is expressed as K P K A A0 .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.

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 × 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.Moreover, the angle between the two lines l 1 and l 2 only is denoted by θ.The vector of each vertex in relation to the origin O = {0, 0} can be stated as follows (noting that cos θ ≤ 0 as θ is an obtuse Parameterization of the hexagonal cell model embedded in an L x by L y box. angle): The perimeter of the shape is p = 2l 1 + 4l 2 , while the area is equal to −2 1 2 (2l 2 sin θ)(l 2 cos θ) + 2l 1 l 2 sin θ = 2l 2 sin θ(l 1 − l 2 cos θ).Accordingly, the dimensionless energy of a single cell (with a unit area) is ( Note that for a regular hexagon with unit area, p = 2l . The above energy depends on three parameters, l 1 , l 2 , and θ.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 compatibleincompatible 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 ≥ 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.

Octagonal Model
In looking towards the three-dimensional model, it will be useful to study an octogonal cell model in two dimen-sions, 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).Assuming the oc- tagon is symmetric along x, y direction about the origin O = {0, 0} and setting θ 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:

B. 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] Here K V represents the cell volume stiffness and K S the cell surface area stiffness.Since L ∼ V with L as the length unit in the three-dimensional scenario, we can express a dimensionless version of tissue energy with where v, s, s 0 are given by v , and . Finally, . 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.A single polyhedron whose copies fill threedimensional space.
2. The polyhedron must have flat, non-curved surfaces.
3. To ultimately compare with earlier work, each vertex in the tiling has four edges attached to it.

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 [5,20].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].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.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 FIG.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.FIG.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.Sec.II A 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 l ′ 2 2 of the hexagon on a plane can be calculated using vertical coordinate z and parameters l 1 , l 2 , ψ if points A, O are known.Therefore, we can use simple parametrization from the octagon shown below.
As defined in Section II A 1, it is possible to define a hexagon using two edge lengths and one angle.The octagonal cell model from Fig. 5b, however, is charac-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 . The dimensionless shape index s is s ≃ 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 x i + g with g = 0 and i is the index of a vertex of a polygon.To express it simply, F x is used to refer to all coordinate vectors of the vertices for n−gons transformed by the matrix 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 δ = 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 θ degree of freedom, for example, which is the relaxed state.We then apply additional deformation by increasing δ.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 δ, which is proportional to the elastic moduli.Since the dimensionless shear modulus G is defined as G = 1 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 II 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 δ 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 , θ, . . . of p 0 and p 0 + ϵ (ϵ ≃ 0.1) are similar.
• Minimization is allowed within the convex region of the shape ( π 2 ≤ θ ≤ π) as long as it stays within the physically reasonable range (l 1 , l 2 ≥ 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 A 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 I.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 A 1, where we have the function e S , defined by parameters l 1 , l 2 , θ, the continuous arc curvature is calculated and displayed in Fig. 8.
We varied C x and p 0 and held the other parameters fixed (l 8 demonstrates that there is effectively no difference between the two approaches, so we move forward using the discrete method.Note that when performing the computations, we typically acquire approximately 10-20 points close to C x = 1 (for instance C x = 1 ± 0.0005).In some instances, there may be distinct solutions for C x = 1 + δ and C x = 1.For some cases, it is related to extreme values (such as θ = π, π/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(∆) is 4 triangle area product of 3 edges , where∆ 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    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 x within a box (L x , L y , L z ) defined by F x The regular truncated octahedron is linear in the z direction, meaning there are just five values (−z, − z 2 , 0, 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 , θ (recall ϕ = 3π 2 − θ) 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 ∼ 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 × 2 results in 24.Hexagons A and B are the same, and we can denote the threedimensional 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 ) = ( r1+r7 2 , r1+r10 2 , r3 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 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 , θ, z relatively easily.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 = √ 1 ± δ, to measure the bulk modulus.
IV. RESULTS

A. 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.

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 (δ), we write e S (δ) = 1 2 (a(δ) − 1) 2 + 1 2 k r (p(δ) − p 0 ) 2 where δ implies C x → 1 + δ, C y → 1/(1 + δ) for the shear modulus, C x → √ 1 + δ, C y → √ 1 + δ for the bulk modulus, and C x → 1 + δ for Young's modulus.Since C x C y = 1 for the shear modulus, a(δ) = 1 if a = 1 (the initial condition).On the other hand, C x C y = 1 + δ for the bulk modulus and the Young's modulus, so a(δ) ̸ = 1.We can calculate the curvature functions for each case as demonstrated below.Specifically, where the primes denote derivatives with respect to δ.
For the hexagon, because a(δ) is same for both bulk and Young's moduli, we arrive at (a(δ) − 1)a ′′ (δ) = 0 and Finally, we obtain Notice that effects of C x → 1+δ and C x → √ 1 + δ, C y → √ 1 + δ act differently for the cell perimeter.These differences will lead to different curves for the bulk and 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 , θ) allows the system to soften even more so.In addition, relaxing both l 2 and θ 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 ≥ 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 , θ 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 nonaffine/affine variables during relaxation and less restrictive geometry constraints.With a fixed boundary, the box size (L x × 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 θ).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.FIG.10: Mechanical response for hexagon.

Octagonal Model
The next step is to compute the mechanical response of the octagonal cell described in Sec.II A 2. 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 , θ) or (l 2 , θ) 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) ≈ 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).12: Young's modulus for the two cellular models for the relaxed case with (C y , θ) as the minimized parameters.The solid red line denotes p * 0 (6), while the dashed red line denotes p * 0 (8).

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.
• 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 : total length in x,y of octagon) • Revise energy function by employing new variables from e S (l 1 , l 2 , l 3 , θ, k r , p 0 ) to e S (w, h 1 , h 2 , θ, 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 θ, h = l 1 + 2l 2 sin θ, has a singularity at θ = 3π/4 (a = (h+w) 2 cot θ+2hw csc θ 2 2(1+cot θ) 2   and cot( 3π 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.
The corresponding angle ϕ is equal to ϕ = 3π/2 − θ.Fig. 14 illustrates the mechanical response of a reparametrized octagon.If we take θ 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 , θ) is minimized in the relaxed case.What is the mechanism behind this twopeaked structure?If we look at Fig. 16, the first case (C y , θ) 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 θ, 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.

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 FIG.16: (a,b) Optimized shapes for relaxed conditions for target p 0 3 − 4.8.(a) For the relaxed case (C y , θ), 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 (θ), 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 θ values for Fig. 15 for the two different relaxed cases.
work with our initial parameterization of the truncated octahedron with l 1 , l 2 , θ,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 twoparameter minimization.If we set l 1 = l 3 = 1 2•2 1/6 and l 2 = √ 2 2•2 1/6 , K V = 10, K S = 1, and θ = 3π 4 , for the constrained case we obtain 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 |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 For the relaxed cases where there is at least a twoparameter 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.Similar to Sec.IV A 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 θ, 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 III 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 π/2.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 θ, Young's modulus remains high, since decreasing the energy requires more than just changing θ.If we measure Young's modulus for h 2 = 0 using (C y , θ) 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.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 compatibleincompatible 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 ≃ 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 ≥ 5.8, we find the cell is more likely to be a staggered hexagonal cell (for h 2 ̸ = 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.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 [20,39].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 ≃ 5.45477 structures with two sides, as demonstrated in Fig. 22, and ≃ 5.58941 for those with four sides.Assuming a reconnection event from one side in Fig. 22, the value of s is ≃ 5.38742 for a transition point.Note that s ≃ 5.38742 ∼ 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 ̸ = 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.FIG.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 ≃ 5.45477).In order to fill space, it is necessary to replace the gray-colored cells with other shapes.
a(0) ∂ 2 es ∂δ 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: Constrained : lim Cx→1 K(e S (fixed vertex coordinates)) Relaxed : lim Cx→1 K(min [e S (vertex coordinates with a given parameterization)]).

FIG. 8 :
FIG.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.
FIG. 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 × L y ) is represented by the shaded area (b).
FIG.12: Young's modulus for the two cellular models for the relaxed case with (C y , θ) as the minimized parameters.The solid red line denotes p * 0 (6), while the dashed red line denotes p * 0 (8).

FIG. 13 :FIG. 14 :
FIG.13: Schematic of the reparameterized octogonal model using h 1 , h 2 , and w.Box (L x × L y ) is represented by the shaded area.

FIG. 15 :
FIG. 15: Young's modulus as a function of p 0 for the reparameterized octagonal model and for two different relaxed cases.
FIG. 18: Transformable structures of the reparamterized truncated octahedron cell model.
Mechanical response of 3d model.

FIG. 19 :
FIG. 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 (θ) and (C y , θ).The first red line is for a value of s 0 = 5.28, and the second one is for s 0 = 5.8.
FIG. 20: (a) Mechanical response of 3d model for relaxed (θ) and (C y , θ).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 .

FIG. 21 :
FIG. 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.
FIG. 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).

Table I and
II summarize the deformation protocol and control parameters.

TABLE II :
Control parameters

TABLE III :
Dimensionless shape index s 0 for various shapes.