A model of 3D confluent tissue behaves as an under-constrained glass

Abstract

The dynamics of glassy materials slows down upon cooling, typically showing either Arrhenius or super-Arrhenius behavior. However, it was recently shown that 2D cell-based models for biological tissues can be continuously tuned between Arrhenius and sub-Arrhenius dynamics. In previous work, using the 2D Voronoi model, we proposed that such atypical dynamical behavior could be a generic feature of the broad class of mechanically under-constrained materials. Yet, our earlier study had left two important points open: (1) many 2D systems are affected by long-wavelength fluctuations and the 2D melting scenario, and (2) the 2D Voronoi model sits exactly at the isostatic point, making it a marginal case rather than a strictly under-constrained one. Both points complicate the interpretation of our 2D Voronoi model results and their generalization to other systems; to remedy this, here we use large-scale simulations to study the glassy behavior of the 3D extension of the Voronoi model. We first show that the structural relaxation time τ_α of the 3D Voronoi model can be tuned between sub-Arrhenius and Arrhenius behavior, like the 2D Voronoi model. We then establish that the four-point susceptibility, the structure factor, and the model's mechanical properties all display trends consistent with the 2D Voronoi model. These results provide strong evidence that sub-Arrhenius glassy dynamics are a generic feature of under-constrained materials across dimensions. Our work thus broadens the class of disordered materials known to have highly unusual glassy phenomenology.

---

1 Introduction

Understanding the origin of glassy dynamics remains one of the central challenges in condensed matter physics. Upon cooling or compression, a wide variety of liquids and soft materials exhibit a dramatic slowdown of relaxation, often accompanied by increasing dynamical heterogeneity and the emergence of an intermediate-time plateau in time correlation functions.^1–4 The fast relaxation processes towards this plateau are called beta relaxation processes, while the much slower relaxation processes away from the plateau are called alpha processes. In conventional glass formers, the behavior of the alpha relaxation time is usually described as either Arrhenius or super-Arrhenius. Arrhenius behavior refers to an alpha relaxation time that scales as τ_α ∼ exp(ΔE/T), which is interpreted as arising from structural relaxation processes that require crossing approximately temperature-independent energy barriers of size ΔE. Such Arrhenius behavior is often associated with strong glass formers. Meanwhile, fragile glass formers display super-Arrhenius behavior, in which a temperature decrease leads to a relaxation time that grows faster than expected from Arrhenius scaling.³ This corresponds to an increase of the effective energy barriers as temperature is lowered. However, the universality of this trend has recently been challenged by the discovery of models that exhibit qualitatively different, sub-Arrhenius dynamics.^5,6 In such systems, the relaxation time grows more slowly than the Arrhenius expectation upon cooling, implying an effective energy barrier that decreases as temperature is lowered. This behavior is anomalous from the perspective of conventional glass formation, and exploring such departures is essential to build a more complete picture of glass formation and of the possible diversity of disordered dynamics.

Often, the glassy behaviors of materials can be closely linked to their mechanical properties. In the zero-temperature limit, the mechanical properties of a system can be classified using the Maxwell constraint-counting criterion, which distinguishes three regimes. Systems that are over-constrained possess more constraints than degrees of freedom and are typically rigid. In contrast, under-constrained systems have fewer constraints than degrees of freedom and are mechanically unstable in the absence of additional interactions. Between these two limits lie isostatic systems, where the number of constraints and degrees of freedom are balanced, leading to marginal stability and weak solid-like behavior. The role of a glass' mechanical properties in determining its dynamics has been the subject of much study,⁷ and the extent to which a system's zero-temperature athermal mechanics dictates its finite-temperature dynamics remains a subject of active debate.

One system where atypical glassy dynamics has been linked to mechanical properties is the cell-based 2D Voronoi model for biological tissues. This model has been shown to exhibit atypical sub-Arrhenius glassy dynamics,⁶ the extent of which can be tuned continuously.⁸ In recent work, we linked this dynamical behavior to this model's unusual, “entropic”, elasticity, i.e. a strong increase of the elastic moduli with temperature.⁹ We further argued that such mechanical behavior is expected for under-constrained systems,^10,11 which include many of the cell-based tissue models.¹² We thus proposed that the tunable sub-Arrhenius behavior may be a generic property of under-constrained systems.⁹

While compelling, the interpretation of these findings is complicated by two aspects. First, the model is two dimensional, and it is well-understood that long wavelength fluctuations play a vital role in two dimensions,¹³ complicating the analysis of quantities like the mean-squared displacement and self-intermediate scattering function.^14,15 Moreover, the two-dimensional melting scenario introduces additional subtleties through the possible emergence of a hexatic phase,^16,17 restricting the parameter regions where disordered dynamics can be studied. Second, the 2D Voronoi model sits precisely at the isostatic constraint-counting transition, with equally balanced degrees of freedom and energetic constraints. Consequently, while its athermal mechanical properties still reflect its proximity to under-constrained tissue models, it does not show an athermal floppy-rigid transition like those other models.¹⁸

