# Connecting anomalous elasticity and sub-Arrhenius structural dynamics in a cell-based model

Chengling Li, Matthias Merkel, Daniel Sussman

► To cite this version:

Chengling Li, Matthias Merkel, Daniel Sussman. Connecting anomalous elasticity and sub-Arrhenius structural dynamics in a cell-based model. Physical Review Letters, 2025, 134 (4), pp.048203. 10.1103/PhysRevLett.134.048203. [hal-04760155v2](https://hal.science/hal-04760155v2)

HAL Id: hal-04760155

https://hal.science/hal-04760155v2

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Connecting anomalous elasticity and sub-Arrhenius structural dynamics in a cell-based model

Chengling Li,¹,* Matthias Merkel,² and Daniel M. Sussman¹,†

¹Department of Physics, Emory University, Atlanta, GA, USA

²Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332),
Turing Centre for Living systems, Marseille, France

(Dated: January 6, 2025)

Understanding the structural dynamics of many-particle glassy systems remains a key challenge in statistical physics. Over the last decade, glassy dynamics has also been reported in biological tissues, but is far from being understood. It was recently shown that vertex models of dense biological tissue exhibit very atypical, sub-Arrhenius dynamics, and here we ask whether such atypical structural dynamics of vertex models are related to unusual elastic properties. It is known that at zero temperature these models have an elasticity controlled by their under-constrained or isostatic nature, but little is known about how their elasticity varies with temperature. To address this question we investigate the 2D Voronoi model and measure the temperature dependence of the intermediate-time plateau shear modulus and the bulk modulus. We find that unlike in conventional glassformers, these moduli increase monotonically with temperature until the system fluidizes. We further show that the structural relaxation time can be quantitatively linked to the plateau shear modulus $G_p$, i.e. $G_p$ modulates the typical energy barrier scale for cell rearrangements. This suggests that the anomalous, structural dynamics of the 2D Voronoi model originates in its unusual elastic properties. Based on our results, we hypothesize that under-constrained systems might more generally give rise to a new class of “ultra-strong” glassformers.

The rheological behavior of disordered materials, in particular in the presence of thermal fluctuations, is still not fully understood. In the athermal (zero-temperature) limit, it depends strongly on the precise disordered geometry and on the Maxwell constraint-counting criterion [1–5], which distinguishes three classes of systems. Roughly, “over-constrained” systems have more constraints than degrees of freedom and are generally rigid. Conversely, “under-constrained” systems have fewer constraints than degrees of freedom and are generally mechanically unstable. Finally, “isostatic” systems have balanced constraints and degrees of freedom, and often behave as weak solids.

Over-constrained systems, which typically show glassy behavior, are the most well-understood. Increasing temperature $T$ above the dynamical glass transition temperature $T_g$ facilitates particle rearrangements, which occur on a structural relaxation time scale $\tau_\alpha$. A simple picture of structural relaxation based on temperature-independent energy barriers of order $\Delta E$ would predict that $\tau_\alpha \sim \exp(\beta\Delta E)$, where $\beta = (k_B T)^{-1}$ is the inverse temperature. Such an exponential scaling of $\tau_\alpha$ with inverse temperature is called Arrhenius scaling, and standard particulate glasses display Arrhenius scaling or super-Arrhenius scaling (in which $\tau_\alpha$ grows faster than exponentially with the inverse temperature, c.f. the red dots in Fig. 1) [7]. A long-standing idea is that since particle rearrangements require local elastic deformations of the material $\Delta E$, and thus $\tau_\alpha$, should be connected to elasticity [8–10]. Above $T_g$ at long times the shear relaxation modulus $G(t)$ vanishes and glasses are fluid, but this is often preceded by an intermediate-time plateau

with value $G_p$. Because transient deformations cost energy, the prediction of the elastic or “shoving” models is that $\Delta E \sim G_p$, i.e. [9, 10]:

$$ \tau_{\alpha} \sim \exp(\beta C G_{p}), \quad (1) $$

where $C$ is a constant. In particulate glasses and other over-constrained systems $G_p$ typically varies extremely modestly with temperature, which would be consistent with Arrhenius behavior. Super-Arrhenius behavior can arise from a slightly larger temperature dependence of

<figure><img src="image_2.png" /><figcaption>FIG. 1. Structural relaxation in the 2D Voronoi model. Scaling of $\tau_\alpha$ of the 2D Voronoi model in the weak solid regime (light blue) and deeper in the solid regime (black), along with data from Ref. [6] for a standard 2D Kob-Andersen glass-former (red dots). Open symbols show the stress relaxation time, $\tau_s$. Both $\tau_\alpha$ and $\tau_s$ were quantified using standard approaches (see End Matter).</figcaption></figure>
$G_p(T)$ which can in turn originate from subtle changes to the inherent states sampled as $T$ changes.

