A geometrically controlled rigidity transition in a model for confluent 3D tissues

The origin of rigidity in disordered materials is an outstanding open problem in statistical physics. Previously, a class of 2D cellular models has been shown to undergo a rigidity transition controlled by a mechanical parameter that specifies cell shapes. Here, we generalize this model to 3D and find a rigidity transition that is similarly controlled by the preferred surface area: the model is solid-like below a dimensionless surface area of $s_0^\ast\approx5.413$, and fluid-like above this value. We demonstrate that, unlike jamming in soft spheres, residual stresses are necessary to create rigidity. These stresses occur precisely when cells are unable to obtain their desired geometry, and we conjecture that there is a well-defined minimal surface area possible for disordered cellular structures. We show that the behavior of this minimal surface induces a linear scaling of the shear modulus with the control parameter at the transition point, which is different from the scaling observed in particulate matter. The existence of such a minimal surface may be relevant for biological tissues and foams, and helps explain why cell shapes are a good structural order parameter for rigidity transitions in biological tissues.

The origin of rigidity in disordered materials is an outstanding open problem in statistical physics. Previously, a class of 2D cellular models has been shown to undergo a rigidity transition controlled by a mechanical parameter that specifies cell shapes. Here, we generalize this model to 3D and find a rigidity transition that is similarly controlled by the preferred surface area: the model is solid-like below a dimensionless surface area of s * 0 ≈ 5.413, and fluid-like above this value. We demonstrate that, unlike jamming in soft spheres, residual stresses are necessary to create rigidity. These stresses occur precisely when cells are unable to obtain their desired geometry, and we conjecture that there is a well-defined minimal surface area possible for disordered cellular structures. We show that the behavior of this minimal surface induces a linear scaling of the shear modulus with the control parameter at the transition point, which is different from the scaling observed in particulate matter. The existence of such a minimal surface may be relevant for biological tissues and foams, and helps explain why cell shapes are a good structural order parameter for rigidity transitions in biological tissues.

I. INTRODUCTION
Many biological tissues are confluent, where there are no gaps or overlaps between cells. These tissues are active, disordered, far-from equilibrium materials, sharing similarities with both fiber networks and disordered particulate or glassy matter. Recent experiments have demonstrated that confluent tissues can exhibit a glassy fluid-to-solid transition [1][2][3][4], which likely plays a role in development and disease [5,6]. Therefore, an interesting open question is how a cell's structure and mechanics influence the rigidity of a confluent tissue.
The onset of rigidity in single cells [7] and groups of cells below confluence [8,9] has been fruitfully studied using simple tensegrity and particle models. To understand rigidity in confluent tissues, we study a simple vertex model [10][11][12][13][14][15][16][17] that incorporates the constraints on cell shape imposed by confluence. The original models focused on 2D monolayers of cells, which are described as networks of cellular polygons that tessellate space, where the degrees of freedom are the vertices of the polygons. Cells pay an energetic penalty when their perimeter differs from a preferred value P 0 and their cross-sectional area differs from a preferred area A 0 . Variations include replacing edges and vertices by a Voronoi tessellation of cell centers [15,16,18,19], and including active dynamics [16,18,20]. Although the dynamical rules for the evolution of cellular potts models are quite different, most also use a similar energy functional based on interfacial energies [21,22].
Recently, it was discovered that such models exhibit rigidity or glass-like transitions controlled by the nondimensionalized preferred perimeter p 0 = P 0 / √ A 0 [11,13,14,18,23]. While the hexagonal ground state becomes linearly unstable if p 0 > 3.72 [11,23], there is * mmanning@syr.edu a bona fide second order rigidity transition in disordered structures with a critical point reported at p * 0 = 3.81 [13,14]. Adding active self-propulsion converts this point into a line of glass-like transitions in vertex models [18] and cellular potts models [24]. A priori model predictions for cell shapes have recently been verified in experiments on glassy 2D monolayers [5]. As in many other systems, the critical point controls global mechanical properties and collective motion, which both contribute to biological function.
But what is the origin of this transition, and why is it so robust -occurring in both simulations (vertex, Voronoi, cellular potts) and experiments? In both particulate granular matter and fiber networks, rigidity is fairly well understood [25,26]. For example, frictionless athermal particles rigidify or jam at a critical value of the packing fraction, and the onset of rigidity is predicted by Maxwell's constraint counting rule [27][28][29]. In models for fiber networks, the bond occupation probability controls the transition point also according to Maxwell's criterion [30][31][32]. Moreover, exploring the effect of dimensionality [33] has helped to discriminate between different theories for the origin of rigidity. Therefore it is natural to study vertex models in different dimensions to investigate whether constraint counting explains rigidity in confluent tissues.
Furthermore, although curved 2D cell sheets embedded in 3D space [34][35][36][37][38][39] are an an active area of research, there is surprisingly little work on confluent bulk 3D tissues [40,41], and that work predates the discovery of the rigidity transition in 2D. While there are currently few experimental observations of glassy behavior in fully 3D tissues [3], rapid advances in microscopy and segmentation algorithms [42,43] are currently generating a host of data on 3D bulk tissues that could be used to test model predictions and make connections to embryonic development and disease in vivo.
Here we generalize an existing 2D Voronoi model for isotropic tissues [18] to three dimensions, and find a rigid- ity transition controlled by the dimensionless preferred shape index with an associated structural order parameter. We demonstrate that in contrast to jammed solids, residual stresses that arise when a cell is unable to attain its preferred geometry rigidify the system. Because the mechanisms for rigidity in vertex and particle-based models are distinct, we can identify several specific differences in response and structure that could be used to distinguish between models in real tissues.
In addition, finite-size scaling indicates that the transition occurs at a precise value of the 3D shape index, leading us to conjecture that there is a well-defined minimal surface for cellular structures that controls this rigidity transition. Our work suggests that the origin of rigidity is purely geometric and does not depend on the details of the vertex energy functional. This provides a possible explanation for why observations of rigidity transitions are so robust in cellular systems, and it strongly suggests that our results hold not only for our specific model but for a much broader class of 3D tissue models.