To test whether the anomalous dynamical properties observed in the 2D Voronoi model are found in truly under-constrained materials, free from the confounding effects of hexatic ordering transitions, we study the 3D Voronoi model, which is extensively under-constrained.¹⁹ As a consequence, like many cell-based tissue models,^12,18,20,21 the athermal 3D Voronoi model exhibits a rigidity transition that is tuned by a target cell surface area s₀: at zero temperature, the model is rigid for s₀ ≲ 5.4 and floppy for s₀ ≳ 5.4.¹⁹

Below, we systematically characterize the 3D Voronoi model at finite temperature, including its dynamical behavior, structural properties, and mechanical response across a range of s₀ values. We show that the relaxation time exhibits sub-Arrhenius behavior whose extent increases with s₀, directly analogous to what has been observed in two dimensions. The self-intermediate scattering function displays universal scaling when rescaled by the α-relaxation time, and the growth of the four-point susceptibility is strongly suppressed at large s₀. The static structure factor reveals a weak but systematic temperature dependence of its first peak, reminiscent of vitrimeric systems in which static correlations are tied to anomalous relaxation. Underpinning these dynamical anomalies, we find that the isotropic tension and shear modulus exhibit a monotonic increase with temperature. The isotropic tension also has a crossover from T^1/2 to T¹ scaling, consistent with theoretical predictions for generic under-constrained systems.^10,11 Together, these results establish that the unconventional glassy dynamics and mechanics of the Voronoi model persist in three dimensions. This further supports the idea that sub-Arrhenius behavior is a generic property of under-constrained systems.

2 Methods

2.1 Model

The 3D Voronoi model describes a biological tissue densely packed with N cells¹⁹ (Fig. 1 inset). A given state of the model is defined by a set of cell center positions r_i in 3D space. Cell shapes are defined by a Voronoi tiling of these centers. Forces on cell positions, f_i = −∂e/∂r_i, are defined using an effective energy e, which reads in dimensionless form:

Fig. 1 The 3D Voronoi model exhibits tunable sub-Arrhenius relaxation dynamics. Angell plot of the α-relaxation time, τ_α, for four values of s₀. Each curve represents a different s₀ from high to low (dark red to light blue, s₀ = {5.5, 5.425, 5.39, 5.37}), and for each curve, T_g is determined separately as the temperature where τ_α = 10⁴. The inset shows an energy-minimized configuration of the 3D Voronoi model at s₀ = 5.5 with N = 8 cells, whose Voronoi centers are shown as gray spheres.

The sum is over all cells i, the variables s_i and v_i denote cell surface areas and volumes, respectively, s₀ is the target surface area, and k_s,i is a surface area modulus. We choose the unit of length to be the cube root of the average cell volume, which we furthermore set to unity. This energy is analogous to how the 2D Voronoi model is defined,²² with harmonic constraints on the cell perimeter and cross-sectional area promoted to constraints on the cell surface area and volume. Like in many cell-based tissue models, the form of the energy is motivated by a low-order expansion of geometric cell properties in terms of target values.^23–25 Note that this energy does not correspond to central force interactions, but rather many-body forces that couple the geometry and topology of the cellular material, which together generate the rigidity of the simulated tissue.^18–20,26 We simulate this model using 3D periodic boundary conditions in a system with total volume N.

Initial work studying the 3D Voronoi model focused on the athermal limit of monodisperse systems,¹⁹ but just as in particulate settings there are model parameter regimes in which monodisperse Voronoi packings have a tendency to crystallize. To maintain our focus on the disordered regime, we use a bidisperse mixture of cells whose surface area stiffness k_s varies with cell type. Here we choose a 80 : 20 mixture of cell types A and B with k_s,A/k_s,B = 1 : 0.1. All cells have the same s₀, and have equal masses (which we take to be unity). We note that we also explored bidisperse mixtures with identical k_s and cell-type-dependent s₀ or cell-type-dependent target cell volume, and we found qualitatively identical results. However, we found that stiffness bidispersity let us explore a somewhat broader range of parameter space without observing crystallization.

To study finite-temperature equilibrium properties of this model, we perform standard NVT simulations,^27,28 with temperature implemented via a standard Nosé–Hoover chain to sample the canonical ensemble. The cell masses are set to m = 1, and the thermostat mass is set to Q = NdT,²⁸ where N is the number of cells, d is the number of dimensions, and T is temperature. The equations of motion are integrated using a Nosé–Hoover chain velocity-Verlet scheme, following the reversible Trotter-splitting algorithm described in Appendix E.2 of ref. 27. The integration time step is dt = 10⁻²τ, where τ is the standard dimensionless simulation unit of time for this model. We use this equilibrium thermostatting scheme because our goal is to characterize the canonical finite-temperature glassy behavior of the 3D Voronoi model, rather than to model the physical overdamped dynamics of real cells. Although Brownian or overdamped dynamics would be more appropriate for modeling cellular motion at low Reynolds number, the equilibrium structural and mechanical quantities studied here are ensemble properties. The Nosé–Hoover dynamics also facilitates direct comparison with prior equilibrium studies of glass-forming models.