Under-constrained systems are less well understood. The most commonly studied under-constrained systems have fixed connectivity, such as spring or fiber networks [12–14]. While generally floppy at $T = 0$, these systems can be rigidified through external strain or by varying model parameters in a way that introduces geometric incompatibility and prestresses. These prestresses in turn govern elastic properties such as the shear modulus. There is relatively little work on the thermal behavior of such systems. A few specific networks have been studied, where the shear modulus was found to scale as a power law with temperature, $G \sim T^\alpha$, for $0.5 \le \alpha \le 1$ [15, 16]. Recent analytical work by one of us on generic under-constrained networks found that for small temperatures there are three regimes; an energetically dominated regime with $\alpha = 0$ close to the athermal rigid regime, an entropically dominated regime with $\alpha = 1$ close to the athermal floppy regime, and a cross-over regime with $\alpha = 0.5$ [17, 18]. Thus, the shear modulus generally increases with temperature. However, because these models have fixed connectivity, no structural relaxation takes place.

There are also under-constrained models whose connectivity is not fixed, for instance “vertex models” of dense biological tissue. These models describe a monolayer of cells as a tiling of polygonal cells [19–21], in which the degrees of freedom are the positions of the polygon corners (i.e., the polygonal vertices). In a standard implementation of the model, the forces on the vertices are given by the negative gradient of an energy functional with dimensionless form [19]:

$$E = \sum_{i=1}^{N} \left[ (a_i - 1)^2 + (p_i - p_0)^2 \right]. \quad (2)$$

Here $N$ is the total number of cells, and $a_i$ and $p_i$ are the area and perimeter of cell $i$. The model parameter $p_0$ corresponds to the target cell perimeter. For simplicity, parameters corresponding to the preferred cell area, and the area and perimeter stiffness of the cells, have been set to unity. This form is often motivated biologically — the area term represents the resistance of the cells to volumetric changes, and the perimeter term represents a competition between active cellular tensions and intercellular adhesion [19, 21] — but it can also be viewed more generically as a low-order expansion in simple geometric quantities [22]. Vertex models have been used to describe many experiments on biological tissues, including observations of cellular jamming transitions and glassy dynamics [23–32].

Like under-constrained systems with fixed connectivity, athermal vertex models have a rigidity transition [12, 19, 33–36]. Specifically, vertex models are floppy with a vanishing shear modulus for $p_0 > p_0^*$ and rigid with a

finite shear modulus for $p_0 < p_0^*$, where the precise value of the transition point (usually quoted as $p_0^* \sim 3.81$) can depend on the network structure [37]. There are also isostatic variants of vertex models, such as the 2D Voronoi model [38–40], where the cell shapes are defined by a Voronoi tessellation, and the degrees of freedom are the Voronoi centers $r_i$ of the cells. There is no floppy regime in this model, but its shear modulus strongly decreases with increasing $p_0$, especially in the weakly solid regime $p_0 \gtrsim 3.8$ [39]. That is, at $T = 0$, the Voronoi model’s behavior seems strongly affected by the transition in the standard vertex model. Importantly, in vertex and Voronoi models structural rearrangements occur as cells are allowed to change their local connectivity, and one of us showed that in the floppy or weak solid regimes, the structural dynamics of these models are of the atypical sub-Arrhenius type (blue solid curve in Fig. 1) [41]. In the athermal solid regime, on the other hand, they show a more normal Arrhenius or super-Arrhenius behavior with temperature (black solid curve in Fig. 1) [42, 43].

Are the atypical structural dynamics of these shape-based models of cells related to the unusual elastic properties of under-constrained systems? Due to its simplicity, we focus on addressing this question for the 2D Voronoi model. We first measure its plateau shear modulus $G_p$ and its bulk modulus $K$. In the weakly solid regime, we find that both moduli monotonically increase with temperature at small $T$. In contrast, deep in the solid regime we find that $G_p$ is $T$-independent for $T \ll 1$ and decreases slightly for larger $T$. We then compare this elastic behavior to the shoving model, Eq. (1), and find data collapse when using a slightly modified form. Our findings suggest that the anomalous structural dynamics of the Voronoi model originate in the unusual temperature dependence of its elastic properties. This also explains why tuning the model from its weak to its deep solid regime leads to a change from sub- to super-Arrhenius structural dynamics [42]. Based on this, we expect to more generally observe sub-Arrhenius structural dynamics in under-constrained systems close to the athermal floppy regime.

We perform standard NVT simulations of the 2D Voronoi model with periodic boundary conditions. The Voronoi centers $r_i$ are subject to forces $f_i = -\partial E/\partial r_i$, where $E$ is the energy functional in Eq. (2). We use the deterministic Nosé-Hoover algorithm to control the temperature; in this setting inertia is present and there is no viscous damping that acts on the cell velocities (see End Matter for details).