II. MODEL AND METHODS
We describe a three-dimensional confluent tissue by a network of N c cells, where each cell i is represented by a position vector R i (Fig. 1a, inset). Cell shapes are described by the Voronoi tessellation of these cell positions. In a straightforward generalization of 2D Voronoi models [18] to 3D, inter-cellular forces are defined as derivatives on an effective energy functional: where V i and S i correspond to cell volumes and surface areas, respectively. The sum is over all cells of the network, K V > 0 is the volume rigidity, V 0 is the preferred volume, K S > 0 is the surface rigidity, and S 0 is the preferred surface area. Our goal is to simulate this model to determine if there is in fact a rigidity transition in 3D, and to shed light on the origin of rigidity. This is a technically challenging question, as we need to identify local energy minima, which correspond to states in mechanical equilibrium, sampling over instantiations of the disorder. We also need to determine precisely if the shear modulus at the minima are distinct from zero over a wide range of model parameters. To perform these tasks, we need very accurate numerical calculations of the (rather complicated) first and second derivatives of the energy functional given by (1). Therefore, in Appendix A, we present details required to develop analytic expressions for these derivatives, which serves two purposes. First, it allows the most accurate computations, and second, it allows us to compare the analytic and numerical derivatives as a consistency check to validate our code.
We use periodic boundary conditions with a fixed cubic box size of L 3 . As in 2D, the preferred volume V 0 only renormalizes the pressure and does not affect the forces between the cells (Appendix A 1) [15,44]. Therefore, without loss of generality we set V 0 = V = L 3 /N c and non-dimensionalize our model with respect to the length unit V 1/3 and the energy unit K S V 4/3 , leading to the dimensionless energy: The remaining three dimensionless parameters are the (b) (a) preferred shape index s 0 = S 0 / V 2/3 , the relative volume rigidity k V = K V V 2/3 /K S , and the number of cells N c .
Extended details of the initialization and minimization are described in Appendix B. Briefly, unless otherwise noted, we fix N c = 512. We vary k V logarithmically and s 0 linearly between 0 and the ideal gas value 5.82, which is the average cell surface area s = ( i s i )/N c for randomly placed cell centers [45] (Appendix B 1). Reported values are averaged over 100 random initial configurations minimized to a local energy minimum using the BFGS algorithm [46], varying the cell positions r i and the simple shear degree of freedom γ (Appendix B 2).

A. Existence of the rigidity transition
For each relaxed state, we calculate the simple shear modulus g = (d 2 e/dγ 2 )/N c using an analytic expression based on the Hessian matrix that describes the second derivatives of the energy functional (Appendices A 5-A 6). We find that there is a continuous transition in g at a preferred relative surface area of s 0 = s * 0 ≈ 5.4 with a solid regime (g > 0) for s 0 below this transition point and a fluid regime (g = 0) above it (Fig. 1a). Neither the transition point s * 0 nor shear modulus g depend strongly on k V . We also studied the precise distribution of transition points p(s * 0 ) (Fig. 1b) (Appendix C 2). We found a finite variance of this distribution, which appears to be only due to finite-size effects (Fig. 1b inset) and identify an extrapolated transition point of s * 0 ≈ 5.413 in the limit N c → ∞ (Appendix C 3).