After verifying that our conclusions are not qualitatively affected by changing the scale of the system, we have focused the bulk of our numerical work on simulations of N = 512 cells. We varied s₀ and temperature T, and obtained 10 equilibrium trajectories for each (s₀,T) point. To this end, for any given value of s₀, we first equilibrated at T_init = 0.1 (for which the characteristic relaxation time τ_α is on the order of 1 − 10τ). We then saved configurations after every 10⁴τ, which is large enough to ensure that the configurations are independent from each other. These then serve as initial configurations for the simulations with the actual temperature T. From these initial configurations, we apply an initial thermalization period at the target temperature T. We use preliminary simulations to measure τ_α at the target temperature, thermalize for at least 10τ_α, and then begin recording trajectory data for another 10τ_α, which we analyze to collect the statistics reported below.

Given the high computational cost of simulating this model, our results primarily focus on four representative s₀ values: s₀ ∈ {5.5, 5.425, 5.39, 5.37}; these values correspond to surface areas well above that of a densely packed crystalline configuration,¹⁹ and are represented by filled symbols in all plots. For each s₀, we select temperatures such that the resulting structural relaxation times span from 10τ to 10⁴τ. In order to further study the s₀-dependence of the mechanical properties, we complement these simulations with additional ones with s₀ ∈ {5.475, 5.45, 5.4, 5.35, 5.32, 5.31, 5.3, 5.29, 5.26, 5.2} but over a more limited temperature range corresponding to τ_α ∼ 10…10³τ; these data – which include surface areas below that of the densest known crystalline packings – are represented by hollow symbols in all plots.

2.2 Dynamical observables

We characterize relaxation using the self-intermediate scattering function,where r_j(t) is the center position of cell j at time t, k = |k|, and the brackets denote an average over initial times, trajectories, and wavevector directions. Unless otherwise stated, we choose k = k_max, where k_max corresponds to the first peak of the static structure factor.

The alpha relaxation time τ_α characterizes the long-time structural relaxation of the system. Operationally, we extract τ_α from the long-time decay of F_s(k_max,t) using the fitting procedure described in Appendix A. The beta relaxation time τ_β characterizes the earlier relaxation process associated with the approach to the intermediate-time plateau. We extract τ_β separately by fitting the early-time decay of F_s(k_max,t), again as described in Appendix A.

We quantify dynamical heterogeneity using the four-point dynamical susceptibility associated with fluctuations of the self-intermediate scattering function,

χ₄(t) = N[〈F_s(k,t)²〉 − 〈F_s(k,t)〉²]. (3)

3 Results

3.1 Dynamical properties

We observe tunable sub-Arrhenius behavior also in the 3D Voronoi model. Fig. 1 shows the dependency of the relaxation time τ_α on temperature T in an Angell plot for several values of the control parameter s₀; these relaxation times were extracted from the behavior of the self-intermediate scattering function, F_s(t), as detailed in Appendix A. In all cases, this representation reveals sub-Arrhenius (i.e. concave) behavior. This behavior depends systematically on s₀: systems with larger s₀ exhibit more pronounced sub-Arrhenius scaling, while decreasing s₀ drives the dynamics closer to an Arrhenius form. This result demonstrates that the 3D Voronoi model inherits the same key feature identified previously in two dimensions: unusual glassy dynamics whose details are tuned continuously via the target shape parameter of the model.^8,29

The top panel of Fig. 2 shows the self-intermediate scattering function, F_s(k,t), evaluated at the wavevector corresponding to the first peak (k = k_max), of the static structure factor for different values of s₀ and T. We choose a representation in which the time axis is scaled by our measured τ_α for each (s₀,T) pair. Remarkably, all curves from different (s₀,T) pairs collapse onto a single master curve in the α-relaxation regime. This collapse demonstrates that, despite quantitative differences in relaxation time, the underlying α-relaxation process is universal across the range of s₀ values and temperatures studied.

Fig. 2 α- and β-relaxation from the self-intermediate scattering function. (top) F_s(t) from the trajectories across all (s₀,T) parameter pairs, plotted with time scaled by τ_α, showing collapse of all curves in the α-relaxation regime. The black line is the fitting result. The inset shows F_s(t) with the time axis scaled by τ_β. (bottom) The measured τ_β as a function of T for several values of s₀; the dotted black line is a guide to the eye with slope −1/2. The inset shows the ratio of the α and β relaxation times as a function of temperature; a clear cross-over is observed from a regime at large s₀ in which the time scales have the same temperature dependence to a more normal regime at low s₀ in which τ_α grows much more rapidly as the temperature is decreased.