We first quantify the plateau shear modulus $G_p$. To this end, we compute $G(t)$ via the shear stress autocorrelation function (see End Matter). Figure 2(a) shows sample $G(t)$ curves across a range of $T$ for $p_0 = 3.825$. We independently verified these results by measuring $G(t)$ from the stress response to a small shear step (Fig. 2(a) inset). We find that the magnitude of $G(t)$ is strongly temperature-dependent at almost all time
<figure><img src="image_3.png" /><figcaption>FIG. 2. Stress relaxation and static structure are strongly temperature-dependent. (a) Shear relaxation modulus as a function of time for $p_0 = 3.825$ and $T = [0.063–0.0005]$, decreasing from dark red to light blue. Dashed lines are stretched exponential fits. Inset: $G(T)$ at $T = 0.0025$ computed via stress autocorrelations ($G_S$) or a small step strain ($G_N$). (b) The isotropic static structure factor $S(q)$ for $p_0 = 3.825$ and the same range of temperatures as in (a). Inset: $S(q)$ for $p_0 = 3.8$ at $T = 0.00385, 0.0031, 0.0028, 0.0025, 0.0022$.</figcaption></figure>

scales (Fig. 2(a)). As $T$ decreases, similar to other glassy systems [10, 44], an intermediate-time plateau gradually emerges, followed by a final decay. We quantify $G_p$ and the duration $\tau_s$ of the intermediate-time plateau by fitting our $G(t)$ data to a stretched exponential form, $G(t) \approx G_p e^{-(t/\tau_s)^\beta}$ at intermediate and late times (dashed curves in Fig. 2(a), see End Matter). We find that $\tau_s$ is typically several times smaller than $\tau_\alpha$ (Fig. 1), suggesting that in this model stress-bearing structures can dissipate before individual cells move a substantial amount.

In Fig. 3, we plot the plateau shear modulus $G_p$ for $p_0 = (3.8, 3.825, 3.85)$ as a function of $T/T_g$, with $T_g$ defined by $\tau_\alpha(T_g) \equiv 10^4$. In contrast to particulate glassformers [10, 44, 45], in the Voronoi model $G_p$ strongly increases upon increasing $T$ until at high enough temperatures we find a pure exponential decay of $G(t)$, which we associate with a fluid phase (Fig. 3). At low tem-

peratures the system enters a hexatic phase (hexagons in Fig. 3), as measured by a peak in the susceptibility of the bond-orientational order parameter [46] – this is consistent with earlier hexatic-phase observations in vertex models [47]. Finally, we note that the $G(t)$ data for $p_0 = 3.85$ at the lowest temperatures was both noisy and poorly fit by a stretched exponential (red dots in Fig. 3), as we begin to hit the floor associated with the statistical error of our autocorrelator. We have ignored this data in our subsequent analysis.

As shown in Fig. 3, for $p_0 \gtrsim 3.8 G_p$ monotonically increases with temperature. In the lower inset we rescale $G_p$ by its value in our highest-$T$ simulations — we find excellent collapse of our data, with scaling that is consistent with $G_p$ crossing over between $T^1$ and $T^{0.5}$ power laws. Although our data cover a numerical range somewhat too limited to be fully conclusive, we note that these results are consistent with observations in various fixed-connectivity spring-network models [15, 16, 48] and those predicted in generic under-constrained systems [17, 18].

Previous work suggests an intuitive explanation for the increase of $G_p$ in under-constrained systems [17, 18]. Specifically, in the athermal floppy regime of the 2D vertex model, the accessible phase space volume increases as the cell perimeters increase above $p_0^*$. This pushes the cell perimeters to higher values for larger temperatures, leading to higher isotropic tensions and, consequentially, a higher plateau shear modulus $G_p$. Hence, the increase of $G_p$ with $T$ is created by a form of entropic elasticity. Indeed, in our isostatic 2D Voronoi simulations we clearly observe such an increase of cell perimeters and isotropic

<figure><img src="image_4.png" /><figcaption>FIG. 3. The plateau shear modulus $G_p$ is strongly temperature-dependent. Here state points corresponding to the hexatic phase are plotted by a hexagon. At the lowest temperatures for $p_0 = 3.85$, the tail of $G(t)$ is both noisy and poorly fit by a stretched exponential; we consider these measurements less reliable and note them with red open circles. Lower Inset: $G_p^* = G_p/G_p(T = 0.039)$ for each $p_0$ as a function of $T$. Dashed lines are guides to the eye of slope 1 and 1/2. Upper Inset: $T$-dependence of a $p_0 = 3.0$ Voronoi model whose structural relaxation is super-Arrhenius [42].</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. The bulk modulus *K* is strongly temperature-dependent. The bulk modulus as calculated from the statis structure factor (open symbols) and the response to isotropic compression (closed symbols). Inset: Contributions to *K* from area (*K*<sub>a</sub>) and perimeter (*K*<sub>p</sub>) terms in the energy.</figcaption></figure>

tension with temperature, consistent with the predicted scaling exponents (Fig. 7, End Matter).

Given that the Voronoi model shows super-Arrhenius behavior deeper in the solid regime, we also measured the plateau shear modulus for the $p_0 = 3.0$ model studied in [42] (Fig. 3, upper inset). We find that $G_p$ is roughly constant at small temperatures and decreases for larger temperatures. That is, in the regime in which the model's dynamics mirror those of a more typical glassformer, so too does its elasticity.

