Short-range correlation of stress chains near solid-to-liquid transition in active monolayers

Using a three-dimensional model of cell monolayers, we study the spatial organization of active stress chains as the monolayer transitions from a solid to a liquid state. The critical exponents that characterize this transition map the isotropic stress percolation onto the two-dimensional random percolation universality class, suggesting short-range stress correlations near this transition. This mapping is achieved via two distinct, independent pathways: (i) cell–cell adhesion and (ii) active traction forces. We unify our findings by linking the nature of this transition to high-stress fluctuations, distinctly linked to each pathway. The results elevate the importance of the transmission of mechanical information in dense active matter and provide a new context for understanding the non-equilibrium statistical physics of phase transition in active systems.

Long-range correlations of mechanical information, i.e. force or stress, prevail in both passive glasses [27,28] and most athermal jammed granular systems [29].In such passive systems, the long-range force or stress correlations are attributed to mechanical equilibrium [30][31][32][33].However, in some non-equilibrium systems, long-range mechanical correlations can also emerge in two seemingly disparate systems: (i) active matter [34], where local energy injection results in the breaking of time-reversal symmetry and detailed balance [35], and (ii) boundary-driven shear of granular systems [36].Whether a broad framework to understand these two systems exists remains to be established [37].
Here, using a three-dimensional model of cell monolayers, we provide evidence for the emergence of active stress chains (figure 1a,b) in the monolayers and show stress percolation criticality near the solid-to-liquid-like transition (figure 1c,d).Importantly, we establish this criticality by driving the transition via two distinct paths by independently increasing (i) cell-cell adhesion strength and (ii) active traction forces, well-established axes on the phase boundary for this transition [38,39].Using finite-size scaling analyses, we further demonstrate that critical exponents from each path correspond to the two-dimensional random percolation universality class near the solid-to-liquid transition.Remarkably, this points to the short-range nature of stress correlations.To explain this, we provide a mechanistic basis by linking the universality class mapping to high fluctuations in stress fields with distinct signatures associated with each path.We further discuss the implications for the short-range nature of stress correlations in active monolayers and its relationship with non-equilibrium glasses, granular and biological systems.

Methods
We consider a cellular monolayer consisting of N = 400 cells on a rigid substrate with its surface normal e → n = e → z = e → x × e → y and periodic boundaries in both e → x and e → y , where e → x , e → y , e → z constitute the global orthonormal basis.Each cell i is represen- ted by a three-dimensional phase field, ϕ i = ϕ i x → , t and initialized with radius R 0 .The phase field ϕ i is evolved through Model A [40] type dynamics with an extra advective term, (2.1) where v → i is the velocity of cell i, Γ represents mobility, ℱ is the free energy functional that stabilizes cell interface and includes mechanical properties such as cell stiffness (E) as well as compressibility (μ), and puts a soft constraint on the cell volume [41][42][43][44] around V 0 = (4/3)πR 0

3
. Additionally, the free energy comprises gradient contributions ( ∇ → ϕ) that account for, and distinguish between, cell-cell (ω cc ) and cell-substrate (ω cw ) adhesions, as introduced recently [45], (2.2) , where κ captures repulsion between cell-cell (subscript cc) and cell-substrate (subscript cw) and ϕ w denotes a static phase field representing the substrate.Based on this free energy function, the interior and exterior of cell i corresponds to ϕ i = 1 and ϕ i = 0 , respectively, connected by a diffuse interface of length λ.To resolve the forces generated at the cellular interfaces, we consider an over-damped dynamics,T , where T → i denotes traction [18,46,47] and contains both active and passive contributions, ξ is substrate friction and F → i act = αF ^i pol represents an active polar force that captures the self-propulsion associated with cell polarity, constantly pushing the system out of equilibrium, with α characterizing the strength of the self-propulsion force.The dynamics of cell polarity are introduced based on contact inhibition of locomotion  royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 21: 20240022 [48,49] by aligning the polarity of the cell to the direction of the total interaction force acting on the cell [50].The polarization dynamics is given by ∂ is the counterclockwise angle of cell polarity measured from e → x , F ^i pol = cos θ i , sin θ i , 0 and η is the Gaussian white noise with zero mean, unit variance, D r is rotational diffusivity, ΔΘ i is the angle between F ^i pol and T → i , with positive constant τ pol setting the alignment time scale between them (see electronic supplementary material for simulation parameters).Previously, this model was tuned directly against experiments [51] and specifically chosen given its success in reproducing stress fields around nematic defects [45,47].
In passive, athermal granular systems, the jamming transition is intimately related to the packing fraction (density) of that system.However, for living cells at confluence, the collective behaviour can change from a liquid to a solid-like behaviour at a constant density.A recent study suggests density has a second-order effect on this transition [12].Consequently, here we keep the cell density constant in our simulations, i.e. no cell proliferation nor extrusion, while ensuring confluency.For each considered pathway, we perform large-scale simulations by incrementally increasing the dimensionless cell-cell to cellsubstrate adhesion ratio ω ~= ω cc /ω cw ∈ [0.1, 0.5] for the first, and the dimensionless traction force α ~= ατ pol /ξR 0 , α ~∈ [0, 0.8], for the second pathway.We consider three distinct realizations for each case.As we show next, the considered ranges for ω ~ and α capture the transition from a solid-like to a liquid state.