B. Mechanism that generates rigidity
To better understand the origin of the rigidity transition, we begin with simple constraint counting which has successfully been used to predict rigidity in many systems, including spring networks [32,47]. Our system has 3N c + 1 degrees of freedom, 3 for each cell position plus the shear degree of freedom. If each of the 2N c terms in Eq. (2) is regarded as a generalized spring and thus contributes a single constraint, the system would be highly under-constrained with at least N + 1 zero modes for all values of s 0 , inconsistent with our observations.
To understand how this seemingly under-constrained system can rigidify, we define the extended HessianD pq : where (z p ) is a (3N c + 1)-dimensional vector of the cell positions and the shear degree of freedom, and the sum is over all cells i. The full analytic expression of the extended Hessian is derived in Appendix A 5. The constraints imposed by the 2N c generalized springs discussed above correspond to the first two terms in Eq. (3), while the last two terms correspond to residual stresses, which are known to rigidify otherwise under-constrained systems [7,28,47,48]. In our system, residual stresses are the 2N c surface tensions 2(s i − s 0 ) and pressures To study whether residual stresses are both necessary and sufficient to generate rigidity, we plot a twodimensional histogram categorizing energy-minimized states with respect to both their shear modulus g and the maximal surface tension magnitude of all cells within the configuration, Fig. 2a. A similar histogram also holds for the pressure (Fig. 7a, Appendix C 5). The fact that the upper left quadrant is largely devoid of configurations shows that residual stresses are necessary for rigidity. The few exceptions are close to both cutoff values and likely due to imperfect minimization.
We independently verified this by computing the overlap between the shear degree of freedom and the infinitesimal zero modes, which is always finite (Appendix C 4). This shows that the system would indeed have no resistance to shear in the absence of residual stresses.
The lower right quadrant in Fig. 2a is also devoid of configurations, except for a handful of systems (about 1 in 1000). These are floppy except for a single rigid cell of a specific polyhedron type with anomalously few neighbors -a 10-sided truncated square trapezohedron (Fig. 7b). Hence, up to these very few exceptions, the rigidity transition is created by the onset of residual stresses, which typically occurs in every cell simultaneously (Appendix C 6).
So what controls the onset of residual stresses in our system? In the floppy regime, when all surface tensions 2(s i − s 0 ) and pressures 2k V (1 − v i ) are zero, cells exactly attain their desired shapes so that the average observed shape index s equals the preferred shape index s 0 , as shown in Fig. 2b and there is zero standard deviation in cell surfaces and volumes (σ /N c and σ v defined analogously, see Fig. 9). We call this geometric compatibility. Because the energy functional Eq. (2) drives residual stresses towards zero, the fact that residual stresses are non-zero in the solid indicates that no geometrically compatible state is reachable by standard energy minimization. This in turn suggests that the transition point s * 0 is determined by a purely geometric criterion: it corresponds to a local minimum in the average surface area s under the constraint that there are no surface and volume fluctuations σ s = σ v = 0 as in the fluid.
Thus, based on our finite-size scaling analysis (Fig. 1b), we conjecture that in the thermodynamic limit there is a minimum possible value of the average surface area s for disordered Voronoi tessellations with σ s = σ v = 0, which is given by s * 0 ≈ 5.413. This conjecture is reminiscent of those for jammed packings of particles, where the distribution of jamming packing fractions approaches a narrowly peaked function in the limit of large system sizes [49][50][51]. It has been suggested that this packing fraction (about 64% in 3D) can be defined as the maximally disordered rigid state of spheres [52], or a divergence in the rate at which accessible disordered states disappear [26].
As in theories for jamming, we must specify that we are restricting ourselves to random tessellations, as s can become smaller than s * 0 for ordered states: the Voronoi packing corresponding to the ordered Kelvin structure has s ≈ 5.315 (Appendix C 1). Identifying which configurations lie on the "disordered" branch is of course difficult. In jamming, it is clear that different protocols for generating packings on the "disordered" branch generate distinct critical packing fractions, although the numbers are all still remarkably close [50]; and we find something similar here. For example, when using a conjugated gradient minimizer instead of BFGS, we have found a disordered transition point slightly shifted down by 0.01. A weak dependence on the minimization protocol was recently also observed for the 2D Voronoi model [19].
C. Universal behavior and scaling of the shear modulus Fig. 2b also shows that in the solid regime, the average surface area s decreases below the transition point s * 0 ≈ 5.413, as the surface and volume fluctuations rise away from zero (Fig. 9). This is not surprising, as Eq. (2) can be rewritten: In the floppy regime all three terms can be zero simultaneously because geometrically compatible states are attainable, but in the solid this is not possible and we find that the minimal surface is a function of the variances: s = s min (σ s , σ v ) [53]. In the vicinity of the transition point, s min (σ s , σ v ) is surprisingly simple and universal: we find numerically that with a s ≈ 2.0 and a v ≈ 6.8, as shown by the data collapse in Fig. 3a. This result is independent of the parameters of the energy functional s 0 and k V . Strikingly, the collapse also shows that the parameters a s and a v are largely independent of the random initial conditions, suggesting that they are universal geometrical properties of disordered 3D Voronoi packings. Possibly, a relation like Eq. (5) with universal coefficients may also hold for general 3D cellular packings. This linear scaling in the minimal surface s min explains the behavior of the shear modulus close to the transition point, which scales linearly with the distance from the transition point, δs 0 = s * 0 − s 0 (Fig. 3b). To understand this, we start from the formula Eq. (A55) for the shear modulus g (Appendix A 7), which we restate here: If and only if the extended HessianD pq has a zero mode with a nonzero shear component, then the shear modulus g is zero. Otherwise, the shear modulus can be computed as: Here the sum is over all strictly positive eigenvaluesω 2 m of the extended HessianD pq , andū m γ are the shear components of the associated eigenvectors. First, we demonstrate that this expression gives the correct answer in the fluid regime, where there are no residual stresses, and the shear modulus should be zero. Without residual stresses there are at least (N c + 1) eigenmodes ofD pq with eigenvalue zero, corresponding to the (N c + 1) zero modes obtained by the naive constraint counting. As shown in Appendix C 4, at least one of these zero modes has a nonzero shear component, which is why the system is floppy with g = 0.
As discussed above, in the solid vicinity of the transition point, the cellular surface tensions 2(s i − s 0 ) and pressures 2k V (1 − v i ) appear and rigidify the system. In particular, while the average pressure is zero because v = 1, Eqs. (4) and (5) imply that the average surface tension scales linearly with the distance to the transition point: 2( s − s 0 ) = 2δs 0 /(1 + a 2 s + a 2 v /k V ) (Appendix C 8). Moreover, it turns out that the surface tensions and pressures of each individual cell are proportional to δs 0 [53]. As a consequence, many of the (N c + 1) eigenvalues that were zero in the fluid increase by an amount proportional to δs 0 , and these small eigenvalues then dominate the formula for the shear modulus g according to Eq. (6). Hence, the shear modulus scales linearly with the distance to the transition point g ∼ δs 0 . We show this numerically by plotting a rescaled shear modulus g = (1 + a 2 s + a 2 v /k V )g over δs 0 . The fact that this quantity collapses for different values for k V demonstrates that the average surface tension 2( s − s 0 ) is a dominant contribution to shear modulus. This is completely different from what is observed in particulate matter, where constraints are added as particles come into contact, generating a shear modulus that scales as the square root of the dimensionless control parameter (the packing fraction).