This representation also makes clear, though, that there are systematic shifts in how the plateau in the two-step relaxation process is approached. We separately measure this initial relaxation time, τ_β, by fitting the initial relaxation of F_s to a stretched exponential decay with the same method as detailed in Appendix A. We show in the inset that this process, too, shares a common functional form across state points. We thus plot the behavior of this faster relaxation process for a range of (s₀,T). As shown in the bottom panel of Fig. 2, we find that in this model the beta relaxation is much more sensitive to temperature than either standard particulate glass formers or the 2D Voronoi model, with a clear dependence across all s₀ values studied. We are not aware of other glassforming models that show such a range of behavior in the competition between the alpha and beta relaxation process scales.

An unscaled representation of representative F_s(k_max,t) curves on a logarithmic time axis, together with the corresponding mean-squared displacement (MSD), is provided in Appendix B. This additional visualization shows the two-step relaxation process directly, without rescaling time by either τ_α or τ_β. The fitting forms used to extract τ_α and τ_β, together with the corresponding stretching/compression exponents, are given in Appendix B.

Fig. 3 shows the four-point dynamical susceptibility χ₄(t) at different temperatures for two representative values of the control parameter, s₀ = 5.5 (bottom) and s₀ = 5.39 (top). The peak height of χ₄(t) is often interpreted as a measure for the magnitude of dynamical heterogeneities in the system.³⁰ For s₀ = 5.5, the peak height of χ₄(t) is essentially independent of temperature, but for s₀ = 5.39 the peak height of χ₄(t) increases systematically upon cooling. This behavior closely parallels previous observations in two dimensions: at higher s₀, dynamical heterogeneity is strongly suppressed, while decreasing s₀ restores the conventional growth of heterogeneous dynamics upon cooling.⁸ The ability of the Voronoi model to tune not only the global relaxation timescale but also the degree of dynamical heterogeneity highlights the unusual character of glass formation in this system.

Fig. 3 The temperature dependence of dynamical heterogeneity is controlled by s₀. Four-point dynamical susceptibility χ₄(t) for s₀ = 5.39 (top) and s₀ = 5.50 (bottom) at different temperatures (temperature decreases from dark red to light blue).

Taken together, the behaviors of structural relaxation and dynamical heterogeneities in the 3D Voronoi model are very close to what we had observed in 2D.^8,9 This suggests that they are robust properties of the Voronoi model, rather than an artifact of dimensionality.

3.2 Structural properties

Fig. 4 shows the static structure factor S(k) for two representative values of the control parameter, s₀ = 5.39 (top) and s₀ = 5.50 (bottom), where the main panels display S(k) over the full range 0 ≤ k ≤ 20. We find that the wavevector k corresponding to the first peak of S(k) is essentially independent of both temperature and s₀, indicating that the underlying length scale of the packing is fixed.

Fig. 4 The temperature dependence of peak height of static structure factor is controlled by s₀. Static structure factor S(k) for s₀ = 5.39 (top) and s₀ = 5.50 (bottom) at different temperatures (temperatures decrease from red to blue). The inset in the top panel magnifies the region corresponding to the first peak in the main panel; the inset in the bottom panel shows the temperature dependence of the first peak height of S(k) for different s₀.

However, the behavior of the height of the first peak, S_m, strongly depends on s₀ (Fig. 4, bottom inset). For s₀ = 5.39, S_m grows systematically as temperature decreases (Fig. 4, both insets), resembling the behavior of conventional glassformers. In contrast, at s₀ = 5.5 the structure factor is nearly temperature independent, with no noticeable change in S_m upon cooling.

This latter trend mirrors observations made in simulated vitrimers, a distinct physical system that also exhibits sub-Arrhenius scaling tuned by density.⁵ In that system, a crossover from sub-Arrhenius to Arrhenius dynamics was similarly linked to the temperature dependence of the first peak in S(k), a connection further supported by schematic mode-coupling theory. The parallels between our results and those vitrimer studies suggest that the temperature dependence of the first peak of S(k) may strongly correlate with the character of glassy dynamics across disparate physical models.

3.3 Mechanical properties

In order to link the observed τ_α behavior (Fig. 1) to mechanical properties, we first tried to quantify an intermediate-time plateau shear modulus G_p, as we had done for the 2D Voronoi model before.⁹ We first calculate the shear stress Σ_xy at each time point, and compute the time-dependent shear modulus G(t) as the stress autocorrelation function, with statistics accumulated using a multiple-tau correlator method.³¹

In Fig. 5 (top inset) we plot G(t) for s₀ = 5.425 and several temperatures. Qualitatively, the behavior matches that of the 2D model, with the expected initial fast relaxation followed by a secondary feature whose magnitude increases as the temperature is raised. This is, however, a notoriously noisy and data-intensive way of extracting mechanical properties.³¹ For the 3D Voronoi model, the data were too noisy to reliably fit the plateau to a stretched exponential form to extract a precise G_p.