Solid-to-liquid transition
To characterize the solid-to-liquid-like transition, we quantify the one-dimensional pair correlation function where ρ denotes number density, shown in figure 2 for time-averaged configurations.For low values of ω ~ and α ~, g r shows a peak near r/R 0 ≈ 2, i.e. the diameter of a single cell, indicating a solid-like state.As ω ãnd α ~ increase, this peak gets smaller while another one, at a shorter distance, appears, indicative of a transition.For the time-averaged configurations, the peak r/R 0 ≈ 0 indicates significant overlaps between cells, relative to their initial positions, during the simulation (see figure 1d).For ω ~= 0.5 and α ~= 0.8, the pair-correlation functions exhibit peaks at a distance much smaller than the cell diameter, suggesting a liquid-like state.Thus, the solid-to-liquid transition takes place near ω ~= 0.45 and α ~= 0.70, mindful of the resolution for ω ~ and α ~ spaces.We also quantify the two-dimensional static structure factor that further confirms the identified transition points (see electronic supplementary material for details).
In the three-dimensional phase field model of active cell monolayers, stress chains reminiscent of force chains in passive, athermal granular systems [52,53] are observed (figure 1).In contrast with passive, athermal, convex granular systems, the stress chains in this active model exhibit both compressive and tensile stress components.This can also be the case in non-convex granular systems with the ability to entangle [54,55].To understand the identified transition between solid-liquid states and its link to the spatial organization of mechanical information, we compute a stress field [56,57] σ = σ( x → , t) that encodes both active and passive contributions on a discretized domain, for node i on the complementary stress lattice, (3.1) where a 0 3 is the volume of the unit cell, N i is the number of neighbours of node i on the original lattice (see electronic supplementary material for details), r →ij = x →i − x →j .Given the definition of active traction for cell i, T → i , the coarse-grained stress field contains contributions from both active and passive forces.Next, we compute the time-averaged isotropic stress field tensor, σ ¯x → = (1/n sim ) ∑ t n sim σ x → , t and focus on the normalized isotropic stress field, σ ¯iso x → = σ ¯iso x → /σ ¯max iso with n sim number of simulation time steps, σ ¯max iso the maximum value of the time-averaged isotropic stress field and σ ¯iso x → = 1/3 trσ ¯x → , providing a measure for expansion (σ iso > 0) and compression (σ iso < 0) in the cell layer.
Visual inspection of the fluctuations in isotropic stress fields and the associated patterns show that for both adhesion parameter ω ~ and traction force α ~, as the control parameter is increased and the system approaches a liquid-like state, a more disordered isotropic stress field emerges.We quantify the spatial organization of emerging stress patterns by lowering a threshold and monitoring the connectivity of the time-averaged isotropic stress field sites larger than the threshold, on a square lattice.Then, we extract the critical exponents for the stress percolation by performing finite-size scaling analyses [58,59].
We begin with quantifying the density of the spanning cluster, P p, ℓ , which is the probability of a site belonging to a spanning cluster as a function of the occupation probability, p, and system size, ℓ.The occupation probability, p ∈ [0, 1], corresponds to sites where σ ¯iso x → is greater than 1 − p × 100th percentile of the isotropic stress distribution, including both compressive and tensile stresses.In the large system size, the 'thermodynamic limit', i.e. ℓ → ∞, we expect the following power-law scaling near the percolation probability, p c , for the density of a spanning cluster P p characterized by the critical .A cluster is a set of connected sites, and two sites are considered connected if they are nearest neighbours (eight-point connectivity).Furthermore, we quantify the average cluster size, S p = ⟨s⟩, where s is the size of a cluster.Near the percolation probability, the power-law scaling of S p quantified by critical exponent γ is S p ∼ | p − p c | −γ and critical exponent ν, which describes the power-law scaling of correlation length ξ, i.e. average distance between two sites in the same cluster, is expected to follow Both the density of the spanning cluster P p, ℓ and the average cluster size S p, ℓ are shown in figure 3 at the onset of the transition, i.e. ω ~ (figure 3a,c) and α ~= 0.7 (figure 3e,g) for various system sizes, ℓ.Additionally, figure 3b,d,f and h shows the collapse of P p, ℓ and S p, ℓ when scaled with critical exponents obtained from the finite-size scaling analyses, accounting for a diverging correlation length ξ near the p c .Remarkably, the critical exponents corresponding to the scaling of the time-averaged isotropic stress field at the onset of the solid-to-liquid transition in the active cell layer (table 1) lead to a reasonable collapse and are in close agreement with those from the two-dimensional random percolation universality class [60].Importantly, this agreement is obtained via two independent pathways, i.e. cell-cell adhesion and active traction force, and points to the short-range nature of stress correlations near this transition.