We also characterize the temperature dependence of the bulk modulus *K* from the low-wavevector behavior of the static structure factor (Fig. 2(b)), which we also verified by measuring the system response to an isotropic compression strain (see End Matter). In Fig. 4, we plot the bulk moduli obtained through both methods, $K_s$ and $K_N$, as a function of $T/T_g$ for $p_0 = (3.8, 3.825, 3.85)$. Again in contrast to typical glassy systems, we find that *K* increases monotonically with *T*: as the system is heated it becomes more resistant to isotropic deformations. As shown in the inset to Fig. 4, this temperature dependence of *K* comes entirely from cell perimeter contributions to the energy, consistent with work on the athermal transition in vertex models [12, 39].

Given the unusual temperature dependence of the Voronoi model's elastic moduli, does a shoving model such as Eq. (1) explain the anomalous structural dynamics? In Fig. 5 we plot $\log(\tau_\alpha)$ vs. $G_p/T$ for the weak solid regime and find neither a linear scaling nor a compelling collapse of the data. However, the inset shows excellent collapse by instead scaling $G_p$ by $T^{1.2}$. This deviation from the typical elastic model suggests the presence of some other very modestly temperature-dependent feature contributes to the energy barrier $\Delta E$.

In summary, we have shown that the atypical sub-Arrhenius structural dynamics in the 2D Voronoi model is directly correlated to the unusual temperature dependence of its elastic moduli. We specifically found

that in the weak solid regime the plateau shear modulus monotonically increases with temperature, in a way which is qualitatively consistent with analytical results on fixed connectivity under-constrained systems [17, 18]. Although the 2D Voronoi model is isostatic rather than under-constrained, we note that at $T=0$ the elastic moduli of the Voronoi model have been shown to be governed by the rigidity transition of the vertex model [39]. Additionally, we note that the structural dynamics of the 2D vertex model is known to be similar to that of the 2D Voronoi model [41], and the finite-temperature elastic behavior is expected to be extremely similar to what we report here for the Voronoi model [17, 18]. Taken together, we expect that the reported sub-Arrhenius structural dynamics of the 2D vertex model is also explained by its elastic properties in the large-$p_0$ regime. In the small-$p_0$ regime, in contrast, the more conventional structural dynamics [43] are likely related to the predicted more conventional elastic properties [17, 18].

More generally, we expect this scenario — unusual elasticity connected to anomalous structural dynamics close to the athermal fluid regime, and more conventional elasticity and structural dynamics in other parts of parameter space — to appear in any under-constrained system that is allowed to dynamically change its connectivity. To our knowledge the only class of systems with glassy behavior similar to the vertex model are low-density vitrimeric polymers [49], where a network of covalently bonded under-constrained monomers can change its connectivity on long time scales via bond-exchange reactions. We suggest that models with similar thermally dependent elasticity may all form a class of “ultra-strong” glasses, complementing the currently classification of “fragile” and “strong” glassformers. If true, our findings may help to develop new materials with such

<figure><img src="image_6.png" /><figcaption>FIG. 5. The mechanics and dynamics are directly connected. In the main plot we show $\log(\tau_\alpha)$ vs $G_p/T$; the inset shows $\log(\tau_\alpha)$ vs $G_p/T^{1.2}$. As in plots above, state points corresponding to hexatic phases or unreliable measurements of $G(t \to \infty)$ are noted.</figcaption></figure>
an ultra-strong glass-forming ability.

Finally, our results have implications for the model-
ing of biological tissues. Some modeling of dense tissue
uses not shape-based under-constrained models but more
typical particle-based models [50, 51]. Given the differ-
ences in elastic and structural properties between partic-
ulate glasses and vertex models, our work shows how this
choice can lead to fundamentally different tissue-scale be-
havior. For instance, under certain conditions, e.g., when
the auto-correlation time of cell motile forces are smaller
than any other relevant time scale of the system, cell
motility may be approximated by thermal noise. For such
systems, our findings suggest that an increase in cellular
motility would not only fluidify the tissue on long times
[38], but may also *increase* the intermediate-time, tissue-
scale rigidity (*G*<sub>p</sub>). It would be interesting to test this
experimentally.

This material is based upon work supported by the
National Science Foundation under Grant No. DMR-
2143815 (DMS and CL). This research used the Delta
advanced computing and data resource which is sup-
ported by the National Science Foundation (award OAC
2005572) and the State of Illinois. Delta is a joint effort
of the University of Illinois Urbana-Champaign and its
National Center for Supercomputing Applications.

<footnote>* chengling.li@ memory.edu</footnote>

<footnote>† daniel.m.sussman@memory.edu</footnote>

<footnote>[1] James Clerk Maxwell. On the calculation of the equi-
librium and stiffness of frames. *Philosophical Magazine
Series 4*, 27(182):294–299, 1864.</footnote>

<footnote>[2] C. R. Calladine. Buckminster Fuller's "Tensegrity" struc-
tures and Clerk Maxwell's rules for the construction of
stiff frames. *International Journal of Solids and Struc-
tures*, 14(2):161–172, 1978.</footnote>

<footnote>[3] TC Lubensky, CL Kane, Xiaoming Mao, Anton Souslov, and Kai Sun. Phonons and elasticity in critically co-ordinated lattices. *Reports on Progress in Physics*, 78(7):073901, 2015.</footnote>

