Origin of yield stress and mechanical plasticity in model biological tissues (2025)

The confluent jamming transition is not unique

To investigate the mechanical behavior of dense epithelial tissues under substantial deformation, we employed a Voronoi-based Vertex model^9,21. The cell centers {r_i} and their geometric configurations are derived from Voronoi tessellation. The biomechanical interactions are captured through a dimensionless mechanical energy functional²⁷ expressed as: \(\varepsilon={\sum }_{i=1}^{N}[{\kappa }_{A}{({a}_{i}-1)}^{2}+{({p}_{i}-{p}_{0})}^{2}],\) where a_i and p_i are the dimensionless area and perimeter of each cell, κ_A is the rescaled area elasticity, and p₀ is the preferred cell shape index illustrating the cells’ homeostatic shape (see “Methods”). To probe the tissue response, we applied quasi-static simple-shear deformation using Lees-Edwards boundary conditions. Strain was incrementally increased, and the FIRE algorithm was used to minimize energy (see “Methods”).

In the absence of shear, it has been demonstrated that the preferred cell shape index p₀ drives a rigidity transition at \({p}_{0}={p}_{0}^{ * }\approx 3.81\), where the linear response shear modulus vanishes²⁰. Recent studies^9,28 have shown that beyond this transition point, in the liquid phase (\({p}_{0} > {p}_{0}^{ * }\)), the modeled monolayer can undergo strain-stiffening, indicating a gain in rigidity upon strain application. In our quasi-static shearing protocol, we explore beyond the linear response shear startup regimes and into the large deformation or steady-shear limit. In this regime, the tissue exhibits plastic flow primarily through cell-cell rearrangements, or T1 transitions. Here, the tissue’s yield stress is given by

Where σ is the xy-component of the stress tensor. The average is taken over all strain values in the steady shear regime. As illustrated in Fig.1a, while the startup shear modulus \({G}_{0}\equiv {\lim }_{\gamma \to 0}\partial \sigma /\partial \gamma\) (see Methods, Eq. (11) for the calculation of shear modulus) vanishes⁹ at the rigidity transition p₀≈3.81, signaling a solid-to-fluid transition, the yield stress σ_yield does not disappear. Instead, it vanishes at a greater cell shape index, p₀≈4. This suggests a marked difference in tissue responses between the transient shear startup and steady-shear regimes.

Under steady shear and at shape indices higher than the rigidity transition associated with shear startup (i.e., at \({p}_{0} > {p}_{0}^{ * }\)), initially fluid-like systems can intermittently exhibit solid-like behavior before reverting to fluid-like states after yielding (Fig.1(c)). Here the solid-like states are characterized by a non-zero shear modulus,(e.g., states c, d, e in Fig.1c), while states that do not resist shear deformation, indicated by having zero shear modulus, are fluid-like (e.g., states a, b, f in Fig.1(c)). To transition from state b to state c, the system gains rigidity under shear, moving from a fluid-like to a solid-like state. Stress then builds up, as indicated by edges experiencing high tension (Fig.1e). Once sufficient stress is stored, the system yields, relaxing this stress via collective rearrangements, commonly referred to as avalanches, and transitions from state c to state d.

In the steady-shear regime, we can compute the instantaneous shear modulus (see “Methods”) to quantify solid vs. fluid-like states. The solid-fluid coexistence is reflected in the bimodal distribution of the shear modulus p.d.f(G) shown in Fig.1b, where the fluid phase is associated with a peak near the numerical noise floor of the shear modulus (~10⁻¹²), while the solid phase corresponds to a finite shear modulus. The shifting behavior of the distributions can be quantified by the fraction of solid states ρ_solid shown in Fig.1b. States below the rigidity transition \({p}_{0}={p}_{0}^{ * }\approx 3.81\) are therefore always in the solid phase, which we term a pure solid. In the range p₀ ∈ [3.81, 4], ρ_solid drops below 1, indicating a solid-fluid coexistence, which we refer to as marginal. For p₀ > 4, the simulated tissue always remains in the fluid phase as it cannot build up stresses in response to shear strain. This is consistent with the yield stress vanishing at p₀ ≃ 4.0. The fact that the material response depends on the application of shear is reminiscent of shear jamming in granular materials, where a state below the isotropic (un-sheared) jamming threshold can be jammed with the application of shear^29,30,31,32. The coexistence of solid and fluid phases also has analogs in dense suspensions near shear jamming³³ and discontinuous shear thickening^10,34.