Fig. 5 The isotropic tension increases with temperature and matches the theoretical prediction for generic under-constrained systems. (top) Isotropic tension as a function of temperature for several s₀, showing a transition from T^1/2 to approximately linear T¹ scaling. Inset: shear relaxation modulus G(t) for s₀ = 5.425, displaying an intermediate-time plateau that decreases monotonically with temperature (from red to blue, temperature decreases). (bottom) When plotting isotropic tension Σ_iso versus isotropic strain ε, a rescaling by leads to a collapse of the curves for the different state points. The dashed lines correspond to eqn (12) with parameter values as given in the main text.

To study a proxy for material rigidity, we use isotropic tension Σ_iso, defined as the trace of the stress tensor. Indeed, in generic under-constrained systems with fixed connectivity and without shear strain, isotropic tension Σ_iso is proportional to the infinite-time shear modulus G.^10,11 It is thus a reliable measure of the system's rigidity while being much more accessible to numerical simulations. The top panel of Fig. 5 shows the temperature dependence of Σ_iso for several values of s₀. As in two dimensions, Σ_iso displays an unusual temperature dependence: rather than decreasing modestly upon heating, it increases monotonically with temperature. At high temperatures, the tension grows approximately as T^1/2, while at low temperatures it follows a nearly linear T¹ dependence. This crossover is consistent with what we previously found in the 2D Voronoi model,⁹ further supporting the robustness of this unusual mechanical response.

The isotropic tensions we measure are quantitatively consistent with a generic theory for under-constrained materials with fixed connectivity.^10,11 Ref. 10 and 11 express the elastic material properties, including Σ_iso, in terms of temperature T and isotropic strain ε as shown in eqn (12) in Appendix C. Moreover, isotropic strain ε can be mapped to the target cell surface area using a length rescaling:¹¹ , where denotes the athermal transition point. The bottom panel of Fig. 5 shows the corresponding scaling plot, which leads to an excellent collapse of curves for different state points (s₀,T).

From a fit of eqn (12) to our simulation data (dashed line in Fig. 5 bottom) we obtain the values of athermal transition point , energetic elasticity κ_E = 13.5 ± 0.2, and entropic elasticity κ_S = 0.22 ± 0.1 (see Appendix C for details). These values of and κ_S lie within the expected ranges: is slightly below the transition point of 5.41 obtained from infinite-temperature quenches,¹⁹ which we attribute to the annealing effect of our finite-temperature simulations. For the same reason we also expect the effective transition point to depend weakly on both s₀ and T. Finally, ref. 10 shows that in dimensionless units, where is the number of first-order zero modes (see the precise definition in ref. 11). Given that the 3D Voronoi model has 3N degrees of freedom (the cell positions) but 2N constraints (area and volume terms in eqn (1)), it is reasonable to assume up to terms of order one, and thus κ_S ≈ 1/6, close to what we obtain here. It is somewhat harder to make clear predictions for κ_E based on past work, since we have used a mix of several cell types here.

Finally, we ask whether these elastic properties alone can explain the observed behavior of τ_α. Specifically, we ask whether structural relaxation can be understood as the crossing of energy barriers of height ∼ΔE, which is expected to scale linearly with a plateau shear modulus G_p (or, similarly, with isotropic tension Σ_iso) in “shoving”-like models.^7,9 In Fig. 6 we plot log τ_α against the scaled isotropic tension Σ_iso/T^1.18. Across these state points, the majority of data points collapse onto a straight line, and we have confirmed that a plot over Σ_iso/T does a much worse job of collapsing the data. There are some systematic deviations at the lowest temperatures for low-s₀ simulations, where τ_α grows more rapidly than the scaling indicated by the data collapse. However, the most severe systematic deviations occur for large s₀. We hypothesize that this may be a result of the measured τ_α being dominated by the unusually large relaxation time of beta processes, as discussed earlier (Fig. 2 top).

Fig. 6 The mechanics and dynamics are directly connected. Log(τ_α) is plotted against scaled isotropic tension Σ_iso/T^1.18 for many different s₀.

At present we do not have a microscopic theory for the numerical value of the temperature exponent of −1.18, which is different from the exponent of −1 expected from the standard activated picture. We therefore treat it as an empirical collapse parameter rather than a universal exponent, although we note that we found a similar exponent in the 2D Voronoi model.⁹ This empirical scaling remains consistent with sub-Arrhenius relaxation, because the increase of Σ_iso with temperature is sufficiently strong. For example, if Σ_iso ∼ T^a with a ≈ 0.5…1 over the relevant temperature range, then Σ_iso/T^1.18 ∼ 1/T^c, where c = 1.18 − a ≈ 0.18…0.68, which varies more weakly with inverse temperature than the Arrhenius form 1/T. Thus, although the exponent −1.18 introduces an additional temperature dependence beyond the simplest shoving-like form, it does not restore conventional Arrhenius behavior. Instead, it suggests that the mechanical rigidity is the dominant, but not the only, temperature-dependent contribution to the effective relaxation barrier.