<footnote>[4] Anaël Lemaître and Craig Maloney. Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature. *Journal of statistical physics*, 123:415–453, 2006.</footnote>

<footnote>[5] Alessio Zaccone and Enzo Scossa-Romano. Approximate analytical description of the nonaffine response of amorphous solids. *Physical Review B—Condensed Matter and Materials Physics*, 83(18):184205, 2011.</footnote>

<footnote>[6] Yan-Wei Li, Chandan K Mishra, Zhao-Yan Sun, Kun Zhao, Thomas G Mason, Rajesh Ganapathy, and Masimo Pica Ciamarra. Long-wavelength fluctuations and anomalous dynamics in 2-dimensional liquids. *Proceedings of the National Academy of Sciences*, 116(46):22977–22982, 2019.</footnote>

<footnote>[7] Mark D Ediger, C Austen Angell, and Sidney R Nagel. Supercooled liquids and glasses. *The journal of physical chemistry*, 100(31):13200–13212, 1996.</footnote>

<footnote>[8] Andrea Cavagna. Supercooled liquids for pedestrians.</footnote>

<footnote>*Physics Reports*, 476(4-6):51–124, 2009.</footnote>

<footnote>[9] Jeppe C Dyre and Wei Hua Wang. The instantaneous shear modulus in the shoving model. *The Journal of chemical physics*, 136(22), 2012.</footnote>

<footnote>[10] F Puosi and Dino Leporini. Communication: Correlation of the instantaneous and the intermediate-time elasticity with the structural relaxation in glassforming systems. *The Journal of chemical physics*, 136(4), 2012.</footnote>

<footnote>[11] Bingyu Cui, Giancarlo Ruocco, and Alessio Zaccone. Theory of elastic constants of athermal amorphous solids with internal stresses. *Granular Matter*, 21(3):69, August 2019.</footnote>

<footnote>[12] Matthias Merkel, Karsten Baumgarten, Brian P Tighe, and M Lisa Manning. A minimal-length approach unifies rigidity in underconstrained materials. *Proceedings of the National Academy of Sciences*, 116(14):6560–6568, 2019.</footnote>

<footnote>[13] Sihan Chen, Tomer Markovich, and Fred C MacKintosh. Field theory for mechanical criticality in disordered fiber networks. *Physical Review Letters*, 133(2):028201, 2024.</footnote>

<footnote>[14] Chase P Broedersz, Xiaoming Mao, Tom C Lubensky, and Frederick C MacKintosh. Criticality and isostaticity in fibre networks. *Nature Physics*, 7(12):983–988, 2011.</footnote>

<footnote>[15] Xiaoming Mao, Anton Souslov, Carlos I Mendoza, and Tom C Lubensky. Mechanical instability at finite temperature. *Nature communications*, 6(1):5968, 2015.</footnote>

<footnote>[16] Leyou Zhang and Xiaoming Mao. Finite-temperature mechanical instability in disordered lattices. *Physical Review E*, 93(2):022110, 2016.</footnote>

<footnote>[17] Cheng-Tai Lee and Matthias Merkel. Generic elasticity of thermal, under-constrained systems. arXiv preprint arXiv:2304.07266, 2023.</footnote>

<footnote>[18] Cheng-Tai Lee and Matthias Merkel. Partition sum of thermal, under-constrained systems. arXiv preprint arXiv:2304.07264, 2023.</footnote>

<footnote>[19] Reza Farhadifar, Jens-Christian Röper, Benoit Aigouy, Suzanne Eaton, and Frank Jülicher. The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. *Current biology*, 17(24):2095–2104, 2007.</footnote>

<footnote>[20] Silvanus Alt, Poulami Ganguly, and Guillaume Salbreux. Vertex models: from cell mechanics to tissue morphogenesis. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 372(1720):20150520, 2017.</footnote>

<footnote>[21] Hisao Honda. Geometrical models for cells in tissues. *International review of cytology*, 81:191–248, 1983.</footnote>

<footnote>[22] Sangwoo Kim, Yiliang Wang, and Sascha Hilgenfeldt. Universal features of metastable state energies in cellular matter. *Physical review letters*, 120(24):248001, 2018.</footnote>

<footnote>[23] Eva-Maria Schoetz, Marcos Lanio, Jared A Talbot, and M Lisa Manning. Glassy dynamics in three-dimensional embryonic tissues. *Journal of The Royal Society Interface*, 10(89):20130726, 2013.</footnote>

<footnote>[24] Thomas E Angelini, Edouard Hannezo, Xavier Trepat, Manuel Marquez, Jeffrey J Fredberg, and David A Weitz. Glass-like dynamics of collective cell migration. *Proceedings of the National Academy of Sciences*, 108(12):4714–4719, 2011.</footnote>

<footnote>[25] Monirosadat Sadati, Nader Taheri Qazvini, Ramaswamy Krishnan, Chan Young Park, and Jeffrey J Fredberg. Collective migration and cell jamming. *Differentiation*, 86(3):121–125, 2013.</footnote>