Predicting the tissue yield stress using a refined Soft Glassy Rheology model

Given the continuous behavior of yield stress across the pure solid - marginal state transition, we aimed to develop a unified model to deepen our qualitative understanding of the steady-shear regime properties using the Soft Glassy Rheology (SGR) framework^35,36. In the SGR model, mesoscopic elements, characterized by an elastic constant k and local strain l, are confined within energy traps E, where they accumulate elastic energy as macroscopic strain increases, approaching a yield point either directly or through an activated “hop" driven by mechanical fluctuations from the yielding of neighboring elements. The material’s dynamics under shear are governed by the probability distribution P(E,l,t) evolving in time t, which follows the Fokker-Planck equation^35,36:

The first term in Eq. (2) represents the motion of the elements driven by the applied shear rate, \(\dot{\gamma }\). The second term describes thermally activated hopping from a trap with an effective depth of E−kl²/2, which corresponds to the remaining distance to yielding. x and g₀ represent the mechanical noise in the system akin to temperature and the hopping rate, respectively. The final term represents the transition to a new trap with energy E drawn from a quenched random distribution ρ(E). The Dirac-delta function δ(l) reflects the assumption that the local strain l is reset to zero after yielding. The total yielding rate at time t,g(t), is explicitly defined in the SI text (seeSI).

In the SGR model, the choice of the functional form of ρ(E) critically influences material behavior³⁶. Direct measurement of energy barrier distributions is challenging, leading prior studies to adopt generic or ad hoc assumptions about ρ(E), such as an exponential distribution^37,38,39,40. In this work, we introduce a novel approach based on the distinct mesoscopic tissue phases observed: (1) fluid elements with zero yielding energy (E=0) and (2) solid elements with finite yielding energy (E>0). Consequently, we propose a refined ρ(E):

Here, f₀ denotes the probability of an element transitioning to a state with E=0, while 1−f₀ corresponds to transitions into states with energy sampled from a k-gamma distribution, parameterized by the mean x₀ and shape factor κ, and Γ(κ) is the Gamma function. This is based on the previous observation that the energy barriers to the T1 transition follow a k-gamma distribution^20,41 with κ≈2. Together, Eqs. (2) and (3) describe three potential transitions in the energy landscape, as depicted in Fig.1: (1) fluid–to–fluid a → b, (2) solid–to–solid c→ d, and (3) solid–to–fluid e → f.

We next examine the steady state behavior of Eq. (2) in the quasi-static limit (\(\dot{\gamma }\to 0\)), with details shown in theSI text. This behavior is governed by three parameters: the dimensionless ratio of mechanical noise to mean yielding energy χ=x/x₀, the probability f₀ of transitioning to a fluid state, and the elastic constant k of solid elements.

An important aspect of the SGR model is that the fluctuations driving rearrangements originate from the mechanical noise generated by other surrounding rearrangement events in the system. These fluctuations are analogous to the energy released during yielding events observed in our simulations. To correlate this mechanical noise with our empirical data, we introduce the following relationship:

Here, 〈ΔE〉 represents the average energy dissipated during yielding events, while 〈E〉 denotes the average energy of cells in the solid state. Next, by analyzing the steady-state solution of Eq. (2), we determine the probability that an element is in the solid phase as a function of f₀ (see SI, Eq.S.9). This precisely corresponds to ρ_solid in our simulations (Fig.1b). Finally, we treat the elastic constant k as independent of the shape index p₀. Given that both \(\langle \Delta E \rangle / \langle E \rangle\) and ρ_solid depend on p₀, the yield stress predicted by the SGR model (seeSI) effectively varies only with p₀. This approach contrasts with previous studies that employed the SGR model^37,38,39,42, in which χ=x/x₀ was often treated as a fitting parameter. In our research, we derive χ directly from simulation data, enhancing the predictive accuracy of our theoretical results and distinguishing our use of the SGR model as predictive rather than merely descriptive. In Fig.1(a), we plot the SGR-predicted yield stress as a function of p₀, demonstrating that the SGR model accurately predicts the vanishing point of the yield stress and its dependence on the cell shape index σ_yield(p₀).

The dual-state SGR model identifies two primary mechanisms responsible for the yield stress transition: (1) As p₀ increases, states with zero yielding energy barriers become more prevalent, leading to frequent yielding under deformation. This behavior is depicted by transitions such as a → b and e → f in Fig.1d; (2) Concurrently, mechanical noise from stress redistributions approaches the scale of the yielding energy barriers, enhancing the likelihood of solid-solid transitions through activated processes induced by neighboring rearrangements, as illustrated by c → d and e → f.