Taken together, like in the 2D Voronoi model,⁹ also in the 3D version, much of the sub-Arrhenius behavior of τ_α can be understood from an increase of mechanical rigidity with temperature.

4 Discussion

In this work, we have demonstrated that the 3D Voronoi model exhibits tunable sub-Arrhenius relaxation dynamics governed by its unusual mechanical properties. By moving from the isostatic 2D limit to an extensively under-constrained 3D system, our results establish that anomalous glassy dynamics are not a dimensional artifact or a consequence of 2D long-wavelength fluctuations, but rather a robust, generic property of under-constrained materials.

The most striking evidence for this is the direct connection between the system's structural relaxation and its mechanical rigidity. Simple elastic models of glass formation posit that the activation energy for structural relaxation ΔE scales with a local, intermediate-time shear modulus G_p.^7,32 Because extracting G_p in 3D is notoriously noisy, we utilized the isotropic tension Σ_iso as a robust proxy for rigidity. Remarkably, plotting log τ_α against the scaled isotropic tension, Σ_iso/T^1.18, yields a collapse of data across a broad range of state points. While a standard shoving model would predict a T⁻¹ scaling, the empirical T^−1.18 dependence confirms that the sub-Arrhenius behavior is primarily driven by a mechanical activation barrier that stiffens at higher temperatures, with the anomalous exponent likely reflecting additional temperature-dependent structural correlations.

Our analytical and structural analyses further illuminate the nature of this under-constrained glassy state. By applying an analytical theory for the elastic properties of thermal, under-constrained systems to our numerical data, we identified a rigidity transition at . This value is lower than the previously reported athermal transition derived from infinite-temperature quenches, likely reflecting the annealing effects of our finite-temperature simulations, though it remains above the crystalline Weaire–Phelan limit (s₀ = 5.288).³³

Moving away from this transition point dramatically alters the dynamical landscape. In the high-s₀ regime, the energy landscape becomes exceptionally shallow; β-relaxation processes become progressively slower and compete with the measured α-relaxation time, while dynamical heterogeneity is largely suppressed. Furthermore, the static correlations in this regime captured by the first peak of the structure factor mirror the atypical temperature independence seen in simulated vitrimers. Conversely, as s₀ decreases toward the rigid regime, the model recovers a more conventional glass-forming phenomenology with pronounced dynamical heterogeneities and growing static structural peaks upon cooling, all while maintaining a universal underlying α-relaxation process. These qualitative features are in excellent agreement with the behavior reported in the 2D Voronoi model.^9,29 In the 2D Voronoi model this qualitative transition has been connected with a change in the underlying fractal structure of the model's energy landscape;^29,34 it remains to be seen whether a similar connection also holds in the 3D model.

Ultimately, the 3D Voronoi model provides a minimal, highly controllable platform for exploring the origins of anomalous glassy dynamics. Our findings broaden the class of disordered materials known to exhibit highly unusual glassy phenomenology. Comparing the mechanics and dynamics of the Voronoi model with other systems (such as vertex models,⁶ vitrimers,⁵ or low-density hard spheres³⁵) may soon allow us to uncover a unifying physical framework for sub-Arrhenius disordered dynamics.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including representative raw trajectories and all of the processed data needed to make the plots, will be made available as a Zenodo repository upon publication.

Appendices

Appendix A: Bayesian inference to extract τ_α

Extracting the relaxation time from stretched exponential fits, A exp[−(t/τ_α)^b], to F_s(t) is notoriously delicate. Because directly fitting this functional form requires a nonlinear least-squares optimization that is highly sensitive to the initial parameter guess, we employed a standard Bayesian inference approach to determine the full posterior distributions of the fitting parameters. Within this framework, we introduce a hierarchical structure that couples parameters across all state points. Crucially, motivated by the excellent collapse of the α-relaxation master curves shown in Fig. 2, we enforce common behavior in A and b while allowing the temperature-dependent relaxation times τ_α to vary independently.

In order to estimate parameters for the relaxation of the self-intermediate scattering function, we performed Bayesian inference using the PyMC probabilistic programming framework and the no-U-turn sampler (NUTS) algorithm.³⁶ The hierarchical model for the stretched exponential decay treated the amplitude A and the stretching exponent b as global parameters, while the relaxation time τ_α was allowed to vary with temperature and s₀. All model parameters were constrained to their physical domains by reparameterization. The stretching exponent b ∈ (0,1) and the amplitude A ∈ (0,1) were represented as sigmoid transforms of unconstrained normal random variables,

q_b ∼ [scr N, script letter N](0,1), b = sigmoid(q_b), (4)

q_A ∼ [scr N, script letter N](0,1), A = sigmoid(q_A). (5)