<footnote>[26] Jin-Ah Park, Lior Atia, Jennifer A Mitchel, Jeffrey J Fredberg, and James P Butler. Collective migration and cell jamming in asthma, cancer and development. *Journal of cell science*, 129(18):3375–3383, 2016.</footnote>
[27] Linda Oswald, Steffen Grosser, David M Smith, and Josef A Käs. Jamming transitions in cancer. *Journal of physics D: Applied physics*, 50(48):483001, 2017.

[28] Jin-Ah Park, Jae Hun Kim, Dapeng Bi, Jennifer A Mitchel, Nader Taheri Qazvini, Kelan Tantisira, Chan Young Park, Maureen McGill, Sae-Hoon Kim, Bomi Gweon, et al. Unjamming and cell shape in the asthmatic airway epithelium. *Nature materials*, 14(10):1040–1048, 2015.

[29] Simon Garcia, Edouard Hannezo, Jens Elgeti, Jean-François Joanny, Pascal Silberzan, and Nir S Gov. Physics of active jamming during collective cellular motion in a monolayer. *Proceedings of the National Academy of Sciences*, 112(50):15314–15319, 2015.

[30] Wenhui Tang, Amit Das, Adrian F Pegoraro, Yu Long Han, Jessie Huang, David A Roberts, Haiqian Yang, Jeffrey J Fredberg, Darrell N Kotton, Dapeng Bi, et al. Collective curvature sensing and fluidity in three-dimensional multicellular systems. *Nature Physics*, 18(11):1371–1378, 2022.

[31] John Devany, Daniel M Sussman, Takaki Yamamoto, M Lisa Manning, and Margaret L Gardel. Cell cycle-dependent active stress drives epithelia remodeling. *Proceedings of the National Academy of Sciences*, 118(10):e1917853118, 2021.

[32] Steffen Grosser, Jürgen Lippoldt, Linda Oswald, Matthias Merkel, Daniel M Sussman, Frédéric Renner, Pablo Gottheil, Erik W Morawetz, Thomas Fuhs, Xiaofan Xie, et al. Cell and nucleus shape as an indicator of tissue fluidity in carcinoma. *Physical Review X*, 11(1):011033, 2021.

[33] Dapeng Bi, JH Lopez, Jennifer M Schwarz, and M Lisa Manning. A density-independent rigidity transition in biological tissues. *Nature Physics*, 11(12):1074–1079, 2015.

[34] Matthias Merkel and M Lisa Manning. A geometrically controlled rigidity transition in a model for confluent 3d tissues. *New Journal of Physics*, 20(2):022002, 2018.

[35] Charlie Duclut, Joris Paijmans, Mandar M Inamdar, Carl D Modes, and Frank Jülicher. Nonlinear rheology of cellular networks. *Cells & development*, 168:203746, 2021.

[36] Sijie Tong, Navreeta K Singh, Rastko Sknepnek, and Andrej Košmrlj. Linear viscoelastic properties of the vertex model for epithelial tissues. *PLoS computational biology*, 18(5):e1010135, 2022.

[37] Xun Wang, Matthias Merkel, Leo B. Sutter, Gonca Erdemci-Tandogan, M. Lisa Manning, and Karen E. Kasza. Anisotropy links cell shapes to tissue flow during convergent extension. *Proceedings of the National Academy of Sciences*, 117(24):13541–13551, June 2020.

[38] Dapeng Bi, Xingbo Yang, M Cristina Marchetti, and M Lisa Manning. Motility-driven glass and jamming transitions in biological tissues. *Physical Review X*, 6(2):021011, 2016.

[39] Daniel M Sussman and Matthias Merkel. No unjamming transition in a voronoi model of biological tissue. *Soft matter*, 14(17):3397–3403, 2018.

[40] Diogo E. P. Pinto, Daniel M. Sussman, Margarida M. Telo da Gama, and Nuno A. M. Araújo. Hierarchical structure of the energy landscape in the voronoi model of dense tissue. *Phys. Rev. Res.*, 4:023187, Jun 2022.

[41] Daniel M Sussman, M Paoluzzi, M Cristina Marchetti, and M Lisa Manning. Anomalous glassy dynamics in simple models of dense biological tissue. *Europhysics Letters*,

121(3):36001, 2018.

[42] Yan-Wei Li, Leon Loh Yeong Wei, Matteo Paoluzzi, and Massimo Pica Ciamarra. Softness, anomalous dynamics, and fractal-like energy landscape in model cell tissues. *Physical Review E*, 103(2):022607, 2021.

[43] Satyam Pandey, Soumitra Kolya, Padmashree Devendran, Souvik Sadhukhan, Tamal Das, and Saroj Kumar Nandi. The structure-dynamics feedback mechanism governs the glassy dynamics in epithelial monolayers. *arXiv preprint arXiv:2306.07250*, 2024.

[44] Elijah Flenner and Grzegorz Szamel. Viscoelastic shear stress relaxation in two-dimensional glass-forming liquids. *Proceedings of the National Academy of Sciences*, 116(6):2015–2020, 2019.