IV. DISCUSSION AND CONCLUSIONS
To our knowledge, this is one of the first systematic numerical studies of the mechanics of confluent bulk biological tissues [40,41]. Our 3D Voronoi-based model, which is a straightforward generalization of successful 2D models for epithelial sheets [18], exhibits a rigidity transition when the dimensionless preferred cell surface area is s * 0 ≈ 5.413. The transition is accompanied by a structural order parameter, the observed average cell surface area s .
In addition, we finally have a numerics-backed conjecture that explains the origin of rigidity in vertex-like models, as well as an explanation for the robustness of the structural order parameter. In contrast to jamming in frictionless spheres, where an unjammed system acquires more constraints as the packing fraction increases and spheres touch, cells in vertex models always have the same number of surface and volume constraints independent of model parameters.
Instead, rigidity is created by residual stresses which are due to geometric incompatibility. Specifically, we conjecture that there is a minimal surface area possible for disordered cellular structures under the constraint that each cell has an identical surface area and volume. These constraints arise naturally at the transition point, because in the fluid phase each cell can exactly attain its desired shape and so there are no fluctuations in those quantities. In the solid, the system would like to have a surface smaller than the minimal one, and so it is stuck there, although small fluctuations in the surface area or volume of each cell can reduce the surface area below the transition value. We believe that this is also the mechanism controlling the transition in 2D vertex models [18,19].
Although we have focused here on the relevance to biological tissues, it is likely that rigidity transitions driven by residual stresses may appear in a wide variety of other models for disparate physical phenomena [47], and the collective nature of the transition could lead to glassy dynamics that is strikingly different from that seen in particulate matter.
Our finding that a purely geometric quantity governs the onset of rigidity in a vertex model may also help explain why the correlation between cell shape and rigidity is so robust [14,18] and even holds in experiments on biological tissues [5]. Our results suggest that any system where cellular structures minimize their surface area and suppress fluctuations to their surface areas and volumes could rigidify when their shape index s drops below 5.413, since there are no available states below this value.
This leads to an immediate and testable prediction for experiments on 3D cell aggregates including embryonic zebrafish cells [3] and human breast cancer cells [6], as well as 3D tissue explants from vertebrate embryos [54]. Our model predicts that in the solid s ≈ 5.4, and s should rise away from that value as tissue fluidizes.
Although a similar prediction was successful in 2D [5], additional mechanical interactions not yet included in our model may be more important in 3D than 2D. Therefore, we think of the model presented here as a useful null hypothesis for establishing a relationship between tissue structure and tissue mechanics in 3D; many different perturbations can and should be studied, as they may alter the transition. For example, we will discuss the influence of cell motility and persistence on the 3D model elsewhere [55], but just as in 2D we find that cell shape is still an excellent predictor of tissue rheology even in the presence of motile forces. Additional useful extensions to the model could account for cell division (which may fluidize the solid phase [56,57]) and polydispersity in the preferred surface areas and volumes.
Another possible perturbation is nuclear rigidity, as an alternate explanation for tissue rigidification is that the nuclei jam. Interestingly, the shape index associated with Voronoi cells for particulate matter at jamming is 5.38 [58], which very close but distinct from the transition we observe here.
Importantly, because the nature of the transitions in vertex and particle models are fundamentally different, we can identify several observables that should distinguish between these two mechanisms. In vertex models, fluctuations in Voronoi volumes and surface areas are minimal in the fluid phase and grow in the solid phase, whereas for particulate matter these fluctuations are large in the fluid phase and get smaller as one approaches the solid phase. As discussed in Section III C we also predict different scaling laws for the shear modulus: in jammed particulate matter the modulus scales with the square root of the control parameter (the pack-ing fraction) [59], while in the vertex model the modulus scales linearly with the control parameter (the target shape index). Finally, there is no bulk modulus in the fluid phase of jammed particles, but always a significant bulk modulus in vertex models, even in the fluid regime.
Although we present numerical evidence for the origin of the transition, it would be interesting to try to develop geometric arguments that predict the value of the dimensionless disordered minimal surface s * 0 ≈ 5.413 as well as analytic arguments that explain why the residual stresses are both necessary and sufficient to create rigidity. Recent work by Moshe and collaborators on lattice-based structures develops a nice framework for studying this problem [60].
It would also be interesting to test the predictions of our model in passive cellular materials like foams or biomimetic cellular materials [61][62][63]. For example, random foams are deep in the solid phase (s 0 → −∞ and k V /|s 0 | → ∞), and since they seek to minimize their surface area under an evolution of cell volumes driven by mean curvature flow, they should also be governed by our minimal surface hypothesis. For example, Ref. [62] reports values for the average shape index s , but not the fluctuations σ s . It would be interesting to revisit these results and determine if the surfaces are related to our minimal family s min (σ s , σ v ).
Finally, our findings could also be used to create artificial cellular materials that may be hyperuniform [64], and can transition between solid-like and fluid-like behavior depending on some microscopic parameter, such as an effective surface tension tuned using chemicals, light, or magnetic fields [63,65], or on macroscopic parameters like the overall volume or the pressure, which change the average cell volume. Here we show that V 0 only offsets the pressure but does not affect intercellular forces and the shear modulus. To this end, we transform Eq. (1) into: Here V tot is the total volume of the system. Note that V 0 only appears in the last term where it only couples to the total system volume. As a consequence, V 0 does not affect any forces except that it offsets the pressure of the system. Hence, we just fixed V 0 = V such that the last term disappears. Note that an analogous argument also holds in 2D [15,44].

Periodic boundary conditions
To implement the periodic boundary conditions in a clean way, we explicitly expressed the total energy e in terms of distance vectors r ij between neighboring cell positions i and j instead of absolute cell positions (see next section, Section A 3). Separately, these distance vectors r ij are expressed in terms of absolute cell positions r i and r j . However, the r ij also depend on the periodic boundary conditions, because cells on opposing sides of the periodic box are neighbors of each other. To implement this, each cell neighbor pair (i, j) is given a integer "periodicity" vector q ij . It is defined such that q ij = 0 whenever one can get from position of cell i to position of cell j without crossing a face of the periodic box. If one needs to cross for instance the upper/lower face of the periodic box once from below in going from i to j, then q z ij = +1, and analogous for the other components of q ij (see also the appendix in [66] for an extensive explanation of the analogous 2D case). The distance vectors r ij = (r x ij , r y ij , r z ij ) then depend as follows on the cell position and the dimensionless periodic box dimensions l x , l y , l z : We are also interested in the effect of simple shear. We thus allow for "skewed" periodic boundary conditions and characterize simple shear by the shear variable γ. It is implemented by modifying the above relations as follows: Note that this relation between cell distance vectors r ij and cell positions r i is the only place where the periodic boundary conditions appear. This allows to completely separate the physics from the boundary conditions in the implementation.

Energy
Here, we express the total energy e in terms of the distance vectors r ij . It is given by the sum of all cell energies e i : with Volume and surface of a cell i are given by: Here, both sums are over all cells n neighboring cell i.
The vector a in denotes the oriented area of the polygonal interface between cells i and n pointing orthogonally towards n. The dot in the volume formula denotes the scalar product and the vertical bars in the surface formula denote the norm of the vector. We thus need to express the oriented area a in of the face between cells i and n in terms of the distance vectors r ij . To this end, we first express a in in terms of the positions h in,m of the vertices that define the polygonal face: being the vertex position relative to the position of cell i. In Eq. (A8), the sum is over all N in vertices of the face, sorted in counter-clockwise order as seen from cell n (right hand rule with the thumb pointing from cell i to cell n). The cross × denotes the vector product and h in,Nin+1 ≡ h in,1 . It remains to compute the position of a vertex that has the four abutting cells i, j, l, p, which we denote by h ijlp , relative to the position of cell i: ∆h ijlp = h ijlp − r i . The relative vertex position in terms of the distance vectors r ij , r il , and r ip is: where That Eqs. (A10)-(A12) yield the right vertex position can be seen as follows. The vertex position h ijlp can be regarded as the intersection of the three planes that respectively orthogonally bisect the lines between the positions of cells i and j, i and l, and i and p. Any point x on each of these planes respectively fulfills