A new method to infer tissue stresses based on instantaneous snapshots

The coexistence of phases observed in the Vertex-based model reinforces the idea that tissues behave as yield stress materials while also exhibit fluid-like behavior. When combined with Soft Glassy Rheology (SGR) theory, this coexistence allows for predictions for tissue yield stress based on the homeostatic cell shape index p₀. Given p₀, the Vertex based model could also provide an estimate of the instantaneous tissue stress using Eqn. (10). However, determining the homeostatic target shape index p₀ remains an experimental challenge. In contrast, segmented cell configurations in tissues are experimentally accessible and have been utilized in several non-invasive stress inference methods such as Bayesian Inference method^43,44 and the image-based Variational method⁴⁵. Although the Bayesian Stress Inference method proposed by Ishihara and Sugimura⁴⁴ has been applied to various systems, including Drosophila notum, retinal ommatidia, germband, and the quail early embryo^46,47, its results strongly depend on a prior distribution of edge tension and cell pressure, which is not necessarily normally distributed as originally proposed. Additionally, the method faces the challenge of having more unknown variables than constraints. In contrast, the Variational Method proposed by Noll et al.⁴⁵ addresses the issue of underconstrained variables but relies on a computationally expensive fitting approach. Here, we propose a fast, non-invasive, image-based tissue stress inference method that is both convenient and accurate. It has been shown that edge length distribution is closely related to the distance to yield stress in monolayers¹⁹. However, the connection between edge length distribution and tissue stress has not been quantified. In this work, we demonstrate that the cumulative distribution of edge length elegantly serves as a robust estimator of tissue stress.

In systems under deformation, the configuration exhibits anisotropy, with preferred orientations dictated by the deformation direction. To analyze this anisotropy and its evolution with strain, we examined the distribution of edge vectors in the system, represented as polar distributions. In these distributions, color encodes the frequency of edges with lengths and orientations defined by radial distance and azimuthal angle, respectively. The polar distribution for an isotropic system therefore looks like a circle with dots scattered randomly inside. This was not observed in our system due to the anisotropy. To systematically study across different systems and time frames, we used the normalized edge length \(l=L/\bar{L}\), obtained by normalizing the edge lengths by the average edge length in the current snapshot. Fig.2a displays the polar distribution of the normalized edge vectors \(\vec{l}\) in a low-stress, stable system far from yielding. Conversely, Fig.2b depicts the edge vector distribution in a high-stress, unstable system close to yielding. In both cases, anisotropy emerges, with longer edges aligning with the shear direction (about 45 degrees with respect to horizontal) and shorter edges perpendicular to it. However, the trend is more pronounced in unstable systems, illustrated by the red color band representing high population of short edges with orientation about 135 degrees with respect to horizontal (Fig.2b), implying a connection between instability and edge vector distribution. Furthermore, these polar plots reveal a strong correlation between edge length and orientation in deformed configurations.

The correlation between edge length and orientation suggests that the edge length distribution can effectively represent the overall edge vector distribution. To investigate how edge configuration influences instability, we compared the edge length distributions of stable and unstable states. As illustrated in Fig.2c, both distributions exhibit a bimodal shape due to anisotropy. However, in the distribution of the unstable system, the separation between the two peaks is larger, reflecting stronger anisotropic effects and an abundance of short edges, which are more prevalent in the unstable configuration. This motivated us to study the evolution of the edge length distribution as strain increases.

By investigating the evolution of the edge length cumulative distribution C(l) in our quasi-static simple shear simulations, we observed a correlation between C(l) and σ, with the correlation level depending on l. As shown in Fig.2d, while the correlation R_σ,C between tissue stress σ and C(l) varies with l, there exists a range of l where C(l) is highly correlated with σ. We denote the l value that maximizes the correlation R_σ,C as l*, and the corresponding cumulative distribution C(l*) as C* (also seeSI text and SI Fig.S1 for details). Interestingly, both l* and R_σ,C remain robust with changes in p₀, with the critical correlation \({R}_{\sigma,{C}^{ * }}\) exceeding 0.9, as shown in Fig.2e. Given the strong correlation between C* and σ, as well as the fact that cell edge lengths can be directly extracted from imaging, C* could serve as a non-invasive metric for inferring tissue-level stress.

