Random traction yielding transition in epithelial tissues

We investigate how randomly oriented cell traction forces lead to fluidisation in a vertex model of epithelial tissues. We find that the fluidisation occurs at a critical value of the traction force magnitude $F_c$. We show that this transition exhibits critical behaviour, similar to the yielding transition of sheared amorphous solids. However, we find that it belongs to a different universality class, even though it satisfies the same scaling relations between critical exponents established in the yielding transition of sheared amorphous solids. Our work provides a fluidisation mechanism through active force generation that could be relevant in biological tissues.

During tissue development, many cells collectively selforganize in dynamic patterns and morphologies.Therefore, a central problem in biophysics of development is understanding the interplay of tissue mechanics and active force generation [1][2][3][4][5][6].Cells in a tissue can generate traction forces through mechanical linkages with a substrate [7][8][9][10] and impairment of this coupling can interrupt the movement of cells as observed, for example, in cancerous spheroid assays of carcinoma and human breast organoids [11,12].The response of biological tissues to mechanical forces is often described as that of viscoelastic active fluids [13][14][15][16].However, recent experimental and theoretical studies have revealed complex mechanical phenomena, including jamming, glass transitions [17][18][19][20][21][22], and yield stress rheology [23,24].These observations suggest that developing biological tissues can behave as active amorphous solids.
Recently there has been an increasing interest in the rheology of active amorphous solids [25,26].In particular, comparing uniform shear to random forcing of particles revealed a very similar non-linear response [27].A hallmark of sheared amorphous solids is a transition from a solid to a plastically flowing state at the yield stress Σ c .The plastic strain rate γ at stress Σ above the yielding transition typically follows the Herschel-Bulkley law γ ∼ (Σ − Σ c ) β , where β ≥ 1 is the flow exponent.Yielding has recently been reported under random forcing in systems of jammed self-propelled particles [28].This raises the question of what is the nature of the yielding transition under random forces and how it is related to the yielding transition under uniform shear.Such random yielding is relevant in the context of biological tissues, allowing them to fluidise through generation of cell traction forces.
Here, we investigate the critical properties of the yielding under random traction forces using a vertex model of epithelial tissues [24,29].Motivated by recent experiments on mouse pancreas spheres which suggest a pres-ence of tissue fluidisation by cell traction forces [30], we consider a vertex model with spherical geometry.This geometry is ubiquitous in multicellular systems such as the early developmental stages of many tissues, including early vertebrate embryos [31,32], and early stages of organoids [30,33,34].We find that randomly oriented traction forces fluidise the cellular network beyond a critical magnitude F c .We call this transition the random yielding transition (RYT).We quantify the critical exponents characterizing overall cell flow, patterns of cell rearrangements, and even the geometry of the cellular network.We compare our results to the properties of the uniform shear yielding transition (YT).Interestingly, some critical exponents differ between the RYT and YT, implying that the transitions belong to different universality classes.Furthermore, we find that RYT critical exponents satisfy the scaling relations between exponents established for the YT [35].These relations imply that the statistical properties of tissue dynamics and cellular geometry are not independent.
Random traction vertex model.We extend the standard vertex model of epithelial tissues [29] to a spherical geometry (Fig. 1a).We represent cells as polygons outlined by straight bonds, and constrain the polygon vertices to move on a sphere of radius R. Geometry of the cellular network evolves following the dynamical equation: where u m is the velocity of vertex m, ζ is the friction coefficient, f a m is the traction force, W is the vertex model energy function, and f n m is the normal force constraining the motion of vertices on the sphere surface (see Supplemental Material (SM) for details).The vertex model energy function that accounts for cell area elasticity and cell bond tension reads: Here, A α is the cell area, L α is the cell perimeter, A 0 is the preferred cell area, K is the area stiffness, and Λ is the perimeter tension magnitude [29,36].We choose units of length, force and velocity to be A We consider a planar cell polarity p α that directs the traction force exerted by cell α on the surrounding matrix (Fig. 1b).We initialize the direction of the cell polarity vectors p α from a uniform distribution, and evolve it following the dynamical equation: where D/Dt denotes a co-rotational time derivative (see SM), and we impose |p α | = 1 at each time.We define the active traction force f a m on a vertex m by uniformly redistributing the cell traction force f α p α of each of the abutting cells with M α number of vertices: Random yielding transition.Random traction forces induce stresses in the vertex model network.The stress magnitude is controlled by the magnitudes of cell traction forces f α .For small magnitudes of traction forces, we find that the elastic forces generated by the vertex model network balance the traction-induced forces, and the network remains solid.However, upon further increasing f α , the network begins to flow through cell rearrangements.To quantitatively explore this transition, we introduce the tissue traction force magnitude F ≡ α f α /N , which in RYT plays the role analogous to the shear stress in the YT.
Application of uniform f α is susceptible to finitesize effects that prevent us from probing the transition.Namely, a finite-size system can by chance reach an unusually stable configuration so that the system does not flow even at high F values.To avoid this issue, we implement a model of traction forces where the attachment of a cell to the substrate moves with speed v along the vector p α and the traction force is transmitted to a spring of stiffness κ that connects the attachment and the cell (Fig. 1b).Therefore, the dynamics of the traction force magnitude for a cell α follows: Here, the term −κp α • u α represents the relaxation of the force in the spring due to motion of the cell with velocity u α .Limits of infinitely soft κ → 0 and infinitely stiff κ → ∞ springs correspond to imposed traction forces and imposed cell center velocities, respectively.In the following, we use κ = 0.01 and vary the imposed spring extension velocity v.
An example of F (t) dynamics as a function of spring displacement vt is shown in Fig. 1c (see also Movie 1).Initially, the cellular network responds elastically, and the traction forces grow linearly with spring displacement vt.As F increases further, the cellular network begins to yield through cell rearrangements, visible as sharp drops of F in Fig. 1c.Finally, in the steady state, the system dynamics consist of periods of elastic loading punctuated by avalanches of cell rearrangements that are visible as sudden drops of F .Ensemble-averaged F (t) for different values of v is shown in Fig. 1d.
The observed behaviour of F is reminiscent of the stress vs strain curve in sheared amorphous solids, such as metallic glasses [37], where sudden drops of stress correspond to avalanches of particle rearrangements [38].In amorphous solids near the YT, the avalanche size, defined as the number of particle rearrangements S in an avalanche, is distributed according to a scaling law P (S) = S −τ f (S/S c ), where S c is the cutoff beyond which P (S) rapidly vanishes.The cutoff is set by the correlation length ξ: S c ∼ ξ d f , where d f is the avalanche fractal dimension [39].However, approaching the YT, ξ diverges and becomes larger than the system size.There- fore, in a finite system of N cells, the cutoff S c is set by the system size S c ∼ N d f /d .Furthermore, the duration of an avalanche T is expected to scale with the size as T ∼ S z/d f , where z is the dynamical exponent [39].
To measure the avalanche size distribution, we measure drops in F in the steady state at the lowest value v = 2 • 10 −4 we used.Then, we estimate the avalanche size corresponding to a force drop ∆F as S N ∆F/κ.We find that the avalanche sizes are indeed power-law distributed, as shown in Fig. 2a, with a system-size dependent cutoff.Moreover, we find τ = 1.35 ± 0.11 [40], which is consistent with the values measured in the YT of 2d elastoplastic models (τ = 1.25 ± 0.05 in Ref. [41] and τ = 1.36 ± 0.03 in Ref. [39]), of a lattice model (τ = 1.342 ± 0.004 [42]) and of a finite element model (τ = 1.25 ± 0.05 [43]).We next estimated the avalanche fractal dimension d f = 0.75 ± 0.15 by finite-size scaling analysis of the avalanche distribution cutoff using S c ∼ S 3 / S 2 , see SM.Finally, we find that the avalanche duration follows a power-law relationship with the avalanche size, see Fig. 2, from which we estimate z/d f = 0.68 ± 0.04.
To test whether the spherical geometry influences the exponent values, we also measured τ , d f and z/d f in a bi-periodic 2d vertex model with identical random traction forces, and we found values of the exponents are consistent with the ones of the spherical model, see SM.
Scaling relations connect cellular dynamics and geometry.Exponents of YT are related through several scaling relations [39].Here we examine two of these relations in the context of the RYT and show that in the vertex model with random traction forces they also provide a relationship between statistics of avalanches of cell rearrangements and cell bond length distribution.
The first scaling relation follows from the fact that in the steady state ∆F = 0 [39], which we now briefly reproduce.Increases of F between avalanches are balanced by decreases during avalanches: The scaling of the average decrease of F with system size can be estimated from the avalanche size distribution as After an avalanche, F will increase until the next T1 transition.Therefore, the increases in F are determined by the network regions closest to a T1 transition.In amorphous solids the density of plastic excitations, defined as local increase in shear stress ∆σ required to trigger a plastic event, exhibits a pseudo-gap P (∆σ) ∼ ∆σ θ , with θ > 0 [35,38].Thus, the average smallest ∆σ in a system of size N scales as ∆σ min ∼ N −1/(1+θ) (see Ref. [39]).Since |∆F | + ∼ ∆σ min it follows that: Using the measured values of τ and d f , this scaling relation predicts θ = 0.32 ± 0.11.This prediction can be tested independently by considering the statistics of the bond length distribution as follows.In a vertex model network, each T1 transition corresponds to a vanishing bond; hence, short bonds anticipate the upcoming T1 transitions.Due to cusps in the vertex model energy landscape at the onset of a T1, it was shown for the planar vertex model [24] that the corresponding ∆σ is proportional to the bond length of disappearing bonds.We show that this relation also holds in the spherical vertex model tissue, by measuring the additional tension ∆f required to shrink a bond of length to 0, see Fig. 3a.In general, local change in shear stress ∆σ will generate a proportional change in the bond tension ∆f .Therefore, observed scaling of imposed ∆f with bond length characterises the scaling of ∆σ.As a consequence, short bonds in the network for F ≤ F c are distributed according to P ( ) ∼ θ .Figure 3b shows the cumulative bond length distribution C( ) = 0 P ( )d obtained in the steady-state simulation at v = 2 • 10 −4 , where we measure bond lengths of networks at time points just after an avalanche.We find that the predicted value of the exponent θ is consistent with the bond length distribution (see also SM).
The second scaling relation reflects that the flow in the  vicinity of the critical point F c is composed of avalanches of spatial extension corresponding to the correlation length ξ ∼ (F − F c ) −ν .Since the average avalanche size scales with S ∼ S c ∼ ξ d f and its duration scales as T ∼ ξ z the contribution of the average avalanche to the overall flow v will scale as v ∼ S/(T ξ d ) ∼ (F −F c ) ν(d−d f +z) [39].This determines the exponent β = ν(z+d−d f ) defined by v ∼ (F −F c ) β .Here we do not directly measure ν and instead we use an additional scaling relation ν = 1/(d−d f ) [39].Therefore, we arrive at the relation: which allows us to estimate β = 1.41±0.098.To test this prediction, we analyze the steady-state flow properties for various magnitudes of loading rate v as shown in Fig. 4 for two sizes N = 100, 200.We find a good agreement between numerical results and the value of β predicted by the scaling relation (7).Discussion.We have described the critical properties of the RYT due to randomly oriented traction forces acting on a spherical epithelium.Our results show that this transition is closely related to the YT of sheared amorphous solids.Furthermore, we find that scaling relations constraining critical exponents of the YT also hold in the RYT, differing from the recent suggestion that one of the relations is violated [28].Furthermore, we independently measure the pseudo-gap exponent θ describing the density of plastic excitations.In our model, this exponent follows directly from the scaling of the distribution of cell bond lengths [24] while it is typically difficult to access in particle models.
We find that the value of fractal dimension d f and pseudo-gap exponent θ are clearly different from the YT values; see Table I.In particular, d f ≈ 1.1 in 2d YT is associated with the one-dimensional shape of avalanches of plastic events, arising from the anisotropy of the Eshelby stress propagator of individual plastic events.In the RYT, the orientation of plastic events is not aligned, which breaks the preference of avalanches to occur along lines, and the value of d f = 0.75 ± 0.15 smaller than 1 shows that the spatial structure of avalanches is sparse.It is interesting to compare RYT to the yielding transition in a mean-field elastoplastic model where the Eshelby stress propagator is randomly redistributed in space, thereby removing all spatial correlations [44], where the pseudo-gap exponent θ = 0.39 ± 0.02 has been reported numerically and supported by analytical calculations.This value is significantly lower than the 2d YT value θ ≈ 0.57 [39].However, since this is consistent with the value θ = 0.32 ± 0.11 that we find in RYT, it would be interesting to test whether RYT is in the mean-field yielding transition universality class by carefully measuring the relevant critical exponents.
To test the influence of spherical geometry on the RYT we have measured the critical exponents τ , z, and d f in flat 2d bi-periodic vertex model simulations (see SM).We found no significant difference in their values, which suggests that spherical geometry does not alter the critical behaviour of the vertex model near the RYT.
The dynamical exponent z describes the dynamics of avalanche propagation T ∼ l z , where l is the linear extension of the avalanche.The value z = 0.51±0.11we find is consistent with reported values in YT in 2d elastoplastic model z = 0.57 ± 0.03 [39] and z 0.5 [45].However, in the thermodynamic limit z < 1 cannot hold due to the finite propagation speed of elastic interactions, which requires z ≥ 2 in overdamped systems.Indeed z ≥ 2 was reported in a large system of disks with overdamped dynamics [28].In finite systems z can be smaller if the elastic interactions propagate through the system faster than the avalanches of cell rearrangements.Note that in elastoplastic models elastic interactions propagate instantaneously.The value of z that we find suggests that, for the biologically relevant system sizes we consider, the elastic interactions in our model propagate much faster than avalanches, effectively behaving as instantaneous.
The fluidisation through the generation of traction forces could allow the biological tissues to transition between a stable solid phase and a malleable fluid phase without the need to alter tissue density [23] or cell mechanical properties [46].We speculate a similar transition could occur in tissues where cells generate randomly oriented active stresses instead of traction forces.
where M α = m∈α M m,α /M α is the average moment of inertia tensor.In the case of a solid-body rotation, δu m,α (t) = 0, such that the angular velocity is simply obtained from Eq. (S.12) as: Duration of an avalanche can be uniquely defined in the quasistatic limit v → 0 as the time period during which F is decreasing due to successive cell rearrangements.After the avalanche no cell rearrangements occur until F in the system increases sufficiently by spring displacement to trigger the next cell rearrangements.Duration of this loading period diverges in the limit v → 0 and avalanches can be precisely identified.
At a finite v, the observed time intervals during which F is decreasing depend on the time resolution δt at which the data is recorded.Each spring moves a distance δs = vδt in a time interval δt (see Eq. 5 in the main text) during which F increases by kδs due to spring movement.Occasionally during an avalanche the F decrease rate can fall below kv and the F slightly increases although the avalanche is still ongoing.If the time resolution δt is very small these slight increases will often be recorded effectively splitting original avalanches into smaller ones.On the other hand choosing too large δt leads to merging of F decrease intervals belonging to different avalanches.
We first determine the extreme limits of low and high δt for which described artifacts leading to splitting and merging of F decrease intervals are clearly visible (see Fig.

VIII. CRITICAL TRACTION FORCE MAGNITUDE
Figure 4 in the main text suggests that the magnitude of critical traction force F c (N ) is system-size dependent and is expected to be of the form F c (N ) − F c ∼ N −1/(dν) , since a system smaller than a correlation length ξ ∼ (F − F c ) −ν cannot be distinguished from the system at F c .Fitting the Herschel-Bulkley law to numerical simulation data (Fig. 4 in the main text) reveals F c (100) ≈ 0.128 and F c (200) ≈ 0.119.The decrease of F c (N ) is qualitatively consistent with our expectation, however, an investigation of F c in a broader range of system sizes would be required to test the scaling prediction and the value of correlation length exponent ν = 0.8 ± 0.096 obtained from the scaling relation ν = 1/(d − d f ).

IX. RANDOM YIELDING TRANSITION IN FLAT BI-PERIODIC VERTEX MODEL
In order to test whether the spherical geometry affects the values of the critical exponents, we have measured the exponents τ , z and d f in a flat two-dimensional vertex model with bi-periodic boundary conditions.We find a power-law avalanche size distribution consistent with the exponent τ ≈ 1.35 we reported for the spherical geometry (

FIG. 1 .
FIG. 1.(a) Spherical vertex model tissue with N = 200 cells with randomly oriented traction forces (red arrows).(b) Traction force is generated by extending a spring of stiffness κ, at speed v in direction of the polarity pα.(c) Example of the tissue traction force magnitude dynamics F as a function of spring displacement vt.(d) Dynamics of ensembleaveraged tissue traction force magnitude.As the spring extension speed v approaches the quasi-static driving limit v → 0, the traction force magnitude averaged over ensemble realisations F converges to its critical value Fc marked by the dashed line (see SM).

FIG. 3 .
FIG. 3. Density of plastic excitations in the tissue.(a) Additional tension ∆f required to collapse the bond as a function of bond length .A linear scaling is observed (solid line).(b) Cumulative bond length distribution C( ) in the steady state for v = 2 • 10 −4 .The predicted value of the exponent θ ≈ 0.32 is indicated by the solid line.At low we observe a linear scaling of C( ) (dashed line), corresponding to a constant bond length distribution, as expected at finite v and for finite system sizes.

FIG. 4 .
FIG. 4. Steady-state flow curve measured in spherical vertex model networks with N = 100 (blue squares) and N = 200 (orange circles).Curves show best fit to v ∼ ( F − Fc) 1.41 for N = 100 (dashed line) and N = 200 (solid line), see SM for discussion of Fc finite size scaling.

1 S
S.2a).We then estimate the uncertainty in our measurements of avalanche size distribution by varying the time resolution δt between these extreme limits and identify the intermediate value of δt for which we test the robustness of the results upon further decreasing v (see Fig. S.2b).Note that varying the time resolution δt does not change avalanche duration scaling with its size (see Fig. S.2c).VI.FRACTAL DIMENSION OF AVALANCHES The m-th moment of avalanche size distribution P (S) reads S m = Sc m P (S)dS, (S.14) where S c is the cut-off size of the avalanches.Considering a power-law normalized avalanche size distribution P (S) ∼ S −τ (Fig. 2a in the main text), it follows S m+1 S m ∼ S c .(S.15) The system size dependent cut-off value in the normalized avalanche size distribution (Fig. 2a in the main text) scales with linear system size as S c ∼ N d f /d , where d f is the fractal dimension of avalanches, and d = 2 is the system dimension.The fit (shown by black lines in Fig. S.3a) yields d f = 0.75 ± 0.15.Taking the estimation of d f ≈ 0.75 and τ ≈ 1.35 (see main text) leads to a collapse of the tail in the avalanche size distribution for various system sizes as is shown in Fig. S.3b.VII.DENSITY OF PLASTIC EXCITATIONSIn our vertex model, we determine the exponent θ by fitting cumulative bond length distribution (see Fig.3bin the main text) for a wide range, and compare the measured values with the range predicted by the scaling relation (main text Eq. 6).Further analysis (Fig. S.4) shows that varying lower and upper limits of the fitting range leads to measurements of θ that vary in a range [0.21, 0.52].We find that both increasing the upper limits while fixing the lower limit (Fig. S.4a) and shifting up the one decade-long fitting interval (Fig. S.4b) lead to an increase in our measurement of θ.To test the quality of each measurement, we quantify the average of squared residuals of the fit r = (1/K) K j=0 (P ( j ) − θ j ) 2 , where sum is over K bonds with lengths in the fitting interval.We find that the confidence of the fit reduces as the measured values goes above the predicted range, marked by dashed lines in Fig. S.4.

FIG. S. 2 .
FIG. S.2.Avalanche statistics.(a) Uncertainty in detecting drops in the average traction force magnitude due to the finite time resolution.(b) Measured avalanche size distribution as a function of time resolution δt allows us to estimate the uncertainty of the avalanche distribution exponent τ = 1.35 ± 0.11.(c) Varying the time resolution does not change avalanche duration scaling with its size.Blacl lines indicate T ∼ S 0.68 .
Fig. S.5a).Moreover, the fractal dimension d f ≈ 0.75 measured in the spherical geometry leads to the collapse of the avalanche size distribution (see Fig. S.5b) and the same value is consistent with the finite-size scaling of the consecutive moments of avalanche size statistics (Fig. S.5d).Finally, we find that avalanche duration scales with avalanche size with an exponent z/d f ≈ 0.68, consistent with the value of exponent z reported in the spherical geometry (Fig. S.5c).

FIG. S. 3 .FIG. S. 4 .
FIG. S.3.Avalanche statistics.(a) Ratio of consecutive moments of avalanche size distribution scales with linear system size √ N as S m+1 / S m ∼ N d f /2 , with d f = 0.75 ± 0.15.(b) Collapse of avalanche size distribution for different system sizes with the exponents τ = 1.35, and d f = 0.75.

FIG. S. 5 .
FIG. S.5.Avalanche statistics in flat bi-periodic vertex model tissues.(a) avalanche size distribution P (S) is a power law with exponent τ ≈ 1.35 indicated by the dashed line.(b) Collapse of P (S) with d f = 0.75 for a wide range of tissue sizes.(c) Avalanche duration vs avalanche size.The solid black line shows T ∼ S z/d f , with z/d f ≈ 0.68.(d) Ratio of consecutive moments of avalanche size distribution scales with linear system size √ N as S m+1 / S m ∼ N d f /2 , with d f = 0.75 marked by the solid black lines.

TABLE I .
[39]critical exponents of RYT on a sphere in comparison with reported values for YT in a 2d elastoplastic model[39].