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

TL;DR

This paper extends a 2D cellular model to 3D, revealing a rigidity transition controlled by cell surface area, with residual stresses playing a key role, and highlights implications for biological tissues and foams.

Contribution

The study introduces a 3D cellular model demonstrating a rigidity transition governed by preferred surface area, with residual stresses necessary for rigidity, differing from particulate matter.

Findings

Rigidity transition occurs at a critical surface area s_0*≈5.413.

Residual stresses are essential for rigidity in the model.

Shear modulus scales linearly with the control parameter at transition.

Abstract

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.

Figures

Equations

$$E=\sum_{i}{\bigg{[}K_{V}(V_{i}-V_{0})^{2}+K_{S}(S_{i}-S_{0})^{2}\bigg{]}}\text{,}$$

$$e=\sum_{i}{\bigg{[}k_{V}(v_{i}-1)^{2}+(s_{i}-s_{0})^{2}\bigg{]}}\text{.}$$

$${\bar{D}}_{pq}$$

$$\displaystyle\qquad\;+(s_{i}-s_{0})\frac{\partial^{2}s_{i}}{\partial z_{p}\partial z_{q}}+k_{V}(v_{i}-1)\frac{\partial^{2}v_{i}}{\partial z_{p}\partial z_{q}}\Bigg{]}\text{,}$$

$$e=N_{c}\Bigl{[}\big{(}\langle s\rangle-s_{0}\big{)}^{2}+\sigma_{s}^{2}+k_{V}\sigma_{v}^{2}\Bigr{]}\text{.}$$

$$s_{\min }({\sigma }_s,{\sigma }_v)=s_0^{\ast }-a_s{\sigma }_s-a_v{\sigma }_v$$

$$g=\frac{1}{N_c}{\left[\sum _m\frac{({\bar{u}}_{\gamma }^m{)}^2}{{\bar{\omega }}_m^2}\right]}^{-1}\text{.}$$

$$E$$

$$+\frac{K_V}{N_c}(V_{\mathrm{tot}}-N_c{V}_0{)}^2\text{.}$$

$$r_{ij}^x$$

$$r_{ij}^y$$

$$r_{ij}^z$$

$$r_{ij}^x$$

$$r_{ij}^y$$

$$r_{ij}^z$$

$$e=\sum _i{e}_i,$$

$$e_i=k_V(v_i-1{)}^2+(s_i-s_0{)}^2.$$

$$v_i$$

$$s_i$$

$${\mathbf{a}}_{in}=\frac{1}{2}\sum _{m=1}^{N_{in}}\Delta {\mathbf{h}}_{in,m}\times \Delta {\mathbf{h}}_{in,m+1}$$

$$\Delta {\mathbf{h}}_{in,m}={\mathbf{h}}_{in,m}-{\mathbf{r}}_i$$

$$\Delta {\mathbf{h}}_{ijlp}=Z_{ijlp}{\mathbf{H}}_{ijlp}\text{,}$$

$$Z_{ijlp}$$

$${\mathbf{H}}_{ijlp}$$

$$+|{\mathbf{r}}_{ip}{|}^2({\mathbf{r}}_{ij}\times {\mathbf{r}}_{il})\text{.}$$

$$0$$

$$\frac{\partial e_{i}}{\partial\bm{r}_{ij}}=2k_{V}\big{(}v_{i}-1\big{)}\frac{\partial v_{i}}{\partial\bm{r}_{ij}}+2\big{(}s_{i}-s_{0}\big{)}\frac{\partial s_{i}}{\partial\bm{r}_{ij}}\text{.}$$

$$\frac{\partial v_i}{\partial r_{ij}^{\alpha }}$$

$$\frac{\partial s_i}{\partial r_{ij}^{\alpha }}$$

$$\frac{\partial a_{in}^{\beta }}{\partial r_{ij}^{\alpha }}$$

$$\displaystyle\qquad\qquad\qquad\qquad+\Delta h_{in,m}^{\gamma}\frac{\partial\Delta h_{in,m+1}^{\delta}}{\partial r_{ij}^{\alpha}}\bigg{)}\text{,}$$

$$\frac{\partial \Delta h_{ijlp}^{\eta }}{\partial r_{ij}^{\alpha }}=\frac{\partial Z_{ijlp}}{\partial r_{ij}^{\alpha }}{H}_{ijlp}^{\eta }+Z_{ijlp}\frac{\partial H_{ijlp}^{\eta }}{\partial r_{ij}^{\alpha }}$$

$$\frac{\partial Z_{ijlp}}{\partial r_{ij}^{\alpha }}$$

$$\frac{\partial H_{ijlp}^{\eta }}{\partial r_{ij}^{\alpha }}$$

$$\frac{\partial e_i}{\partial {\mathbf{r}}_j}=\frac{\partial e_i}{\partial {\mathbf{r}}_{ij}}$$

Full text

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 Angelini et al. (2011); Sadati et al. (2013); Schotz et al. (2013); Nnetu et al. (2012), which likely plays a role in development and disease Park et al. (2015); Pawlizak et al. (2015). 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 Ingber (1997) and groups of cells below confluence Belmonte et al. (2008); Henkes et al. (2011) has been fruitfully studied using simple tensegrity and particle models. To understand rigidity in confluent tissues, we study a simple vertex model Honda et al. (1984); Farhadifar et al. (2007); Fletcher et al. (2014); Bi et al. (2014, 2015); Su and Lan (2016); Barton et al. (2016); Alt et al. (2017) 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 Bi et al. (2016); Su and Lan (2016); Barton et al. (2016); Sussman and Merkel (2017), and including active dynamics Sato et al. (2015); Bi et al. (2016); Barton et al. (2016). 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 Graner and Glazier (1992); Kabla (2012).