While l*≈0.61 is consistent across various p₀ values in our quasi-static simulations, its practical value may depend on system properties and the shearing direction. Thus, developing a metric to experimentally determine l* is crucial for the stress inference method. Recognizing the strong relationship between edge length and orientation in anisotropic systems, we investigated tissue stress and edge orientation correlations. To quantify the edge orientation distribution, we used the proportion of edges a given orientation ω:

Here, N(ω−Δω,ω+Δω) represent the number of edges in a cone with the open angle of 2Δω at the ω direction and N_edges is the total number of edges in the system. In this section, all angular values are in degrees. Tracking the evolution of F(ω) at different ω as strain increased, we observed a critical orientation ω*≈120°, where F(ω) was highly correlated with tissue stress, with a correlation coefficient \({R}_{\sigma,F({\omega }^{ * })} > 0.9\) (Fig.2e inset). Although R_σ,F depends on the cone wideness Δω, which we chose as Δω=18°=π/10 to optimize R_σ,F, the critical orientation ω* is independent of Δω. Additionally, ω* was robust across variations in p₀ (Fig.2e inset). Given the high correlations R_σ,C and R_σ,F, the critical edge length l* could potentially be determined from the correlation between edge length and edge orientation, R_C,F. As suggested by Fig.2e, we only focus on the range l<1 in which C(l) is known to positively correlates with tissue stress σ. By computing the correlation R_C(l),F(ω) at different (l,ω), we identified a region where the correlation is exceptionally high (greater than 0.93), as indicated by the dashed circle in Fig.2f. The l values in this region align well with the critical edge length l*, suggesting that l* can be estimated as the normalized edge length value l maximizing R_C(l),F(ω). With l* in hand, ones could used C*=C(l*) to extract an estimation of the tissue stress as it evolves during the deformation process. Our stress inference method could be conveniently validated by experiments such as the monolayer stretching protocol³. Although C* could be used directly as the inference for tissue stress, various functional relationships can be established between C* and stress (seeSI text).

Dynamics of tissue plasticity

So far, we know that for a system in the coexistence phase to transition from solid to fluid, a collective rearrangement event is required to significantly remodel the configuration. However, the mechanism that governs the occurrences of these events and the evolution of the system during the events is still not fully understood. To explore the yielding behavior of biological tissues, it is essential to describe the dynamics of plastic events during avalanches. In this study, we examine the plasticity dynamics by investigating the spatiotemporal evolution of the plastic rearrangements created by other rearrangements during an avalanche^48,49. Fig.3a displays a space-time plot of the occurrences of T1 transitions during an avalanche, with the cells labeled in red indicating participation in the T1 transitions.

As depicted in Fig.3a, an avalanche involving numerous plastic rearrangements can originate from a single T1 transition, which we refer to as the initial trigger. Starting from the initial trigger, the stress redistribution from each event can also stimulate surrounding cells to become unstable and undergo T1 transitions. These soft cells, whose mechanical characteristics are distinct from the surrounding cells, experience larger deformation than the neighboring cells and have been observed in multicellular tumor spheroids using 3D light microscopy⁵⁰. From an amorphous solid point of view, these cells are referred to as soft spots^23,24,25,26 or Shear Transformation Zones^51,52. This cascade of cellular rearrangements can therefore lead to an avalanche, which continues until the population of soft spots is sufficiently depleted. In Fig.3b, we show the number of accumulated T1 transitions and the tissue shear stress during a typical avalanche. The stress relaxation due to an avalanche is therefore the origin of the discontinuous yielding of stress during quasistatic shear. In Fig.3a, the location of rearrangements over time suggests a preferred direction for avalanche propagation. In order to quantify this and establish a causal relationship in time, we define a two-point, two-time correlation function:

where P(r,t) is a binary field, representing the occurrence of a T1 transition (1 if a T1 transition occurs at r and time t, and 0 otherwise). 〈...〉 denotes spatial and ensemble averaging. With this definition, ϕ represents the conditional probability of observing a T1 transition at location r₀+r and time t₀+Δt, given that a transition has already occurred at (r₀,t₀).

In Fig.4a, we illustrate the evolution of the field ϕ as Δt increases. The field ϕ resembles a wave that propagates away from the causal rearrangement and is primarily in the x-direction, aligned with the direction of the external shear force. The evolution of the field ϕ reflects the interplay between the stress redistribution from a plastic rearrangement and the population of soft spots which could rearrange under the effect of the strain field.