Mechanism
Our results thus far do not explain why short-range correlations that are characteristic of random percolation emerge in the solid-to-liquid transition of the studied active monolayers.To explore this, we analysed the stress field statistics extensively and identified two mechanical signatures associated with each path to transition.For the first path, increasing cellcell adhesion leads to a solid-like to liquid transition.Although this counterintuitive behaviour, called the adhesion paradox, was observed experimentally and captured by a two-dimensional vertex model previously [61], the mechanism responsible for it is yet to be explained.Our analyses suggest increasing cell-cell adhesion results in higher global mean and field fluctuations, quantified by susceptibility, χ( . )= n ⟨ . 2 ⟩ − ⟨ .⟩ 2 with n denoting field dimensions, in the out-of-plane component of the time-averaged stress tensor, σ ¯zz x → = e → z ⋅ σ ¯x → ⋅ e → z .This is shown in figure 4a for σ ¯zz x → = σ ¯zz /σ ¯zz max .This high mean out-of-plane stress, neither accessible in two-dimensional models nor experiments, can lead to local buckling owing to cell shape changes caused by adhesion increase [62] providing a slightly more free volume needed for a transition to a liquid state [63].Moreover, increased fluctuations are consistent with the emergence of short-range stress correlations and random percolation at the onset of the transition.For the second path to transition, increasing active polarity and hence in-plane active traction forces lead to an increase in both the global mean and field fluctuations in the maximum in-plane shear field, ).This is shown in figure 4b for the time-averaged field, σ ¯τ x → , normalized by its maximum, σ ¯τ max , σ ¯τ x → = σ ¯τ/σ ¯τ max .Thus, the transition driven by increasing active traction leads to an increase in both mean and fluctuations of local, maximum shear, from  which emerges stress percolation, near the transition.We note that an important observation here is the encoding of different types of mechanical information in the isotropic stress field.In this vein, our approach has the potential to unify disparate observations in biological systems, such as percolations based on (i) cell connectivity [25] and (ii) edge tension network [64,65] through understanding the transmission of mechanical information.