[45] Zhibin Lu, Jiangong Li, Hang Shao, H Gleiter, and Xia Ni. The correlation between shear elastic modulus and glass transition temperature of bulk metallic glasses. *Applied Physics Letters*, 94(9), 2009.

[46] BI Halperin and David R Nelson. Theory of two-dimensional melting. *Physical Review Letters*, 41(2):121, 1978.

[47] Yan-Wei Li and Massimo Pica Ciamarra. Role of cell deformability in the two-dimensional melting of biological tissues. *Physical Review Materials*, 2(4):045602, 2018.

[48] Sadjad Arzash, Anupama Gannavarapu, and Fred C MacKintosh. Mechanical criticality of fiber networks at a finite temperature. *Physical Review E*, 108(5):054403, 2023.

[49] Simone Ciarella, Rutger A Biezemans, and Liesbeth MC Janssen. Understanding, predicting, and tuning the fragility of vitrimeric polymers. *Proceedings of the National Academy of Sciences*, 116(50):25013–25022, 2019.

[50] Philipp Germann, Miquel Marin-Riera, and James Sharpe. Ya|a: GPU-Powered Spheroid Models for Meshchyme and Epithelium. *Cell Systems*, 8(3):261–266.e3, 2019.

[51] Silke Henkes, Kaja Kostanjevec, J. Martin Collinson, Rastko Sknepnek, and Eric Bertin. Dense active matter model of motion patterns in confluent cell monolayers. *Nature Communications*, 11(1):1–9, 2020.

[52] Glenn J Martyna, Mark E Tuckerman, Douglas J Tobias, and Michael L Klein. Explicit reversible integrators for extended systems dynamics. *Molecular Physics*, 87(5):1117–1157, 1996.

[53] Daan Frenkel and Berend Smit. *Understanding molecular simulation: from algorithms to applications*, volume 2. Academic Press San Diego, 2002.

[54] Daniel M Sussman. cellgpu: Massively parallel simulations of dynamic vertex models. *Computer Physics Communications*, 219:400–406, 2017.

[55] Skanda Vivek, Colm P Kelleher, Paul M Chaikin, and Eric R Weeks. Long-wavelength fluctuations and the glass transition in two dimensions and three dimensions. *Proceedings of the National Academy of Sciences*, 114(8):1850–1855, 2017.

[56] Bernd Illing, Sebastian Fritschi, Herbert Kaiser, Christian L Klix, Georg Maret, and Peter Keim. Mermin–wagner fluctuations in 2d amorphous solids. *Proceedings of the National Academy of Sciences*, 114(8):1856–1861, 2017.

[57] Jorge Ramírez, Sathish K Sukumaran, Bart Vorselaars, and Alexei E Likhtman. Efficient on the fly calculation of time correlation functions in computer simulations. *The Journal of chemical physics*, 133(15), 2010.
[58] Jean-Pierre Hansen and Ian Ranald McDonald. *Theory of simple liquids: with applications to soft matter*. Academic press, 2013.

[59] Vitaliy Zhuravlyov, J Goree, Paolo Elvati, and Angela Violi. Finite-size effects in the static structure factor s (k) and s (0) for a two-dimensional yukawa liquid. *Physical Review E*, 108(3):035211, 2023.

## End Matter

## Simulation details

We perform standard NVT simulations [52, 53] of the 2D Voronoi model under periodic boundary conditions using the open-source *cellGPU* package [54]. Temperature is controlled by coupling the system to a Nosé-Hoover thermostat with fictitious degrees of freedom $\zeta_i$. The equations of motion of the standard Nosé-Hoover thermostat are

$$m_i \ddot{\mathbf{r}}_i = \mathbf{f}_i - \zeta m_i \mathbf{v}_i \quad (3)$$

$$\dot{\zeta} = Q^{-1} \sum_{i}^{N} m_i v_i^2 - N_f k_B T, \quad (4)$$

where $m_i$ and $v_i$ are the mass and velocity of particle i, $N_f$ is the number of degrees of freedom in the system, T is the target temperature, and Q is the thermostat parameter that determines its relaxational dynamics. This non-Hamiltonian scheme is one of the most accurate and efficient ways of simulating a constant-temperature ensemble, and we further make use of a formulation with a short chain of fictitious thermostat variables [52]. We then employ for each timestep a standard update of the equations of motion based on a Trotter factorization [53].

We measure the characteristic relaxation time, $\tau_\alpha$, by measuring the decay of the cage-relative self-intermediate scattering function $F_s^{CR}$ [55, 56]. We ensure that all simulations have been equilibrated for at least 10 times longer than the measured $\tau_\alpha$. We report results for systems of N = 4096 cells, and have checked that increasing this to N = 32768 does not change our results.

## Quantification of the plateau shear modulus $G_p$

At each simulation step we analytically calculate the shear stress $\sigma_{xy} \equiv (dE/d\gamma)/A$, where $\gamma$ is the shear strain deforming the periodic box, and A is the total system area. We then compute $G(t)$ via the stress autocorrelation function