This choice places weakly informative priors on b and A, centered around 0.5 with broad support over (0,1).

For the temperature-dependent relaxation times τ_α, we worked with logarithmic variables,

Log τ_α ∼ [scr N, script letter N](μ_τ,σ_τ²), (6)

with μ_τ = log(10) and σ_τ = 5, providing wide support over several decades in τ. The exponential transform ensures τ_α > 0. The likelihood is constructed using a Gaussian probability with the mean and variance computed directly from the distribution of the 100 independent F_s trajectories at each state point.

Posterior sampling was carried out using four independent chains, with 2000 warm-up iterations and 2000 sampling iterations per chain, yielding approximately 8000 posterior draws. Convergence was assessed using the potential scale reduction statistic R̂ ≤ 1.01 and by inspecting effective sample sizes. Posterior predictive checks confirmed that the fitted model reproduced both the decay shape and the variance structure of the measured F_s(t). From the resulting posterior distributions we extracted credible intervals for τ_α(T). This hierarchical Bayesian approach produces stable and statistically consistent fits, quantifies uncertainty in a principled way, and exploits shared physical information across the entire dataset.

Appendix B: Additional visualization of two-step relaxation and fitting parameters

In the main text, the self-intermediate scattering functions are shown after rescaling time by either the alpha relaxation time τ_α or the beta relaxation time τ_β. Those rescaled plots are useful for demonstrating the collapse of the relaxation curves across different state points. However, to more directly visualize the two-step relaxation process, we also show representative unscaled curves in Fig. 7.

Fig. 7 Unscaled self-intermediate scattering function and mean-squared displacement. Main panel: Self-intermediate scattering function F_s(k_max,t) for s₀ = 5.425 at different temperatures, plotted with a logarithmic time axis. Temperature decreases from red to blue. The labels β and α indicate the beta- and alpha-relaxation regimes, respectively. Black dashed curves show the corresponding fits. Inset: mean-squared displacement for the highest and lowest temperatures shown in the main panel. The blue dashed guide lines indicate the short-time ballistic regime, MSD ∼ t², and the long-time diffusive regime, MSD ∼ t. Both F_s(k_max,t) and the MSD show a clear two-step relaxation process with an intermediate plateau-like regime.

The main panel of Fig. 7 shows F_s(k_max,t) for s₀ = 5.425 over a range of temperatures. The curves display a clear two-step decay. At early times, F_s(k_max,t) decreases from its initial value toward an intermediate-time plateau; we identify this part of the relaxation with the beta-relaxation regime. At longer times, the system escapes this plateau and F_s(k_max,t) decays to zero; we identify this final decay with the alpha-relaxation regime. As the temperature decreases, both relaxation processes shift to longer times, while the separation between the early and late relaxation regimes becomes more apparent. The black dashed curves show the fits used to extract the corresponding beta and alpha relaxation times.

The inset of Fig. 7 shows the mean-squared displacement,

for the highest and lowest temperatures shown in the main panel. The dashed guide lines indicate the short-time ballistic regime, 〈Δr²(t)〉 ∼ t², and the long-time diffusive regime, 〈Δr²(t)〉 ∼ t. Between these two limits, the MSD develops a plateau-like regime, corresponding to transient caging. This behavior is consistent with the intermediate-time plateau observed in F_s(k_max,t) and provides an independent visualization of the two-step relaxation dynamics.

For the alpha-relaxation regime, we fit the long-time decay of the self-intermediate scattering function using

where τ_α is the alpha relaxation time and b_α is the stretching or compression exponent. This fit is applied only to the second, long-time stage of the decay.

We fit the beta-relaxation regime separately using the early-time decay toward the plateau. Specifically, we use

where τ_β characterizes the beta relaxation time, and b_β is the corresponding stretching or compression exponent. We use the notation b_α and b_β for these exponents to avoid confusion with the beta relaxation process itself.

For the all state points we studied, the fitting parameters are

A_α = 0.6, A_β = 1.0, (10)

b_α = 0.82, b_β = 1.74. (11)

Values of b < 1 correspond to stretched exponential relaxation, while values of b > 1 correspond to compressed exponential relaxation. The fitted curves shown in Fig. 7 demonstrate that the two-step relaxation is well captured by fitting the beta and alpha regimes separately.

Appendix C: Prediction of elastic properties of 3D Voronoi tissue

In generic under-constrained systems without shear strain, the isotropic tension can be expressed as:^10,11where κ_E is the energetic isotropic bulk modulus and κ_S is the entropic isotropic bulk modulus. Moreover, ε is the isotropic strain, measured with respect to the athermal transition point; note that athermal under-constrained systems can be driven across the floppy-rigid transition by applying isotropic strain.^10–12,37 In other words, varying isotropic strain ε is equivalent to varying the target cell surface s₀.^11,12,37 According to the non-dimensionalization in ref. 11:with being the athermal transition point.