Recently, it was discovered that such models exhibit rigidity or glass-like transitions controlled by the non-dimensionalized preferred perimeter $p_{0}=P_{0}/\sqrt{A_{0}}$ Farhadifar et al. (2007); Staple et al. (2010); Bi et al. (2014, 2015, 2016). While the hexagonal ground state becomes linearly unstable if $p_{0}>3.72$ Farhadifar et al. (2007); Staple et al. (2010), there is a bona fide second order rigidity transition in disordered structures with a critical point reported at $p_{0}^{\ast}=3.81$ Bi et al. (2014, 2015). Adding active self-propulsion converts this point into a line of glass-like transitions in vertex models Bi et al. (2016) and cellular potts models Chiang and Marenduzzo (2016). A priori model predictions for cell shapes have recently been verified in experiments on glassy 2D monolayers Park et al. (2015). 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 Sarkar et al. (2013); Kamien and Liu (2007). 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 Maxwell (1864); Calladine (1978); Liu and Nagel (2010a). In models for fiber networks, the bond occupation probability controls the transition point also according to Maxwell’s criterion Jacobs and Thorpe (1995); Ellenbroek et al. (2015); Lubensky et al. (2015). Moreover, exploring the effect of dimensionality Charbonneau et al. (2017) 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 Honda et al. (2008); Hannezo et al. (2014); Monier et al. (2015); Okuda et al. (2015); Bielmeier et al. (2016); Misra et al. (2016) are an an active area of research, there is surprisingly little work on confluent bulk 3D tissues Honda et al. (2004); Viens and Brodland (2007), 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 Schotz et al. (2013), rapid advances in microscopy and segmentation algorithms Keller (2013); Stegmaier et al. (2016) 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 Bi et al. (2016) to three dimensions, and find a rigidity 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 $\bm{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 Bi et al. (2016) to 3D, inter-cellular forces are defined as derivatives on an effective energy functional:

[TABLE]

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) Su and Lan (2016); Yang et al. (2017). Therefore, without loss of generality we set $V_{0}=\langle V\rangle=L^{3}/N_{c}$ and non-dimensionalize our model with respect to the length unit $\langle V\rangle^{1/3}$ and the energy unit $K_{S}\langle V\rangle^{4/3}$, leading to the dimensionless energy:

[TABLE]

The remaining three dimensionless parameters are the preferred shape index $s_{0}=S_{0}/\langle V\rangle^{2/3}$, the relative volume rigidity $k_{V}=K_{V}\langle V\rangle^{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 [math] and the ideal gas value $5.82$, which is the average cell surface area $\langle s\rangle=(\sum_{i}{s_{i}})/N_{c}$ for randomly placed cell centers Meijering (1953) (Appendix B.1). Reported values are averaged over 100 random initial configurations minimized to a local energy minimum using the BFGS algorithm Fletcher (2000), varying the cell positions $\bm{r}_{i}$ and the simple shear degree of freedom $\gamma$ (Appendix LABEL:sec:energyMinimization).

III Results

III.1 Existence of the rigidity transition

For each relaxed state, we calculate the simple shear modulus $g=(\mathrm{d}^{2}e/\mathrm{d}\gamma^{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}^{\ast}\approx 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}^{\ast}$ nor shear modulus $g$ depend strongly on $k_{V}$. We also studied the precise distribution of transition points $p(s_{0}^{\ast})$ (Fig. 1b) (Appendix LABEL:sec:transitionPoint). 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}^{\ast}\approx 5.413$ in the limit $N_{c}\rightarrow\infty$ (Appendix LABEL:sec:fss).

III.2 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 Lubensky et al. (2015); Vermeulen et al. (2017). 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 Hessian $\bar{D}_{pq}$:

[TABLE]

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 Calladine (1978); Ingber (1997); Alexander (1998); Vermeulen et al. (2017). In our system, residual stresses are the $2N_{c}$ surface tensions $2(s_{i}-s_{0})$ and pressures $2k_{V}(1-v_{i})$.