We first focus on the angular dependence of ϕ, which shows an anisotropic four-fold pattern. This anisotropy is consistent with the stress redistribution field due to a rearrangement in an elastic medium as predicted by elastoplastic models^19,53. However, it differs from the isotropic probability field observed in ductile, soft disk systems⁵⁴. This contradiction likely arises from the difference in shape anisotropy between soft disks and cells in our system. While soft disk systems exhibit minimal particle shape anisotropy, cells in our system can sustain large deformations and have highly anisotropic shapes. Consequently, deviatoric strain triggers rearrangements in soft disk systems, whereas simple shear strain is responsible for triggering rearrangements in our system.

To better understand how the rearrangement probability field propagates, we looked at \({\phi }_{x}=\frac{\phi (x,y=0)}{{\sum }_{x}\phi (x,y=0)}\) and \({\phi }_{y}=\frac{\phi (x=0,y)}{{\sum }_{y}\phi (x=0,y)}\) separately. ϕ_x, as observed in Fig.4b, is a bimodal distribution that evolves in time such that the distance between the two peaks d_peaks increases as time progresses. The diffusing bimodal distribution suggests that the propagation results from a combination of convection and diffusion, and the drift velocity could be captured by the rate at which the peaks’ separation increases. Conversely, ϕ_y is a bell-like shape distribution that gets broaden over time (Fig.4c), indicating that the propagation in the y-direction is similar to a purely diffusion process with the diffusivity can be captured by the evolution of the FWHM. Fig.4d shows that d_peaks increases faster than the FWHM-ϕ_y. Since the rearrangement probability field is the result of the shear stress redistribution and the population of soft spots, the propagating mechanism of the field also should agree with the propagation of the stress redistribution.

The tendency of shear stress redistribution to propagate in the direction of shear has also been observed in particle-based systems governed by inverse-power-law pairwise potentials⁴⁹. This behavior bears a resemblance to the propagation of elastic waves. In the theory of elasticity, longitudinal waves, characterized by displacement in the direction of propagation, outpace transverse waves⁵⁵. Moreover, longitudinal elastic waves involve changes in local density⁵⁵, akin to the x-propagating excitation wave’s modulation of local density via T1 transitions. Conversely, transverse elastic waves do not induce density changes, resembling the infrequent involvement of T1 transitions in y-propagating excitation waves. Notably, the mechanism driving stress redistribution to preferentially propagate in the shear direction appears universal, independent of p₀. However, since p₀ governs the elastic response in our system, with higher values corresponding to a less elastic state, there is a negative correlation between the stress redistribution wave’s speed and p₀.

Statistics of tissue avalanches

In addition to the universal propagation mechanism, we wondered if the statistics of tissue yielding events also exhibit universality. In Fig.5a, both the average yielding size \(\bar{S}\), denoting the total number of T1 transitions after a yielding event, and the average stress drop, representing the amount of stress relaxed by the event, exhibit the same dependence on p₀. In the solid regime, while the stress decreases with increasing p₀, the average yielding size shows minimal variation. This trend of \(\bar{S}\) versus p₀ is akin to that observed in Fig.1b for the proportion of the solid state.

However, in the marginal phase, there are different types of yielding events as discussed previously (Fig.1c). In the yielding events that occur while the system is fluid-like, illustrated by the a → b transition in Fig.1c (Type I), the tissue lacks rigidity and therefore is unable to transmit stress to initiate a cascade of rearrangements. Conversely, yielding events following a solid state, illustrated by c → d and e → f transitions in Fig.1c (Type II), tend to be cascading as the rigid tissue is capable of propagating the stress redistribution. It is this latter type that we refer to as tissue avalanches from now on. Since the avalanches growing mechanism is universal, we expect their statistics to be independent of p₀ as well.

Excluding yielding events of type I from the analysis and specifically analyzing only the avalanches, we indeed find that the average avalanche size \(\overline{{S}_{s}}\) does not vary significantly with p₀ (Fig.5a inset), suggesting universal avalanches size statistics. To rigorously assess this universality, we examine the distribution of avalanche sizes across various p₀ values (see Fig.5b). We observe a consistent power-law distribution, reminiscent of the Gutenberg-Richter law observed in earthquakes^56,57, with an exponent τ=−1.36, which agrees with the analogous exponent observed in overdamped elastoplastic models under shear^58,59 and in vertex model on spherical surface⁶⁰. This shared characteristic suggests a parallel between biological and seismic avalanches and supports the argument that the vertex model and epithelial tissues belong to the universality class of plastic amorphous systems.