Forces
To compute the first derivative of the total energy e, we successively use the chain rule of differentiation. We first compute the derivative of the cell energy e i with respect to a distance vector r ij . We have: The derivatives of volume and surface are: Here, both sums are over all neighbors k of cell i. Moreover, here and in the following, Greek letters represent dimension indices: α, β, · · · ∈ {x, y, z} and we use Einstein convention, i.e. summation over same indices is implied. The derivative of the oriented area is: where ε βγδ denotes the Levi-Civita Symbol. It remains to compute the derivative of the position of a vertex abutting cells i, j, l, p relative to the position of cell i, ∆h ijlp , with respect to a distance vector r ij : We now know the derivative of a cell energy e i with respect to the distance to any of its neighbors j, ∂e i /∂r ij . Using the chain rule, the derivative of e i with respect to neighbor j's position is and the derivative of the energy of cell i with respect to its own position is Here, the sum is over all neighbors j of cell i. The total force on a cell i can be written as f i = − k ∂e k /∂r i , where the sum is over all cells k of the network. Inserting Eqs. (A23)-(A24), this becomes where both sums are over all cells k of the network and the second sum is over all neighbors j of cell i.

Hessian matrix
We start by deriving the derivative of the energy of a cell i with respect to the cell distance vectors r ij where cell j is among the neighbors of cell i.
(A26) The first derivatives of volume and surface are documented in the previous section. The second derivatives are: The second derivative of the relative vertex position ∆h ijlp is: where cell k is one of j, l, p. For k = j, we obtain for the second derivatives of Z ijlp and H ijlp : where δ αβ is the Kronecker symbol. For k = l, the derivatives are: Finally, for k = p, they are: The second derivative of the total energy e with respect to the absolute cell positions r i is the Hessian matrix of the system: Here, the first sum is over all cells i in the network, the second sum is over all neighbors l of cell j, the third sum is over all neighbors m of cell k, and the fourth double sum is over all neighbors l, m of cell j.
To compute the shear modulus, we also need the second derivatives of the energy with respect to the simple shear variable γ. We obtain the derivatives involving simple shear using the chain rule and Eqs. (A3)-(A4): where the sum runs over all cells i and all combination of neighbors j, k and l y is the dimensionless length of the box in y direction. For the mixed derivative involving the simple shear, we obtain: where the first sum is over all cells i of the network and all neighbors k of i, and the second sum is over all combination of neighbors l, k of cell j.

Computation of the shear modulus
In Section A 3 we expressed the energy e of the system only in terms of the distance vectors of neighboring cells, r ij , and in Section A 2 we expressed these distance vectors in terms of the absolute position vectors r i , and the simple shear variable γ. Thus, the energy can be expressed as e = e({r i }, γ).
To introduce the long-time simple shear modulus g, we define e min (γ) as the minimum of e = e({r i }, γ) for fixed γ: We denote the set of cell positions with minimal e for given γ by r min i (γ). Then we define the simple shear modulus as: where N c is the number of cells and thus the volume of the periodic box in dimensionless units.
The simple shear modulus is directly related to the Hessian matrix of the system. To derive this relation, we evaluate the derivative in Eq. (A42): The second term in the parenthesis can be transformed as follows: To simplify the sum over m in Eq. (A44), we distinguish between zero modes, for which ω m = 0 and non-zero modes with ω m > 0. For zero modes m, the second term in Eq. (A47) has to be zero. Because this is the same term as the first factor in the sum over m in Eq. (A44), zero modes do not contribute to this sum. For the nonzero modes, we can substitute Eq. (A47) into Eq. (A44) and obtain from Eq. (A43): In this equation, the sum over m excludes zero modes of the Hessian D jα,kβ .

Alternative computation of the shear modulus using the extended Hessian
Alternatively to Eq. (A48), the shear modulus can also be computed from directly the eigen spectrum of the extended HessianD where (z p ) is a (3N c + 1)-dimensional vector comprising all cell positions and the shear degree of freedom γ: (z p ) = (r 1 , . . . , r Nc , γ). The eigenvalues ofD pq areω 2 m and the normalized eigenvectors areū m p : Note that for shear-stabilized states, the extended Hes-sianD pq is non-negative.
To derive an alternative formula for the shear modulus, we combine Eqs. (A43) and (A46) into: where we sum again only over non-zero modes m, becausė z min γ = 1. Insertion of Eq. (A52) finally yields: Here, the sum is over all non-zero modes m.    [67]). Cell surface areas and volumes, forces, and the shear modulus were then computed as described in Sections A 2-A 6. In particular, shear moduli g were computed using a cutoff value of 10 −14 below which eigenvalues of the Hessian were regarded as zero modes and thus disregarded for the sum to compute g. We checked that g was largely independent of this cutoff over a range of cutoff values. To diagonalize the Hessian, we used the Eigen3 library (version 3.2.7, [68]).
We studied the cases k V = 10 −3 , 10 −2 , . . . , 10 3 and s 0 in the range [0, 5.9] in steps of 0.1 and in the range [5.35, 5.45] in steps of 10 −3 . For each parameter pair (s 0 , k V ), we ran 100 minimizations each of which was initialized with a different set of random cell positions. The system size was N c = 512 unless stated otherwise.

Numerical energy minimization
To minimize the energy of the system, we used the BFGS2 multidimensional minimization routine of the GNU Scientific Library (GSL, version 2.1, [69]) [46]. The parameters used for one GSL minimization are listed in Table I, where N dof is the number of degrees of freedom varied during the minimization. We tested that different values for the individual parameters did not improve the minimization, i.e. the norm of the total force vector after the minimization was not smaller for different parameter values.
Often, the GSL library could not further minimize the energy and did not reach the total force cutoff listed in Table I. In these cases, we tested whether the total force norm was at least below a cutoff of C = 10 −7 √ N dof . If that was not the case, we started another GSL minimization starting with the last set of cell positions. We repeated GSL minimizations until the cutoff C was reached or 10 GSL minimizations had been performed.
To obtain shear-stabilized force-balanced states, we first ran up to 10 GSL minimizations varying all N c cell positions r i , i.e. N dof = 3N c . Afterwards, to shearstabilize the system, we included the shear degree of freedom γ into the minimization, simultaneously varying N dof = 3N c + 1 degrees of freedom, running again up to 10 GSL minimizations. We discarded all simulation runs that had a total force norm larger than 10 −4 after the minimization procedure or a negative shear modulus smaller than g < −10 −5 . An exception are Figs. 1a and 2b, where we needed to increase the force cutoff to 10 −3 . This is because for large k V deep in the solid regime, we couldn't minimize the total force below 10 −4 in many cases.