Discussion
The critical exponents obtained by finite-size scaling analyses of time-averaged isotropic stress fields map the solid-to-liquid transition in active monolayers to a two-dimensional random percolation universality class.Thus, we reveal the short-range nature of stress correlations near this transition, attuned to differing critical exponents for percolation in systems with longrange correlated disorder [66].In this vein, our findings lend further support to the viewpoint that in passive amorphous solids, it is the mechanical equilibrium condition that gives rise to long-range stress correlations [30][31][32][33].Importantly, we establish the criticality of the active solid-to-liquid transition via two distinct paths, which are both well established experimentally [23,61].
A particularly appealing feature of our approach is the implicit manifestation of physical interactions as mechanical information in isotropic stress fields, given the two distinct drives considered in our study.This is important since direct measurement of traction forces and mechanical stress fields is experimentally accessible [46,67] and could be used to confirm the proposed critical behaviour.Our results are particularly illuminating since recent experiments highlight the dominant role of traction forces in deriving the solid-to-liquid transition in epithelial tissues [23].Activity in our study is owing to active polarity.The particular model we use is tuned directly by experiments [51] and specifically chosen given its success in representing cell motility by reproducing stress patterns around nematic defects [45,47].In this vein, a recent study suggests that the active glassy dynamics are mute to the details of such a model [68].
The critical exponents were estimated by extracting squared subsystems of characteristic length ℓ ∈ [20, 100], spanning almost a decade of length scales, with dimensions (ℓ + 1) 2 given the lattice spacing a 0 = 1 ≪ R 0 , from the time-averaged isotropic stress fields with total dimensions of (L x + 1) × (L y + 1) = 103 041 (see electronic supplementary material for finite-size scaling analyses).While our analyses can definitely benefit from more statistics, we are confident in its robustness since we obtain similar critical exponents from two distinct and independent paths towards solid-to-liquid transition.Additionally, we used site percolation on a square lattice, given the nature of our simulations, but it is well known that critical exponents for random percolation universality class are robust to model details [60,66,69], e.g. bond versus site percolation and lattice type.Thus, we do not expect our results to change based on lattice types, e.g.triangular, excluding the Bethe lattice.The finite-size scaling analyses were performed on two-dimensional isotropic stress fields that embed the out-of-plane stress component σ zz .Furthermore, the three-dimensional nature of our modelling approach presents other unique advantages including the ability to tune independently and explicitly for cell-cell and cell-substrate interactions.Most importantly, we can ensure that, within the range of studied parameters, the solid-to-liquid transition is not owing to an extrusion event, altering the confluent state.This cannot be done with two-dimensional modelling approaches.
We also examined potential links to jamming transition in passive granular systems [70,71].Interestingly, a recent study suggests short-range force correlation in jammed granular systems [72], in agreement with our findings.However, this contrasts with other studies that point to long-range correlations in such systems [29,53].This may be owing to the dependence of critical parameters associated with jamming transition on the nature of the global load, e.g.shear versus isotropic compression, as well as the loading rate [73,74].Moreover, our results differ from the critical exponents associated with standard rigidity percolation [75,76], indicating distinct universality classes.In this vein, efforts to link the jamming transition to rigidity percolation remain elusive [71,77].Additionally, contact connectivity percolation [78,79] has critical exponents that are different from those for random percolation [80,81].Our results also contrast with non-equilibrium, boundary-driven sheared granular systems exhibiting long-range force correlation [36] and delineate the difference between an active matter system, driven by local subsystems with characteristic length ℓ and their collapse after performing a finite size scaling analysis for P(p, ℓ) (b, f) and S(p, ℓ) (d, h).First row (a-d) corresponds to increasing dimensionless cell-cell adhesion ω ~= ω cc /ω cw and the second row (e-h) corresponds to the increasing strength of active traction forces α ~.
royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 21: 20240022 injection of energy and boundary-driven non-equilibrium systems.In this vein, our study provides a new context to explore potential links between dense active matter and glasses [82][83][84].Exploring other pathways to transition, such as intercellular friction [85] and the role of mechanochemical waves [86], are exciting avenues for near-future studies.Moreover, progress is made towards direct measurement of traction forces and mechanical stress fields in three dimensions [87] which can be used to probe the proposed universality in cellular monolayers in experiments.Furthermore, understanding the spatial organization of stresses and their correlations in active p-atic liquid crystal models beyond onefold (p=1) rotational symmetry, i.e. active polarity, and coupling with active poroelasticity [88] are the natural next steps for this study.

Figure 1 .
Figure 1.An example showing stress chains in an active monolayer.Snapshots of isotropic stress fluctuations, δσ iso = σ iso − ⟨σ iso ⟩, normalized by maximum compression for (a) solid and (b) liquid states outlining both compressive (blue) and tensile (red) chains.The trajectories of individual cells during total simulation time steps (n sim ) corresponding to (c) a solid-like state and (d) a liquid-like state.The colour bar indicates simulation time.(c) and (d) correspond to simulations shown in (a) and (b), respectively.

Figure 2 .
Figure 2. Characterization of solid-to-liquid transition via one-dimensional pair-correlation function, corresponding to time-averaged configurations, for (a)-(d) cell-cell adhesion drive and (e)-(h) active traction forces.

Figure 3 .
Figure3.Finite-size scaling analyses of isotropic stress criticality.The density of spanning cluster P(p, ℓ) (a, e) and the average cluster size S(p, ℓ) (c, g) for different subsystems with characteristic length ℓ and their collapse after performing a finite size scaling analysis for P(p, ℓ) (b, f) and S(p, ℓ) (d, h).First row (a-d) corresponds to increasing dimensionless cell-cell adhesion ω ~= ω cc /ω cw and the second row (e-h) corresponds to the increasing strength of active traction forces α ~.

Figure 4 .
Figure 4. Global field statistics for time-averaged, normalized, mechanical fields.(a) The mean and the susceptibility of time-averaged out-of-plane stress fields for the cell-cell adhesion drive.(b) The mean and the susceptibility of the time-averaged maximum in-plane shear for the active traction drive.All values are normalized by those corresponding to ω ~= 0.2 for (a) and α ~= 0.2 for (b).

Table 1 .
Critical exponents obtained from finite-size scaling analyses compared with those from two-dimensional random percolation (RP).