Elasticity, stability, and quasi-oscillations of cell-cell junctions in solid confluent epithelia
EElasticity, stability, and quasi-oscillations of cell-cell junctionsin solid confluent epithelia
Cl´ement Zankoc and Matej Krajnc ∗ Joˇzef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
Macroscopic properties and shapes of biological tissues depend on the remodelling of cell-celljunctions at the microscopic scale. We propose a theoretical framework that couples a vertexmodel of solid confluent tissues with the dynamics describing generation of local force dipoles inthe junctional actomyosin. Depending on the myosin-turnover rate, junctions either preserve stablelength or collapse to initiate cell rearrangements. We find that noise can amplify and sustaintransient oscillations to the fixed point, giving rise to quasi-periodic junctional dynamics. We alsodiscover that junctional stability is affected by cell arrangements and junctional rest tensions, whichmay explain junctional collapse during convergence and extension in embryos.
INTRODUCTION
In confluent tissues, the adjacent cells adhere to oneanother through narrow joints known as the adherensjunctions. These regions are rich with protein com-plexes, which govern cell-cell adhesion and couple cells’cytoskeletons [1, 2]. Forces transmitted along the junc-tions play a major role in morphogenesis. In particu-lar, contractile tensions generated in the actomyosin [3–5] drive various types of cell deformations and tissue-scale movements [6, 7]. For instance, the actomyosin net-work collapses during cell ingression and oscillates duringdorsal closure in
Drosophila [8–12], directed junctionalcollapse drives convergence and extension in the early
Drosophila embryo [13–16], whereas junctional fluctua-tions establish the arrangement of cells in
Drosophila pu-pal notum [17] and fluidize the tissue during vertebratebody axis elongation [18].Recent measurements of junctional dynamics duringtissue remodelling revealed that the rate of junctionalcollapse often increases with contraction [16, 17]. Thissuggests a positive feedback loop between junctional con-traction and generation of junctional tension. Apart fromdriving collapse of the actomyosin, such a feedback loopcould establish oscillatory dynamics as previously pro-posed within a generic theoretical framework of activecontractile elements [19].To explore these dynamics in the context of confluenttissues, we develop a vertex model with a feedback loopbetween cell-scale junctional contractions and generationof force dipoles at the level of the junctional actomyosin.We find that the nonlinear elastic response of solid tis-sues to local force dipoles does not satisfy conditions toyield a limit cycle of junctional oscillations. Neverthe-less, the variety of dynamical regimes remains rich andincludes junctions that can either sustain stable length orcollapse. One of our main findings shows that junctionalnoise establishes quasi-periodic junctional dynamics. Wealso discover that junctional stability depends on cell ar-rangements as well as on distribution of junctional resttensions, which could be relevant for the active junctional remodelling during morphogenesis [13, 14].
MATERIALS AND METHODSThe model
The tissue is represented by a planar polygonal net-work of cell-cell junctions parametrized by the positionsof cell vertices r i = ( x i , y i ) [20–22]. Forces on vertices F i are assumed conservative such that F i = −∇ i W , where W = κ P k ( A k − A ) + P ij γ ij l ij is the total potentialenergy of the system and ∇ i = ( ∂/∂x i , ∂/∂y i ). The firstsum in the energy goes over all cells and describes cell-area elasticity ( A k and A being the actual and the pre-ferred cell areas, respectively), whereas the second sumgoes over all junctions and describes the line energy dueto adhesion and cortical tension ( l ij being the junctionallength); indices i and j denote head and tail verticesof the junctions. For simplicity, we disregard the usualcell-perimeter-squared energy term [23, 24], which is notneeded to describe solid tissues studied here. Due tostrong friction, the motion of vertices is overdamped: η ˙ r i = F i = − κ X k ( A k − A ) ∇ i A k − X j γ ij ∇ i l ij , (1)where η is the friction coefficient. The tension γ ij ( t ) = γ + ∆ γ ij ( t ), where γ and ∆ γ ij ( t ) are the average (rest)tension and the time-depending part, respectively. To en-sure stability, γ is assumed positive, which implies thatthe cortical tension dominates over the adhesion [17, 25].From now on we use dimensionless quantities by choosing √ A , τ = ( ηγ ) − , and γ as the units of length, time,and tension, respectively, and rescaling the modulus as κA / /γ → κ .Tissues with constant line tensions ( γ ij = const . )behave as passive elastic materials capable of sustain-ing shear stresses [23, 26]. To examine their responseto a locally applied force dipole, we assume an excessline tension ∆ γ on a single junction. This deforms thejunction away from its rest length l , and the elastic a r X i v : . [ c ond - m a t . s o f t ] S e p ED l eng t h time t m odu l u s κ turnover rate 1/ τ m A f o r c e f ( ) length -1 0.90.6-302-21 1.2 C τ m onoff Δ γ f ( ) actinmyosin junctioncontractiontension increasereducerelax < . τ m TOS local model S C B linear termhysteresis ( ) ≠ . static ( ) = . CTOS > . FIG. 1. (A) Static force-extension relation (red curve) andhysteresis (grey curve) for a single junction in the honey-comb cell tiling. (B) Schematic of the interplay between theelastic restoring force f ( l ) and active tension ∆ γ generatedby the junctional myosin. (C) The feedback loop betweenjunctional contraction and tension. (D) Junction length vs.time in the local model. (E) Phase diagram of local modelexhibits stable (S), collapse (C), and stable dynamics withtransient oscillations to the fixed point (TO). Red, blue, andblack curves (panel D) and circles (panels E) correspond to(1 /τ m , κ ) = (1 , . , , deformation propagates into the bulk (Supporting Ma-terials and Methods, Fig. S1A). Force balance implies0 = − γ − f ( l ), where f ( l ) is the elastic restoring force,which is well described by f ( l ) = P m =1 α m ( l − l ) m with α /α = 0 .
33 and α /α = − .
10 for a regular hon-eycomb cell tiling (Fig. 1A and Supporting Materials andMethods, Fig. S1B,C). The significant contribution of thesecond-order term yields softening and stiffening behav-iors for compression and stretch, respectively (Fig. 1A).Furthermore, due to friction, f ( l ) displays hysteresis atnon-zero rates of change of junctional length (Fig. 1A).Tensions change on time scales associated with the dy-namics of the underlying actomyosin (Fig. 1B). We as-sume a linear relation between junctional tension andmyosin concentration [defined by a number of molecularmotors per junction length c ij ( t ) = N ij /l ij ]: c ij ( t ) = αγ ij ( t ), where α is a constant proportionality factor.The total rate of change of the myosin concentration˙ c ij = ˙ N ij /l ij − N ij ˙ l ij /l ij , where the first and the secondterm describe changes of the number of motors at a fixedjunction length and changes of the junction length at afixed number of motors, respectively. In particular, themotor-actin binding and undbinding contribute a timescale τ m and are described by ˙ N ij = ( − N ij + c l ij ) /τ m ,where c = αγ is the ambient myosin concentration.In turn, ˙ l ij can be explicitly written in terms of the forces acting at the vertices as ˙ l ij = r ij · F ij /l ij , where r ij = r j − r i and F ij = F j − F i . Overall, this yields adeterministic dynamic equation for tension fluctuations:˙∆ γ ij = − τ m ∆ γ ij − γ ij r ij · F ij l ij , (2)where the first term describes tension relaxation due tomyosin turnover, whereas the second term describes acoupling between tension dynamics and mechanics of thevertex model. This coupling gives rise to a feedbackloop between generation of tension and junctional con-tractions (Fig. 1C). In contrast to some previous studies,which describe similar feedback mechanisms using spe-cific chemomechanical models [11, 15, 27], here the feed-back follows directly from the relation between junctionaltension and myosin concentration. RESULTSLocal model
First, we apply the model on a quartet of hexagonalcells (inset to Fig. 1D). We derive the equation of statefor the central junction, which reads f ( l ) = − κl (cid:18) ll − (cid:19) − − l/l − p ( l/l ) − l/l + 7 , (3)where l = 3 − / √ ≈ .
62 is the side length of a regularhexagon with unit area. Interestingly, despite its simplic-ity, the local model exhibits similar elasticity than thefull tissue–including softening/stiffening behavior uponjunction compression/stretch (Supporting Materials andMethods, Sec. I and Fig. S2A,B). Next, we simulate thedynamics described by Eqs. (1) and (2), starting in aconfiguration away from the fixed point ( l, γ ) = ( l , l ( t ), depends on the initial conditions.A linear stability analysis (Supporting Materials andMethods, Sec. II) reveals the analytical condition for theHopf bifurcation: κ ∗ = 13 l − l τ m . (4)In particular, the junction is stable for κ > κ ∗ and un-stable for κ < κ ∗ (Fig. 1E). Transient oscillations ap-pear within the stable regime at κ − < κ < κ + , where κ ± = 1 / (3 l ) + 1 / (3 l τ m ) ± (1 / l ) p / ( l τ m ) (Fig. 1E).Importantly, the first Lyapunov coefficient shows that theHopf bifurcation is always subcritical, meaning that noparameter values yield a stable limit cycle that would de-scribe periodic oscillations like those previously observedin similar dynamical models [11, 12, 19] (Supporting Ma-terials and Methods, Sec. III).Next, to search for the various types of junctional be-haviors in the full-tissue setting, we apply our model tothe honeycomb cell tiling and numerically preform thelinear stability analysis (Supporting Materials and Meth-ods, Sec. IV). We identify the same dynamical regimesas in the local model (Fig. 2A). However, in contrast tothe local model, the collapse regime extends to large κ values (Fig. 2A). This is because collective cell deforma-tions in the full-tissue model allow preservation of cellareas even upon junction collapse. FIG. 2. (A) Phase diagram of the full tissue exhibits sta-ble (S), collapse (C), and transient oscillations to the fixedpoint (TO). Eigenvalues of the Jacobian were evaluated at(1 /τ m , κ )-pairs on an equidistant grid of 100 ×
100 points.(B,C) Junctional length (B) and tension (C) vs. time for ex-amples with quasi-oscillations (grey and orange curves) andstochastic fluctuations (blue). (D) Amplitude of length fluctu-ations vs. myosin turnover rate 1 /τ m for σ = 0 . − . /τ m , κ ) = (3 . , , , . σ = 0 . Quasi-oscillations
Binding and undbinding of myosin are stochastic pro-cesses, which cause stochastic fluctuations of junctionaltensions [17, 25, 28, 29]. To explore the role of noise, weadd an extra term to the equation for tension, which now reads ˙ γ ij = ∆ ˙ γ ij + p σ /τ m ξ ij . Here the first term obeysEq. (2) like before, whereas the second term describes thewhite noise with long-time variance σ and h ξ ij ( t ) i = 0, h ξ ij ( t ) ξ kl ( t ) i = δ ik δ jl δ ( t − t ). Note that in the absenceof the coupling term in Eq. (2), our stochastic tension dy-namics reduces to the classical Ornstein-Uhlenbeck (OU)process, which was previously used to describe stochasticjunctional fluctuations [17, 25].We simulate the stochastic dynamics in a honeycombcell tiling and find that while in the stable (S) regime,length (and tension) fluctuations are noisy (blue curvesin Fig. 2B and C), the movements become more regu-lar in the regime of transient oscillations (TO) and, incontrast to the deterministic case, manage to sustain awell-defined amplitude (orange and grey curves in Fig. 2Band C). These movements in fact correspond to the tran-sient oscillations, which eventually die out in the purelydeterministic case (Fig. 1D), but get amplified (Fig. 2D)and become sustained indefinitely in the presence of thenoise. This gives rise to a quasi-periodic trajectory orthe so-called quasi-cycle, which arises in dynamical sys-tems in the presence of the noise when the Jacobian ma-trix has complex eigenvalues with strictly negative realparts [30, 32 ? ]. We refer to these junctional movementsas quasi-oscillations.Next, we examine the observed dynamics in the Fourierspace. The quasi-oscillations give rise to peaks in powerspectral densities (PSD) of tension- and length fluctu-ations, PSD γ and PSD l , respectively. These peaks getdominated by the noise when moving away from the bi-furcation point (Fig. 2E and Supporting Materials andMethods, Sec. V). In fact, in the limit 1 /τ m → ∞ wherethe relative contribution of the coupling term in the ten-sion dynamics [Eq. (2)] becomes negligible compared tothe relaxation term, the PSDs agree with the model inwhich tensions obey a pure OU process (blue curves inFig. 2E).In confluent tissues, junctions are interconnected andso need to synchronize their quasi-oscillations. Sincethree junctions meet at each vertex, junction networksare geometrically frustrated, which can lead to nontrivialspatial patterns of cell deformations (Supporting Mate-rials and Methods, Sec. VI). Inhomogeneous tissues
The stability condition of our local model [Eq. (4)]suggests that disorder of cell packing may affect junc-tional stability. Indeed, it can be recast as l > l ∗ , wherethe l ∗ , is the critical rest length for collapse. Since therest lengths are distributed in disordered tissues (inset toFig. 3A), there might exist a fraction of junctions thatare shorter than l ∗ . These junctions would necessarilycollapse and trigger cell rearrangements. To test thispossibility, we examine an ensemble of 100 disorderedtissues (Supporting Materials and Methods, Sec. VII).Starting close to the fixed point (i.e., all junctions attheir rest lengths and all tensions equal 1), we simulatethe dynamics [Eqs. (1) and (2)] at fixed κ = 15 andrecord the rest lengths of the first 10 collapsing junctions.Here, each junctional collapse initiates a T1 transition,which is performed as soon as the junction length dropsbelow 0 .
01. The length of the newly created junctionis set to 0.001, whereas its tension is reset to γ . Un-like in the honeycomb lattice where junctions at κ = 15collapse only for 1 /τ m < . /τ m values (Fig. 3A). This confirms that the distributionof rest lengths in disordered tissues importantly affectslocal junctional stability. We note in passing that fre-quent collapses of short junctions drive partial orderingof disordered tissues (Supporting Materials and Methods,Fig. S6 and Videos. S1 and S2). A rest length P D F c oun t s τ m rest tension γ h i gh - c o ll ap s e γ f r a c t i on o f high γ BC first 10 junctions first 10 junctions FIG. 3. (A) Probability distribution function (PDF) ofthe rest length of first 10 collapsing junctions for 1 /τ m =0 . , , , , . γ ≥ γ collapsed junctions versus γ . Finally, the dimensional form of the stability condi-tion, κ ∗ > γ / (3 l ) − η/ (3 l τ m ), suggests that the restjunctional tension γ may also affect junctional stabil-ity. In particular, junctions with higher γ are expectedto have an increased critical rest length l ∗ , possibly re-sulting in their collapse. We test this prediction by ana-lyzing a tissue that contains a supracellular cable of en-riched junctional myosin, which increases the rest ten-sions on the corresponding junctions compared to otherjunctions (Fig. 3B). Without the contraction-tension cou-pling [Eq. (2)], this cable would be stable despite beingunder higher tension, since every vertex within the cableis acted upon by a pair of equal but opposite forces. How-ever, as predicted by the stability condition, the couplingindeed affects the stability of the junctions that are underhigher rest tension, making them more succesptible forcollapse. To show this, we record the first 10 collapsingjunctions and find that the fraction of those with high γ increases with γ (Fig. 3C). DISCUSSION
We studied a mechanical model of tissues with a feed-back loop between junctional contractions and the dy-namics of junctional tensions (Fig. 1). In particular, weused a previously proposed description of force genera-tion at the level of the actomyosin [19] and combined itwith the vertex model of solid confluent tissues, whichprovided a faithful representation of the elasticity under-lying the response of solid tissues to local force dipolesat the junctions. While nonlinearities in this system donot meet the conditions to yield a stable limit cycle ofjunctional oscillations (Fig. 1A), we discovered that junc-tional noise can amplify and sustain quasi-periodic junc-tional dynamics at biologically relevant myosin-turnoverrates [Fig. 2 and Ref. [17]]. Importantly, this dynami-cal regime does not even require nonlinear elasticity andmay thus be more common than the limit cycle of pe-riodic junctional oscillations. Another important resultof our work highlights the role of cell arrangements andthe distribution of rest tensions within the tissue forthe junctional stability (Fig. 3). Both effects may bepresent during convergence and extension in
Drosophila embryo, where the so-called parasegmental boundaries,enriched with the junctional myosin, frequently collapsetheir junctions [13, 14, 33].An interesting future direction would be to employ theArea- and Perimeter-Elasticity vertex model [23, 24] anduse it to explore tissues that are closer to the solid-fluidtransition. These tissues are associated with highly non-linear elasticity [24, 26, 34], which could, in contrast toour model, yield a stable limit cycle of junctional oscilla-tions [19]. In turn, this could lead to spontaneous orga-nization of junctions into groups of locally synchronizedoscillators. Furthermore, the role of correlations betweenjunctions needs to be further investigated. In particular,these correlations can appear because individual junc-tions are under tension exerted by the actomyosins from two adjacent cells and because cell membranes belong-ing to the same cell share a common pool of molecularmotors. These effects would provide additional sourcesof coupling between junctions, which could significantlyaffect their correlated movements. Finally, an alternativemodel could also assume that junctional tension is pro-portional to the number of myosin motors rather than totheir concentration.
ACKNOWLEDGMENTS
We thank Primoˇz Ziherl, Jan Rozman, Tomer Stern,Guillaume Salbreux, and Fabio Staniscia for critical read-ing of the manuscript. We acknowledge the financialsupport from the Slovenian Research Agency (researchproject No. Z1-1851 and research core funding No. P1-0055). ∗ [email protected][1] T. J. C. Harris and U. Tepass, Nat. Rev. Mol. Cell Biol. , 502 (2010).[2] C. G. Vasquez and A. C. Martin, Dev. Dyn. , 361(2016).[3] T. Lecuit and P.-F. Lenne, Nat. Rev. Mol. Cell Biol. ,633 (2007).[4] T. Lecuit, P.-F. Lenne, and E. Munro, Annu. Rev. CellDev. Biol. , 157 (2011).[5] M. Murrell, P. W. Oakes, M. Lenz, and M. L. Gardel,Nat. Rev. Mol. Cell Biol. , 486 (2015).[6] R. Fernandez-Gonzalez, S. de Matos Simoes, J.-C. R¨oper,S. Eaton, and J. A. Zallen, Dev. Cell , 736 (2009).[7] C.-P. Heisenberg and Y. Bella¨ıche, Cell , 948 (2013).[8] S. Simes, Y. Oh, M. F. Wang, R. Fernandez-Gonzalez,and U. Tepass, J. Cell Biol. , 1387 (2017).[9] A. Jayasinghe, S. Crews, D. Mashburn, and M. Hutson,Biophys. J. , 255 (2013).[10] A. C. Martin, M. Kaschube, and E. F. Wieschaus, Nature , 495 (2009).[11] S.-Z. Lin, B. Li, G. Lan, and X.-Q. Feng, Proc. Natl.Acad. Sci. USA , 8157 (2017).[12] W.-C. Lo, C. Madrak, D. P. Kiehart, and G. S. Edwards,Phys. Rev. E , 062414 (2018).[13] C. Bertet, L. Sulak, and T. Lecuit, Nature , 667(2004).[14] M. Rauzi, P.-F. Lenne, and T. Lecuit, Nature , 1110(2010).[15] M. F. Staddon, K. E. Cavanaugh, E. M. Munro,M. L. Gardel, and S. Banerjee, Biophys. J. , 1739(2019).[16] T. Stern, S. Shvartsman, and E. Wieschaus, PLOSComp. Biol. , 1 (2020).[17] S. Curran, C. Strandkvist, J. Bathmann, M. de Gennes,A. Kabla, G. Salbreux, and B. Baum, Dev. Cell , 480 (2017).[18] A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shel-ton, D. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio,J. Giammona, and O. Camp`as, Nature , 401 (2018).[19] K. Dierkes, A. Sumi, J. Solon, and G. Salbreux, Phys.Rev. Lett , 148102 (2014).[20] A. Fletcher, M. Osterfield, R. Baker, and S. Shvartsman,Biophys. J. 106, 2291 (2014).[21] S. Alt, P. Ganguly, and G. Salbreux, Philos. Trans. RoyalSoc. B , 20150520 (2017).[22] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek,PLoS Comp. Biol. , e1005569 (2017).[23] R. Farhadifar, J.-C. R¨oper, B. Aigouy, S. Eaton, andF. J¨ulicher, Curr. Biol. , 2095 (2007).[24] D. Bi, J. Lopez, J. Schwarz, and M. L. Manning, Nat.Phys. , 1074 (2015).[25] M. Krajnc, Soft Matter , 3209 (2020).[26] D. B. Staple, R. Farhadifar, J.-C. R¨oper, B. Aigouy,S. Eaton, and F. J¨ulicher, Eur. Phys. J. E , 117 (2010).[27] L. C. Siang, R. Fernandez-Gonzalez, and J. J. Feng,Phys. Biol. , 066008 (2018).[28] M. Krajnc, S. Dasgupta, P. Ziherl, and J. Prost, Phys.Rev. E , 022409 (2018).[29] S. Okuda, E. Kuranaga, and K. Sato, Biophys. J. ,1159 (2019).[30] R. P. Boland, T. Galla, and A. J. McKane, J. Stat. Mech.:Theory Exp. (2008), P09001.[31] C. A. Lugo and A. J. McKane, Phys. Rev. E 78, 051911(2008).[32] C. Zankoc, D. Fanelli, F. Ginelli, and R. Livi, Phys. Rev.E 96, 022308 (2017).[33] M. Rauzi, Philos. Trans. R. Soc. , 20190552 (2020).[34] P. Sahu, J. Kang, G. Erdemci-Tandogan, and M. L. Man-ning, Soft Matter , 1850 (2020). lasticity, stability and quasi-oscillations of cell-cell junctionsin solid confluent epithelia: Supplemental Material Cl´ement Zankoc and Matej Krajnc
FIG. S1. (A) Deformation field due to a local force dipole (redarrows) in a honeycomb cell tiling. (B) Coefficients α i vs. restlength l from force-extension relations f ( l ) = P m =1 α m ( l − l ) m for junctions extracted from disordered tissues. (C) Thedistribution of second- and third-order non-linear coefficientsin disordered tissues. The green circle in panel B denotesthe pair of nonlinear coefficients corresponding to the stateequation of an edge in the honeycomb lattice of cells. I. EQUATION OF STATE
We treat a quartet of hexagonal cells in which thecentral pair of vertices with positions r = (0 , y ) and r = (0 , y ) move symmetrically with respect to the x -axis (i.e., y = − y ), whereas all the other verticesare fixed (Fig. S2A). While such vertex displacementspreserve the total area of this four-cell neighborhood,the areas of individual cells A A , A B , A C , and A D dochange with vertex displacements. Nevertheless, the sys-tem possesses two axes of mirror symmetry, meaning that A A = A B , A C = A D , whereas the length of all four edgesadjacent to the central edge, λ ( l ) = vuut √ l ! + (cid:18) l − l (cid:19) . (S1)We refer to these edges as peripheral edges and to pre-serve the symmetry, we assume that they are under the y x AB CD y y γγ p γ p γ p γ p λ ( l ) l fixed vertexmobile vertex A B f o r c e f ( ) / f ( ) length -0.5 0.90.6-1.500.51.5 local model ( κ = κ = -1.51.0 1.2 C -8 0.1 1 10 1000.01-60 e i gen v a l ue s Λ turnover rate 1/ τ m realimaginary collapse oscill. stable L y apuno v c oe ff . τ m D FIG. S2. (A) Schematic of the local model. Vertices denotedby green crosses are fixed, whereas the blue vertices are al-lowed to move along the y direction. Variables of the modelare the length of the central edge l , and tensions γ and γ p onthe central and 4 peripheral edges, respectively. The lengthsof peripheral edges λ ( l ) can be explicitly calculated from l .(B) Force-extension curves of the full vertex model (blackcurve) and the local model (purple curve). Dashed line showsthe linear part of the elastic restoring force. (C) Eigenvaluesof the Jacobian versus turnover rate 1 /τ m at fixed κ = 0 . /τ m . same tension γ p = γ p ( t ); tension on the central edge isdiffrerent from that at peripheral edges and is denotedby γ .The total energy of the system is given by W = κ X α ( A α − A ) + X ij γ ij l ij . (S2)We are interested in the dynamics of the central edgelength l = y − y . Therefore, we need to first derive thedynamics equations for y and y , which obey η d y / d t = − ∂W/∂y and η d y / d t = − ∂W/∂y , respectively. Inturn, this gives η d y d t = − κ X α ( A α − A ) ∂A α ∂y − γ ∂l∂y − γ p ∂λ ( l ) ∂y (S3) a r X i v : . [ c ond - m a t . s o f t ] S e p and η d y d t = − κ X α ( A α − A ) ∂A α ∂y − γ ∂l∂y − γ p ∂λ ( l ) ∂y . (S4)Next, we express cell areas with respect to l as follows: A A = A B = √ l (5 l + l ) (S5)and A C = A D = √ l (7 l − l ) . (S6)The derivatives of the areas with respect to y and y read ∂A A ∂y = ∂A B ∂y = − √ l , (S7) ∂A A ∂y = ∂A B ∂y = √ l , (S8) ∂A C ∂y = ∂A D ∂y = √ l , (S9)and ∂A C ∂y = ∂A D ∂y = − √ l , (S10)whereas the derivatives of the lengths are ∂l/∂y = − ∂l/∂y = 1, and ∂λ ( l ) ∂y = − ∂λ ( l ) ∂y = − l − l λ ( l ) . (S11)Together, this yields the differential equation for thelength l : η d l d t = − κl (cid:18) ll − (cid:19) − γ − γ p l − l λ ( l ) . (S12)Finally, setting γ = γ p = γ yields the equation ofstate f ( l ) = − κl (cid:18) ll − (cid:19) − γ − γ ( l − l ) p l − l l + 7 l , (S13)which describes the elastic restoring force acting on thejunction that is subject to a local force dipole.The equation of state of our local model agrees wellwith the one describing junctional elasticity in the fulltissue. This is somewhat surprising, since any junctiondeformation in the local model neccassarily changes theareas of all 4 cells involved (due to limited number of de-grees of freedom), whereas multiple-cells deformations inthe full tissue model allow preserving cell areas. Never-theless, provided that the local model allows significant area deformations ( κ = 1 in Fig. S2B), the equations ofstate are qualitatively the same in both systems.The complete dynamics of the local model is given bythe following set of differential equations: η d l d t = − κl (cid:18) ll − (cid:19) − γ − γ p l − l λ ( l ) , (S14)d γ d t = − τ m ( γ − γ ) − γl d l d t , (S15)d γ p d t = − τ m ( γ p − γ ) − γ p λ ( l ) d λ d l d l d t . (S16) II. LINEAR STABILITY ANALYSISA. Hopf bifurcation
We assume a small perturbation from the fixed pointby setting l = l + δl , γ = γ + δγ , and γ p = γ + δγ p ,where δl, δγ and δγ p (cid:28)
1. The linearized system canbe written in a matrix form as d ρ / d t = J ρ , where ρ =( δl, δγ, δγ p ) and the Jacobian matrix J = − ( κl + γ ) ηl − η η γ ( κl + γ ) ηl γ ηl − τ m − γ ηl − γ ( κl + γ ) ηl − γ ηl γ ηl − τ m . (S17)The Jacobian has three eigenvalues: The first one, Λ = − /τ m , is always negative, whereas the remaining tworead Λ ± = b a − ± r − acb ! , (S18)where a = 2 ηl τ m > b = 2 (cid:0) ηl + 3 κl τ m − γ τ m (cid:1) ,and c = 3 (cid:0) κl + γ (cid:1) >
0. The eigenvalues are plottedin Fig. S2C.The system is at the critical point when Λ + = 0; i.e.,Λ + < + > b = 0. This gives thecritical area-compressibility modulus, which reads κ ∗ = γ l − η l τ m . (S19) B. Transient oscillations
By considering the stable system, i.e., κ > κ ∗ , we de-rive the condition for the existence of transient oscilla-tions, which imply that the eiganvalues have a non-zeroimaginary part. This is true when 1 − ac/b < κ < γ l + η l τ m + 13 l r ηγ l τ m . (S20)The frequency of transient oscillations is given by ω =(1 / p ac/b − ω = s ηl τ m (2 κl + γ )( ηl + 3 κl τ m − γ τ m ) − . (S21) III. FIRST LYAPUNOV COEFFICIENT
To examine the stability of the limit cycle, we calculatethe first Lyapunov exponent l at the critical point, i.e.,for κ = γ / (3 l ) − η/ (3 l τ m ). First, we compute eigen-vectors and eigenvalues of the Jacobian J at the criticalpoint and we denote by q the eigenvector associated tothe eigenvalue iω where ω = p − ( b − ac ). We dothe same for the Jacobian transpose J T and denote by p its eigenvector associated with eigenvalue the − iω , suchthat J T p = − iω p . We then Taylor-expand our systemaround its fixed point to third order: d ρ dt = F ( ρ ) = J ρ + 12 b ( ρ , ρ ) + 16 c ( ρ , ρ , ρ ) , (S22)where the components of b and c read B j ( u, v ) = X k,l =1 ∂ F j ( ρ ) ∂ρ k ∂ρ l (cid:12)(cid:12)(cid:12)(cid:12) F P u k v l (S23)and C j ( u, v, w ) = X k,l,m =1 ∂ F j ( ρ ) ∂ρ k ∂ρ l ∂ρ m (cid:12)(cid:12)(cid:12)(cid:12) F P u k v l w m . (S24) The first Lyapunov coefficient l is then given by l = 12 ω < [ h p , c ( q , q , ¯ q ) i − h p , b ( q , J − b ( q , ¯ q )) i + h p , b ( ¯ q , (2 iω − J ) − b ( q , q )) i ] (S25)and it turns out strictly positive for 1 /τ m ∈ (0 , /l ]and therefore the system never displays a stable limitcycle (Fig. S2D). IV. LINEAR STABILITY ANALYSIS OF FULLTISSUE
To study the stability of the full tissue, we need towrite down the system of equations in its full form. Inparticular, given N v vertices and N e edges, we have 2 N v differential equations for positions of the vertices and N e differential equations for tensions. In particular, the forceexerted on vertex i reads F i = η d r i d t = − X j γ ∇ i l ij − K X k ( A k − A ) ∇ i A k , (S26)whereas the rate of change of tension at junction con-necting vertices i and j readsd γ ij d t = − τ m ( γ ij − γ ) − γ k ( F j − F i ) · ( r j − r i ) l ij η . (S27)We perturb the system around the fixed point, wherethe vertices assume their positions in a regular hexagonallattice, whereas tensions on all junctions equal 1. Welinearize the system, which then rereads δ F i = η d δ r i d t = X j δ Λ lij + X k δ Λ Aik , (S28)and d δγ ij d t = − τ m δγ ij − δ Γ ij . (S29)Here δ Λ lij = γ l ij ( y i − y j ) [( y i − y j ) ( − δx i + δx j ) + ( x i − x j )( δy i − δy j )]( x i − x j ) [( y i − y j ) ( δx i − δx j ) + ( x i − x j )( − δy i + δy j )] + − ( x i − x j ) δγ ij /l ij − ( y i − y j ) δγ ij /l ij , (S30) δ Λ Aik = K (cid:0) y ν k ( i )+1 − y ν k ( i ) − (cid:1) (cid:2) − P ν ∈ k ( y ν +1 − y ν − ) δx ν + P ν ∈ k ( x ν +1 − x ν − ) δy ν (cid:3)(cid:0) x ν k ( i )+1 − x ν k ( i ) − (cid:1) (cid:2)P ν ∈ k ( y ν +1 − y ν − ) δx ν − P ν ∈ k ( x ν +1 − x ν − ) δy ν (cid:3) , (S31)and δ Γ ij = γ ( δ F j − δ F i ) · ( r j − r i ) l ij η = γ l ij η (cid:2) δF xj ( x j − x i ) + δF yj ( y j − y i ) − δF xi ( x j − x i ) − δF yi ( y j − y i ) (cid:3) . (S32)Here r i = ( x i , y i ) and l ij = p ( x i − x j ) + ( y i − y j ) arepositions of vertices and edge lengths, respectively, atthe fixed point. Figure S3 shows the dependence of themaximal eigenvalue of the Jacobian as a function of 1 /τ m at fixed κ = 0 . m a x . e i gen v a l ue Λ τ m
10 100
Re Im κ = collapse oscill. stable FIG. S3. Maximal real and imaginary parts of the eigenvaluesof the Jacobian versus 1 /τ m at fixed κ = 0 . V. NOISE
We rewrite the system of differential equations, addingthe noise to the tension dynamics: η d l d t = − Kl (cid:18) ll − (cid:19) − γ − γ p l − l λ ( l ) , (S33)d γ d t = − τ m ( γ − γ ) − γl d l d t + s σ τ m ξ ( t ) , (S34)d γ p d t = − τ m ( γ p − γ ) − γ p λ ( l ) d λ d l d l d t + s σ τ m ξ p ( t ) , (S35)where ξ ( t ) is a Gaussian random variable with zero av-erage and variance equal to 1. Its correlator reads h η ( t ) η ( t ) i = h η p ( t ) η p ( t ) i = δ ( t − t ) (S36)and h η ( t ) η p ( t ) i = 0 . (S37) The noise increases the probability for the system to es-cape the basin of attraction and thus expands the collapseregime in the (1 /τ m , κ ) parameter space. Close to theboundary between the collapse and stable regimes, thesystem exhibits quasi-oscillations. These are instanceswhere the noise amplifies the transient oscillations andsustains them for an indefinite time period.To compute the power spectrum of the quasi-oscillations, we again linearize the system around thefixed point: l = l + δl , γ = γ + δγ , and γ p = γ + δγ p .The linearized system of stochastic differential equationscan then be written in a matrix form asdd t v = J v + η . (S38)Next, we write v in the Fourier space as ˜ v ( ω ) = R v ( t ) exp( iωt ) to obtain ( iω − J ) ˜ v ( ω ) = ˜ η ( ω ). De-noting Φ( ω ) = iω − J and solving for ˜ v ( ω ) gives˜ v ( ω ) = Φ − ( ω ) ˜ η ( ω ) . (S39)Now we can calculate the power spectral density ma-trix (PSDM) with the elements given by P ij ( ω ) = h ˜ v i ( ω )˜ v † j ( ω ) i = X l =1 3 X m =1 Φ − il ( ω ) D lm (cid:2) Φ − mj (cid:3) † ( ω ) , (S40)where D = 2 σ τ m . (S41)The diagonal entries of the PSDM are real and coin-cide with the power spectra for the fluctuations, asso-ciated with the three species l , γ , and γ p . In particular,the power spectrum for the junctional length fluctuationsreads P l ( ω ) = 16 σ τ m g ( ω ) (cid:18) τ m + ω (cid:19) , (S42)where g ( ω ) = (cid:12)(cid:12)(cid:12)(cid:12) − iω + (cid:18) γ l − κl − τ m (cid:19) ω + i (cid:18) τ m + γ l τ + 6 κl τ m (cid:19) ω + 3 γ l τ m + 3 κl τ m (cid:12)(cid:12)(cid:12)(cid:12) . (S43)Note that the spectrum of the Ornstein-Uhlenbeck pro-cess is given by P ( ω ) = 2 (cid:0) σ /τ m (cid:1) / (cid:0) ω + 1 /τ m (cid:1) . VI. CORRELATION OF QUASI-OSCILLATIONS
In tissues, cell-cell junctions are physically coupled andare therefore forced to synchronize their movements. For
FIG. S4. Amplitude of length fluctuations vs. myosinturnover rate 1 /τ m for σ = 0 . − . example, junctions can lengthen in expense of the short-ening of their immediate neighbors. In polygonal net-works, three or more junctions meet at a vertex, whichprevents this pairwise relation from being simultane-ously satisfied for all pairs of adjacent junctions. To ex-plore what patterns of junctional movements this geomet-ric frustration establishes, we simulate collective quasi-oscillations within a regular honeycomb cell tiling. Wemeasure the Pearson correlation coefficient, p , betweenjunctional lengths l ij ( t ) and the length of the referencejunction l ref ( t ). Not surprisingly, we find that overall, thecorrelation strength | p | decreases with the topological dis- d A a v e r age P ea r s on c o rr e l a t i on c oe ff . | p | B - p / | p | C -6 + P D +1+1 - - - - μ -1-0.500.51 p FIG. S5. (A) Pearson correlation coefficient | p | vs. topologicaldistance d averaged over all junctions with equal d . (B) Colormap of correlation between junctions and the reference (black)junction. Black arrows point to cell chains with increased cor-relation strength of their parallel junctions with the referencejunction. (C) Color map of the sign of correlation betweenjunctions and the reference (black) junction. (D) Color mapof cells’ cumulative correlation of their adjacent junctions P .Results in all panels are calculated at (1 /τ m , κ ) = (3 . , tance from the reference junction d (Fig. S5A); the topo-logical distance d is defined as the integer shortest pathbetween a junction and the reference junction; d = 1 forthe nearest neighbors, d = 2 for the next-nearest neigh-bors, etc. However, | p | ( d ) is non-monotonic, which sug-gests a nontrivial structure of the correlation map withinthe cell-junctions network. Indeed, the correlation pat-tern exhibits two clear features: (i) correlation is strongerbetween junctions that are mutually parallel (Fig. S5B)and (ii) chains of cells running perpendicularly to thereference junction alternate in the strength of correla-tion (black arrows in Fig. S5B).It is also interesting to observe how junctions synchro-nize their quasi-oscillations within individual cells. Inparticular, we focus on the sign of the correlation coef-ficient p/ | p | , which tells whether the pair of junctionsis correlated ( p/ | p | = +1) or anti-correlated ( p/ | p | = −
1) (Fig. S5C). We define cumulative correlation of ad-jacent junctions in a given cell by P = P µ ∈ cell ( p µ / | p µ | ) · ( p µ +1 / | p µ +1 | ), where the index µ = 1 , ...,
6. As a resultof geometric frustration at the vertices, not all cells canhave P = − P alternates between − f r a c t i on o f he x agon s t BA f r a c t i on o f he x agon s t
40 500.80.9 10 20 30
FIG. S6. Fraction of hexagons versus time for a disordered ini-tial configuration at 1 /τ m = 5 . κ = 10, and σ = 0 .
01 (panelA) and an ordered initial configuration at 1 /τ m = 1, κ = 10,and σ = 0 .
01 (panel B). Tissue simulations for both cases areshown in Videos S1 and S2.
VII. DISORDERED TISSUE SAMPLES
To prepare disordered tissue samples, with distributedjunction rest lengths, we use a pure Ornstein-Uhlenbeckscheme to describe tension dynamics:˙∆ γ ij = − τ m ∆ γ ij + s σ τ m ξ ij . (S44)In particular, we start with a regular honeycomb lattice,which we fluidize by running simulations at σ = 0 . τ m = 1 for t = 1000. Next, we quench the system instan-taneously by setting σ to 0 and simulate the model for tt