Probing the solid vicinity of the transition point
We performed dedicated simulations to explore the solid vicinity of the transition point. To this end, we first created an energy-minimized at its transition point and then decreased s 0 using exponentially increasing steps.
To tune a configuration right at its transition point, we used bisection on the s 0 parameter with the initial left and right bracket values of s 0 = 5.38 and s 0 = 5.44. In each bisection step, the s 0 value is set to the average of the current bracket values and the energy is minimized as described in the previous section. If the new state is solid, it is kept for the next minimization, and the next left bracket value is set to the current s 0 value. However, if the state is fluid, the system is reverted to the last solid state (corresponding to the left bracket value) and the next right bracket value is set to the current s 0 value. We choose to keep the solid but not the fluid states to reduce the probability to switch the "inherent state" during the bisection. As another measure to avoid switching the "inherent state", we included a check verifying that whenever the energy-minimized state at the current s 0 value is solid, the system at the right bracket is still fluid if we start the minimization from the energy-minimized state at the current s 0 value. We performed 13 such bisection steps, and a configuration was deemed solid if its shear modulus was larger than 10 −7 .
Once we obtained an energy-minimized state at its rigidity transition, we explored the solid regime by iteratively reducing s 0 in exponentially growing steps, each time minimizing the energy as described in the previous section. As a final measure to exclude simulation runs where the "inherent state" changed, we computed the mean-squared deviations R 2 of the cell positions r i at some point s 0 in the solid regime from the cell positions at the transition point s * 0 (corrected by any overall translation). To exclude simulations where the "inhered state" changed, we made sure R 2 was not too large in our simulations. More precisely, we excluded simulations with R 2 > 100δs 2 0 for any of the s 0 probed when exploring the solid vicinity of s * 0 .
Type of cellular polyhedron Position Appendix C: Additional numerical analysis of the 3D Voronoi model

Periodic Voronoi packings
The average cell surface areas s for Kelvin ( s ≈ 5.306) and Weaire-Phelan packings ( s ≈ 5.288) are well-known. However, these figures correspond to configurations with curved cell outlines.
To numerically obtain the s values for the Voronoi packings corresponding to the Kelvin and Weaire-Phelan structures, we proceeded as follows. For the Kelvin Voronoi structure, we prepared a cubic periodic box with side length 2 1/3 with two cells at r 1 = (0, 0, 0) and r 2 = 2 −2/3 (1, 1, 1) and minimized the total surface area. This initial configuration already corresponded to a minimal surface area with s ≈ 5.315 and equal volume of both cells.
For the Weaire-Phelan Voronoi structure, we prepared a state with a cubic periodic box with side length of 2 and 8 cells at the positions listed in Table II. Without constraining the volume, the surface was already minimal for these initial positions with s ≈ 5.295. Note however that the volumes of the cells were slightly different. All 14-faced cells had v i ≈ 1.0078 while the two 12-faced cells had v i ≈ 0.9765. We tried to constrain all cell volumes to be equal but did not succeed. We thus compare our result for the disordered transition point only to the Kelvin structure in the main text.

Determination of the precise transition point
In order to determine the precise transition point, we first quantified the fraction F of rigid networks for different values of the preferred surface area s 0 . In the limit of an infinite number of simulation runs, the function F (s 0 ) can be regarded as the integrated probability distribution of transition points P (s * 0 ): and (s * 0 ) 2 = ∞ −∞ (s * 0 ) 2 P (s * 0 ) ds * 0 can be evaluated using partial integration. Assuming that the transition point is always non-negative, we obtain: where σ(s * 0 ) denotes the standard deviation of the distribution of transition points.
We defined an energy-minimized network as rigid whenever its shear modulus g was larger than a cutoff value. To determine this cutoff value for a given k V , we plotted the shear modulus g for each simulation run over s 0 (Fig. 4a). In such plots, we find a clear separation of at least a decade between two clusters of networks. We interpret all networks belonging to the respective lower cluster as non-rigid and those belonging to the upper cluster as rigid. For k V = 10, our chosen cutoff is indicated by a magenta dashed line in Fig. 4a. For the cases k V = 0.001 and k V = 0.01, both clusters were not clearly enough separated to define a sensible cutoff value. We have thus ignored these parameter values throughout this article.
The resulting values obtained for average s * 0 and standard deviation σ(s * 0 ) of the transition points are listed in Fig. 4b for the different values of k V . The values for the average transition point range from s * 0 ≈ 5.414 to s * 0 ≈ 5.416. When respectively varying the cutoff on g between the two clusters, the value of s * 0 varied by up to ∼ 10 −3 . Thus, the average transition points for the different values of k V are not significantly different from each other.

Finite-size scaling
We studied the behavior of the transition point distribution for varying system size. However, to compute the shear modulus for very large systems, we would have needed to diagonalize very big Hessians of size N dof × N dof , which quickly exhausted the memory of our machines. We thus chose an alternative approach, where we first realized that the shear modulus correlated very well with the median surface tension in the system (Fig. 5a). Note in particular that this is the case even for the exceptions visible in Fig. 2a and discussed in Section C 5 a below. The cutoff for the median, 10 −6 , was chosen such that for each system size N c it clearly separated the two clusters that appear in the histograms showing the distribution of medians depending on s 0 (Fig. 5b).
Using this median cutoff, we extracted the functions F (s 0 ) for each N c for k V = 10. Using Eqs. (C2) and (C3), we then computed the average (Fig. 5c) and variance (Fig. 1b inset) of the transition point distribution. For the average, we find a scaling exponent of −0.9 ± 0.4 and a limit value of s * 0 (N c → ∞) = 5.413 ± 0.001. For the variance, we find a scaling exponent of −0.90 ± 0.04.
In Fig. 1b, the transition point distributions have been computed based on the F (s 0 ) function. To numerically compute the derivative while suppressing noise, for each system size N c , we convoluted F with the derivative of a Gaussian with a standard deviation of σ(s * 0 )/2.