Furthermore, we find that the same power law applies to the distribution of average stress drops during avalanches when scaled by the average (Fig.5b inset). This collapse after rescaling implies that the stress relaxation mechanism via avalanches is independent of the shape index p₀, and p₀ only affects the average stress relaxed by governing tissue overall stiffness. Moreover, the similarity in the stress drop distribution and avalanches size distribution indicates that each plastic rearrangement, on average, releases a similar amount of stress that depends only on p₀. The convergence of these distributions suggests that the growth of avalanches remains unaffected by changes in p₀, providing additional evidence for the universal propagation mechanism discussed earlier.

Predicting tissue avalanches based on static structural information

While the first cells to undergo a T1 transition trigger the avalanche, in order for the avalanche to grow, it is necessary to have soft spots in the system that are susceptible to undergo T1s. In the framework of the elastoplastic model, it has been established that the distance to yield x, which represents the additional stress required to trigger a yielding event, follows a power-law distribution, p(x) ∝ x^{θ 61,62,63}. The exponent θ has been suggested as a measure of the system’s instability, with a higher value indicating a more stable state.

In the vertex-based model family, it has been proposed that the distance to yield x exhibits a linear relationship with the length of cell edges L, and that the distribution of edge lengths should follow the same power-law behavior as ρ(x)¹⁹. While this argument establishes a connection between system configuration and instability, the efficacy of using the distribution of short edges to describe instability remains uncertain.

To address this ambiguity, we investigate this concept within our Voronoi model and observed an intriguing correlation between the exponent θ and system instability (seeSI text and SI Fig.S2). However, θ is not a reliable metric because it is derived from a power-law fit that heavily depends on the range of fitting¹⁹. In practice, the cumulative distribution function (c.d.f) of L exhibits power-law behavior only within a specific range, which varies from sample to sample. Therefore, we propose using C* as the parameter of instability to avoid the uncertainties and biases associated with fitting.

Using C* as the parameter of instability, with higher C* corresponding to a more unstable state, we observed a relationship between instability and avalanche size. To quantify the relation between instability and avalanche size, we focus on the failure states, i.e states that are about to yield, indicated by having a stress drop in the next strain step and a corresponding avalanche. We computed C* for those states and then grouped the avalanche sizes based on C* (Fig.6a). At low C*, if the tissue yields, the size of the yielding event is likely to be small, indicated by a rapid increase to 1 in the cumulative distribution function (c.d.f.) of S. As C* increases, the likelihood of larger yield events grows. The probability of observing an avalanche of size 25 or greater is summarized in Fig.6b.

While C* can provide predictions about avalanche size, it cannot determine when an avalanche will occur. Therefore, we require an additional tool to forecast yielding events. In amorphous solids, the locations of plastic rearrangements during an avalanche largely depend on the material’s structural configuration, with areas more likely to experience plastic events called soft spots. Various frameworks have been proposed to link structure and plasticity, such as the Shear Transformation Zone theory^{25,51,52,64,65} and lattice-based models. The most promising theoretical approach for predicting the locations of soft spots involves identifying these areas based on soft vibrational modes^24,66,67,68. As a system approaches a plastic rearrangement, at least one normal mode is supposed to approach zero frequency^24,69,70. However, the vibrational mode analysis is not applicable to the vertex-based model family due to the cuspiness of the energy landscape^20,41. In such systems, the energy cusp at a plastic event prevents the corresponding low-frequency mode from vanishing as it would in systems with smooth, analytic energy landscape. As shown in the SI Fig.S3, in our system, no low-frequency mode approaches zero frequency except at the onset of the plastic event. Hence, an alternative approach is necessary.

In our model, the deformation of edges immediately following shear strain is deterministic. Through simple geometry, we deduce the existence of a range of orientations wherein edges are prone to shortening upon shearing, rendering them more susceptible. In the vertex model, under the condition \(\dot{\gamma }\ll 1\), the change in edge length δL due to shear is approximated as \(\delta L\approx \dot{\gamma }L\sin (2\phi )\), where L and ϕ represent the edge length and orientation, respectively. Consequently, in the vertex model, the most susceptible orientation is \(\frac{3\pi }{4}\). If an edge is sufficiently weak (or short) and happens to align with this susceptible orientation, it may yield under the influence of shearing, potentially triggering further rearrangements in its vicinity. We refer to these susceptible edges as triggers. Since triggers are local elements, their presence is not captured by the cumulative distribution function of edge lengths, c.d.f(L), thereby explaining why C* alone cannot predict imminent system failure. In summary, the presence of a trigger is the necessary condition, while a high C* in the current state is the sufficient condition for a large avalanche in a tissue monolayer.