Strictly speaking, the analytical results in ref. 10 and 11 were derived for systems lacking a state of self-stress (SSS) in the athermal floppy regime. In contrast, cellular models like the 3D Voronoi model strictly possess a SSS, a direct mathematical consequence of the total system volume being fixed (i.e., the sum of all cell volumes is invariant). However, based on previous analysis of the 3D Voronoi model,¹² the macroscopic constraint restricting the N volume degrees of freedom largely decouples from the surface area constraints. Consequently, the “trivial” bulk modulus arising from the total volume constraint drops out of the dominant low-energy mechanics, allowing us to safely apply the formalism such that κ_E in eqn (12) corresponds to the isotropic bulk modulus contributed exclusively by the cell surface tension terms.

Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant No. DMR-2143815. This research used the Delta advanced computing and data resource which is supported 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.

Notes and references

1. A. Cavagna, Phys. Rep., 2009, 476, 51–124.
2. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259–267.
3. M. D. Ediger, C. A. Angell and S. R. Nagel, J. Phys. Chem., 1996, 100, 13200–13212.
4. W. Kob and H. C. Andersen, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 4626.
5. S. Ciarella, R. A. Biezemans and L. M. Janssen, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 25013–25022.
6. D. M. Sussman, M. Paoluzzi, M. C. Marchetti and M. L. Manning, Europhys. Lett., 2018, 121, 36001.
7. J. C. Dyre and W. H. Wang, J. Chem. Phys., 2012, 136(22), 224108.
8. H. S. Ansell, C. Li and D. M. Sussman, Phys. Rev. E, 2025, 112, L013403.
9. C. Li, M. Merkel and D. M. Sussman, Phys. Rev. Lett., 2025, 134, 048203.
10. C.-T. Lee and M. Merkel, Phys. Rev. Lett., 2024, 133, 268201.
11. C.-T. Lee and M. Merkel, Phys. Rev. E, 2024, 110, 064147.
12. M. Merkel, K. Baumgarten, B. P. Tighe and M. L. Manning, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 6560–6568.
13. E. Flenner and G. Szamel, Nat. Commun., 2015, 6, 7392.
14. S. Vivek, C. P. Kelleher, P. M. Chaikin and E. R. Weeks, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 1850–1855.
15. B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret and P. Keim, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 1856–1861.
16. B. Halperin and D. R. Nelson, Phys. Rev. Lett., 1978, 41, 121.
17. Y.-W. Li and M. P. Ciamarra, Phys. Rev. Mater., 2018, 2, 045602.
18. D. M. Sussman and M. Merkel, Soft Matter, 2018, 14, 3397–3403.
19. M. Merkel and M. L. Manning, New J. Phys., 2018, 20, 022002.
20. D. Bi, J. Lopez, J. M. Schwarz and M. L. Manning, Nat. Phys., 2015, 11, 1074–1079.
21. K. Kim, T. Zhang and J. M. Schwarz, New J. Phys., 2024, 26, 043009.
22. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
23. R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Curr. Biol., 2007, 17, 2095–2104.
24. D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Eur. Phys. J. E:Soft Matter Biol. Phys., 2010, 33, 117–127.
25. S. Kim, Y. Wang and S. Hilgenfeldt, Phys. Rev. Lett., 2018, 120, 248001.
26. D. M. Sussman, J. Schwarz, M. C. Marchetti and M. L. Manning, Phys. Rev. Lett., 2018, 120, 058001.
27. D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Elsevier, 2023.
28. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys., 1996, 87, 1117–1157.
29. Y.-W. Li, L. L. Y. Wei, M. Paoluzzi and M. P. Ciamarra, Phys. Rev. E, 2021, 103, 022607.
30. L. Berthier, G. Biroli, J.-P. Bouchaud and R. L. Jack, Dynamical heterogeneities in glasses, colloids, and granular media, 2011, vol. 150, p. 68.
31. J. Ramrez, S. K. Sukumaran, B. Vorselaars and A. E. Likhtman, J. Chem. Phys., 2010, 133(15), 154103.
32. F. Puosi and D. Leporini, J. Chem. Phys., 2012, 136(4), 041104.
33. D. L. Weaire and S. Hutzler, The Physics of Foams, Clarendon Press, 1999.
34. D. E. Pinto, D. M. Sussman, M. M. Telo da Gama and N. A. Araújo, Phys. Rev. Res., 2022, 4, 023187.
35. L. Berthier and T. A. Witten, Europhys. Lett., 2009, 86, 10001.
36. O. Abril-Pla, V. Andreani, C. Carroll, L. Dong, C. J. Fonnesbeck, M. Kochurov, R. Kumar, J. Lao, C. C. Luhmann, O. A. Martin, M. Osthege, R. Vieira, T. Wiecki and R. Zinkov, PeerJ Comput. Sci., 2023, 9, e1516.
37. C.-T. Lee and M. Merkel, Soft Matter, 2022, 18, 5410–5425.

---