Effect of cellular rearrangement time delays on the rheology of vertex models for confluent tissues
EEffect of cellular rearrangement time delays on the rheology of vertex models forconfluent tissues
Gonca Erdemci-Tandogan
1, 2, ∗ and M. Lisa Manning
1, 3, † Department of Physics, Syracuse University, Syracuse, NY 13244, USA Institute of Biomedical Engineering, University of Toronto, Toronto, ON M5S 3G9, Canada BioInspired Institute, Syracuse University, Syracuse, NY 13244, USA
Large-scale tissue deformation during biological processes such as morphogenesis requires cellularrearrangements. The simplest rearrangement in confluent cellular monolayers involves neighborexchanges among four cells, called a T1 transition, in analogy to foams. But unlike foams, cellsmust execute a sequence of molecular processes, such as endocytosis of adhesion molecules, tocomplete a T1 transition. Such processes could take a long time compared to other timescales inthe tissue. In this work, we incorporate this idea by augmenting vertex models to require a fixed,finite time for T1 transitions, which we call the “T1 delay time”. We study how variations in T1delay time affect tissue mechanics, by quantifying the relaxation time of tissues in the presence ofT1 delays and comparing that to the cell-shape based timescale that characterizes fluidity in theabsence of any T1 delays. We show that the molecular-scale T1 delay timescale dominates over thecell shape-scale collective response timescale when the T1 delay time is the larger of the two. Weextend this analysis to tissues that become anisotropic under convergent extension, finding similarresults. Moreover, we find that increasing the T1 delay time increases the percentage of higher-foldcoordinated vertices and rosettes, and decreases the overall number of successful T1s, contributingto a more elastic-like – and less fluid-like – tissue response. Our work suggests that molecularmechanisms that act as a brake on T1 transitions could stiffen global tissue mechanics and enhancerosette formation during morphogenesis.
INTRODUCTION
In processes such as development and wound healing,biological tissues must generate large-scale changes to theglobal shape of the tissue [1]. In confluent tissues, wherethere are no gaps or overlaps between cells, such globalchanges necessarily correspond to either changes in in-dividual cell shape or cell rearrangements [2], and large-scale deformation almost always requires a large numberof rearrangements.For epithelial monolayers, the geometry of most suchrearrangements is quite simple: viewing the apical sideof the layer, four cells come together at a single four-foldvertex, which subsequently resolves into two three-foldvertices where cells have exchanged neighbors. This pro-cess is called a T1 transition, adopted from the literatureon foams [3]. In some tissues, it is also common to ob-serve higher-fold vertices called rosettes [4, 5].Historically, there have been different perspectives onhow to understand and quantify such rearrangements.At the molecular scale , a concerted sequence of processesmust occur to allow such a change, including localizationof non-muscle myosin and actin to shorten interfaces [6–10], unbinding of adhesion molecules and trafficking awayfrom the membrane via endocytosis [11], exocytosis ofadhesion molecules to newly formed interfaces and newhomotypic binding, and reorganization of the cytoskele-ton to stabilize the new edges. Moreover, molecules suchas tricellulin [12, 13] are known to localize at tricellularjunctions and must be reorganized [14, 15]. In additionto all this, there is recent evidence that some cell types possess mechanosensitive machinery that will only trig-ger this molecular rearrangement cascade if tension onthe interface is sufficiently large [16].A complementary perspective has focused on the col-lective behavior of cells in a tissue. Specifically, a combi-nation of theoretical modeling [17, 18] and experimentaldata [19–21] has suggested that the collective mechanicsof a tissue has a huge impact on the rate of cell rearrange-ments, and that the collective mechanics are dominatedby a simple observable, the cell shape.Shapes of cells in a confluent tissue are generated by abalance between contractility generated by the cytoskele-ton and adhesion generated by molecules such as cad-herins, as well as active force generation by cells [21, 22].This suggests that cell rearrangement rates are governedby cell-scale features, such as overall expression levelsof adhesion and cytoskeletal machinery. Moreover, itis possible to identify a “collective response” timescale τ α [23] that describes the typical timescale over whichcells change neighbors. In both isotropic and anisotropictissues, this timescale depends on cell shape and align-ment [19, 21, 22].Given the strong experimental support for cell-scaleand molecular-scale perspectives, we hypothesize thatboth kinds of mechanisms must be working in concertto drive cell rearrangement rates in confluent tissues.While vertex models have been successful in predictingrearrangement rates in many cases, standard versions ofvertex models are missing the fact that specific molecularcascades and triggers are needed to allow rearrangements.Specifically, standard vertex model formulations requirethat cell neighbor exchanges proceed instantaneously af- a r X i v : . [ phy s i c s . b i o - ph ] F e b ter the creation of a higher-fold coordinated vertex. Insystems where molecular mechanisms delay cell neighborexchanges after the creation of higher-fold vertices, ver-tex models will make incorrect predictions for the globaltissue response.An extreme example is the amnioserosa tissuethat is required for proper germband extension in Drosophila [24]; although the cell shapes become ex-tremely elongated in the tissue [25], and standard ver-tex models would predict high numbers of cell rearrange-ments in response, experiments demonstrate that cellsdo not change their neighbor relationships at all. This isimportant for amniocerosa function: an elastic responsegenerated when cells maintain neighbor relationships al-lows the amniocerosa to pull strongly on the germbandtissue and elongate it [24, 26, 27].While experimental work is ongoing to understandthe precise molecular mechanisms that prevent cell re-arrangements in the amnioserosa, there have been recentattempts to augment vertex models by incorporating var-ious molecular mechanisms that affect how cells exchangeneighbors. One model already highlighted above stud-ies how a stress or strain threshold for T1 transitionsaffects cell shapes and rates of cell rearrangement [16].Other models incorporate strongly fluctuating line ten-sions [28, 29] that can trap edges so they cannot executeT1 transitions.In this manuscript, we instead augment vertex modelswith a model parameter we term the “
T1 delay time ”,which is intended to incorporate a broad range of molec-ular mechanisms that act as a brake on T1 transitions.For every situation where the cell-scale dynamics gener-ate a four-fold or higher coordinated vertex, the vertexis prevented from undergoing a T1 transition for a du-ration we call the T1 delay time. A similar mechanismhas also been studied in independent concurrent workby Das et al [30], which focuses on how controlled T1timescales generate intermittency and streaming statesin glassy isotropic tissues. Here, we study how such amechanism alters the global response of a tissue, bothin isotropic tissues and in tissues where there is a globalanisotropic change to tissues shape, such as during con-vergent extension in development.We find that the “molecular-scale” T1 delay timescaledominates over the “cell-scale” collective responsetimescale when the T1 delay timescale is the larger ofthe two, slowing down tissue dynamics and solidifyingthe tissue. In addition, we find that increasing T1 delaytime enhances the rosette formation in anisotropic sys-tems. This suggests that organisms might utilize specificmolecular processes that act as a brake on the resolutionof four-fold and higher coordinated vertices in order tocontrol the global tissue response in processes such aswound healing and convergent extension.
MODEL
Vertex model with T1 rearrangement time.
Weintroduce a new model parameter, T1 delay time, both inisotropic and anisotropic vertex models. A vertex modeldefines an epithelial tissue as confluent tiling of N cellswith an energy functional for the preferred geometries ofthe cells [17, 22, 31] E i = N (cid:88) i K A ( A i − A i ) + K P ( P i − P i ) . (1)In this definition, both the area and perimeter of a cellact like an effective spring. Here, A i and A i are theactual and preferred areas of cell i while P i and P i arethe actual and preferred perimeters. Using open-sourcecellGPU software [32], we simulate over-damped Brown-ian dynamics, where the positions of the cell vertices areupdated at each time step according to∆ r αi = µF αi ∆ t + η αi (2)using a simple forward Euler method. Here F αi = −∇ i E is the force on vertex i in α direction, µ the inversefriction, ∆ t the integration time step and η αi is a nor-mally distributed random force with zero mean and (cid:104) η αi ( t ) η βj ( t (cid:48) ) (cid:105) = 2 µT ∆ tδ ij δ αβ . The temperature T setsa thermal noise on vertices of each cell. The integrationtime step is set to ∆ t = 0 . τ where τ is the natural timeunit of the simulations: τ = 1 / ( µK A A ).As is standard in the literature, for instantaneous cell-neighbor exchanges, a T1 transition proceeds wheneverthe distance between two vertices is less than a thresh-old value, l c = 0 .
04 in natural simulation units. Wehave checked that our results are not sensitive to theprecise value of this cutoff. We set K A = 1, K P = 1, µ = 1 and all cells to be identical so that A i = A and P i = P . We nondimensionalize the length by the nat-ural unit length of the simulations l = √ A . This yieldsa target shape index or a preferred perimeter/area ratio p = P / √ A .The T1 delay time, t T , is a finite T1 rearrangementtime that acts as a brake on T1 transitions, which couldarise from a broad range of molecular mechanisms as dis-cussed in the introduction. We implement this feature inour simulations so that when the cell-scale dynamics gen-erate a four-fold or many-fold vertex (cellular junctionsthat satisfy l < l c criteria), the vertex is prevented fromundergoing a T1 transition for the duration of the T1delay time, t T (Fig. 1(A,B)). While the edge waits for a t T time, the configuration is not on hold – the systemevolves according to Eq. 2 and edges can still lengthenand shorten. In particular, every edge of a cell is associ-ated with two timers (one for each connected vertex) tokeep track of the time delays. When the timer reaches t T , if an edge l is still less than a critical length l c , theassociated cell undergoes the T1 process. Anisotropic Vertex Model.
We introduceanisotropy in our model using two different sets of simu-lation methods, aniosotropic line tensions and shear, sim-ilar to the protocols described in [21]. For the first setof simulations, we introduce an additional line tensionto nearly vertically-oriented edges. This type of pertur-bation was first developed to model dynamic anisotropicmyosin distribution due to planar cell polarity pathwaysin the germband extension in
Drosophila [9, 33]. This isimplemented via a vertex model energy functional: E i = N (cid:88) i K A ( A i − A i ) + K P ( P i − P i ) + (cid:88)
0, thenwe increase it by 5 × − increments at every simulationstep. After each shear, we minimize the system by updat-ing the vertex positions but keeping the box dimensionsfixed. RESULTS
Relaxation time of the tissue.
To quantify theglobal mechanical properties of the tissue, we character-ize a relaxation time ( τ Sα ) as a function of T1 delay time t T and target shape index p using the decay of a self-overlap function. The self-overlap function is a standardcorrelation function used to quantify glassy dynamics inmolecular and colloidal materials [34]. It represents thefraction of particles (vertices) that have been displacedby more than a characteristic distance a in time t, Q s ( t ) = 1 N N (cid:88) i =1 w ( (cid:12)(cid:12) r i ( t ) − r i (0) (cid:12)(cid:12) ) (5)where r i is the position of vertex i and the function w , w ( r ≤ a ) = 1 and w ( r > a ) = 0. The characteristic relaxation time of the system, τ Sα , is the time which mostof the vertices are displaced less than a characteristicdistance 1 /e : Q s ( τ Sα ) = 1 /e .We run simulations across a three-dimensional parame-ter space: T1 delay time t T , target shape index parame-ter p = P / √ A , and temperature T , with 100 indepen-dent simulations initialized from different configurationsfor each point in parameter space. All simulations arethermalized at their target temperature for 10 τ beforerecording the data. We then run simulations for addi-tional 3 × τ to ensure the system reaches a steadystate.We first use the self-overlap function to compute thecharacteristic relaxation time of the tissue, τ Sα , in theabsence of any T1 delays. This is simply the time atwhich the self-overlap function decays to 1 /e of its orig-inal value for simulations where the T1 rearrangementis instantaneous, t T = 0. Therefore, τ Sα is the typicalcollective response timescale that depends on cell shapeand alignment in vertex models. Fig. 1(C) shows thistimescale as a function of the target shape parameter p with T = 0 .
02. As expected, τ Sα decreases monotonicallyas p increases, demonstrating that lower values of p areassociated with glassy behavior and increasing relaxationtimes.Next, we study the behavior of the self-overlapfunction as a function of the T1 delay time.Fig. 1(D) shows that the self-overlap function changesa function of the T1 delay time, with t T =0 , . , . , . , . , . , . , . τ (darkgreen to yellow) for fixed p = 3 .
74 and T = 0 .
02. Fromthis data and additional simulations at other values of p , we extract the characteristic characteristic relaxationtime in the presence of T1 delays, τ Sα . The inset toFig. 1(E) shows the behavior of τ Sα as a function of t T for different values of p = 3 . , . , . ... . T = 0 .
02 and N = 256.As τ Sα represents the inherent relaxation timescale ofthe tissue controlled by p in the absence of T1 delays,we attempt to collapse this data by rescaling both the T1delay timescale t T and the observed relaxation timescale τ Sα by the inherent timescale for each value of p . Thedata collapses, showing that the mechanical properties ofthe tissue remain unchanged for any t T delay time be-low ∼ τ Sα . For t T > ∼ τ Sα , the relaxation timeincreases significantly, with a slope approximately equalto unity, indicating τ Sα ≈ t T when t T > τ Sα . Thebest linear fit (Fig. S1) to this region has a coefficient m = 1 . ∼
1, indicating that the relaxation time is justthe T1 delay time: τ α = t T . This suggests that whenthe “molecular-scale” T1 delay timescale is larger thanthe cell-scale collective timescale τ Sα , it dominates theresponse and solidifies the tissue. Moreover, the data col-lapse suggests that the cellular rearrangement timescale τ Sα is a good proxy for the mechanical response of thetissue over a wide range of model parameters. Figure 1. Vertex model with T1 delay time
T1 delay time A l l A) Example cell configuration in an isotropic vertex model at a finite temperature T = 0 . 02 and fixed system size N = 256.B) Schematic of a cellular rearrangement. An edge length of l shrinks to a length less than a critical length l c , forming afour-fold vertex. The edge is prevented from undergoing the T1 transition for a t T delay time as described in the main text.C) The characteristic relaxation time in the absence of T1 delays, defined by the self-overlap function, for various p values.The tissue becomes more viscous as p decreases at fixed temperature. D) Self-overlap function for T1 rearrangement delaytime of t T = 0 , . , . , . , . , . , . , . τ (darker green to yellow) for p = 3 . 74. The dotted lines indicatewhere Q s ( τ Sα ) = 1 /e in the absence of a T1 delay. E) Log-log plot showing collapse of the characteristic relaxation time τ Sα as a function of T1 delay time normalized by the collective response timescale τ Sα without a T1 delay. Colors correspondto different values of p = 3 . , . , . ... . T = 0 . 02, and N = 256. The inset shows thecharacteristic relaxation time τ Sα as a function of T1 delay time without any normalization, for the same values of p . Convergent extension rate is disrupted by T1rearrangement time. As many biological processes re-quire large-scale, anisotropic changes in tissue shape, wenext focus on the role of T1 delays in models that areanisotropic.We first study the dynamics of a vertex model withanisotropic line tensions, described by Eq. 3 and 4,where we initialize the simulations on a square do-main. Fig. 2(A) shows snapshots of the evolution of ananisotropic tissue in our simulations.For the anisotropic tissue, we first note that the stan-dard overlap function Q s defined in Eq. 5 is not a goodmetric for the rheology of a tissue with global shapechanges or tissue flow. This is because cells may stop overlapping their initial positions due to the macroscopicflow instead of due to local neighbor exchanges that areimportant for rheology. Therefore, we use a differentneighbors-overlap function Q n [35] to capture the rheol-ogy, which represents the fraction of cells that have losttwo or more neighbors in time t (see Supporting Infor-mation (SI) Section 1 for further details). In isotropicsystems, Q n and Q s are very similar, but Q n is muchbetter at identifying rheological changes in anisotropicsystems. Formally, it is defined as Q n ( t ) = 1 N N (cid:88) i =1 w (6)where w = 0 if a cell has lost two or more neighbors and w = 1 otherwise. Then the characteristic relaxation timeof the system measured by the nearest neighbors overlapfunction, τ Nα , is the time when Q n ( τ Nα ) = 1 /e .As in the isotropic case, we run simulations across arange of t T , p , and temperature T , with 100 indepen- dent simulations at each point in parameter space, andwhere all simulations are thermalized at their target tem-perature for 10 τ . We then run simulations for additionalsimulation time until the box reaches to a height of aboutfour cells to avoid numerical instabilities due to periodicboundary conditions. ax i s o f h i gh t e n s i on ABD C E Figure 2. Anisotropic vertex model with T1 delay - - t T1 / (cid:1) (cid:1) N E l ong a t i on R a t e ( (cid:1) ) A s p ec t R a t i o - - t T1 / (cid:1) (cid:1) t AR F (cid:1) (cid:1) (cid:1) N t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N FIG. 2. Anisotropic vertex model with T1 delay time A) Simulations of an anisotropic tissue. An anisotropic linetension on vertical edges is introduced to obtain global anisotropic changes to tissue shape. B) The collective response timescale for various γ values, the anisotropic line tension amplitude. C) Data collapse for p = 3 . . , . ... . T = 0 . N = 256 and γ = 1 . 0. The characteristic relaxation time, τ Nα as a function of T1 rearrangement delaytime normalized by the collective response timescale τ Nα ( t T = 0). The dotted line is a slope of 1. D) The aspect ratio of thesimulation box over time for T1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and 1000 τ (dark green to yellow), p = 3 . T = 0 . N = 256 and γ = 1 . 0. E) The time ( t AR ) at which the system first goes above the plateau value as afunction of t T for each aspect ratio curve in (D). F) The rate of elongation obtained from the aspect ratio curves in (D) as afunction of t T delay time. (D), (E) and (F) are from 10 independent simulation runs and the rate values are average ± onestandard error. First, we find that the characteristic relaxation τ Nα iscontrolled not only by p but also by the magnitude ofthe applied anisotropic line tension in Eq. 4, γ . Thisdata is shown for fixed p = 3 . 74 in Fig. 2(B).Figure 2(C) illustrates the characteristic relaxationtime, τ Nα , for the same values of p and T as shown inFig. 1(E), but with a fixed anisotropic line tension am-plitude of γ = 1 . 0. Both axes are normalized by thecollective response timescale τ Nα which corresponds to the case where the T1 rearrangement is instantaneous, t T = 0. The relation between the molecular-scale T1delay timescale and collective response timescale is sim-ilar to that of the isotropic tissue, where the molecular-scale T1 delay timescale dominates over the cell basedcollective response timescale when the t T delay time islarger than the collective response timescale. Support-ing Fig. S2 shows that this result is independent of themagnitude of the line tension in the anisotropic model.Next we analyze the rate of convergent extension for afixed p = 3 . 74 and γ = 1 . t T . We see that for t T < ∼ τ Nα ,there is a smooth elongation process until the simulationends, consistent with a fluid-like response (or like thebehavior of a yield-stress solid above the yield stress.)In this regime, increasing t T increases the rate of elon-gation slightly. In contrast, for t T > ∼ τ Nα , the systemfirst plateaus at a specific aspect ratio (which is aboutthree for the parameter values shown here), and only be-gins to elongate beyond that value for timescales greaterthan t T . The plateau value occurs in the absence of anyrearrangements, so it is entirely due to changes in indi-vidual cell aspect ratios. It is therefore governed by abalance between γ and k P P , and can be predicted ana-lytically, as described in the SI Section 3. These featuresare not strongly dependent on system size, as shown in(Fig. S3(D)).The change in behavior at t T ∼ τ Nα is highlighted inthe inset of Fig. 2(E), where we plot the time ( t AR ) atwhich the system first goes above the plateau value asa function of t T . Similar features can be seen in a plotof the elongation rate as a function of T1 delay time,shown in Fig. 2(F). We calculated the rate of elongation(Fig. 2(F)) as the growth constant of an exponential fit tothe aspect ratio over time (Fig. 2(D)) for each T1 delaytime value.To see if these observations are specific to systemswhere the anisotropy is generated by internal line ten-sions, or instead a generic feature of anisotropic sys-tems, we study a vertex model in the presence of anexternally applied pure shear strain (Fig. 3(A)). Fig. 3shows the collective response of the tissue to pure shear.Fig. 3(B) illustrates the inherent relaxation timescale ex-tracted from Q n as a function of p , and the fact that isvery similar to that in Fig. 1(C) suggests both that Q n and Q s are providing similar information and that thetissue rheology is robust across different perturbations(fluctuations vs. shear). Fig. 3(C) is similar to bothFig. 2(C) and Fig. 1(E), confirming that our observationthat T1 delays do not affect the tissue mechanics until t T ∼ τ α , and beyond that value the tissue me-chanics is dominated by that timescale, is robust acrossall perturbations studied. T1 rearrangement time delays and tissueanisotropies contribute to rosette formation. While T1 transitions are the simplest type of rearrange-ments in epithelial monolayers, it is common for verticesthat connect more than 4 cell edges to appear during de-velopmental processes. These are termed “rosettes” andthey appear often in tissue morphogenesis [4, 5] and col-lective cell migration [36].Although our model only allows three-fold coordinatedvertices, previous work by some of us has shown that A Figure 3. Anisotropic vertex model with T1 time delay (Pure Shear) B t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N t T1 ( (cid:1) ) (cid:1) (cid:1) N I n c r eas i ng ex t e r n a l pu r e s h ea r s t r a i n p (cid:1) (cid:1) N C FIG. 3. Anisotropic vertex model with external pureshear and T1 delay time A) We apply an external pureshear on the simulation box which is initially a squareddomain of L x = L y = L . We shear the box such that L x = e (cid:15) L and L y = e − (cid:15) L . Snapshots are from simulationswith p = 3 . T = 0 . N = 256 and t T = 1000 τ . B) Thecollective response time scale τ Nα ( t T = 0) for various p val-ues. C) Data collapse for p = 3 . . , . ... . T = 0 . 02 and N = 256. The characteristic relax-ation time, τ Nα as a function of T1 rearrangement delay timenormalized by the collective response timescale τ Nα ( t T = 0).The dotted line is a slope of 1. Inset shows τ Nα the charac-teristic relaxation time as a function of T1 delay time beforenormalization for p0 = 3.74, 3.76, 3.78...3.9 (darker to lightblue). vertices connected by short edges can be considered as aproxy for higher-order coordinated vertices [29]. In thatwork, a cutoff of 0 . √ A was used to threshold veryshort edges as a proxy for multi-fold coordination. Some-thing similar is also explicitly the case in experiments,where due to microscope resolution it is not possible todistinguish between very short edges and multi-fold co-ordinated vertices [21]. In that work, a cutoff of 0 . √ A was imposed by microscope resolution, and also adoptedin analysis of vertex models. Moreover, many-fold ver-tices are shown to be stable at heterotypic interfaces [37].In our simulations of tissues with anisotropic line ten-sions, either the additional anisotropic tension on inter-faces, or the edges that are prevented from undergoingT1 transitions, or both, could generate an increase inthe number of observed rosettes. Therefore, to study therole of T1 time delays on higher-order vertex formation,we analyze the number of very short edges per cell overtime in our simulations. For figures in the main text weadopt the larger cutoff of 0 . √ A used in [21], but inthe supporting information (Fig. S4) we show that theresults remain qualitatively similar for the smaller cutoffof 0 . √ A used in [29], although of course the over-all number of very short edges is smaller with the lowerthreshold. Figure 4. Many-fold vertex formation (note: per cell, p0=3.74,\gamma0=1,anisotropic) tT1=21.54 t=0.1,1,100 tau Snapshots of configurations at (A) t =0 . 1, (B) t = 1, and (C) t = 100 τ for t T = 77 . τ for a tissuewith p = 3 . T = 0 . N = 256 and γ = 1 . 0. D) Numberof very short edges per cell ξ for an anisotropic tissue as afunction of time. Shaded lines represent different T1 delaytimes t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and1000 τ (dark green to yellow), for a tissue with τ Nα = 0 . τ – p = 3 . T = 0 . N = 256 and γ = 1 . 0. E) Ensemble-averaged maximum value of ξ over a simulation timecourse( ξ max ) vs. the T1 delay time t T normalized by τ Nα . Theaverage is taken over 10 independent simulations, and errorbars correspond to one standard error. Figure 4 panels (A-C) highlights snapshots of typ-ical cellular structures at different timepoints in ananisotropic simulation with intermediate T1 delay time.Initially (panel A), cells are isotropic, and after about onenatural time unit (panel B) vertically oriented edges un-der anisotropic tension have shrunk to near zero length,resulting in a significant number of 4-fold coordinatedvertices and elongated rectangular shapes with an as-pect ratio set by a balance of anisotropic ( γ ) andisotropic( κ P P ) tensions, as discussed in SI Section 3.At timescales larger than the T1 delay time (panel C),T1 transitions allow some short edges to resolve and relax the structure.Figure 4(D) shows the number of short edges (SE)per cell ( ξ = 2 ∗ N avgSE /N cell , where the factor of tworeflects that edges are shared by two cells) over timefor T1 delay of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5,77.4, 278.2 and 1000 τ for an anisotropic tissue gener-ated using anisotropic internal tensions. For all valuesof the T1 delay time, there is an initial sharp rise to amaximum value, followed by a decrease. This decreasesuggests that, although there is a significant populationof higher-fold vertices that remain unresolved when theT1 delay is large, some many-fold vertices resolve onlonger timescales. This aligns well with experimental ob-servations of germ-band extension in Drosophila , wherecell junctions with dorsal-ventral orientation collapse toform higher order rosette structures and the rosettes areresolved by the extension of new junctions in anterior-posterior orientation [38].In addition, Fig. 4(E) shows that as the T1 delay timeincreases, the maximum number of short edges increases.Again, we see that there is a change in behavior around t T ∼ τ Nα . For t T < ∼ τ Nα , it is a monotonically increas-ing function, while for t T > ∼ τ Nα , it plateaus at the samelarge value of about 1 edge per cell. This is consistentwith our previous discussion of the mechanisms drivingelongation, where we noted that for large t T there wereno cellular rearrangements until t = t T , and so prior tothat timepoint the cells individually deform until theyform a nearly rectangular lattice. A perfect rectangularlattice would have ξ max = 2, so that two edges of thehexagon have shrunk to zero length, whereas in our dis-ordered systems we find approximately one short edgeper cell. Nevertheless, the maximum in the number ofshort edges is associated with these maximally deformedcell shapes.We also find that in isotropic tissues, increasing the T1delay times increases the number of many-fold vertices(Fig. S4(B,C)), although the numbers per cell are muchsmaller than that in the anisotropic tissue, as expected.All together, our results suggest that anisotropic linetension can collapse cell-cell junctions, resulting highernumber of many-fold structures. The number of suchstructures increases as the T1 rearrangement time delayincreases. Therefore, cellular rearrangement time couldbe a mechanism to regulate multicellular rosette struc-tures during morphogenesis. Number of T1 rearrangements. To study the im-pact of T1 delays on number of T1 rearrangements, wecalculate the number of successful T1 transitions per cell( η = N avgT irr ) /N cell ) over time for a T1 delay of t T =0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and 1000 τ (Fig. 5). In particular, in our model if an edge length l isless than the critical length l c and the associated T1 de-lay timer reaches to zero, the edge orientation is flipped.However, the same edge can flip back and forth to itsoriginal orientation easily until its final steady state con-dition is obtained. Hence, we define successful T1 tran-sitions as the arrangements that the cells rearrange andstay in their new configurations. In other words, the ar-rangements are irreversible (see SI Section 2 for details). Figure 5. Number of T1 transitions - - t T1 / (cid:1) (cid:1) N (cid:1) m ax BA ( (cid:1) ) (cid:1) I n c r eas i ng t T FIG. 5. Number of successful T1 transitions A) Num-ber of successful (irreversible) T1 transitions per cell η as afunction of simulation time t for an anisotropic tissue withT1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4,278.2 and 1000 τ (dark green to yellow), with τ Nα = 0 . τ – p = 3 . T = 0 . γ = 1 . N = 256, averagedover 10 independent realizations. B) Number of successful(irreversible) T1 transitions at the maximum averaged over10 realizations. Error bars represent one standard error. The data for the the number of successful (irreversible)T1 transitions share some similarities with the analysisof short edges in 4. Specifically Fig. 5(A) shows that thenumber of successful T1s grows towards a maximum andthen decays, which is consistent with a picture that cellsfirst become deformed and then execute T1 transitionsat longer timescales to facilitate large-scale tissue defor-mation. Fig. 5(B) shows that this maximum decreasesrapidly with increasing T1 delays, highlighting that thetissue response is much more elastic-like in the limit oflarge T1 delays. Again, we see a crossover in behavioraround t T ∼ τ Nα . Similarly, we study the effect of T1 de-lays on the number of T1 transitions for an isotropic setof simulations. The number of successful T1 transitionsdecreases as T1 delays increase (Fig. S5(E)), althoughthe effect is much weaker in isotropic tissues as there areno large-scale deformations driving T1 transitions in thatcase. Viscoelastic behavior of cellular junctions. Al-though the previous sections focus on global tissue re-sponse, we wanted to briefly investigate the impact of T1delays on localized cellular junction dynamics. We applya contractile tension on a cellular junction in a simula-tion with anisotropic internal tension. Specifically, westart from a configuration in the final steady state for agiven γ and p , and we choose p = 4 . t T = 278 . τ ; in this case the tissue ex-hibits viscoelastic features and recoils and recovers backto %70 − %90 of its initial length in high (Fig. 6(A) or-ange curve) and low (Fig. 6(B) orange curve) stress casesrespectively. Even though the cell shapes are similar inboth cases, the local viscous response of the cellular junc-tions change depending on the T1 delays in the tissue.Again, even at this smaller scale, systems with a largerT1 delay time are more elastic while those with a smallerT1 delay are more viscous. C C’B Figure 6. Viscous response of the junctions. Tissue becomes viscoelastic as tT1 increases. tT1=0 and p0=4 (green line) and tT1=280\tau and p0=4 (orange). A ** l ( t ) / l *** * l ( t ) / l FIG. 6. Viscous response of the cellular junctions (A,B) Junction length, l , normalized by the initial length of thejunction l over time in units of simulation time steps, duringand after a stress application on an edge in an anisotropicsimulation. Grey regions indicate the time period of an highapplied stress which shrinks the edge by %50 (A) and a lowapplied stress which shrinks the edge by %20 (B). The bluecurves are from the simulations without a T1 delay, t T = 0and the orange curves correspond to the simulations with a T1delay time of t T = 278 . τ . C) Snapshot from the simulationsillustrating the edge dynamics before, during and after thehigh stress application (asterisks in (A) indicate the exacttime points of the snapshots in (C) and (C’)) with t T = 0 (C)or with t T = 278 . τ (C’) T1 delay time. Other parametersare γ = 1 . p = 4 . T = 0 . N = 256. DISCUSSION AND CONCLUSIONS While standard vertex models for confluent tissues as-sume that T1 transitions proceed immediately after theconfiguration attains a multi-fold vertex, it is clear thatsome molecular processes may act as a brake on suchtransitions, generating a delay in the time required toresolve a higher-order vertex. In this work, we demon-strate that such T1 delays affect the tissue mechani-cal response in similar ways in isotropic, anisotropicallysheared, and internally anisotropic tissues. Specifically,we demonstrate that the relaxation timescale associatedwith neighbor exchanges in the absence of T1 delays, τ α ,is an excellent metric for glassy tissue response in thesedisparate systems. Moreover, in systems with T1 delays,the observed relaxation timescale τ α is related to τ α in aremarkably simple manner. For t T < ∼ τ α , the T1 delaysdo not strongly affect the system and τ α ∼ τ α , whilefor t T > ∼ τ α the T1 delay dominates the macroscopicdynamics and τ α = t T . In a related observation, we findthat the number of successful T1 transitions, where cellsneighbor exchange occurs and does not reverse at a latertime, decreases significantly for t T > ∼ τ α .This suggests that in tissues where molecular mech-anisms generate large T1 delays, the standard vertexmodel picture – where tissue fluidity is correlated withcell shapes, adhesion, and cortical tension – breaks down.While such molecular mechanisms are not able to speedup rearrangements, they can generically slow them down,provided that the T1 delays are larger that the inherentrelaxation timescale of the tissue. As an example, in pro-cesses such as convergent extension of the body axis in Drosophila , the aspect ratio changes by a factor of twoin about 30 minutes [8]. Given that cells in a hexago-nal vertex model would normally change neighbors aftera change of about two in the aspect ratio [39, 40], ourresults suggest that in Drosophila germband extension,molecular processes that require on the order of tens ofminutes or more to complete would be effective at inter-fering with global extension rates.Interestingly, the behavior of T1 transitions over timefor large T1 delays (Fig. 5(A) exhibits similar featuresto those observed in the germband of Drosophila snailtwist and bnt mutant embryos [21, 41]. In such embryosthe rearrangement rate does decrease significantly or dis-appear altogether, despite the fact that their cell shapeswould suggest a high rearrangement rate in a standardvertex model [21]. This is consistent with the hypothe-sis that molecular mechanisms in these mutants act as abrake on T1 transitions across all of developmental time.It is interesting to speculate that even in wild typeembryos, such molecular brakes could be deployed atdifferent stages of development to “freeze in” structuressculpted previously while the tissue was a fluid-like phase.For example, after the initial rapid elongation of the bodyaxis in fruit fly described in the previous paragraph, thecellular rearrangement rates decrease fairly precipitouslyabout 20 minutes after the elongation process initiates,even though the cell shapes are elongated and becomeeven more so [21].An additional observation is that in anisotropic sys-tems, increased T1 delay times are associated with in-creased persistence of higher-fold coordinated vertices,which we track by identifying very short edges in ourcomputer model. Specifically, for t T > ∼ τ α we find thatthe number of very short edges per cell increases dra-matically and remains high throughout the simulation.Again, this is consistent with the observation of a signif-icant number of rosettes in the later stages of Drosophila body axes elongation [4].In concurrent work, Das et al. [30] have studied a sim-ilar mechanism with an embargo on cell neighbor ex-change time, for a constant target shape index param-eter and in isotropic tissues. They discover interest-ing streaming glassy states where cells migrate in in-termittent coherent streams, similar to what is seen inspheroid/ECM experiments. Our work is complemen-tary, as we study both isotropic and anisotropic tissuesover a range of shape index parameters. This allowsus to emphasize the importance of the competition be-tween the collective response timescale driven by cell-scale properties and the T1 delay timescale driven bymolecular scale proerties at vertices. In addition, our fo-cus on global anisotropic changes to tissue shape allowsthis work to serve as a starting point for understandinghow T1 delays impact developmental processes such asthe convergent extension and rosette formation duringbody axis elongation.Taken together, our work suggests that moving for-ward it is really important to design experiments thatinvestigate which types of molecular processes are act-ing as brakes on T1 transitions. Obvious candidates areplayers in the cooperative disassembly and reassemblyof complex adhesive cell-cell junctions such as adherensjunctions [42] and/or desmosomes [43], or dynamics ofmolecules such as tricellulin that localize to three-foldcoordinated vertices [12].In recent work, Finegan et al. [14] study sdk mu-tants that lack the adhesion molecule Sidekick(Sdk)which localizes at tricellular vertices. They show that Drosophila sdk mutants exhibit a 1-minute delay in cellrearrangement timescales compared to wt embryos dur-ing Drosophila axis extension, accompanied by moreelongated cell shapes during the extension. In addition,they develop a vertex model that explicitly allows theformation of rosettes (higher fold-coordinated vertices)and additionally specifies that rosettes take longer to re-solve than simple 4-fold coordinated vertices. The modelrecapitulates cell shapes and global tissue deformationsseen in experiments. The spirit of the vertex model inthat work is very similar to the one we report here, ex-cept that we do not require any special rules for rosettes.In our model they form naturally in systems where T1delays occur, and they take longer to resolve simply be-cause the individual vertices that comprise them eachtake longer to resolve. Therefore, it would be interest-ing to see if sdk mutants are quantitatively consistentwith the model presented here, and whether one couldestimate the T1 delay timescale by fitting to the model,and then look for molecular processes at the vertices thatoccur on that same timescale that might be driving thedelay.In related recent work, Yu and Zallen [44] study Canoe,a different tricellular junctional protein. They find thatrecruitment of Canoe to tricellular junctions is correlated0with myosin localization during Drosophila convergentextension, and cells are arrested at four-fold vertex con-figurations in embryos that express vertex-trapped Ca-noe. The arrested cell rearrangements are 4 min longercompared to the rearrangements in wt embryos. In com-bination with our work, this suggests Canoe might alsoregulate T1 delays, and that tissues with perturbed Ca-noe dynamics or expression might another good systemfor testing our predictions about the relationship betweenT1 delays and global tissue mechanical properties.Lastly, one particularly intriguing avenue suggested byrecent work [16] is whether mechano-sensitive moleculesmay generate a stress-dependence for the T1 delaytimescale. In other words, our work here focuses on theeffects of a fixed T1 delay timescale that is independentof any local mechanical features of the cells. However, itis possible to introduce a feedback loop where T1 delaysare longer or shorter depending on the magnitude of thestresses on nearby edges, mimicking the behavior of well-known catch or slip bonds except now at a cell- or tissue-scale. Such feedback loops could lead to interesting pat-terning and dynamical behavior. ACKNOWLEDGEMENTS We thank Daniel Sussman, Matthias Merkel, Pe-ter Morse, Karen Kasza, and Xun Wang for fruit-ful discussions. This work was supported by NIHR01GM117598 and the Simons Foundation (grant ∗ [email protected] † [email protected][1] K. E. Kasza and J. A. Zallen, Dynamics and regulationof contractile actin-myosin networks in morphogenesis,Curr. Opin. Cell Biol. , 30 (2011).[2] M. Merkel, R. Etournay, M. Popovic, G. Salbreux,S. Eaton, and F. Julicher, Triangles bridge the scales:Quantifying cellular contributions to tissue deformation,Phys. Rev. E , 32401 (2017).[3] D. L. Weaire and S. Hutzler, The Physics of Foams (Ox-ford University Press, 2000).[4] J. T. Blankenship, S. T. Backovic, J. Sanny, O. Weitz,and J. A. Zallen, Multicellular rosette formation linksplanar cell polarity to tissue morphogenesis, Developmen-tal Cell , 459 (2006).[5] M. Harding, H. F. McGraw, and N. A., The roles andregulation of multicellular rosette structures during mor-phogenesis, Development , 2549 (2014).[6] Z. J. A., Planar polarity and tissue morphogenesis, Cell , 1051 (2007).[7] L. V. Goodrich and D. Strutt, Principles of planar po-larity in animal development, Development , 1877(2011). [8] K. E. Kasza, D. L. Farrell, and J. A. Zallen, Spatiotempo-ral control of epithelial remodeling by regulated myosinphosphorylation, PNAS , 11732 (2014).[9] R. J. Tetley, G. B. Blanchard, A. G. Fletcher, R. J.Adams, and B. Sanson, Unipolar distributions of junc-tional myosin ii identify cell stripe boundaries that drivecell intercalation throughout Drosophila axis extension,eLife , e12094 (2016).[10] M. Rauzi, P. F. Lenne, and T. Lecuit, Planar polarizedactomyosin contractile flows control epithelial junctionremodelling, Nature , 1110–1114 (2010).[11] H. Katsuno-Kambe and A. S. Yap, Endocytosis, cad-herins and tissue dynamics, Traffic , 268 (2020).[12] J. Ikenouchi, M. Furuse, K. Furuse, H. Sasaki, S. Tsukita,and S. Tsukita, Tricellulin constitutes a novel barrier attricellular contacts of epithelial cells , Journal of Cell Bi-ology , 939 (2005).[13] Y. Oda, T. Otani, J. Ikenouchi, and M. Furuse, Tricel-lulin regulates junctional tension of epithelial cells at tri-cellular contacts through cdc42, Journal of Cell Science , 4201 (2014).[14] T. M. Finegan, N. Hervieux, A. Nestor-Bergmann, A. G.Fletcher, G. B. Blanchard, and B. Sanson, The tricellularvertex-specific adhesion molecule sidekick facilitates po-larised cell intercalation during drosophila axis extension,PLOS Biology , 1 (2019).[15] H. Uechi and E. Kuranaga, The tricellular junction pro-tein sidekick regulates vertex dynamics to promote bi-cellular junction extension, Developmental Cell , 327(2019).[16] M. F. Staddon, K. E. Cavanaugh, E. M. Munro, M. L.Gardel, and S. Banerjee, Mechanosensitive junction re-modeling promotes robust epithelial morphogenesis, Bio-physical Journal , 1739 (2019).[17] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, Adensity-independent glass transition in biological tissues,Nature Physics , 1074–1079 (2015).[18] M. Chiang and D. Marenduzzo, Glass transitions in thecellular potts model, EPL , 28009 (2016).[19] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T.Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim,B. Gweon, J. Notbohm, R. Steward Jr, S. Burger, S. H.Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore,E. Israel, D. A. Weitz, D. J. Tschumperlin, E. Henske,S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen,and J. J. Fredberg, Unjamming and cell shape in theasthmatic airway epithelium, Nature Materials , 1040(2015).[20] J. Devany, D. M. Sussman, M. L. Manning, and M. L.Gardel, Cell cycle-dependent active stress drives epitheliaremodeling, bioRxiv 10.1101/804294 (2020).[21] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan,M. L. Manning, and K. E. Kasza, Anisotropy links cellshapes to tissue flow during convergent extension, PNAS , 13541 (2020).[22] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Motility-Driven Glass and Jamming Transitions in Bio-logical Tissues, Physical Review X , 021011 (2016).[23] D. M. Sussman, M. Paoluzzi, M. Cristina Marchetti, andM. Lisa Manning, Anomalous glassy dynamics in simplemodels of dense biological tissue, EPL , 1 (2018).[24] W. T. McCleery, J. Veldhuis, M. E. Bennett, H. E.Lynch, X. Ma, G. W. Brodland, and M. S. Hutson, Elon-gated Cells Drive Morphogenesis in a Surface-Wrapped Finite-Element Model of Germband Retraction, Biophys-ical Journal , 157 (2019).[25] F. Schock and N. Perrimon, Cellular processes associatedwith germ band retraction in drosophila, DevelopmentalBiology , 29 (2002).[26] D. P. Kiehart, J. M. Crawford, A. Aristotelous, S. Ve-nakides, and G. S. Edwards, Cell sheet morphogenesis:Dorsal closure in drosophila melanogaster as a model sys-tem, Annu. Rev. Cell Dev. Biol. , 169 (2017).[27] M. E. Lacy and M. S. Hutson, Amnioserosa developmentand function in drosophila embryogenesis: Critical me-chanical roles for an extraembryonic tissue., Developmen-tal Dynamics , 558–568 (2016).[28] M. F. Staddon, K. E. Cavanaugh, E. M. Munro, M. L.Gardel, and S. Banerjee, Solid–fluid transition and cellsorting in epithelia with junctional tension fluctuations,Soft Matter , 3209 (2020).[29] T. Yamamoto, D. M. Sussman, T. Shibata, and M. L.Manning, Non-monotonic fluidization generated by fluc-tuating edge tensions in confluent tissues, arXiv (2020).[30] A. Das, S. Sastry, and D. Bi, Controlled neighbor ex-changes drive glassy behavior, intermittency and cellstreaming in epithelial tissues, arXiv (2020).[31] R. Farhadifar, J. C. R¨oper, B. Aigouy, S. Eaton, andF. J¨ulicher, The Influence of Cell Mechanics, Cell-Cell In-teractions, and Proliferation on Epithelial Packing, Cur-rent Biology , 2095 (2007).[32] D. M. Sussman, cellGPU: Massively parallel simulationsof dynamic vertex models, Computer Physics Communi-cations , 400 (2017).[33] M. Rauzi, P. Verant, T. Lecuit, and P. Lenne, Natureand anisotropy of cortical forces orienting drosophila tis-sue morphogenesis, Nature Cell Biology , 1401–1410(2008).[34] T. Castellani and A. Cavagna, Spin-glass theory forpedestrians, J. Stat. Mech. , P05012 (2005). [35] F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita,M. L. Manning, R. Cerbino, and M. C. Marchetti, Flock-ing transitions in confluent tissues, Soft Matter , 3471(2018).[36] G. Trichas, A. M. Smith, N. White, V. Wilkins,T. Watanabe, A. Moore, B. Joyce, J. Sugnaseelan, T. A.Rodriguez, D. Kay, R. E. Baker, P. K. Maini, and S. Srini-vas, Multi-cellular rosettes in the mouse visceral endo-derm facilitate the ordered migration of anterior visceralendoderm cells., PLOS Biology , e1001256 (2012).[37] D. M. Sussman, J. M. Schwarz, M. C. Marchetti, andM. L. Manning, Soft yet Sharp Interfaces in a VertexModel of Confluent Tissue, Physical Review Letters ,1 (2018).[38] D. Kong, F. Wolf, and J. Großhans, Forces directinggerm-band extension in Drosophila embryos., Mecha-nisms of Development , 11 (2017).[39] S. A. Khan, Foam rheology: Relation between exten-sional and shear deformations in high gas fraction foams,Rheologica Acta , 78 (1987).[40] D. Weaire and R. Phelan, The physics of foam, Journalof Physics: Condensed Matter , 9519 (1996).[41] D. L. Farrell, O. Weitz, M. O. Magnasco, and J. A. Za-llen, Segga: a toolset for rapid automated analysis ofepithelial cell polarity and dynamics, Development ,1725 (2017).[42] C. D’Souza-Schorey, Disassembling adherens junctions:breaking up is hard to do, Trends in Cell Biology , 19(2005).[43] Y. Kitajima, Mechanisms of desmosome assembly anddisassembly, Clinical and Experimental Dermatology ,684 (2002).[44] H. H. Yu and J. A. Zallen, Abl and canoe/afadin medi-ate mechanotransduction at tricellular junctions, Science , 1 (2020). SUPPORTING INFORMATION1. The neighbors-overlap function We define a neighbors-overlap function Q n , as described in the text, which represents the fraction of cells that havelost two or more neighbors in time t . We decided upon the cutoff value of two or more neighbors based on the followingobservation. We compare the results obtained from the neighbors-overlap function with the results obtained from theself-overlap function. The two should give similar results for an isotropic tissue. The characteristic relaxation timeobtained using the standard self-overlap function is similar to the one obtained from the neighbors-overlap function(Fig. S6(A)). On the contrary, a definition based on a cutoff value of losing three or more neighbors show a differencebetween the results obtained from the self-overlap function and that of the neighbors-overlap function (Fig. S6(B)). 2. Irreversible T1 transitions In our model, cells go through a T1 transition whenever the edge l between four neighboring cells becomes lessthan a critical length l c , and the T1 delay count is reached. The same group of four cells can go back and forthbetween their original configuration and their after T1 configuration until the final steady state condition is obtained(Fig. S5(A). Hence, these flipping events can cause overcounting of the number of true T1 transitions. To avoid this,we calculate the number of successful T1 events where cells rearrange and stay in their new configurations. In otherwords, we calculate the number of irreversible T1 events.To calculate the irreversible events ( N irr ), we first characterize the reversible events ( N rev ) in our simulations.When four cells go through a T1 transition, the two cells share an edge lose a vertex and the two cells that arenot neighbor before the T1 transition gain a vertex (Fig. S5(A)). We record the cell configurations that go throughthe T1 transition and scan through the rest of the simulation steps to check whether the cell configurations areever back to their “before T1” configurations in the next τ + t T steps. Here, τ is the natural time unit of thesimulations. If the configuration is repeated, we count the event as reversible, and we call the time between the twosame configurations reversibility time, t R . As we are interested in instantaneous successful T1 events over time, wesearch for the reversibility in short timescales, namely in τ steps. Fig. S5(B) is a histogram of reversibility time for t T = 0 which fits to a power law (Fig. S5(B) inset) with a long tail which indicates that most of the flipping T1events happen in very short timescales.We calculate the fraction of irreversible events, f = N irr /N total . For an anisotropic tissue, 70% of the T1 eventsare irreversible for t T = 0 while almost all of the events are irreversible for t T > t T = 0 (Fig. S5(D)). The irreversibility increases as the t T delay time increases (Fig. S5(D)). To calculate the successful T1 events (main text Fig. 5(A) and (B)), we multiplythe total number of T1 transitions by the fraction of irreversible events f (Fig. S5(C)) at each t T value. Similarly, foran isotropic tissue (Fig. S5(E)), we multiply the total number of T1 transitions by the fraction of irreversible events f (Fig. S5(D)) at each t T value. 3. Plateau value for the tissue aspect ratios in the absence of T1 transitions As discussed in the main text and seen in Fig. 2(D), in the absence of T1 transitions ( i.e. , at simulation timepointsbefore the T1 delay timescale) tissues under anisotropic tension will increase their aspect ratio up to a plateau valueand get stuck there. Numerical simulations confirm that this plateau value varies with both the isotropic tension κ P P in the bare vertex model and the additional anisotropic line tension γ (Fig. S7).To predict this plateau value analytically, we hypothesize that the plateau occurs when the shapes are at an energyminimum under the constraint there are no T1 transitions. We propose that the cells undergo a two-step process tominimize their energy, which is easiest to see starting from an ordered isotropic hexagonal packing, shown in SupportingFig.S8(A). First, the large line tension on the vertical edges will cause those to shrink to zero length, generating adiamond pattern with 4-fold coordinated vertices, as shown in Fig.S8(B). This process should be independent of γ and P for values of γ that are sufficiently large, and one can show it results in a change of aspect ratio by a factorof 1.58: AR = 1 . l and horizontal length l , andassume the area of each cell is fixed to unity so that l = 1 /l . Taking the derivative of the part of the vertex energyfunctional related to line tensions, E ∼ γ l + κ P ( P − P ) , and then setting the derivative with respect to l equalto zero results in a 4th order polynomial equation: l + (cid:18) γ κ P − P (cid:19) l + P l − . (7)For a given values of κ P , P , and γ , we can solve this equation for its positive roots and identify the energy minimizingvalue of l min . The addition change to the global aspect ratio of the tissue allowed by these rectangles is simply ratioof l = 1 /l to l : AR = 1 / ( l min ) .Finally, the total change to the aspect ratio under both processes is simply AR tot = AR × AR = 1 . / ( l min ) .This analytic prediction is illustrated by the blue squares in Fig.S8(D). The observed plateau values are given bythe red circles, and they are in fairly good agreement. The analytic prediction overestimates the aspect ratio for thelowest value of γ , likely because the tension isn’t sufficiently large to shrink vertical edges to zero in the first process.4 4. Supporting figures Fig S1 t T1 / (cid:1) (cid:1) S (cid:1) (cid:1) S / (cid:1) (cid:1) S FIG. S1. The characteristic relaxation time τ Sα on linear scale (log-log scale as in Fig. 1) as a function of T1 delay timenormalized by the collective response timescale τ Sα without T1 delay. The dashed line is the best linear fit to high T1 delayregion with a slope of m = 1 . 13. Colors correspond to different values of p = 3 . , . , . ... . T = 0 . 02, and N = 256. Figure S2. Anisotropic vertex model with T1 rearrangement time (note: a,a’ (pure shear), b) different \gamma (compared to Fig.2 ), c) different T C B t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N t T1 ( (cid:1) ) (cid:1) (cid:1) N t T1 ( (cid:1) ) (cid:1) (cid:1) N t T1 ( (cid:1) ) (cid:1) (cid:1) N t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N D t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N A t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N FIG. S2. A) The characteristic relaxation time τ Nα as a function of T1 rearrangement delay time for anisotropic tissuesimulations with anisotropic line tension amplitude of γ = 0 . 01 (A), γ = 0 . γ = 0 . p = 3 . , . , . ... . 9, and otherparameters are T = 0 . 02, and N = 256. Each data set are normalized by the corresponding characteristic relaxation time τ Nα where the T1 rearrangement is instantaneous, t T = 0. Figure S3. Anisotropic vertex model with T1 rearrangement time. N=1024 AB t T1 / (cid:1) (cid:1) N (cid:1) (cid:1) N / (cid:1) (cid:1) N ax i s o f h i gh t e n s i on - - t T1 / (cid:1) (cid:1) N t AR ( (cid:1) ) A s p ec t R a t i o - - t T1 / (cid:1) (cid:1) N E l ong a t i on R a t e CD E F ( (cid:1) ) (cid:1) - - t T1 / (cid:1) (cid:1) N (cid:1) m ax FIG. S3. A) Simulations of an anisotropic tissue with a larger size. An anisotropic line tension on vertical edges isintroduced to obtain global anisotropic changes to tissue shape. B) Data collapse for p = 3 . . , . ... . T = 0 . N = 1024 and γ = 1 . 0. The characteristic relaxation time, τ Nα as a function of T1 rearrangement delay timenormalized by the collective response timescale τ α ( t T = 0). The dotted line is a slope of 1. C) The number of T1 transitionsper cell over time and at the maximum averaged over 10 realizations (inset) for T1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99,21.5, 77.4, 278.2 and 1000 τ (dark green to yellow), p = 3 . T = 0 . N = 1024 and γ = 1 . 0. D) The aspect ratio of thesimulation box over time for T1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and 1000 τ (dark green to yellow), p = 3 . T = 0 . N = 1024 and γ = 1 . 0. E) The time ( t AR ) at which the system first goes above the plateau value as afunction of t T for each aspect ratio curve in (D). F) The rate of elongation obtained from the aspect ratio curves in (D) as afunction of t T delay time. Both (C), (D), (E) and (F) are from 10 realizations. Error bars represent one standard error. Top panel cutoff=0.11Bottom panel cutoff=0.04Figure S4. Many-fold vertex formation (note: per cell, p0=3.74) isotropic) B’ ( (cid:1) ) (cid:1) - - - - t T1 / (cid:1) (cid:1) N (cid:1) m ax - - - - t T1 / (cid:1) (cid:1) N (cid:1) m ax B ( (cid:1) ) (cid:1) C C’ - - t T1 / (cid:1) (cid:1) N (cid:1) m ax A’ ( (cid:1) ) (cid:1) A FIG. S4. Number of very short edges per cell ξ for an anisotropic tissue as a function of time (A). Shaded lines representdifferent T1 delay times t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and 1000 τ (dark green to yellow), for a tissue with τ Nα = 0 . τ – p = 3 . T = 0 . N = 256 and γ = 1 . 0. The cutoff value to threshold very short edges as a proxy formulti-fold coordination is 0 . √ A . (A’) Ensemble-averaged maximum value of ξ over a simulation timecourse vs. the T1 delaytime t T normalized by τ Nα . The average is taken over 10 independent simulations, and error bars correspond to one standarderror. Number of very short edges per cell ξ over time (B, C) and the ensemble-averaged maximum value of ξ over a simulationtimecourse (B’, C’) for an isotropic tissue simulations of T1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4, 278.2 and1000 τ (dark green to yellow), p = 3 . T = 0 . 02 and N = 256. The cutoff value to threshold very short edges as a proxy formulti-fold coordination is 0 . √ A in (B, B’) and 0 . √ A in (C, C’). The average is taken over 10 independent simulations,and error bars correspond to one standard error. A 21 3456 2 671 3 45 t R ( (cid:1) ) R eve r s i b l eeve n t s C BD Figure S5 - - - t T1 / (cid:1) (cid:1) (cid:1) e q - - - t T1 / (cid:1) (cid:1) N f - t T1 / (cid:1) (cid:1) N f anisotropic E FIG. S5. A) Schematic of a reversible T1 transition. The two cells that are neighbor before the T1 transition lose a vertex(top and bottom cells) and the cells that are not neighbors before the T1 transition gain a vertex and become neighbors (rightand left cells). Vertices of the cell on the left are labeled as 1,2,...6. The same cell has a configuration of vertices as 1,2,...7 afterthe T1 transition. If the T1 transition is reversible, then the cell goes back to its original configuration of vertices 1,2,...6. Wetrack such cell configuration changes in our simulations to determine if a T1 event is reversed. B) Histogram of reversibilitytime t R which fits to a power law (inset black points) for an anisotropic tissue of p = 3 . T = 0 . N = 256, t T = 0 and γ = 1 . 0. C) Fraction of irreversible events as a function of T1 delay time for an anisotropic tissue of p = 3 . T = 0 . N = 256 and γ = 1 . 0. D) Fraction of irreversible events as a function of T1 delay time for an isotropic tissue of p = 3 . T = 0 . 02 and N = 256. E) Number of successful (irreversible) T1 transitions at the equilibrium averaged over 10 realizationsfor an isotropic tissue of p = 3 . T = 0 . 02 and N = 256. Error bars represent one standard error. Figure S6. New overlap function comparison A t T1 / (cid:1) (cid:1) (cid:1) (cid:1) / (cid:1) (cid:1) t T1 / (cid:1) (cid:1) (cid:1) (cid:1) / (cid:1) (cid:1) B FIG. S6. The characteristic relaxation time, τ α , obtained using neighbors-overlap function (based on losing 2 or more neighbors(A) and losing 3 or more neighbors (B)), as a function of T1 rearrangement time for p = 3 . , . , . ... . T = 0 . 02 for an isotropic tissue. Both (A) and (B) are plotted together with the characteristic relaxation timeobtained using the self-overlap function for comparison for p = 3 . , . , . ... . τ α where the T1 rearrangement is instantaneous, t T = 0. ( (cid:1) ) A s p ec t R a t i o ( (cid:1) ) A s p ec t R a t i o Figure S7 A B FIG. S7. The aspect ratio of the simulation box over time for T1 delay time of t T = 0, 0.13, 0.46, 1.67, 5.99, 21.5, 77.4,278.2 and 1000 τ (dark green to yellow), p = 3 . T = 0 . 02 and N = 256. The anisotropic line tension amplitude γ = 0 . γ = 0 . ∼ . ∼ . FIG. S8. A) A hexagonal unit box with a hexagon side of a . B) The cells meet at a 4-fold coordinated vertex after ahigh tension is applied on the vertical edges of the hexagons in (A). Cells form a diamond pattern with a side length of a .Both in (A) and (B), L x and L y are the sides of the box along the horizontal and vertical direction. C) A rectangular gridof cells with sides l along the vertical and l along the horizontal direction. D) Analytical prediction of the plateau valuesobserved in simulations with various γ0