To study whether residual stresses are both necessary and sufficient to generate rigidity, we plot a two-dimensional 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. 4a, Appendix B.3). 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 B.2). 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. 4b). 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 B.4).

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 $\langle s\rangle$ equals the preferred shape index $s_{0}$, as shown in Fig. 2b and there is zero standard deviation in cell surfaces and volumes ($\sigma_{s}=\sigma_{v}=0$ with $\sigma_{s}^{2}=(\sum_{i}{[s_{i}-\langle s\rangle]^{2}})/N_{c}$ and $\sigma_{v}$ defined analogously, see Fig. 6). 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}^{\ast}$ is determined by a purely geometric criterion: it corresponds to a local minimum in the average surface area $\langle s\rangle$ under the constraint that there are no surface and volume fluctuations $\sigma_{s}=\sigma_{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 $\langle s\rangle$ for disordered Voronoi tessellations with $\sigma_{s}=\sigma_{v}=0$, which is given by $s_{0}^{\ast}\approx 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 Xu et al. (2005); Chaudhuri et al. (2010); Charbonneau et al. (2012). It has been suggested that this packing fraction (about 64% in 3D) can be defined as the maximally disordered rigid state of spheres Torquato et al. (2000), or a divergence in the rate at which accessible disordered states disappear Kamien and Liu (2007).

As in theories for jamming, we must specify that we are restricting ourselves to random tessellations, as $\langle s\rangle$ can become smaller than $s_{0}^{\ast}$ for ordered states: the Voronoi packing corresponding to the ordered Kelvin structure has $\langle s\rangle\approx 5.315$ (Appendix LABEL:sec:periodicPackings). 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 Chaudhuri et al. (2010); 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 $\lesssim 0.01$. A weak dependence on the minimization protocol was recently also observed for the 2D Voronoi model Sussman and Merkel (2017).

III.3 Universal behavior and scaling of the shear modulus

Fig. 2b also shows that in the solid regime, the average surface area $\langle s\rangle$ decreases below the transition point $s_{0}^{\ast}\approx 5.413$, as the surface and volume fluctuations rise away from zero (Fig. 6). This is not surprising, as Eq. (2) can be rewritten:

[TABLE]

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: $\langle s\rangle=s_{\mathrm{min}}(\sigma_{s},\sigma_{v})$ Merkel and Manning (2017a). In the vicinity of the transition point, $s_{\mathrm{min}}(\sigma_{s},\sigma_{v})$ is surprisingly simple and universal: we find numerically that

[TABLE]

with $a_{s}\approx 2.0$ and $a_{v}\approx 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_{\mathrm{min}}$ explains the behavior of the shear modulus close to the transition point, which scales linearly with the distance from the transition point, $\delta s_{0}=s_{0}^{\ast}-s_{0}$ (Fig. 3b). To understand this, we start from the formula Eq. (61) for the shear modulus $g$ (Appendix A.7), which we restate here: If and only if the extended Hessian $\bar{D}_{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:

[TABLE]

Here the sum is over all strictly positive eigenvalues $\bar{\omega}_{m}^{2}$ of the extended Hessian $\bar{D}_{pq}$, and $\bar{u}_{\gamma}^{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 of $\bar{D}_{pq}$ with eigenvalue zero, corresponding to the $(N_{c}+1)$ zero modes obtained by the naive constraint counting. As shown in Appendix B.2, 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 $\langle v\rangle=1$, Eqs. (4) and (5) imply that the average surface tension scales linearly with the distance to the transition point: $2(\langle s\rangle-s_{0})=2\delta s_{0}/(1+a_{s}^{2}+a_{v}^{2}/k_{V})$ (Appendix B.6). Moreover, it turns out that the surface tensions and pressures of each individual cell are proportional to $\delta s_{0}$ Merkel and Manning (2017a). As a consequence, many of the $(N_{c}+1)$ eigenvalues that were zero in the fluid increase by an amount proportional to $\delta 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\sim\delta s_{0}$. We show this numerically by plotting a rescaled shear modulus $g^{\prime}=(1+a_{s}^{2}+a_{v}^{2}/k_{V})g$ over $\delta s_{0}$. The fact that this quantity collapses for different values for $k_{V}$ demonstrates that the average surface tension $2(\langle s\rangle-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 Honda et al. (2004); Viens and Brodland (2007). Our 3D Voronoi-based model, which is a straightforward generalization of successful 2D models for epithelial sheets Bi et al. (2016), exhibits a rigidity transition when the dimensionless preferred cell surface area is $s_{0}^{\ast}\approx 5.413$. The transition is accompanied by a structural order parameter, the observed average cell surface area $\langle s\rangle$.

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 Bi et al. (2016); Sussman and Merkel (2017).

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 Vermeulen et al. (2017), 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 Bi et al. (2015, 2016) and even holds in experiments on biological tissues Park et al. (2015). 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 $\langle s\rangle$ 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 Schotz et al. (2013) and human breast cancer cells Pawlizak et al. (2015), as well as 3D tissue explants from vertebrate embryos Hopyan (2017). Our model predicts that in the solid $\langle s\rangle\approx 5.4$, and $\langle s\rangle$ should rise away from that value as tissue fluidizes.

Although a similar prediction was successful in 2D Park et al. (2015), 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 Merkel and Manning (2017b), 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 Puliafito et al. (2012); Matoz-Fernandez et al. (2017)) 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 Morse and Corwin (2016), 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.3 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 packing fraction) Liu and Nagel (2010b), 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}^{\ast}\approx 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 Moshe et al. (2017).

It would also be interesting to test the predictions of our model in passive cellular materials like foams or biomimetic cellular materials Weaire and Hutzler (2001); Kraynik et al. (2003); Jorjadze et al. (2011). For example, random foams are deep in the solid phase ($s_{0}\rightarrow-\infty$ and $k_{V}/|s_{0}|\rightarrow\infty$), 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. Kraynik et al. (2003) reports values for the average shape index $\langle s\rangle$, but not the fluctuations $\sigma_{s}$. It would be interesting to revisit these results and determine if the surfaces are related to our minimal family $s_{\mathrm{min}}(\sigma_{s},\sigma_{v})$.

Finally, our findings could also be used to create artificial cellular materials that may be hyperuniform Bi (2016), 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 Brown et al. (2013); Jorjadze et al. (2011), or on macroscopic parameters like the overall volume or the pressure, which change the average cell volume.

Acknowledgements.

MM and MLM thank Peter Morse, Daniel M. Sussman, and Michael Moshe, as well as Steffen Grosser, Jürgen Lippoldt, and Josef Käs for fruitful discussions. Both authors also acknowledge funding from the Alfred P. Sloan Foundation, the Gordon and Betty Moore Foundation, the Research Corporation for Scientific Advancement, and computational support through NSF ACI-1541396. MLM also acknowledges support from the Simons Foundation under grant numbers 446222 and 454947, and NSF-DMR-1352184 and NSF-PHY-1607416.

Appendix A Analytic results that allow efficient computation and analysis of the 3D Voronoi model

A.1 Effect of preferred 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:

[TABLE]

Here $V_{\mathrm{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}=\langle V\rangle$ such that the last term disappears. Note that an analogous argument also holds in 2D Su and Lan (2016); Yang et al. (2017).

A.2 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 $\bm{r}_{ij}$ between neighboring cell positions $i$ and $j$ instead of absolute cell positions (see next section, Section A.3). Separately, these distance vectors $\bm{r}_{ij}$ are expressed in terms of absolute cell positions $\bm{r}_{i}$ and $\bm{r}_{j}$. However, the $\bm{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 $\bm{q}_{ij}$. It is defined such that $\bm{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_{ij}^{z}=+1$, and analogous for the other components of $\bm{q}_{ij}$ (see also the appendix in Merkel (2014) for an extensive explanation of the analogous 2D case). The distance vectors $\bm{r}_{ij}=(r_{ij}^{x},r_{ij}^{y},r_{ij}^{z})$ then depend as follows on the cell position and the dimensionless periodic box dimensions $l_{x},l_{y},l_{z}$:

[TABLE]

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 $\gamma$. It is implemented by modifying the above relations as follows:

[TABLE]

Note that this relation between cell distance vectors $\bm{r}_{ij}$ and cell positions $\bm{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.

A.3 Energy

Here, we express the total energy $e$ in terms of the distance vectors $\bm{r}_{ij}$. It is given by the sum of all cell energies $e_{i}$:

[TABLE]

with

[TABLE]

Volume and surface of a cell $i$ are given by:

[TABLE]

Here, both sums are over all cells $n$ neighboring cell $i$. The vector $\bm{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 $\bm{a}_{in}$ of the face between cells $i$ and $n$ in terms of the distance vectors $\bm{r}_{ij}$. To this end, we first express $\bm{a}_{in}$ in terms of the positions $\bm{h}_{in,m}$ of the vertices that define the polygonal face:

[TABLE]

with

[TABLE]

being the vertex position relative to the position of cell $i$. In Eq. (14), 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 $\times$ denotes the vector product and $\bm{h}_{in,N_{in}+1}\equiv\bm{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 $\bm{h}_{ijlp}$, relative to the position of cell $i$: $\Delta\bm{h}_{ijlp}=\bm{h}_{ijlp}-\bm{r}_{i}$. The relative vertex position in terms of the distance vectors $\bm{r}_{ij}$, $\bm{r}_{il}$, and $\bm{r}_{ip}$ is:

[TABLE]

where

[TABLE]

That Eqs. (16)–(18) yield the right vertex position can be seen as follows. The vertex position $\bm{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 $\bm{x}$ on each of these planes respectively fulfills

[TABLE]

Insertion verifies that $\bm{x}=\bm{h}_{ijlp}$ as defined by Eqs. (16)–(18) fulfills all three equations. Note that Eqs. (10)–(18) express the system energy $e$ only in terms of distance vectors $\bm{r}_{ij}$, where $i$ and $j$ are neighboring cells.

A.4 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 $\bm{r}_{ij}$. We have:

[TABLE]

The derivatives of volume and surface are:

[TABLE]

Here, both sums are over all neighbors $k$ of cell $i$. Moreover, here and in the following, Greek letters represent dimension indices: $\alpha,\beta,\dots\in\{x,y,z\}$ and we use Einstein convention, i.e. summation over same indices is implied.

The derivative of the oriented area is:

[TABLE]

where $\varepsilon^{\beta\gamma\delta}$ 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$, $\Delta\bm{h}_{ijlp}$, with respect to a distance vector $\bm{r}_{ij}$:

[TABLE]

with

[TABLE]

We now know the derivative of a cell energy $e_{i}$ with respect to the distance to any of its neighbors $j$, $\partial e_{i}/\partial\bm{r}_{ij}$. Using the chain rule, the derivative of $e_{i}$ with respect to neighbor $j$’s position is

[TABLE]

and the derivative of the energy of cell $i$ with respect to its own position is

[TABLE]

Here, the sum is over all neighbors $j$ of cell $i$. The total force on a cell $i$ can be written as $\bm{f}_{i}=-\sum_{k}{\partial e_{k}/\partial\bm{r}_{i}}$, where the sum is over all cells $k$ of the network. Inserting Eqs. (29)–(30), this becomes

[TABLE]

where both sums are over all cells $k$ of the network and the second sum is over all neighbors $j$ of cell $i$.

A.5 Hessian matrix

We start by deriving the derivative of the energy of a cell $i$ with respect to the cell distance vectors $\bm{r}_{ij}$ where cell $j$ is among the neighbors of cell $i$.

[TABLE]

The first derivatives of volume and surface are documented in the previous section. The second derivatives are:

[TABLE]

The second derivative of the oriented area is:

[TABLE]

The second derivative of the relative vertex position $\Delta\bm{h}_{ijlp}$ is:

[TABLE]

where cell $k$ is one of $j,l,p$. For $k=j$, we obtain for the second derivatives of $Z_{ijlp}$ and $\bm{H}_{ijlp}$:

[TABLE]

where $\delta^{\alpha\beta}$ is the Kronecker symbol. For $k=l$, the derivatives are:

[TABLE]

Finally, for $k=p$, they are:

[TABLE]

The second derivative of the total energy $e$ with respect to the absolute cell positions $\bm{r}_{i}$ is the Hessian matrix of the system:

[TABLE]

To derive an expression for it, we can proceed analogously to Eqs. (29)–(31). Based on the derivatives of the cell energies $e_{i}$ with respect to relative cell positions $\bm{r}_{ij}$, we obtain:

[TABLE]

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 $\gamma$. We obtain the derivatives involving simple shear using the chain rule and Eqs. (9)–(10):

[TABLE]

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:

[TABLE]

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

A.6 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, $\bm{r}_{ij}$, and in Section A.2 we expressed these distance vectors in terms of the absolute position vectors $\bm{r}_{i}$, and the simple shear variable $\gamma$. Thus, the energy can be expressed as $e=e(\{\bm{r}_{i}\},\gamma)$.

To introduce the long-time simple shear modulus $g$, we define $e^{\mathrm{min}}(\gamma)$ as the minimum of $e=e(\{\bm{r}_{i}\},\gamma)$ for fixed $\gamma$:

[TABLE]

We denote the set of cell positions with minimal $e$ for given $\gamma$ by $\bm{r}_{i}^{\mathrm{min}}(\gamma)$. Then we define the simple shear modulus as:

[TABLE]

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. (48):

[TABLE]

The second term in the parenthesis can be transformed as follows:

[TABLE]

where the sum is over all eigenvalues $\omega_{m}^{2}$ of the Hessian, and $u_{j\alpha}^{m}$ are the corresponding normalized eigenvectors:

[TABLE]

Note that since we minimize with respect to the $\bm{r}^{\mathrm{min}}$, the Hessian $D_{j\alpha,k\beta}$ is non-negative. To further transform Eq. (50), we take the total derivative of the force balance condition $0=\partial e(\{\bm{r}_{i}^{\mathrm{min}}(\gamma)\},\gamma)/\partial r_{j}^{\alpha}$ with respect to $\gamma$, which yields:

[TABLE]

Here, $\dot{r}_{k}^{\mathrm{min},\beta}=\mathrm{d}r_{k}^{\mathrm{min},\beta}/\mathrm{d}\gamma$. Inserting Eq. (51) and computing the $3N_{c}$-dimensional scalar product with $u_{j\alpha}^{m}$ yields:

[TABLE]

To simplify the sum over $m$ in Eq. (50), we distinguish between zero modes, for which $\omega_{m}=0$ and non-zero modes with $\omega_{m}>0$. For zero modes $m$, the second term in Eq. (53) has to be zero. Because this is the same term as the first factor in the sum over $m$ in Eq. (50), zero modes do not contribute to this sum. For the non-zero modes, we can substitute Eq. (53) into Eq. (50) and obtain from Eq. (49):

[TABLE]

In this equation, the sum over $m$ excludes zero modes of the Hessian $D_{j\alpha,k\beta}$.

A.7 Alternative computation of the shear modulus using the extended Hessian

Alternatively to Eq. (54), the shear modulus can also be computed from directly the eigen spectrum of the extended Hessian

[TABLE]

where $(z_{p})$ is a $(3N_{c}+1)$-dimensional vector comprising all cell positions and the shear degree of freedom $\gamma$: $(z_{p})=(\bm{r}_{1},\dots,\bm{r}_{N_{c}},\gamma)$. The eigenvalues of $\bar{D}_{pq}$ are $\bar{\omega}_{m}^{2}$ and the normalized eigenvectors are $\bar{u}_{p}^{m}$:

[TABLE]

Note that for shear-stabilized states, the extended Hessian $\bar{D}_{pq}$ is non-negative.

To derive an alternative formula for the shear modulus, we combine Eqs. (49) and (52) into:

[TABLE]

with $\delta$ being the Kronecker delta and $(\dot{z}_{q}^{\mathrm{min}})=(\dot{\bm{r}}_{1}^{\mathrm{min}},\dots,\dot{\bm{r}}_{N_{c}}^{\mathrm{min}},1)$. The $(3N_{c}+1)$-dimensional scalar product with an eigenvector $\bar{u}_{p}^{m}$ yields:

[TABLE]

for each eigenvalue $m$.

To evaluate the shear modulus, we distinguish two cases. First, whenever there exists a zero mode $m$ (i.e. $\bar{\omega}_{m}=0$) with nonzero shear component $\bar{u}_{\gamma}^{m}$, then Eq. (58) implies that the shear modulus $g$ has to vanish.

Second, all zero modes $m$ have vanishing shear components: $\bar{u}_{\gamma}^{m}=0$. In this case, the following relation

[TABLE]

holds also when only summing over non-zero modes $m$. As a consequence, the following relation holds:

[TABLE]

where we sum again only over non-zero modes $m$, because $\dot{z}_{\gamma}^{\mathrm{min}}=1$. Insertion of Eq. (58) finally yields:

[TABLE]

Here, the sum is over all non-zero modes $m$.

Appendix B Numerical implementation of the 3D Voronoi model

B.1 Initial conditions and parameter values

For all simulations, the cells are initially assigned independent random positions $\bm{r}_{i}$ drawn from a uniform distribution. To generate the Voronoi tessellations for a given set of cell positions $\{\bm{r}_{i}\}$, we used the oro++ library by Chris Rycroft (ersion 0.4.6, Rycroft (2009)). 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 igen3 library (version 3.2.7, \footnote{Available at: \protect\url{http://eigen.tuxfamily.org}}). We studied the cases $k_V=10^{-3},10^{-2},\dots,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. \subsection{Numerical energy minimization} \label{sec:energyMinimization} \begin{table} \begin{tabular}{p{6cm}|r} \bf Quantity & \bf Value \\hline Initial step size, \verb step_size & $0.01$ \ Line minimization tolerance, \verb tol & $0.01$ \ Total force norm cutoff, \verb epsabs & $10^{-12}\sqrt{N_\mathrm{dof}}$ \ Maximal number of iteration steps per minimization & $100N_\mathrm{dof}$ \end{tabular} \caption{Parameter values used for the multidimensional energy minimization. $N_\mathrm{dof}$ is the number of degrees of freedom varied during the minimization. For any minimization, all cell positions are included. Thus, depending on whether the shear degree of freedom is also included, $N_\mathrm{dof}=3N_c$ or $N_\mathrm{dof}=3N_c+1$.\label{tab:minParameters}} \end{table} To minimize the energy of the system, we used the BFGS2 multidimensional minimization routine of the GNU Scientific Library (GSL, version 2.1, \footnote{Available at: \protect\url{https://www.gnu.org/software/gsl/}}) \cite{Fletcher2000}. The parameters used for one GSL minimization are listed in \tref{tab:minParameters}, where $N_\mathrm{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 \tref{tab:minParameters}. In these cases, we tested whether the total force norm was at least below a cutoff of $C=10^{-7}\sqrt{N_\mathrm{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 $\vec{r}_i$, i.e.\ $N_\mathrm{dof}=3N_c$. Afterwards, to shear-stabilize the system, we included the shear degree of freedom $\gamma$ into the minimization, simultaneously varying $N_\mathrm{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 \fsref{fig:modelAndTransition}a and \sfref{fig:residualstressesAndMinSurface}b, 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. \subsection{Probing the solid vicinity of the transition point} \label{sec:solidVicinity} 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 $\vec{r}_i$ at some point $s_0$ in the solid regime from the cell positions at the transition point $s_0^\ast$ (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\delta s_0^2$ for any of the $s_0$ probed when exploring the solid vicinity of $s_0^\ast$. \section{Additional numerical analysis of the 3D Voronoi model} \subsection{Periodic Voronoi packings} \label{sec:periodicPackings} The average cell surface areas $\langle s\rangle$ for Kelvin ($\langle s\rangle\approx5.306$) and Weaire-Phelan packings ($\langle s\rangle\approx5.288$) are well-known. However, these figures correspond to configurations with curved cell outlines. To numerically obtain the $\langle s\rangle$ 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 $\vec{r}_1=(0,0,0)$ and $\vec{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 $\langle s\rangle\approx5.315$ and equal volume of both cells. \begin{table} \begin{tabular}{p{4.7cm}|c} \bf Type of cellular polyhedron & \bf Position \\hline \multirow{6}{3cm}{6 cells with 14 faces} & $(0.0, 1.5, 1.5)$ \ & $(1.0, 1.5, 1.5)$ \ & $(0.5, 0.0, 0.5)$ \ & $(0.5, 1.0, 0.5)$ \ & $(1.5, 0.5, 0.0)$ \ & $(1.5, 0.5, 1.0)$ \\hline \multirow{2}{3cm}{2 cells with 12 faces} & $(0.5, 0.5, 1.5)$ \ & $(1.5, 1.5, 0.5)$ \end{tabular} \caption{Initial cell positions to obtain the surface area corresponding to the Weaire-Phelan Voronoi structure. These cells are put into a in a cubic periodic box with side length 2.\label{tab:WeairePhelanVoronoi}} \end{table} 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 \tref{tab:WeairePhelanVoronoi}. Without constraining the volume, the surface was already minimal for these initial positions with $\langle s\rangle\approx5.295$. Note however that the volumes of the cells were slightly different. All 14-faced cells had $v_i\approx1.0078$ while the two 12-faced cells had $v_i\approx0.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. \begin{figure} \includegraphics{figure-suppl-preciseTransitionPoint.pdf} \caption{Determination of the precise transition point $s_0^\ast$. \textbf{(a)} Scatter plot of the shear modulus $g$ over the preferred surface area $s_0$ for $k_V=10$ in order to determine a cutoff value on $g$ (magenta dashed line) to define a network as rigid. Note that we chose different cutoff values for different values of $k_V$ (plots are similar for other values of $k_V$). The color of the dots indicate the norm of the total force vector at the end of the minimization, $\lvert f\rvert$. \textbf{(b)} Fraction of rigid networks depending on the preferred surface area $s_0$. The shaded regions indicate the respective uncertainties, computed as twice the standard deviation of the corresponding binomial distribution. The legend lists the resulting values for average and standard deviation of the transition point. System size $N_c=512$.\label{fig:preciseTransitionPoint}} \end{figure} \subsection{Determination of the precise transition point} \label{sec:transitionPoint} 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^\ast)$: \begin{equation} F(s_0) = 1 - \int_{-\infty}^{s_0}{P(s_0^\ast),\d s_0^\ast}\text{,} \end{equation} or conversely, $P(s_0^\ast)=-F’(s_0^\ast)$. The first two moments of the transition point $\langle s_0^\ast\rangle=\int_{-\infty}^\infty{s_0^\ast P(s_0^\ast)\,\d s_0^\ast}$ and $\langle (s_0^\ast)^2\rangle=\int_{-\infty}^\infty{(s_0^\ast)^2 P(s_0^\ast)\,\d s_0^\ast}$ can be evaluated using partial integration. Assuming that the transition point is always non-negative, we obtain: \begin{align} \langle s_0^\ast\rangle &= \int_0^\infty{F(s_0),\d s_0} \label{eq:tpAvg}\ \sigma^2(s_0^\ast) &= 2\int_0^\infty{s_0F(s_0),\d s_0} - \langle s_0^\ast\rangle^2 \label{eq:tpVar}\text{,} \end{align} where $\sigma(s_0^\ast)$ 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$ (\fref{fig:preciseTransitionPoint}a). 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 \fref{fig:preciseTransitionPoint}a. 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 $\langle s_0^\ast\rangle$ and standard deviation $\sigma(s_0^\ast)$ of the transition points are listed in \fref{fig:preciseTransitionPoint}b for the different values of $k_V$. The values for the average transition point range from $\langle s_0^\ast\rangle\approx5.414$ to $\langle s_0^\ast\rangle\approx5.416$. When respectively varying the cutoff on $g$ between the two clusters, the value of $\langle s_0^\ast\rangle$ varied by up to $\sim10^{-3}$. Thus, the average transition points for the different values of $k_V$ are not significantly different from each other. \begin{figure} \includegraphics{figure-suppl-median.pdf} \caption{Finite-size scaling of the transition point distribution using the median of $s_i-s_0$. \textbf{(a)} The median of $s_i-s_0$ is a very good indicator for rigidity as measured by the shear modulus $g$. Magenta lines indicate cutoff values. \textbf{(b)} Determination of the cutoff value for the median based on histograms of the median of $s_i-s_0$ depending on $s_0$ with $k_V=10$. \textbf{(c)} Average of the transition point distribution $\langle s_0^\ast\rangle$ depending on system size $N_c$. The black dashed line indicates a power law fit with offset $s_0^\ast\approx5.413$ and exponent $-0.9\pm0.4$, where we only used $N_c\geq512$, because a common fit through all data points was not consistent with the $N_c=256$ data point. \label{fig:median}} \end{figure} \subsection{Finite-size scaling} \label{sec:fss} 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_\mathrm{dof}\times N_\mathrm{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 (\fref{fig:median}a). Note in particular that this is the case even for the exceptions visible in \fref{fig:residualstressesAndMinSurface}a and discussed in \sref{sec:exceptions} 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$ (\fref{fig:median}b). Using this median cutoff, we extracted the functions $F(s_0)$ for each $N_c$ for $k_V=10$. Using \esref{eq:tpAvg} and \seref{eq:tpVar}, we then computed the average (\fref{fig:median}c) and variance (\fref{fig:modelAndTransition}b inset) of the transition point distribution. For the average, we find a scaling exponent of $-0.9\pm0.4$ and a limit value of $s_0^\ast(N_c\rightarrow\infty)=5.413\pm 0.001$. For the variance, we find a scaling exponent of $-0.90\pm0.04$. In \fref{fig:modelAndTransition}b, 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 $\sigma(s_0^\ast)/2$. \begin{figure} \centering \includegraphics{figure-suppl-overlapInfZeroModesWithShear.pdf} \caption{Overlap of the infinitesimal zero modes with the shear degree of freedom for different values of $s_0$ and $k_V=10$ (plots for other values of $k_V$ are similar). ach dot represents a single energy-minimized network configuration that is (a) shear stabilized or (b) not shear stabilized. In the shear-stabilized case, the overlap in the solid regime is finite, which indicates that residual stresses are necessary to create rigidity. 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 $(3N_{c}+1)$-dimensional (a) or $(3N_{c})$-dimensional (b) total force vector after the minimization.

B.2 Residual stresses were necessary to create rigidity.

In the main text, we discuss the full extended Hessian $\bar{D}_{pq}$ and find that the system is rigid below the transition point $s_{0}^{\ast}$. To numerically check whether residual stresses are necessary to create such rigidity, we study an unstressed extended Hessian $\bar{D}_{pq}^{\prime}$ where the terms that depend on the residual stresses have been removed:

[TABLE]

Here $(z_{p})$ is a $(3N_{c}+1)$-dimensional vector of the cell positions and the shear degree of freedom. $\bar{D}_{pq}^{\prime}$ 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) $\bar{D}_{pq}^{\prime}$ 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 that $\bar{D}_{pq}^{\prime}$ has such zero modes in the solid phase, we can conclude that residual stresses are necessary to rigidify the original system $\bar{D}_{pq}$. The zero modes of $\bar{D}_{pq}^{\prime}$ correspond to the infinitesimal zero modes of the system. Infinitesimal zero modes are $(3N_{c}+1)$-dimensional vectors $\bm{w}$ in the kernel of the so-called compatibility matrix $\mathrm{C}$ Lubensky et al. (2015), i.e. $\mathrm{C}\cdot\bm{w}=0$. For our system, $\mathrm{C}$ is a $(2N_{c})\times(3N_{c}+1)$ matrix defined as:

[TABLE]

Note that $\mathrm{C}$ is the same for both original and modified system. Hence, an infinitesimal zero mode corresponds to a collective change of cell center positions $\{\bm{r}_{i}\}$ and shear variable $\gamma$ 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 $\bm{w}^{q}$ with $q=1,\dots,N_{0}$ and $\mathrm{C}\cdot\bm{w}^{q}=0$, which are orthonormal: $\bm{w}^{p}\cdot\bm{w}^{q}=\delta_{pq}$. We extracted such an orthonormal set of zero modes $\bm{w}^{q}$ using singular value decomposition of $\mathrm{C}$. As a tool to study the shear components $w^{q}_{\gamma}$ of all infinitesimal zero modes $\bm{w}^{q}$ at once, we define the overlap $r$ of the kernel of $\mathrm{C}$ with the shear degree of freedom by:

[TABLE]

By definition, $0\leq r\leq 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. B.1a, 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\sim 10^{-3}$ both in the fluid regime ($s_{0}>s_{0}^{\ast}$) and in the solid regime ($s_{0}<s_{0}^{\ast}$), 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}^{\ast}\approx 5.4$, the overlap is occasionally smaller than $10^{-3}$. However, this occurs only in a fraction of the cases at $s_{0}^{\ast}$. Moreover, the overlap was always larger than $10^{-5}$. To ensure that the observed values of $r\in[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. B.1b. We find that the measured overlap in the solid phase is always below $r\sim 10^{-7}$ and is typically on the order of $r\sim 10^{-14}$, which is significantly smaller than $r\in[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 $\sigma_{xy}$, we first consider the force balance condition combined with the definition for the shear stress:

[TABLE]

Here, $\delta$ 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:

[TABLE]

Here, $\bm{\hat{\gamma}}$ is the $(3N_{c}+1)$-dimensional vector $\bm{\hat{\gamma}}=(0,\dots,0,1)$. For $\sigma_{xy}\neq 0$, we obtain an expression for $\bm{\hat{\gamma}}$ from Eq. (66), and insertion into $w^{q}_{\gamma}=\bm{\hat{\gamma}}\cdot\bm{w}^{q}$ yields indeed $w^{q}_{\gamma}=0$, because $\mathrm{C}\cdot\bm{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.

B.3 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 $\mathrm{max}\,2\lvert s_{i}-s_{0}\rvert$ (Fig. 2a in main text) and maximal pressure magnitude $\mathrm{max}\,2k_{V}\lvert v_{i}-1\rvert$ (Fig. 4a). In addition, we marked in each case the cutoff defined in Section LABEL:sec:transitionPoint below which we interpret a configuration as fluid by a horizontal magenta dashed line (compare Fig. LABEL:fig:preciseTransitionPoint). 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 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 $\Delta s_{0}\sim 10^{-1}$ away from the transition point, we varied $s_{0}$ in steps of $0.1$ (compare Section LABEL:sec:energyMinimization). 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 $\mathrm{max}\,2\lvert s_{i}-s_{0}\rvert$ or $\mathrm{max}\,2k_{V}\lvert v_{i}-1\rvert$. 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 4a), we verified using analogous plots that our conclusions are unchanged for other values of $k_{V}$ between $0.1$ and $1000$.

B.3.1 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\dots 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 4a, 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. 4b). 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_{i}^{2/3}\approx 5.44364$. We tested whether $5.44364$ is the minimal surface-to-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_{\mathrm{TST}}=5.44363528(9)$. This suggests the following picture for the vast majority of cases that are fluid while containing residual stresses. Because $s_{\mathrm{TST}}$ is larger than the transition point $s_{0}^{\ast}$, 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_{i}^{2/3}=s_{\mathrm{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_{\mathrm{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.

B.4 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. 5, we show that if one residual stress was non-vanishing, all were non-vanishing. First we show that whenever there was at least one non-vanishing cell surface tension, then there was a non-vanishing cell pressure, and vice versa. To this end, we correlate maximum surface tension with maximum pressure (Fig. 5a). The horizontal and vertical magenta dashed lines represent again the respective cutoffs shown in Figs. 4 and defined in Section B.3. 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 non-vanishing cell surface tension, all cells had finite surface tensions. We correlate maximal with minimal cell surface tension in Fig. 5b. The vertical dashed lines are again taken from Section B.3. 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 lower-right 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. 5c. 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 simultaneously 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. 5 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$ ($\sim 1\%$ of all configurations), we found that in all these cases a TST-shaped 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.

B.5 Surface and volume standard deviations

Here, we display the plots for the surface standard deviation $\sigma_{s}$ and volume standard deviation $\sigma_{v}$ (Fig. 6). 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.

B.6 Scaling of the average cell surface tension in the solid vicinity of the transition point

Here, we derive the scaling of the average surface tension $2(\langle s\rangle-s_{0})$ as a function of the distance from the transition point $s_{0}^{\ast}$. Because we study energy-minimized states and the energy can be written in terms of $\sigma_{s}$ and $\sigma_{v}$ (Eq. (4) with $\langle s\rangle=s_{\mathrm{min}}(\sigma_{s},\sigma_{v})$), we can have:

[TABLE]

Insertion of Eqs. (4) and (5) yields:

[TABLE]

Together with Eq. (5) and $\langle s\rangle=s_{\mathrm{min}}(\sigma_{s},\sigma_{v})$, we obtain indeed that:

[TABLE]