Controlled neighbor exchanges drive glassy behavior, intermittency and cell streaming in epithelial tissues

Cell neighbor exchanges are integral to tissue rearrangements in biology, including development and repair. Often these processes occur via topological T1 transitions analogous to those observed in foams, grains and colloids. However, in contrast to in non-living materials the T1 transitions in biological tissues are rate-limited and cannot occur instantaneously due to the finite time required to remodel complex structures at cell-cell junctions. Here we study how this rate-limiting process affects the mechanics and collective behavior of cells in a tissue by introducing this important biological constraint in a theoretical vertex-based model as an intrinsic single-cell property. We report in the absence of this time constraint, the tissue undergoes a motility-driven glass transition characterized by a sharp increase in the intermittency of cell-cell rearrangements. Remarkably, this glass transition disappears as T1 transitions are temporally limited. As a unique consequence of limited rearrangements, we also find that the tissue develops spatially correlated streams of fast and slow cells, in which the fast cells organize into stream-like patterns with leader-follower interactions, and maintain optimally stable cell-cell contacts. The predictions of this work is compared with existing in-vivo experiments in Drosophila pupal development.

Cell neighbor exchange is fundamental to a host of active biological processes, from embryonic development [1,2] to tissue repair [3,4]. Often referred to as cell intercalations [5], it is the leading mode of rearrangement in a confluent cellpacking, and in the simplest form described by a topological T1 process observed in foams [6]. During development, T1 rearrangements lead to diverse reorganization patterns [7][8][9][10][11], based on polarized global cues [12,13], or apparently random localized active fluctuations [14,15]. Disruption of neighbor exchanges lead to defects in developing embryos [12,13], while an increase has been shown to alleviate disease conditions [7,16]. It is, therefore, important that the underlying biological programs of neighbor exchanges are regulated in a tissue.
The steps associated with a T1 event, namely the shrinkage and restoration of a cell-cell junction, are driven actively by myosin II motors, present inside the cell and at the junctions [13,15,[17][18][19][20][21]. Different global and local signals drive this activity [13,15,22], and perturbing them can disrupt or enhance neighbor exchanges [12,13,15]. However, not much is known about how these signals determine the rate of neighbor exchanges at the level of a single cell and how these effects influence multicellular behavior.
While there is a scarcity of quantitative measurement of cell neighbor exchange rates [15], recent evidences suggest that limiting the rate of neighbor exchanges influences many tissue level properties. The fluidity of an epithelium [23] is affected by differential rates of neighbor exchanges in development [15,24]. Neighbor exchanges and dynamic remodeling of cell-cell contacts [25] are extremely crucial in numerous examples of collective cell migration through dense environments, examples include the migration of border cells in developing Drosophila oocyte [26], the formation of dorsal branches in Drosophila tracheae [27], the formation of Zebrafish posterior lateral line primordium [28] and the remodeling of lung airway epithelium under asthmatic conditions [29].
Invading cancer cells often form multicellular streams to migrate through narrow spaces in rigid tissues [30,31] that depends crucially on neighbor exchanges, as suggested recently [32].
Despite these evidences, it is not clear how the neighbor exchanges and dynamic remodeling of cell-cell contacts are related, and what roles they play in an emergent multicellular behavior. In this work we try to understand this relation within a theoretical approach. We extend the well-known 2D vertex model [33,34] which has already provided useful insights on morphogenesis [35,36], epithelial maturation [29], unjamming [37], and wound-healing [38]. In all these cases, the T1 events have been considered instantaneous or implicit. Furthermore, the interactions of active cell motility and the rate of T1 events have not been considered before to the best of our knowledge.
The focus of this work is to study the effect of cell motility and rate of T1 events regulated at the level of individual cells. We find that changing the intrinsic cell-level persistence time for T1 events induces a gradual slowing of the cellular dynamics. Though this is reminiscent of the dynamical arrest in glassy systems, the nature of the glassy state that emerges is distinct and not previously described. The interplay of cell motility and intrinsic persistence of T1 events gives rise to out-of-equilibrium states that behave largely like a glassy system, yet lacks conventional dynamic heterogeneity and surprisingly contains a population of mobile cells that can migrate via an unusual coordinated stream-like motion.

RESULTS
Model -Our appropriately modified "dynamic vertex model" (DVM) [37] retains most of the classic features [39,40]. Cells in a 2D epithelial monolayer are described by irregular polygons, defined using vertices, which constitute the (1) Due to the biomechanical interactions cells resist changes to their shapes, described by the tissue mechanical energy [34,39], where N is the number of cells, A i and P i are the area and the perimeter of cell i, respectively. A 0 and P 0 are the equilibrium cell area and perimeter, respectively, considered uniform across the tissue [41]. K A and K P are the elastic moduli associated with deformations of the area and perimeter, respectively. The second term in the right hand side of Eq. 2 yields a dimensionless target cell shape index p 0 = P 0 / √ A 0 [29,37]. Here we choose a constant p 0 = 3.6, which describes a solid-like tissue when cell motility is low [40]. The force on any vertex due to cell shape fluctuations is Each cell in the DVM behaves like a self-propelled particle with motility force v 0 [40,42] acting on the geometric cell center. The total active force on vertex i is whereñ i is the direction of average active force on vertex i ( Fig. S1 and Methods). The most important ingredient in our model is that we allow T1 processes with a time delay τ T 1 between successive T1 events (Fig. 1a). For each cell, a timer δt c is kept which records the time elapsed since the last T1 transition involving the cell. Then all edges adjacent to the cell can undergo a T1 if and only if the length of the edge is shorter than a threshold l min = 0.1 and δt c > τ T 1 . Initially, cells are seeded with a random δt c chosen from a uniform distribution [0, τ T 1 ]. This rule introduces a persistent memory governing intercalations of cells.
Dynamical slow down in cellular motion due to persistence in T1 events -According to previous results using similar models [40,42], solid-fluid transition in a motile tissue at p 0 = 3.6 happens for v 0 > 0.6. In order to understand the effects of single cell persistence time of T1 events, we start with a high motility fluid at this p 0 . To quantify the cell dynamics we compute the mean square displacements (MSD), ∆r(t) 2 of the geometric centers of the cells (Fig. 1b and S2a) as well as the self-intermediate scattering function, F s (q,t) (Fig. 1c), a standard measure in glassy physics to quantify structural relaxation (see Methods for calculation details). When τ T 1 is small, the cell motion is nearly ballis- tic at short times and diffusive at long times ( ∆r(t) 2 ∝ t) as expected in a typical viscoelastic fluid. This situation corresponds to nearly instantaneous T1 events typically considered previously [29,34,35,38]. Here the F s (q,t) decays sharply with a single timescale of relaxation indicating the tissue fluidizes quickly due to motility (Fig. 1c). Interestingly, as τ T 1 is increased the MSD starts to develop a 'plateau' after the early ballistic regime and requires increasingly longer time to become diffusive. The development of this plateau is indicative of the onset of kinetic arrest and shows up as a two-step relaxation in terms of F s (q,t) (Fig. 1c). In glassy physics [43][44][45][46], the origin of the two-step relaxation has been attributed to βrelaxation, taking place at intermediate times when each cell just jiggles inside the cage formed by its neighbors, which is also responsible for the plateau in MSD; and α-relaxation, taking place at later times when a cell is uncaged and undergoes large scale motion to make the MSD rise after the plateau. We extract τ β , the timescale of β-relaxation [47] by locating the point of inflection in MSD curve on a log-log plot (see Methods and Fig. S2b). The time scale τ α associated with αrelaxation can be extracted from F s (q,t) according to convention (Methods). In Fig. 1d we show the dependence of τ β and τ α on τ T 1 . For low τ T 1 (fast T1's), the two timescales coincide indicating τ T 1 has no discernible effect on the dynamics but starting at τ T 1 ≈ 20 and onward the difference between the two timescales grow drastically. This gives a characteristic timescale where the hindrance of T1's causes effective caged motion and kinetic arrest in cells. As these behaviors are remi-niscent of the onset of glassy behavior in super-cooled liquids, the natural question is whether introducing a T1-delay merely provides another route leading to a conventional glassy state, similar to e.g. lowering the temperature? To answer this question, we quantify the dynamical heterogeneity in states suffering kinetic arrest due to T1-delay and compare with more conventional glassy states obtained by lowering v 0 at short τ T 1 .
We measure the four-point susceptibility χ 4 (t) which is a conventional measure of dynamical heterogeneity in glassy systems [47,48]. For nearly instantaneous T1's ( Fig. 2a), the tissue behaves like a conventional super-cooled glass as v 0 is decreased: χ 4 (t) exhibits a peak which shifts towards larger times with decreasing v 0 . Together with the increase in the peak magnitude of χ 4 (t), these results indicate the lengthscale and timescale associated with dynamic heterogeneity become increasingly larger due to lowering v 0 , which plays the role of an effective temperature [40]. Using the peak value χ * 4 ( Fig. 2c), the glass transition can be located at v 0 ≈ 0.8. In contrast, the behavior of χ 4 (t) for large large T1 delay is very different. Here, χ 4 (t) develops a plateau over decades in time (Fig. 2b), which is significantly broader than in the low τ T 1 regime. Further, hindering T1's here has also reduced the dependence of dynamic heterogeneity on v 0 and the peak value of χ * 4 ( Fig. 2c) no longer exhibits any peaks as function of v 0 . This suggests that v 0 has been replaced as the ratelimiting factor determining cell rearrangements, which are instead dominated by τ T 1 . This analysis thus reveals that while the high τ T 1 regime becomes kinetically arrested similar to a super-cooled liquid, it does so in a manner distinct from lowering the temperature. Here, the effects of effective temperature (v 0 ) on dynamics get 'washed away' and the characteristic timescale is mostly set by τ T 1 . Therefore, we essentially have a glassy state (with a τ α that can be directly controlled using τ T 1 ) with a low degree of dynamic heterogeneity.
Intermittency in T1 events points to a dynamic regime distinct from a conventional glass -Since the origin of a growing dynamic heterogeneity in glasses is attributed to highly intermittent motion of individuals [46], we explicitly investigate how intermittency of T1 rearrangements depends on v 0 and τ T 1 . By tracking all neighbor exchanges we maintain a time dependent list {τ i w } illustrated in Fig. 1a. We define an instantaneous waiting timeτ inst w (t) (Fig. 2d-e) by averaging the entries in this list at a given time t. We also consider distributions of the waiting times, P(τ w ) (Fig. S3).
When T1 rearrangements are nearly unconstrained and motility is high,τ inst w (t) exhibits steady fluctuations about a mean value that is slightly larger compared to the lower bound set by τ T 1 . As motility is lowered T1 occurrences become less frequent andτ inst w (t) increases overall. Near the glass transition, the dynamics becomes highly intermittent, and consequently, T1 events are separated by a broad distribution of waiting times (Fig. S3) and multiple T1 events can take place simultaneously akin to avalanches [48]. Deeper in the glass phase (v 0 < 0.8) we seeτ inst w (t) increases slowly over time, reminiscent of the aging behavior observed in glassy materials [46,48]. The intermittent fluctuations are dampened significantly as T1 delay becomes very large (Fig. 2e).
To further characterize the nature of fluctuations and the intermittency, we compute the kurtosis κ of the τ w distributions. In the small τ T 1 regime, the distribution P(τ w ) decays with power-law tails for v 0 < 1 (Fig. S3), and this is also the location where κ vs v 0 exhibits a very pronounced peak (Fig. 2f) at the point of glass transition (v 0 ≈ 0.8). Away from the glass transition (v 0 > 0.8), κ ≈ 3 corresponding to Gaussian fluctuations (Fig. 2d,e and Fig. S3). However, when τ T 1 is increased, the intermittency fades away and the peak in κ disappears.
Universal scaling separates fast and slow regimes -Next we analyze the interplay between the T1 delay and the observed mean waiting times τ w (Fig. 3a). When τ T 1 is small, the ratio τ w /τ T 1 remains always larger than unity and varies by orders of magnitude depending on v 0 . As the delay increases, τ w /τ T 1 approaches unity and depends weakly on v 0 . This behavior suggests the following universal scaling ansatz In Eq. (5), f (x) is the dynamical crossover scaling function with x = τ T 1 v z 0 . Re-plotting the data using Eq. (5) we find good scaling collapse with exponent z = 4 ± 0.1 (Fig. 3a-inset). This uncovers two distinct regimes. For small x, f (x) ∼ x −1 , which implies τ w ∝ v −z 0 . We refer to this as the fast rearrangement regime. For large x, f (x) → 1 indicating a slow rearrangement regime where τ w ∝ τ T 1 and independent of v 0 . The transition between the two regimes occurs at x = x * ≈ 100, corresponding to a scaling relation of τ T 1 = x * v −z 0 that constitutes the boundary separating these two rearrangement regimes.
The dynamics in the fast regime is not hindered by τ T 1 and is largely driven by the effective temperature T e f f ∝ v 2 0 /D r [40]. For high values of v 0 in this regime, the motile forces are sufficient to overcome the energy barrier associated with T1 transitions; however as v 0 drops below v * 0 ≈ 0.8, the tissue enters into a glassy solid state where cell motion become caged.
In the slow regime, τ T 1 dominates the dynamics as it is the longest timescale in the system and therefore the ultimate bottleneck to rearrangements. This leads to τ w to depend linearly on τ T 1 while being insensitive to motility. In addition to glassy-behavior which is a source of non-equilibrium fluctuations, here τ T 1 constitutes another possible route that can take the system out-of-equilibrium, effectively slowing down the dynamics.
Here, the uncaging timescales (τ α ) grows with τ T 1 (Fig. 1d) but remain quite finite, as in a highly viscous or glassy fluid. Phase diagram of cellular dynamics -A summary of these results suggest that the phase diagram on the v 0 − τ T 1 plane (Fig. 3b) can be categorized into three phases: • glassy solid, • fluid, and an unusual • active streaming glassy state (ASGS), to be discussed in depth below. The solid-fluid phase boundary is given by peak-positions in κ (Fig. 2f) and χ * 4 ( Fig. 2c), and effective diffusivity, D e f f (Methods, Fig. 3b, and Fig. S2) [40]. The boundary between fluid and the ASGS phase is given by the crossover scaling (Eq. 5) between the slow and fast regimes.
The ultra-low D e f f values in the ASGS phase resemble a solid more than a fluid. However the finite τ α (Fig. 1d) and the nature of waiting times (Fig. S2) indicate that these states are only solid-like for timescales less than τ T 1 . Furthermore, it appears that the eventual transition to diffusive motion (at times > τ T 1 ) occurs in a highly heterogeneous manner and the motion is dominated by a small population of fast moving cells. To understand this peculiar nature of fluidity we quantify the fraction of relatively fast moving or mobile cells, given by f mob (Methods). As expected, f mob is large and unity for fluid and zero for glassy solid. However, in the ASGS phase, although f mob decreases slowly as τ T 1 increases, curiously enough, it remains finite. These evidences suggest that the ASGS phase always includes states with a distribution of fast and slow cells. This could be signature of another kind of growing heterogeneity which we try to understand next. Efficient fast cells organize into cellular streams -We observe that depending on the phase they belong to the fast cells exhibit different degrees of neighbor exchanges. To quantify this, we define a single cell intercalation efficiency I (Fig. 4a and Methods). For fast cells I increase with increasing τ T 1 at a fixed motility (Fig. 4a). At the same time, the spatial distribution of the fast cells become increasingly heterogeneous, and stream-like. The distributions of I (Fig. 4b) also reflect this change, as they become wider and heavy-tailed, with 2- 5 fold increase in the mean efficiencyĪ as τ T 1 goes up. The fast cells also exhibit mutual spatiotemporal alignment (Methods), which grows rapidly with τ T 1 (Fig. 4c). The overall alignment probability φ a (Methods, Fig. 4d) can be used as an order parameter for streaming behavior. At small motilities, φ a is insensitive to τ T 1 and large due to collective vibrations in a solid-like tissue [40]. At higher motilities, φ a becomes increasingly dependent on τ T 1 . Taken together, these analyses pinpoint the necessary conditions for stream-like behavior: increasing T1 delay alone can induce alignments, but higher motilities are also crucial.
To visualize this streaming behavior we follow cell trajectories ( Fig. 4e-top), which are uniformly distributed and randomly oriented for low τ T 1 , but gradually become sparse and grouped into stream-like collectives as τ T 1 increases. These collectives are highly correlated and persistent. Interestingly, here we have used a time interval much smaller than the corresponding β-relaxation timescales to compute the displacement vectors of cell centers. Therefore, this type of collective behavior occurs even before the cells uncage.
To quantify the spatial correlations arising during streaming we adapt a quantity V [49], computed in the co-moving frame of a given cell, which represents the average velocities at different locations around it. Collective migratory behavior shows up as vectors of similar length and direction near the cell, whereas solid or fluid-like behavior results in isotropic organization of vectors of uniform sizes. In Fig. 4e-bottom, signatures of collective motions are absent when τ T 1 is small, but gradually appears with increase in τ T 1 as the sizes and orientations of the vectors become more correlated along the direction of motion and remain uncorrelated along the direction perpendicular to it. Such anisotropic vector fields are hallmarks of cellular streaming [49] that involves high frontback correlations (Fig. 4f-top), and low left and right corre-lations ( Fig. 4f-bottom). This also provides the size of a typical stream, about 8 cells long and 4 cells wide. We further find that these streams are driven by leader-follower interactions [50,51] that increase with τ T 1 (Fig. S4). In this light, after quantifying the streaming behavior, we consider its possible underlying mechanisms. One big candidate is dynamic remodeling of cell-cell contacts.

Delayed T1 events result in effective cell-cell cohesion -
The delays between successive T1 events in our model essentially provide a dynamical way to control stability of cell-cell contacts. Two neighboring cells can maintain an effective adhesion for a time-period equal to τ T 1 or longer, depending on the time-evolution of the shared junction. In simulations, we have direct access to the characteristic time of this effective adhesion, τ ad , which can be obtained from the cell-cell adjacency information. Alternatively, a more experimentally accessible measurement can be performed using cell tracking. Fig. 5a illustrates quantification of cell-cell cohesion in our simulations using both techniques. When τ T 1 is small, cellcell contacts are frequent but short lived compared to the large τ T 1 case for which the contacts become long-lived. The distributions of τ ad 's ( Fig. 5b) are very narrow for low motilities, with a single peak at τ ad ≈ τ total , the total simulation time, as expected in the glassy solid state. Increase in motility widens the distributions and introduces a new peak at τ ad ≈ τ w . At even higher motility, the peak at τ total disappears. The new peak shifts to smaller values but occurs at τ ad > τ T 1 in the fluid phase, and at τ ad = τ T 1 in the ASGS. The mean adhesion times τ ad also reflect these trends (Fig. 5c) and the two methods of extracting τ ad 's agree nicely. An important feature is that the variance of τ ad is always large for high motility states, indicating a heterogeneity in the remodeling of cell-cell contacts. In the ASGS, τ ad has a regime of strong dependence on τ T 1 (Fig. 5c, bottom), indicating tunable window of The emergence of cell streaming from the effective cellcell cohesion is an amplified tissue level response to the persistence memory introduced by the T1 delay. This can be shown by inspecting the consequences of a single T1 event (Fig. 5d). It not only resets the timers on the two cells swapping the shared junction, but also resets timers on nearby cells that now become neighbors after the swap. These events introduce a cohesive 4-cell unit that is stable for a time period of at least τ T 1 . This time-restriction also sets a natural lower bound to the τ ad 's, as reflected in Fig. 5b. A consequence of this effective cohesion is that these 4 cells can now only move coherently for a period τ T 1 or more (Fig. 5e). The only way this cluster can move is via T1 transitions involving junctions at the periphery, as illustrated in Fig. 5e-i. This also allows this cluster to grow into a stream only when τ T 1 is optimal, because very frequent or rare T1 events would not be useful. That is why we see a drop in f mob (Fig. 3b) and saturation of τ ad (Fig. 5c) as τ T 1 become very high. The movements of such streams are always unidirectional at any given timepoint which automatically introduces a leader-follower interaction among these cells, characterized in Fig. S4.
The above analysis also reveals that different glassy regimes (Fig. 3b) can be distinguished based on measuring Prediction of T1 rates and comparison with experimental data on Drosophila pupa development. Mean rate of T1 events and mean tension, calculated per junction from our simulations at v 0 = 1.6 as function of τ T 1 . The dashed lines represent the experimentally measured rates of neighbor exchange in Drosophila pupal notum under three situations from Ref. [15]: the wild type notum (red), with over-expression of a constitutively active Rho-kinase (Rok-CAT, green) that enhances junctional myosin-II activity and hence increases tension, and with a Rho-kinase RNAi (Rok-RNAi, blue) that reduces junctional myosin-II activity and hence decreases tension. The observed changes in T1 rate due to these two perturbations are consistent with our predicted relationship between junctional tension and T1 rates.
the characteristic time of cell-cell cohesion and its statistical distribution P(τ ad ). For example, consider two different glassy states, one with a large τ T 1 and one with small motility. In terms of conventional measures both states would exhibit glassy features (vanishing D e f f , large χ 4 ), however, our work predicts that the state with a large τ T 1 should have a broadly distributed P(τ ad ) compared to the motility driven glass transition. Predictions for Drosophila pupa development -Cell intercalations in our model are spatially uncorrelated with no directional polarity. Such unpolarized cell intercalations have been discovered recently in Drosophila notum in the early pupal stage [15]. Evidences indicate that random fluctuations in junction lengths, strongly correlated with activity of junctional myosin-II, drive such unpolarized intercalations. Rate of intercalation drops as the pupa ages, and at the same time there is an increase in junctional tensions. We observe a similar decrease in the rate of T1 events in the model as τ T 1 increases, accompanied by an increase in junctional tension (Fig. 6) (Methods), consistent with the experimental observation. We can also compare the observed trends of experimental perturbations to the wild-type tissue with our results and predict that increasing or decreasing τ T 1 in our model influences T1 rates in a manner similar to perturbing junctional tension by up-regulating or inhibiting junctional myosin-II ac-tivity, respectively.

DISCUSSION AND CONCLUSION
To model collective dynamics of confluent epithelial cells, we introduced an inherent timescale (τ T 1 ) for cells to undergo rearrangements based on the important observation that T1 events in real tissues do not occur instantaneously. When τ T 1 is short compared to other timescales in the model, we recover a motility-driven glass transition, occurring at motilities large enough to overcome the energetic barriers that cause a cell to become caged. Near this glass transition, we observe highly intermittent cell motion, as well as neighbor exchanges events concomitant with growing dynamical heterogeneity.
However, when τ T 1 grows we discover a rich dynamical regime where the system appears glassy on timescales governed by τ T 1 but become fluidized at longer times. Compared to the motility-driven glass transition, this regime has a completely different kind of heterogenous glassy behavior characterized by disappearance of conventional glass-fluid boundary and appearance of spatially distributed pockets of fast and slow cells. The origin of this new glassy behavior stems from the effective cell-cell cohesion caused by the inability to undergo local cellular rearrangements. Surprisingly, this local frustration actually serves to enhance collective migration by facilitating stream-like patterns, reminiscent of leader-follower behavior, albeit without any explicit alignment interactions between cells. Interestingly, this connection diminishes as mean effective adhesion times saturate when persistence of rearrangements becomes too low or too high ( Figure 5c). Therefore, the cell-cell contacts need to be dynamic, but optimally stable to maintain the streaming mode. Note that the stream-like motion we observe occurs in the absence of any dynamic heterogeneity and this sets it apart from string-like cooperative motion observed in a 3D supercooled liquid [52] which is a direct consequence of dynamic heterogeneity.
This result, consistent with the existing biological mechanisms [31,53], sheds new light on the nature of cell-cell adhesion associated with cell streaming in cancerous tissues [31].
The effective adhesion picture that emerges from our study is also consistent with the recent observations regarding polarized intercalations in extending germ-band of Drosophila embryo [18]. Oscillations in adherens junction protein Ecadherin has been deemed necessary for successful intercalations, and inhibition of this oscillation led to a decrease in successful T1 events. Our predictions reverse engineer this effect by limiting the rate of T1 events which leads to an increase in effective cell-cell adhesion and stability of junctions.
Increasing the E-cadherin levels at the junctions would be the biochemical way of achieving this as shown before [18].
The unusual nature of the active streaming glassy state has multitude of implications. The slow but finite structural relaxation gives the material a tunable viscosity which can be extremely useful for preparation of biology-inspired sheet-like objects of controllable stiffness [54,55]. The control of T1 rate can be potentially translated to gene-level control of the signaling [56] associated with developmental events and disease conditions that strongly depend on cell intercalations, e.g. body-axis extension and kidney-cyst formation [2], respectively. This might even allow design of organisms with programmable development [57], or one with controlled diseasespreading rates.

ACKNOWLEDGMENTS
We acknowledge the support of the Northeastern University Discovery Cluster and the Indo-US Virtual Networked Joint Center project titled "Emergence and Re-modeling of force chains in soft and Biological Matter No. IUSSTF/JC-026/2016.

Active cell motility in DVM -
The explicit form of the active force on any vertex i depends on the motility force contributions from the adjacent cells: where a c = l c /(2z i ∑ c↔i l c ) is the weight associated with cell c adjacent to the vertex i, l c is the total length of the edges shared by vertex i and cell c, z i is the connectivity at vertex i, and the factor 2 takes care of double contributions from the same cell. This averaging scheme is different than the recent approaches using a flat average of active forces on adjacent cells [29,58,59] and ensures that the active force on vertex i is dominated by the cell sharing the longest edges attached to i. The polarization of the self-propulsion on cell−c is given by n c = (cos θ c , sin θ c ), where the polarization angle θ c is perturbed only by a white-noise [60][61][62][63][64]: where the ζ θ is a Gaussian white noise with zero mean and variance 2D r which sets the repolarization timescale for the cells in our model given by 1/D r . We use a fixed D r for all cells throughout the present study. Simulation details -Our simulations are overdamped dynamics of 256 cells periodic boundaries along both x and y directions. We use √ A 0 as the unit of length, K P A 0 as the unit of energy and ζ/K P as the unit to measure time t in our simulations. The vertex positions are updated by solving Eq. 1 using Euler's scheme. Our dynamical simulations are initialized from random Voronoi configurations which have been subjected to energy minimization using the conjugategradient algorithm. All simulations have been done in the Surface Evolver program [65] with a fixed equilibrium cell area A 0 =Ā = 1 (Ā, the mean cell area), time step of integration ∆t = 0.04. We run each of our simulations for ∼ 10 6 steps and collect data for subsequent analyses after the tissue properties like the mechanical energy reaches steady state. We scan the following parameter space: v 0 ∈ [0.2, 2] and τ T 1 ∈ [0. 4,2000] at a fixed D r = 0.5. For each combination of v 0 and τ T 1 we perform 10-20 independent simulations. Measuring junctional tension -Tensions arise in the vertex model due to mismatch between the actual cell perimeters and the equilibrium perimeter. It can be defined in terms of the preferred scaled perimeter or target cell shape index p 0 . We calculate tension on any junction shared between cells i and j by where p i and p j are respective scaled cell perimeters, measured during the simulation. Determination of τ β -We extract the β-relaxation timescale from any given MSD vs t plot by locating the minimum in the time derivative of MSD given by: dln( ∆r 2 (t) ) dln(t) [66]. Plots of this time derivative for MSDs shown in Fig. 1b are shown in Fig. S2b.

Analysis of self-intermediate scattering function F s (q,t) -
We have used the following definition of self-intermediate scattering function: F s (q,t) = e iq·∆r(t) where q is the wave vector corresponding to our length scale of choice and the angular brackets represent ensemble average and averages over angles made by q and ∆r(t), cell displacement vectors for time-delay t. Conventionally, emphasis is given to the behavior at q = 2π/σ, where σ is the inter-particle particle separation for configurations with particles just touching. However, we choose to focus on a even more restrictive wave vectorq = π/σ, corresponding to a length scale of 2 cell diameters. The reason for this is a technical one: since the degrees of freedom in the DVM are the vertices rather than the cell centers, the cell centers (calculated at every step based on vertices) can exhibit unusual fluctuations even when cells are completely caged. The current choice ofq allows us to consider relaxations where these artificial fluctuations contribute much less. To eliminate contributions from any local or global drift we consider temporal changes in nearest-neighbor separations as ∆r(t), instead of pure displacements of cell centers.
We define the α-relaxation timescale τ α as follows: F s (q, t = τ α ) = 0.2, following the definition used recently by another cell-based model study on similar tissues [67]. Definition of quantities associated with mobility -We have used several order parameters to describe different regimes of cell dynamics in our model tissue. Below we define them one by one. We define effective diffusivity by: D e f f = D s /D 0 where D s = lim t→∞ ∆r(t) 2 /4t, the long time self-diffusion coefficient and D 0 = v 2 0 /2D r , the free diffusion coefficient of a single isolated cell. D s has been computed using the value of mean-square displacement at the maximum time delay allowed in our simulation.
To define the fraction of mobile cells f mob , we follow individual cell MSD and define the net displacement of a cell: d ∞ = lim t→∞ ∆r 2 (t) . Then we find out the number of cells N mob that have d ∞ 2 cell diameters. Finally, f mob = N mob /N. This definition is consistent with the definition of mobile particles used in Ref. [68].
We define individual cell intercalation efficiency as I = d ∞ /n T 1 where n T 1 is the net T1 count for the given cell in the entire simulation. Analysis of orientation alignment in cell trajectories -To capture the orientational order and spatial organization of fast cells, we concentrate on the cells with intercalation efficiency I ≥Ī, the mean intercalation efficiency. This gives us a list of fast cells. Then we consider the whole simulation trajectory and calculate the probabilities of the angle, θ a between the instantaneous displacement vectors of any cell pair chosen from our list of fast cells, where the cell-cell separation is less or equal to 2 cell diameters. We consider all possible cell pairs satisfying this criterion and pool all the θ a to generate the probability density function P(θ a ). We define the net alignment probability φ a , by the following: where t c = 30 • . We follow a recent analysis [51] that probes how the motion of any two cells are correlated as a function of both their distance and the direction from one cell, chosen as reference, to another. We find that this directional correlation, shown above for different τ T 1 values at v 0 = 1.6, is low and nearly uniform in all directions, and decays very sharply within 1-2 cell diameters for low τ T 1 . As τ T 1 increases we see the correlation contours getting longer-ranged, polarized and elongated along the direction of motion of the reference cell. These features are classical signatures of predominant leader-follower behavior. For τ T 1 ≥ 100 we see very strong, anisotropic directional correlations that remain significant even at 6 cell diameters away from the reference cell. This approximate correlation length is consistent with that found from vector field analysis of cell streaming. These results show that the fast cells play major roles driving the cellular streaming observed for very large T1 delays.