$$G(t) = \frac{A}{k_B T} \langle \sigma_{xy}(t_0) \sigma_{xy}(t_0 + t) \rangle, \quad (5)$$

where $\langle ... \rangle$ represents an average over $t_0$ for each simulation, with statistics accumulated using a multiple-tau correlator method [57].

We verified these results by independently measuring $G(t)$ from the response to a small shear step $\Delta\gamma$ at $t = 0$ via $G(t) = \Delta\sigma_{xy}(t)/\Delta\gamma$. Specifically, we imposed a pure shear deformation on the orthonormal shape of the primary unit cell, numerically measured the derivative of energy with respect to this pure shear deformation, and converted the pure shear modulus to the simple shear modulus, G.

We obtain $G_p$ and $\tau_s$ by fitting $G(t)$ to a stretched exponential form, $G(t) \approx G_p e^{-(t/\tau_s)^\beta}$ at intermediate and late times. More precisely, we fit data from 10 simulations for $t \in [t_i, t_f]$ for each set of parameters ($p_0$, T); $t_i$ is chosen as the earliest time that an approximate plateau appears. If there is no such plateau, we vary $t_i$ across a range to ensure a qualitatively good fit encompassing as many data points as possible. We choose $t_f$ as the time at which the magnitude of $G(t)$ is comparable to either the thermal fluctuations or the statistical error of our autocorrelation method [57].

## Unusual cage dynamics at short times

We note that the shear relaxation modulus $G(t)$ shown in Fig. 2 has substantially more structure at short times t than is usually seen in glassforming systems. This is a direct reflection of the unusually large cage fluctuations, in which cells undergo surprisingly large displacements without exchanging neighbors, that characterize single-particle motion in these models. This can also be seen in the unusual behavior of the scattering function itself. As shown in Fig. 6, one sees that the initial plateau of $F_s^{CR}$ is itself substantially lower than in particulate systems, indicating that a large number of cells move substantially but not so far that they escape their cage and keep diffusing. Furthermore, there is evidence of a secondary re-

<figure><img src="image_7.png" /><figcaption>FIG. 6. The cage-relative scattering function $F_s^{CR}$ shows unusual cage dynamics. Curves correspond to the same systems ($p_0$ and T) as in Fig. 2, and the scattering function is measured at a wavevector corresponding to the first peak of the structure factor.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Scaling of perimeter-related contributions to the energy with temperature. The main plot shows the behavior of $X_t^*$ (the isotropic cell tensions) and $X_p^*$ (the average deviation of cell perimeters from their preferred value) at different values of $p_0$. The inset shows the probability density function of standardized cell perimeters (i.e., $z = (p - \mu)/\sigma$) — as shown, these distributions are extremely close to Gaussian.</figcaption></figure>

relaxation process (the minima and non-monotonicity after the initial decay) indicative of “cage-rebounding” dynamics. This process itself has a characteristic temperature dependence.

## Quantification of the bulk modulus K

We quantify the bulk modulus $K$ by measuring the low-wavevector behavior of the static structure factor, $S(\mathbf{q}) \equiv \frac{1}{N} \sum_{j=1}^{N} \sum_{k=1}^{N} e^{-i\mathbf{q} \cdot (\mathbf{r}_j - \mathbf{r}_k)}$, where $j, k$ index the cells. The bulk modulus is $K_s = Nk_B T/AS(q \to 0)$ [58, 59]. Fig. 2(b) shows representative plots of $S(q)$ across a range of $T$ for $p_0 = 3.825$. For most of our sys-

tems there is an unambiguous low-$q$ plateau, although for the lowest temperatures for $p_0 = 3.8$ there is a non-monotonicity at small $q$ (Fig. 2(b)). We thus use the minimum value attained by the structure factor, $S(q_{min})$, as a lower bound for $S(q \to 0)$ at those temperatures.

We verified these results by directly measuring $K$ from the response to an isotropic compression strain of $10^{-3}$, computing the bulk modulus as $K_N = -AdP/dA$, for $P = -(\sigma_{xx} + \sigma_{yy})/2$. The greatest deviations between these two measures are in the regime where $S(q)$ has these low-wavevector non-monotonicities and the structure begins hexatically ordering.

## Temperature dependence of cell perimeters and isotropic tension

Intuitively, previous work by one of us suggests that the “entropic elasticity” behavior of the shear modulus stems from changes in the accessible volume of phase space as the perimeter increases above $p_0^*$. In this regime, higher temperatures favor larger mean values of the actual cell perimeters, which in turn imply monotonically increasing perimeter tensions. These tensions are in turn what rigidifies the system. We present data for this in Fig. 7, which shows the scaling of both $X_p^* = \langle p \rangle - p_0$ and the isotropic cell tensions, $X_t^* = (\sigma_{xx} + \sigma_{yy})/2$ as a function of temperature. We find good agreement with the theoretical expectations based on generic underconstrained systems [17, 18], and a high correlation between this measurement and the data on $G_p(T)$ itself presented in Fig. 3. Although the plateau shear modulus is dominated by perimeter contributions, we note that we do measure non-zero contributions from the area elasticity to the stress autocorrelation function.