4.
Residual stresses were necessary to create rigidity.
In the main text, we discuss the full extended Hessian D pq and find that the system is rigid below the transition point s * 0 . To numerically check whether residual stresses are necessary to create such rigidity, we study an unstressed extended HessianD pq where the terms that depend on the residual stresses have been removed: To distinguish whether the overlaps that we measure in (a) are finite because of the physics or because of numerical noise, we also plotted the not shear-stabilized case in (b). We expect the overlap in this case to be zero in the solid regime and indeed, we measure values that are several decades smaller than in the shear-stabilized case (a). The color of the dots represents the norm of the (3Nc +1)-dimensional (a) or (3Nc)dimensional (b) total force vector after the minimization.
Here (z p ) is a (3N c + 1)-dimensional vector of the cell positions and the shear degree of freedom. D pq corresponds to the Hessian of a system with a modified energy functional where we set the preferred surface area of each cell to its actual surface area and the preferred volume of each cell to its actual volume. This modified system thus corresponds to the original system without residual stresses (cf. Eq. (3) in the main text).
If (and only if)D pq has a zero mode with a nonzero shear component, then the shear modulus of the modified system is zero (Appendix A 7). Therefore, if we can show thatD pq has such zero modes in the solid phase, we can conclude that residual stresses are necessary to rigidify the original systemD pq .
The zero modes ofD pq correspond to the infinitesimal zero modes of the system. Infinitesimal zero modes are (3N c + 1)-dimensional vectors w in the kernel of the socalled compatibility matrix C [32], i.e. C · w = 0. For our system, C is a (2N c ) × (3N c + 1) matrix defined as: (C5) Note that C is the same for both original and modified system. Hence, an infinitesimal zero mode corresponds to a collective change of cell center positions {r i } and shear variable γ that leaves all cell surfaces and volumes constant to linear order.
The compatibility matrix has a given number N 0 of independent infinitesimal zero modes. More precisely, this means that there is a set of N 0 infinitesimal zero modes w q with q = 1, . . . , N 0 and C · w q = 0, which are orthonormal: w p · w q = δ pq . We extracted such an orthonormal set of zero modes w q using singular value decomposition of C.
As a tool to study the shear components w q γ of all infinitesimal zero modes w q at once, we define the overlap r of the kernel of C with the shear degree of freedom by: By definition, 0 ≤ r ≤ 1. In particular, r is nonzero if and only if there is an infinitesimal zero mode that has a finite shear component. In other words, r is nonzero if and only if the shear modulus of the modified system is zero. In Fig. 6a, we plot the overlap r depending on the preferred cell surface area s 0 for k V = 10, where each dot represents an energy-minimized and shear-stabilized state. We find that the overlap is finite and on the order of r ∼ 10 −3 both in the fluid regime (s 0 > s * 0 ) and in the solid regime (s 0 < s * 0 ), which holds independent of the value of k V (data not shown). This shows that the modified system always has zero shear modulus and thus residual stresses are indeed necessary to rigidify the original system.
Note that close to the transition point s * 0 ≈ 5.4, the overlap is occasionally smaller than 10 −3 . However, this occurs only in a fraction of the cases at s * 0 . Moreover, the overlap was always larger than 10 −5 . To ensure that the observed values of r ∈ [10 −5 , 10 −3 ] in shear-stabilized systems are demonstrably different from zero, we measure the numerical noise in a system where we know the overlap r should be zero. We expect that configurations that have not been shear-stabilized will generically have a finite shear stress in the solid phase, and as we discuss below, a finite shear stress induces a zero overlap r. Therefore, we plot the overlap for non-shear-stabilized configurations with k V = 10, as shown in Fig. 6b. We find that the measured overlap in the solid phase is always below r ∼ 10 −7 and is typically on the order of r ∼ 10 −14 , which is significantly smaller than r ∈ [10 −5 , 10 −3 ]. This supports our interpretation that the overlap is indeed finite in the shear-stabilized systems.
To see why the overlap is zero for non-shear-stabilized states with finite shear stress σ xy , we first consider the force balance condition combined with the definition for the shear stress: Here, δ denotes the Kronecker delta and we use dimensionless units so that N c corresponds to the system volume. Using chain rule on the left-hand side and using vector notation, this transforms into: Here,γ is the (3N c + 1)-dimensional vectorγ = (0, . . . , 0, 1). For σ xy = 0, we obtain an expression for γ from Eq. (C8), and insertion into w q γ =γ · w q yields indeed w q γ = 0, because C · w q = 0. Thus, the overlap is generally zero in the not shear-stabilized case: r = 0. Intuitively, in order for the system to support a finite shear stress, there can be no infinitesimal shear mode that involves the shear degree of freedom, i.e. there can be no collective displacement that leaves all surfaces and volumes constant to linear order while also shearing the system. 5. Residual stresses were in most cases sufficient to create rigidity.
In the previous section, we have shown that residual stresses were necessary to create rigidity in all our simulations. Here, we provide evidence that, at least in the vast majority of simulations, they were also sufficient to create rigidity.
To relate residual stresses to rigidity, we sorted all energy-minimized configurations into two-dimensional histograms with respect to the shear modulus g and with respect to the maximal surface tension magnitude max 2|s i − s 0 | (Fig. 2a in main text) and maximal pressure magnitude max 2k V |v i − 1| (Fig. 7a). In addition, we marked in each case the cutoff defined in Section C 2 below which we interpret a configuration as fluid by a horizontal magenta dashed line (compare Fig. 4).
In both plots, we find one cluster that clearly lies in the fluid regime and several clusters appearing in the solid regime. We think that the appearance of more than  7. (a) Relation between rigidity and the existence of residual stresses: Two-dimensional histogram of shear modulus g and maximal magnitude of cell pressure max 2kV |vi − 1| for kV = 1 (cf. Fig. 2a; plots for different values of kV are similar). Note that all values of g that are smaller than 10 −10 are mapped onto the g = 10 −10 line, and similarly for max 2kV |vi − 1|. (b) Fluid configurations in which residual stresses occur typically contain a truncated square trapezohedron (TST). A TST consists of two opposing squares, each of which is connected to four irregular pentagons, totaling to eight pentagons with minimal surface-to-volume ratio sTST ≈ 5.4436. one solid cluster is merely due to our sampling of the s 0 values: Close to the transition point, we varied s 0 in small steps of 10 −3 while some distance ∆s 0 ∼ 10 −1 away from the transition point, we varied s 0 in steps of 0.1 (compare Section B 2).
We realized that the fluid cluster in both plots can be separated from the solid clusters not only by the horizontal line, but also by a vertical line, i.e. by a cutoff on max 2|s i − s 0 | or max 2k V |v i − 1|. We interpret configurations right of this line as having residual stresses of the respective kind, while configurations left of this line have no residual stresses. The fact that the lower-right of the quadrants formed by the two lines in both plots is mostly devoid of configurations indicates that residual stresses were in the vast majority of cases sufficient to create rigidity in our simulations. Finally, although we present here only the plots for the case k V = 1 (Figs. 2a and 7a), we verified using analogous plots that our conclusions are unchanged for other values of k V between 0.1 and 1000.
a. Rare cases of fluid configurations with residual stresses To study the rare exception cases of fluid configurations with residual stresses, we ran dedicated sets of simulations for k V = 1 and s 0 = 5.35 . . . 5.45 in steps of 10 −3 where we additionally stored all cell positions of the resulting shear-stabilized energy-minimized states. Like in Figs. 2a and 7a, we observed again few simulation runs that were in the fluid regime while containing significant residual stresses. Analyzing the individual cells in these states, we found that in most cases, only one cell deviated significantly from preferred surface and volume, while the other cells had no residual stresses (i.e. their surface and volume deviations were below the respective cutoffs). Interestingly, in almost all of these cases, the cell with residual stresses always had the shape of a so-called truncated square trapezohedron (TST) (Fig. 7b). In our simulations, these TSTs always had a surface larger than s 0 and a volume smaller than 1. Moreover, they always had a surface-to-volume ratio of s i /v 2/3 i ≈ 5.44364. We tested whether 5.44364 is the minimal surfaceto-volume ratio that a TST can assume. To this end, we derived analytical expressions for surface S and volume V of a TST, assuming 4-fold symmetry around the axis connecting the midpoints of the two squares and mirror-rotational symmetry perpendicular to this axis. Our expressions for S and V contain three free parameters, one of which is a linear scaling parameter. Numerical minimization of s = S/V 2/3 with respect to the other two parameters yielded indeed a minimum at s TST = 5.44363528 (9).
This suggests the following picture for the vast majority of cases that are fluid while containing residual stresses. Because s TST is larger than the transition point s * 0 , for parameters s 0 in between both values, all cells can attain their preferred surfaces and volumes except for TST-shaped cells. If a TST-shaped cell is contained in a configuration in this parameter regime, it will adjust surface area and volume to minimizes its own energy and will thus attain surface s i and volume v i such that s i /v 2/3 i = s TST . Its surface area s i will thereby be stretched above s 0 while its volume v i is compressed below 1. The actual value of v i depends on k V . For small k V the deviation v i − 1 will be larger than for large k V .
Note that theoretically, the fact that for s 0 < s TST TST-shaped cells attain a volume smaller than 1 creates a slight shift in the transition point for these configurations. This is because the non-TST cells now have to occupy a slightly larger volume. Thus, they are also forced to attain a larger surface area. This induces an increase of the rigidity transition point by a very small offset. Also note that very rarely (3 times out of ca. 10,000), we have also encountered fluid configurations with residual stresses that were not explained by any of the reasons given in this section. In these cases, many cells had surface and volume deviations above the cutoff in these cases. Due to numerical limitations, it was not possible to distinguish whether these cases represent real physics or just regions in the energy functional with very shallow gradients. 6. The onset of residual stresses occurred in all cells at once.
Here we show that the rigidity transition was in our simulations a truly collective transition in the sense that the onset of residual stresses occurred in all cells at once. To this end, using the two-dimensional histograms in Fig. 8, we show that if one residual stress was nonvanishing, all were non-vanishing.
First we show that whenever there was at least one non-vanishing cell surface tension, then there was a nonvanishing cell pressure, and vice versa. To this end, we correlate maximum surface tension with maximum pressure (Fig. 8a). The horizontal and vertical magenta dashed lines represent again the respective cutoffs shown in Figs. 7 and defined in Section C 5. Indeed, in this plot the upper-left and lower-right quadrants are completely empty, indicating that there are finite surface tensions if and only if there are finite pressures.
We next show that, whenever there was one nonvanishing cell surface tension, all cells had finite surface tensions. We correlate maximal with minimal cell surface tension in Fig. 8b. The vertical dashed lines are again taken from Section C 5. We find that the fluid clusters can be clearly separated not only by a vertical, but also by a horizontal cutoff line, separating configurations where at least one cell has no surface tension from simulations where all cells have a finite surface tension. The fact that only very few configurations are in the lowerright quadrant indicate that mostly whenever one cell deviated from its preferred surface, all cells did. There are a few exceptions to these observations, appearing in the lower-right quadrant, which are typically due to the appearance of a TST-shaped cell (see previous section).
Finally, we show that, whenever there was one finite cell pressure, all cell pressures were finite. We correlate maximal with minimal cell pressure in Fig. 8c. We find again that the fluid clusters can be clearly separated by a horizontal cutoff line, separating configurations where at least one cell has no pressure from simulations where all cells have a finite pressure. The fact that only few configurations appear in the lower-right quadrant indicates that except for a handful of exceptions, whenever one cell deviated from its preferred volume, all cells did.
These findings indicate that the rigidity transition in the 3D Voronoi model is indeed a truly collective one, where in the vast majority of cases, all cells simultane- ously acquire both nonzero surface tension and nonzero pressure at the onset of rigidity. This was true for all tested k V values between 0.1 and 1000. However, it is in principle possible that these findings are due to finite-size effects. We thus created plots similar to those in Fig. 8 for N c = 4096 cells. Although we found that there were more configurations where only some of the cells had a finite surface tension that for N c = 512 (∼ 1% of all configurations), we found that in all these cases a TSTshaped cell with its minimal possible surface-to-volume ratio appeared. Thus, with the exception of the appearance of TSTs, even for N c = 4096 the onset of residual stresses is collective.

Surface and volume standard deviations
Here, we display the plots for the surface standard deviation σ s and volume standard deviation σ v (Fig. 9). Note that both standard deviations vanish in the floppy regime, become nonzero at the transition point, and increase when going deeper into the solid regime. Here, we derive the scaling of the average surface tension 2( s − s 0 ) as a function of the distance from the transition point s * 0 . Because we study energy-minimized states and the energy can be written in terms of σ s and σ v (Eq. (4) with s = s min (σ s , σ v )), we can have: Together with Eq. (5) and s = s min (σ s , σ v ), we obtain indeed that: