How does the extracellular matrix affect the rigidity of an embedded spheroid?
Amanda Parker, M. Cristina Marchetti, M. Lisa Manning, J. M. Schwarz
HHow does the extracellular matrix affect the rigidity of an embedded spheroid?
Amanda Parker , M. Cristina Marchetti , M. Lisa Manning , J. M. Schwarz , Department of Physics and BioInspired Syracuse, Syracuse University,Syracuse, NY 13244, USA Department of Physics,University of California at Santa Barbara, Santa Barbara, CA,93106, USA Indian Creek Farm, Ithaca, NY, 14850, USA (Dated: June 30, 2020)Cellularized tissue and polymer networks can both transition from floppy to rigid as a functionof their control parameters, and, yet, the two systems often mechanically interact, which may affecttheir respective rigidities. To study this interaction, we consider a vertex model with interfacialtension (a spheroid) embedded in a spring network in two dimensions. We identify two regimeswith different global spheroid shapes, governed by the pressure resulting from competition betweeninterfacial tension and tension in the network. In the first regime, the tissue remains compact,while in the second, a cavitation-like instability leads to the emergence of gaps at the tissue-networkinterface. Intriguingly, compression of the tissue promotes fluidization, while tension promotescellular alignment and rigidification, with the mechanisms driving rigidification differing on eitherside of the instability.
Introduction.
Cellularized tissue– groups of cells ex-hibiting collective behavior– is held together by adhe-sive cell-cell junctions, which regulate an incredible rangeof functions, including wound healing and embryogene-sis [1–5]. In addition, there is typically a non-cellularcomponent present, the extracellular matrix (ECM)– anetwork often largely composed of cross-linked collagenfibers– to which cells mechanically couple via cell-ECMadhesions. Through these adhesions, the ECM imposesforces on cells to provide structural support and regulatecell-cell interactions, driving changes in cell behavior and,ultimately, tissue mechanics [6–10]. Maintenance of theseinteractions is crucial for healthy tissue, and, in fact, dis-ruptions in forces exerted by the ECM, or in the cells’ability to sense or respond to mechanical signals fromthe ECM, are hallmarks of disease [11–15].A key in vitro model system for studying the com-bined effects of cell-cell and cell-ECM interactions ontissue behavior is a tissue spheroid embedded in a col-lagen network. It has been shown that tensile forcesand stiffer collagen fibers in the ECM facilitate tumorinvasion [16, 17], and that as various ECM proteins aremodified, cell shapes may be a good predictor of invasionpotential [18].Despite the recent spate of experimental work on mul-ticellular systems embedded in biopolymer networks, nu-merical models for cells interacting with ECM have fo-cused on single cells, potentially missing important fea-tures associated with collective cell behavior [19–21]. Yet,models for bulk tissue and collagen network have sepa-rately been quite successful. One such model is a cell-based vertex model [22–24], which predicts a density-independent rigidity transition in disordered confluenttissues and micro-demixing in tissue mixtures [24–26].Both phenomena have been verified experimentally [26–28]. Meanwhile, the strain-stiffening behavior observedin biopolymer networks has been captured by under- constrained spring network models [29, 30]. Interest-ingly, in the absence of energetic penalties for bending,the spring network’s floppy-to-rigid transition at finiteshear strain is yet another density-independent rigiditytransition in the same universality class as that of thevertex model [31].Here we couple a tissue-based vertex model to a sur-rounding spring network in two dimensions with line ten-sion at the interface, and study the rheology and mor-phology of the tissue. Our model is minimal as our goalis to investigate the effect of mechanical coupling aloneon the collective behavior of pre-migratory tissue cells.However, even for this minimal bi-material, without feed-back to dynamically up- or down-regulate model param-eters, nontrivial changes in cell shape, tissue phase, andoverall tissue structure can be induced. Using mean-fieldapproximations, we analytically quantify the state of thetissue in two distinct regimes, separated by an instabilityin the tissue boundary, to predict the fluidity and geom-etry of the tissue cells, using arguments that are gener-alizable to any biological system where the competitionbetween external forces and interfacial tension plays animportant role.
Model.
Our model is composed of a tissue of N cells cells,each with area, a α , and perimeter, p α , embedded in anetwork of springs, each with equilibrium spring length, l , and stiffness, k sp , with line tension, γ , at the inter-face (see Fig. 1 (a), detail). The dimensionless energyfunctional is e total = X cells,α ( a α − + k p X cells,α ( p α − p ) + 12 k sp X sp h ij i ( l sp h ij i − l ) + γ X int h ij i l int h ij i , (1)with the first two terms representing the standard vertexmodel for the tissue, the third term the elasticity of thespring network with edge lengths l sp h ij i , and the fourth a r X i v : . [ q - b i o . CB ] J un (d)(a) (e) ECM edges (springs)Interface edgesTissue cell edges (b)(c)
FIG. 1.
Tissue configurations, average cell shape, and fluid-ity, for fixed p = 3 . . (a)-(c) Minimal energy configurationsfor fixed equilibrium spring length, l = 0 .
56, and varying in-terfacial tension, γ = 1.0, 0.22, and 0.0, including a detailshowing springs in orange, interface edges in black, and celledges in green. (d) The average tissue cell area, h a i , versus γ for a series of l values equally-spaced from 0 .
45 (purple)to 0 . D eff , of tissue cells as a function of γ and l . Theline indicates the analytical prediction for the cavitation in-stability, γ MFc , while Xs denote the onset of the instability, γ c ,measured from simulations. The black symbols correspond tothe same three sets of parameters throughout. term the interfacial tension, with edge lengths l int h ij i .The control parameter for rigidity in the standard vertexmodel is the preferred cell shape index, p , with values of p below a threshold inducing rigidity in the cellularizedtissue (see Supplementary Material (SM) Sec. A). In ahomogeneous, under-constrained spring network system,the control parameter for rigidity (tension) is the restlength of the springs, l , with a rest length below somethreshold resulting in a rigid network (see SM Sec. B).For simplicity, we study the spring network in the limit ofzero bending energy where the rigid-to-floppy transitionis well-defined.Using zero-temperature energy-minimization and low-temperature Brownian dynamics (see SM Sec. D), weexplore minimal-energy configurations of this model, fo-cusing on the emergent morphology and rheology of thetissue. In doing so, the rules for tissue cell T1 transitionsnear the boundary must be carefully considered (see SMSec. C). Results.
We find that both the phase and morphology ofthe tissue depend on the balance of interfacial tension, γ , and the equilibrium spring length, l , of the surroundingnetwork, for a fixed p of the tissue cells. As l decreases,the mean spring length, h l sp i , decreases to minimize thetotal energy. This decrease necessitates the shrinking ofthe area taken up by the springs, either by increasing thesize of the tissue or by opening gaps at the tissue-ECMinterface. At the same time, changes to the global area ofthe tissue spheroid cost energy due to a quadratic penaltyfor deviations from the preferred cell area, while changesto the tissue-ECM interface are penalized by the inter-facial tension and induced deviations from the preferredcell perimeter.The competition between these effects is illustrated inFig. 1 for p = 3 .
91. The configurations in (a)-(c) de-pict the effect on the overall tissue shape of decreasing γ ( γ ( a ) = 1 . γ ( b ) = 0 . γ ( c ) = 0 .
0) for fixed l . Weobserve roughly circular shapes for the tissue in the firsttwo cases, while at γ = 0 cavities appear, yielding a dif-ferent global shape. Fig. 1 (d) shows the mean cell area, h a i , as a function of γ , for varying l values (see caption).Fig. 1 (e) shows the effects of γ and l on the tissue fluid-ity, quantified by the effective diffusion coefficient, D eff (see SM Sec. D2), where higher values correspond tomore fluid-like tissues. Points with the same parametersas in (a)-(c) are marked with corresponding symbols in(d) and (e).For a fixed l corresponding to rigid ECM, increasing γ results in a smooth increase in D eff , as shown in (e),while h a i changes non-monotonically, as shown in (d). Aswe increase γ beyond some critical value, γ c , defined bythe maximum of each curve in (d), h a i decreases, whilethe overall tissue geometry remains roughly circular. Aswe decrease γ below γ c , h a i also decreases, but ECMtension induces cavities. This behavior is indicative ofan instability, whose location we predict analytically andplot as a solid line in (e), while the Xs denote the loca-tions of γ c determined numerically. Cavitation Instability.
To analytically predict the on-set of the instability occurring at γ c , we devise a mean-field model in which a circular tissue containing a circularcavity is embedded in a regular, hexagonal spring net-work. We assume that all cells remain at their preferredarea and perimeter, leaving contributions only from thetissue boundary and spring network. This assumptionallows us to write the total energy as a function of thecavity radius, r , only (see SM Sec. F). For a fixed l value,we find two energy minima— one at r = 0, and anotherat r > γ which we denote γ MFc .For γ > γ
MFc , the energetically-favorable cavity radiusis zero, and the system prefers circular, compact tissuegeometries. For γ < γ
MFc , the system prefers a cavityof finite radius, which grows as γ decreases. The com-puted γ MFc agrees well with the measured γ c identifiedby the location of the peak in the h a i vs. γ curve forthe corresponding l (see Fig. 1 (e) and SM Sec. F). As l decreases, γ MFc increases, since larger interfacial ten- (b)(c)(d)
Cell shape, s (a)
FIG. 2.
Compact tissue regime results for fixed p = 3 . . (a) Minimum-energy configurations for a series of interfacialtension, γ , values. From top to bottom, γ =0.0, 0.53, 1.47, and1.79. Each cell center is colored by the cell’s shape, s = p/ √ a ,where p and a are the cell’s perimeter and area. (b) Theeffective diffusion coefficient, D eff , vs. γ . Below a criticalvalue of γ , D eff ≈ D eff > γ value at which theeffective target shape index, ˜ p , is equal to 3.81 is indicated bya vertical dashed line. (c) The effective target shape index,˜ p , as a function of γ . The horizontal dashed line denotes˜ p = 3 .
81, and the corresponding γ value is denoted with avertical dashed line. (d) The observed mean cell shape, h s i ,vs. ˜ p . Simulation data for varying γ is plotted with circles,while data for a bulk vertex model with varying p is plottedwith triangles for comparison. sion is required to overcome the increased tension in thenetwork and prevent the formation of cavities. Compact tissue regime.
When γ > γ c , the tissueis roughly circular, or “compact,” with no cavities. Tounderstand this regime, we first focus on the limit of zeroECM tension. Thus, increasing the interfacial tension, γ ,drives the tissue boundary to shrink, reducing the tissueradius and compressing the tissue, as shown in Fig. 2 (a).To study the effect of this compression on tissue fluidity,we compute the cells’ effective diffusion coefficient, D eff ,as a function of interfacial tension, γ . Fig. 2 (b) showsresults for a tissue with p = 3 .
68. At γ = 0, D eff ≈
0, and the tissue is solid-like. As γ increases and thetissue cells are compressed, D eff becomes non-zero at aparticular value of γ . This compression-induced fluidization can be under-stood via a mapping of the embedded tissue to an ef-fective bulk vertex model. By approximating the tis-sue geometry as circular, ignoring fluctuations in the cellperimeters, cell areas, and interface edge lengths, and as-suming that cell areas depend on the pressure balancegenerated by interfacial tension and spring tension, wegenerate an effective bulk model (see SM Sec. E):˜ e tot = ( ˜ h a i − + ˜ k p ( ˜ h p i − ˜ p ) . (2)The effective model parameters are denoted by tildesand are functions of the original model parameters, in-cluding γ , and the number of tissue cells, N cells . A plotof the effective parameter ˜ p is shown as a function of γ in Fig. 2 (c). Given the form of this effective energyfunctional, we predict that a system with a finite bound-ary and non-zero interfacial line tension will behave as abulk vertex model with effective parameters, ˜ p and ˜ k p .Although the critical p , p ,c , in the bulk vertex modelhas been shown to depend on temperature and simula-tion protocol, the zero-temperature value of p ,c = 3 . p ( p , γ, N cells ) ≈ .
81 (see Refs.[32–34] and SM Sec. A). Indeed we find that the valueof γ at which ˜ p = 3 .
81 corresponds well to the pointat which D eff becomes non-zero, as noted by the dashedlines in Fig. 2 (b) and (c). This mapping explains our ob-servation of fluid-like behavior with increasing interfacialtension; even for p < .
81, we expect fluid-like behaviorfor associated ˜ p > .
81. Moreover, since the relation-ship between p and mean observed cell shape, h s i , hasbeen shown for the bulk vertex model, our mapping alsopredicts the behavior of h s i with increasing γ , or, equiv-alently, ˜ p . This relationship, along with that for a bulkvertex model alone with varying p , is illustrated in Fig.2 (d).Although these results focus on the limit of floppyECM, adding tension to the network, by decreasing l ,simply shifts the mean cell area, as shown in Fig. 1 (d).By identifying another mean-field approximation relat-ing the mean spring length to the mean cell area, oureffective model parameters become functions of both γ and l (see SM Sec. E). Irregular tissue regime. At γ c , the tissue boundarybecomes unstable, and for γ < γ c , we observe irregulartissue boundaries facilitated by the opening of cavities–empty spaces along which cells are not coupled to ECMsprings. To study this regime, we first explore the limitof increasing tension in the spring network with no tis-sue interfacial tension. In this limit, there is no cost toopening gaps that increase the length of the spheroid-ECM interface, and so, as l decreases below the criticalthreshold and the tension in the ECM increases, gapsopen. Minimum energy configurations are illustrated inFig. 3 (a). As the spring network tension increases, thetissue boundary becomes increasingly irregular and theaverage cavity size grows.To study the phase of the tissue, we again run low-temperature simulations and extract the effective diffu-sion constant, D eff , in this case for a range of l values.For a fixed p and decreasing l , D eff decreases. Ex-amples of this are shown in Fig. 3 (b). To determineif the observed cell shape is correlated with the decreasein diffusivity, we ignore fluctuations in cell perimeters,areas, and spring lengths and treat each cell as border-ing a cavity and being stretched along the cavity’s circu-lar circumference. We also assume the cell areas do notincrease, which is substantiated by numerics. As eachcell becomes highly elongated in the tangential directionand compressed in the radial direction, the relationshipbetween the cell perimeter and the radius of the cavityapproaches a linear form (see SM Sec. G). Therefore, theperimeter of the cell scales with the square root of thearea of the cavity, and, since we assume the cell area isfixed, so does the cell shape. Assuming the total cav-ity area, A cav , is a multiple of the single cavity area, weobtain h s i = B p A cav + B . (3)In Fig. 3 (c), simulation results for p = 3 .
70 and afit of the form of Eq. 3 are shown. Since we antic-ipate bulk tissue behavior when A cav = 0, we expect B ≈ .
81 (see SM Sec. A). We find that B = 0 . B = 3 . com-pact tissue regime , a mean cell shape above p ,c does notnecessarily correspond to fluid-like tissue.Prior work has shown that bulk tissue underanisotropic strain can rigidify due to cell-cell alignmentfor h s i > p ,c [34]. To determine if this phenomenonis relevant here, we quantify cell-cell alignment with theparameter Q [34], which is similar to a nematic order pa-rameter except that it also contains information aboutthe shape of the cells (see SM Sec. H). As illustratedin Fig. 3 (d), the decrease in D eff is indeed correlatedwith an increase in cell shape and an increase in cell-cell alignment. Zero-temperature energy barrier mea-surements also support this finding (see SM Sec. I).Note that although there is a clear decrease in D eff as l decreases, a difference remains across p values. Thissuggests that the rigidification may not necessarily be abulk phenomenon but become enhanced toward the pe-riphery of the spheroid. Energy barrier measurements asa function of distance from the periphery suggest this.Finite-size and temperature studies are in SM Secs. Jand K. Discussion.
We have generalized the bulk vertex modelfor cellularized tissue to include a coupling to a surround-ing spring network representing the ECM. By studyingthis bi-material, we gain insight into the competition be-tween cell-cell and cell-environment interactions. Our (a) (c)(b)(d)
Spring tension
Cell shape
FIG. 3.
Irregular tissue regime. (a) Minimium-energy con-figurations for tissues with p = 3 .
70 and increasing ECMtension (decreasing l ) for γ = 0. From top to bottom, l =0.63, 0.61, 0.60, and 0.59. Each spring is colored by its tension( T h ij i = k sp ( l h ij i − l )), and each cell center is colored by itsshape, s = p/ √ a , where p and a are the cell’s perimeter andarea. (b) The effective diffusion coefficient, D eff , vs. l of thesurrounding spring network, for a range of p values for thetissue cells. (c) The mean cell shape, h s i , vs. the total cavityarea, A cav , for a tissue of p = 3 .
70 (points), and a fit of theform in Eq. 3 (dashed). (d) A heat map colored by D eff , forvarying mean cell shape, h s i , and cell-cell alignment, Q . findings represent the essential first step in understandinga model that can now easily be augmented with specificfeatures. When the tissue remains compact, isotropiccompression– driven by increasing interfacial tension–leads to increased cell shapes and fluidization of the tis-sue, as predicted by a mean-field mapping to a bulk tis-sue. Moreover, a transition from a compact tissue ge-ometry to an irregular one happens suddenly, due to ageneric instability in the balance of negative pressure in-duced by the surrounding network and positive pressurecaused by the interfacial tension. Once gaps open be-tween the tissue and the ECM, the rheology of the tissuedepends both on cell shape and alignment, as in the caseof tissues under shear deformations [34]. We note thatthe tissue deformation in our system is non-uniform, dueto the system’s disorder and heterogeneity, resulting inheterogeneity in alignment and shape among the tissuecells. It has been shown that heterogeneity itself canpromote rigidity in tissues [35], which perhaps also con-tributes to the induced rigidity of the system.Our results suggest that although stiff networks havebeen shown to induce tissue breakup and cell migration,this may not mean they induce tissue fluidity. In fact,stiff networks may induce stiffer tissues, due to emer-gent cellular alignment, despite irregular tumor bound-aries and elongated cell shapes. In other words, describ-ing cancer invasion as cellular unjamming may be over-simplifying the phenomenon, as some experiments haverecently indicated [36]. The pre-invasion stage of the tu-mor may involve solidification due to cellular alignment,ultimately facilitating cellular streaming upon enhancedcellular activity, with the tissue becoming an active fluid,flowing with respect to the ECM. Measurements of cel-lular alignment are needed to test this idea. Experi-ments searching for compression-induced fluidization incompact tumors or looking for cavity formation at thetumor-ECM interface would also test the model, since,given the generic nature of our arguments, we expect thephenomena to survive in three-dimensions. Finally, ourmodel lays the groundwork for understanding more com-plex scenarios. For example, the ability for cells to formor degrade cell-ECM contacts, cell self-propulsion, ECMreorganization, and feedback between cell-cell and cell-ECM interactions are all fascinating directions for futurestudy.We acknowledge nsf-pols 1607416 for financial support.We would like to thank Daniel Sussman and Preeti Sahufor stimulating discussions and guidance on use of thecellGPU codebase. We acknowledge Mingming Wu forfeedback on early results. [1] G. Cooper, The Cell: A Molecular Approach. 2nd edition (Sunderland (MA): Sinauer Associates, 2000).[2] C. M. Niessen, D. Leckband, and A. S. Yap, Physiol.Rev. (2011).[3] L. Li, Y. He, M. Zhao, and J. Jiang, Burns & trauma ,21 (2013).[4] M. A. Garcia, W. J. Nelson, and N. Chavez, Cold SpringHarbor Perspect. Biol. (2018).[5] B. Bruckner and A. Janshoff, Sci. Reps. (2018).[6] C. G. Galbraith and M. P. Sheetz, Curr. Opin. Cell Biol. (1998).[7] C. M. Lo, H. B. Wang, M. Dembo, and Y. L. Wang,Biophys. J. (2000).[8] A. Liu, O. Chaudhuri, and S. Parekh, Int. J. Integr. Biol. (2017).[9] J. Humphrey, E. Dufresne, and M. Schwartz, Nat. Rev.Mol. Cell Biol. (2014).[10] C. Frantz, K. Stewart, and V. Weaver, J. Cell Sci. (2010).[11] C. Walker, E. Mojares, and A. del Ro Hernndez, Int. J.Mol. Sci. (2018). [12] M. Giussani, G. Merlino, V. Cappelletti, E. Tagliabue,and M. G. Daidone, Semin. Cancer Biol. (2015).[13] I. Acerbi, L. Cassereau, I. Dean, Q. Shi, A. Au, C. Park,Y. Y. Chen, J. Liphardt, E. S. Hwang, and V. M.Weaver, Int. Biol. (2015).[14] J. L. Leight, A. P. Drain, and V. M. Weaver, Annu. Rev.Cancer Biol. (2017).[15] A. J. McKenzie, S. R. Hicks, K. V. Svec, H. Naughton,Z. L. Edmunds, and A. K. Howe, Sci. Rep. (2018).[16] K. S. Kopanska, Y. Alcheikh, R. Staneva, D. Vignjevic,and T. Betz, PLOS ONE (2016).[17] Y. J. Suh, M. S. Hall, Y. L. Huang, S. Y. Moon, W. Song,M. Ma, L. J. Bonassar, J. E. Segall, and M. Wu, Integr.Biol. (2019).[18] J. P. Baskaran, A. Weldy, J. Guarin, G. Munoz, M. Kot-lik, N. Subbiah, A. Wishart, Y. Peng, M. A. Miller,L. Cowen, and M. J. Oudin, bioRxiv (2020).[19] Y. Zheng, H. Nan, Y. Liu, Q. Fan, X. Wang, R. Liu,L. Liu, F. Ye, B. Sun, and Y. Jiao, Phys. Rev. E (2019).[20] M.-C. Kim, Y. R. Silberberg, R. Abeyaratne, R. D.Kamm, and H. H. Asada, PNAS (2018).[21] Y. L. Han, P. Ronceray, G. Xu, A. Malandrino, R. D.Kamm, M. Lenz, C. P. Broedersz, and M. Guo, Proc.Nat. Acad. Sci. (2018).[22] T. Nagai and H. Honda, Philos. Mag. B - Physics of Con-densed Matter, Statistical Mechanics, Electronic, Opti-cal, and Magnetic Properties (2001).[23] R. Farhadifar, J.-C. Rper, B. Aigouy, S. Eaton, andF. Jlicher, Curr. Biol. (2007).[24] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning,Nat. Phys. (2015).[25] M. Merkel and M. L. Manning, New J. of Phys. (2018).[26] P. Sahu, D. M. Sussman, M. Rbsam, A. F. Mertz,V. Horsley, E. R. Dufresne, C. M. Niessen, M. C.Marchetti, M. L. Manning, and J. M. Schwarz, Soft Mat-ter (2020).[27] 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, and et. al., Nat. Mater. (2015).[28] C. Malinverno, S. Corallino, M. Giavazzi, Fabioand Berg-ert, Q. Li, M. Leoni, A. Disanza, E. Frittoli, A. Oldani,E. Martini, T. Lendenmann, and et. al., Nat. Mater. (2017).[29] A. J. Licup, S. M¨unster, A. Sharma, M. Sheinman, L. M.Jawerth, B. Fabry, D. A. Weitz, and F. C. MacKintosh,PNAS (2015).[30] K. A. Jansen, A. J. Licup, A. Sharma, R. Rens, F. C.MacKintosh, and G. H. Koenderink, Biophys. J. (2018).[31] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L.Manning, PNAS (2019).[32] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Phys. Rev. X (2016).[33] D. M. Sussman and M. Merkel, Soft Matter (2018).[34] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan,M. L. Manning, and K. E. Kasza, PNAS (2020).[35] X. Li, A. Das, and D. Bi, Phys. Rev. Lett. (2019).[36] W. Kang, J. Ferruzzi, C.-P. Spatarelu, Y. L. Han,Y. Sharma, S. A. Koehler, J. P. Butler, D. Roblyer, M. H.Zaman, M. Guo, Z. Chen, A. F. Pegoraro, and J. J. Fred-berg, bioRxiv (2020). SUPPLEMENTARY MATERIALA. Vertex models for tissues
Vertex models of epithelial tissues represent groups ofcells as space-tiling polygons. A configuration is definedby a set of vertex positions and, for each vertex, a setof neighbor vertices. The connections between verticesand their neighbors form the edges of the polygonal tis-sue cells. The degrees of freedom in the system are thecoordinates of the vertices, which evolve in a simulation,enabling changes in the shapes of the tissue cells overtime. For a 2D vertex model, a typical energy functionalis E V M = N cells X α =1 K A,α ( A α − A ,α ) + N cells X α =1 K P,α P α + 12 N edges X h ij i =1 Λ h ij i l h ij i (S1)where A α and P α are the area and perimeter of cell α , l h ij i is the length of edge h ij i , and A ,α , K A,α , K P,α and Λ h ij i are the target area, area stiffness, contractil-ity strength, and line tension of cell α , respectively. Theterms in the first sum in Eq. S1 capture the 2D incom-pressibility of the cells. The second term captures activecontraction in acto-myosin cables, and the third termcaptures cell-cell adhesion and cortical tension generatedby acto-myosin [S1, S2]. If we assume that all constantsare the same for each cell and edge, we can pull the con-stants outside the sums, and write A ,α as simply A .The sum in the last term can now be written in terms ofcell perimeters, since it is a sum over all cell edge lengths.We then have E V M = K A N cells X α =1 ( A α − A ) + K P N cells X α =1 P α + Λ N cells X α =1 P α . (S2)Now, for each cell, we combine the last two terms and,dropping constants that simply shift the total energy, weobtain˜ E V M = K A N cells X α =1 ( A α − A ) + K P N cells X α =1 ( P α − P ) , (S3)where P = − Λ2 K P is the target cell perimeter, commonto all cells. Non-dimensionalizing Eq. S3, using K A A asour unit of energy, reduces the number of free parametersto leave us with e V M = N cells X α =1 ( a α − + k p N cells X α =1 ( p α − p ) , (S4)where we use lower-case letters to indicate dimensionlessparameters (and drop the tilde on the left-hand side), k p = K P K A A and p = P √ A . The parameter p is the dimensionless target perime-ter or the target, or preferred, “shape index” of the cells.This follows from defining the observed “shape” of a cellas s = P/ √ A , where P and A are the actual perimeterand area of the cell, respectively. This quantity, s , cap-tures how circular a cell is, with lower values being morecircular and higher values being less circular (which inmany cases means more elliptical or elongated). Since p is defined as P √ A , it therefore acts as a “target” shapefor the cells.It has been shown that the preferred shape, p , can actas a tuning parameter for the bulk 2D vertex model. Oneway to measure the phase of a tissue is to calculate aneffective diffusion coefficient, D eff , for the cells, which isapproximately zero for p < p ,c and greater than zero for p > p ,c [S2]. In order for cells to diffuse in a system withno gaps between them (a “confluent” tissue), they mustundergo T1 transitions, which allow them to exchangeneighbors even with the confluency constraint (see Sec.C for more information).This transition also manifests in the relationship be-tween the observed mean cell shape, h s i , and preferredcell shape, p . In the 2D, disordered vertex model, for p < p ,c , cells become frustrated as they cannot achievetheir preferred shape, while for p > p ,c , h s i = p . Thiscomes from the fact that for a fixed area, the perimeterof a polygon has a lower bound, but no upper bound.Previous work has shown that the value of p ,c dependson temperature (the magnitude of thermal fluctuations)and simulation protocol, and the often-cited p ,c = 3 . p = 3 .
81, but donot expect it to match exactly.
B. Spring network models
A spring network model is a system, like the vertexmodel for tissues, consisting of vertices and vertex-vertexneighbors, with the positions of the vertices as the de-grees of freedom. A simple Hamiltonian for a springnetwork model assigns a spring-like energy cost to a de-viation of a vertex-vertex neighbor distance, L h ij i , fromsome “preferred” or “target” distance, L : E SNM = 12 N sp X sp h ij i K sp h ij i ( L sp h ij i − L ) . (S5)This model does not include bending penalties that aretypically included in fiber network models. However, re-cent work suggests that the mechanics of such networks a r X i v : . [ q - b i o . CB ] J un is largely controlled by the critical point in simple springnetwork models [S6].This model does not have particles or cells that cantransition from a diffusive phase to a caged phase. Onecan, however, quantify the difference between a fluid andsolid phase by considering the response to an imposedstrain. A floppy (fluid-like) system has no mechanical re-sponse to strain, while a solid system resists strain in theform of an energy cost. Spring networks, even withoutmore physically-realistic bending energies, also exhibitthis behavior.The method of Maxwell constraint counting demon-strates that determining whether a network will be floppyor rigid at small strain is equivalent to determiningwhether there are more degrees of freedom or constraintsin the network [S7]. If the number of degrees of freedomis greater than the number of constraints, the systemis referred to as “under-constrained,” and the system isfloppy. However, even a floppy network can become rigidabove a critical strain value, at which point the networkbecomes geometrically frustrated. In other words, due tothe connectivity of the network, there is a point at whichthe new imposed deformation of the network cannot beaccommodated without an energy cost, and the systemis rigid [S8, S9].In order to impose a strain on a spring network, one canexternally deform the network’s boundaries. In simula-tions, this is often done by shearing or uniformly expand-ing the boundaries. Alternatively, for fixed boundaries,one can impose a new L value for the springs. The effectof this is similar to that of imposing an external uniformexpansion of the network, without having to introducea new strain parameter. Indeed, we see that a simpletwo-body spring network will transition from a floppy torigid phase as a function of L , as shown in Fig. S1.It is not a coincidence that both the spring networkand the vertex model have a rigid-to-floppy transitionthat is tuned by a preferred length parameter ( L and p , respectively). In fact, it has been shown that theseparameters are deeply related, and that the transitionsin both models are of the same universality class [S9]. C. T1 transitions
In a confluent tissue, there are no gaps between cells.In a fluid-like state, however, cells move through the tis-sue, flowing like a liquid under external forces. This ispossible due to a process called a T1 transition. Dur-ing a T1 transition, an edge in the network shrinks to apoint and a new edge then expands from that point, inthe direction perpendicular to the original edge (see Fig.S2). In this way, cells can exchange neighbors and trans-late, while staying in the plane of the tissue and withoutneeding empty space to move into.In our bi-material, any edge between two tissue-like
FIG. S1. The average spring network tension as a functionof spring rest length, l , for a system of 500 polygons, whereeach edge is assigned a spring energy. In this example, thesystem is purely a spring network, with no embedded tissuecells. At a critical value of l , the tension in the networkincreases by about 10 orders of magnitude. FIG. S2. An illustration of a T1 transition in bulk tissue.Tissue cells are labeled with letters and vertices defining theedge undergoing the T1 are labeled with numbers. Edge 1-2 shrinks to a point (not shown) and then elongates in thedirection perpendicular to its initial configuration. Noticethat cells B and D are no longer neighbors after the T1, whilecells A and C have become neighbors. This demonstrates thatcells can exchange neighbors and therefore move through atissue, even when there are no spaces between cells. vertices is allowed to undergo a T1. This refers to edgesin the tissue bulk. In the spring network, the opposite istrue: no T1 transitions are allowed along spring edges.Edges that are between two interface vertices or an in-terface vertex and a tissue vertex, however, have to beconsidered separately. In describing the cases of inter-est, we will refer to Fig. S3, where interface edges arerepresented as thick lines, tissue edges are representedas double lines, and spring edges are represented as aspring icons. Tissue cells are labeled with capital lettersand the edge that is undergoing the T1 has vertices la-beled 1 and 2. Polygons that are not tissue cells but areinstead empty space in the spring network are given nolabel.The first case to consider is illustrated in Fig. S3 (a).In this case, vertex 1 is also a neighbor to a spring-typevertex and another interface-type vertex, and vertex 2 (a) (b)(c) (d)
FIG. S3. Illustrations of four special cases of T1 transitionsinvolving edges along the tissue-ECM boundary. The ver-tices of the edge undergoing a T1 are labeled 1 and 2 in eachcase. Tissue cells are labeled with capital letters, while emptyspace in the ECM network is unlabeled. Spring edges arerepresented using a coil icon, edges on the interface are rep-resented with solid lines, and edges between two tissue cellsare represented with double lines. is also a neighbor to a tissue-type vertex and anotherinterface-type vertex. In this case, a T1 simply pushesthe cell that shared edge 1-2, cell B, along the boundary.The edge remains an interface-type edge. The reverse ofthis case (moving from the right-hand-side of the figureto the left) is also allowed and the procedure is identicalto the forward case, except that the vertex labels 1 and2 are switched.The second case, Fig. S3 (b), concerns an interfaceedge in which both vertices 1 and 2 have a spring-typeneighbor vertex and another interface-type neighbor ver-tex. Here, a T1 pushes the cell that shares edge 1-2, cellA, away from the boundary, into the bulk. Edge 1-2 isthen shared by two spring-network-type polygons, andshould be a spring-type edge. Therefore it automaticallygets relabeled as a spring. The reverse of this is not al-lowed, as T1s are not allowed along spring-type edges. Inother words, an interface edge that goes through this T1and becomes a spring cannot revert back to an interfaceedge. Note that an interface edge is lost in this T1.The third case, Fig. S3 (c), is analogous to the sec-ond, except that instead of springs, vertices 1 and 2 eachhas a tissue-type neighbor in addition to their interface-type neighbor. Here the T1 also forces the adjacent cell,cell B, into the bulk, but the edge now changes from aninterface-type to a tissue-type due to its new cell neigh-bors. The reverse of this procedure is allowed and is de-scribed in case 4. Case 3 results in a loss of an interfaceedge, while the reverse leads to interface edge creation.Case four, Fig. S3 (d), concerns a T1 between aninterface-type vertex and a tissue-type vertex. After aT1, this originally tissue-type edge is now on the bound-ary, as a bulk cell, cell C, has moved from the bulk tothe boundary. Therefore this edge gets updated to aninterface-type edge, i.e. a new interface edge emerges. The reverse of this process is described in case three.While both interface edge creation and removal areallowed, interface edge creation via the cavitation insta-bility is much more common than interface edge removal.
D. Simulation details: algorithms and parameters
To numerically study the bi-material, we implementa modified version of the open-source cellGPU soft-ware [S10]. This molecular dynamics software is writ-ten specifically for vertex-like models of tissues and hasthe ability to run on GPUs for especially efficient perfor-mance, although it can also be run on CPUs.As is true for most molecular dynamics software, cell-GPU uses a force field-plus-equation of motion approach,where one can essentially “mix-and-match” potentialenergy functions and equations of motion. Therefore,adding the force field (potential energy and correspond-ing force functions) for our coupled tissue-ECM model,and all necessary data structures not already included(like spring moduli, interface edge moduli, etc.), allowsus to then interface with the already-existing algorithms(described below) for zero-temperature energy minimiza-tion and low-temperature Brownian dynamics.For all simulations, we begin with N total = 1000 poly-gons tiling a periodic box. We then define a circle whosecenter is the box center and whose area is 20% of the totalbox area. Any polygon whose center is within this circleis labeled a tissue-type cell, and any polygon outside thiscircle is labeled an ECM-type polygon. The edges sharedby different types of polygons are labeled interface-type.The number of tissue-type cells varies slightly from con-figuration to configuration, but is always very close to1000 × .
20 = 200 cells. We have studied additional ra-tios of the number of cells to springs (see Sec. J).
1. Energy minimization simulation procedure
All results except those for the effective diffusion co-efficient are obtained using energy minimization simula-tions. Specifically, we use the FIRE minimization algo-rithm [S11], where the equation of motion is given by ~v = ~v + ~F ∆ t + α ( | ~v | ˆ F − ~v ) . For a given particle’s (vertex’s) current position and ve-locity, the algorithm calculates the current forces on theparticle ( ~F = −∇ E ( ~x )) and the power, P = ~F · ~v . If P is negative, it means that the forces would move usback “uphill,” so we decrease our timestep, dt , by somefraction, f dec , such that dt = dt × f dec , set the velocityto zero, and set α back to its initial value. If instead P ispositive, and the number of timesteps since P was neg-ative is greater than some chosen minimum value, n min , Parameter Value dt start × − α start . dt max × dt start f inc . f dec . α dec . n min F cut − off × − TABLE SI. Typical parameter values used during FIRE min-imization simulations.Parameter Value µ . T × − dt . k B . we increase our timestep by some fraction, f inc (unless weare already at our maximum allowed timestep, dt max , inwhich case we do not change the timestep) and decrease α by some fraction f α . The simulation ends when themaximum force in the system is less than or equal to thechosen force cut-off, F cut − off . Typical values for theseparameters used in our simulations are given in table SI.
2. Brownian dynamics simulation procedure and calculationof mean squared displacements and effective diffusioncoefficients
In cases where we are not just interested in the staticenergy-minimum of a system, we use Brownian dynamics(over-damped Langevin dynamics) to simulate the mo-tion of vertices over time. This equation of motion isgiven by ~x = ~x + µ ~F dt + p µk B T dt ~R ( t )where ~F = −∇ E ( ~x ), µ is an inverse damping coeffi-cient, T is the temperature, and ~R ( t ) is a vector of delta-correlated random numbers with a Gaussian distributionand zero mean. Typical values used for the parametersin our simulations are given in table SII.For each simulation, we run at least 1 × initialtimesteps at which point the system has come to a steadystate. We then run for an additional 1 × timesteps FIG. S4. Mean squared displacement (MSD) for fixed p =3 .
68 and varying interfacial tension, γ , values (and thereforevarying ˜ p values), vs. simulation time. The final data pointsare used to calculate the effective diffusion coefficient for each˜ p , which are reported in the main text. during which data on the state of the system is peri-odically saved. We use this data to calculate the meansquared displacement (MSD) of the tissue cell centers.This is defined asMSD( τ ) = h ( ~r ( t + τ ) − ~r ( t )) i where t is the time of some reference state, in units of ourtimestep. This average, for a single simulation, is overcells and reference times. The MSDs for all simulationsare then averaged and this mean and standard deviationis used as our result.The effective diffusion coefficient, D eff , for the cells isdefined as D eff = D s T = 1 T lim τ →∞ MSD( τ )4 τ . where, instead of taking τ to infinity, we use τ = 4 × (in units of the simulation timestep length) as our long-time estimate.An example of a set of MSDs are shown in Fig. S4.The D eff values calculated from this particular exampleare used to produce Fig. 2 (b) in the main text. E. Analytical model: compact tissue regime
To model the full compact tissue regime , our goal isto express the sum over interfacial edges and sum oversprings in terms of the mean cell area, in order to con-struct an effective theory of only cell areas and perime-ters. To do so, we first simplify the energy function (Eq.1 of the main text) by ignoring fluctuations from themean cell area, h a i , mean cell perimeter, h p i , mean in-terface edge length, h l int i , and mean spring length, h l sp i .This gives e total = N cells ( h a i − + k p N cells ( h p i − p ) + γN int h l int i + 12 k sp N sp ( h l sp i − l ) , (S6)where N sp is the total number of springs, N cells is thetotal number of cells, and N int is the total number ofinterface edges.The tissue’s total perimeter is equal to the sum overthe lengths of the interface edges. Ignoring fluctuationsgives P tissue = X int h ij i l int h ij i = N int h l int i . (S7)At the same time, since the tissue in this regime is highlycircular, the perimeter is also related to the tissue radius, R tissue : P tissue = 2 πR tissue . (S8)Therefore, the R tissue is given by R tissue = N int h l int i π . (S9)Since the tissue is circular, its area is also a function ofits radius, and we can write the total tissue area in termsof Eq. S9: A tissue = πR tissue = N int h l int i π . (S10)The area of the tissue must also be the sum of the con-stituent cell areas, which we approximate as each equalto the mean cell area. This gives A tissue = X cell,α a α = N cells h a i . (S11)Combining Eq. S10 with Eq. S11 gives the average cellarea as h a i = N int h l int i πN cells (S12)Rearranging Eq. S12, we write the mean interface edgelength as a function of the mean cell area: h l int i = s πN cells N int p h a i (S13)Now, to write h l sp i in terms of h a i , we approximatethe spring network representing the ECM as a regular,hexagonal lattice. In that case, the number of hexagonsin the network is N hex = 16 ( N springs )(2 hexagonsspring ) = 13 N springs , (S14) since each spring is shared by two hexagons. The factorof 1 / A sp.net = N hex A hex = 13 N sp ( 3 √ h l sp i ) , (S15)where A hex is the area of a regular hexagon with edgelength equal to h l sp i . The total area taken up by thespring network is also equal to the total box area minusthe total tissue area, which is itself a sum of the cell areas.Ignoring fluctuations in cell areas as before, this gives A sp.net = A box − A tissue = A box − N cells h a i . (S16)Combining Eqs. S15 and S16 and solving for h l sp i gives h l sp i = s √ s A box − N cells h a i N sp . (S17)Now let h a i = 1 + (cid:15) , where | (cid:15) | < (cid:15) can be positiveor negative. Then Eq. S13 becomes h l int i = s πN cells N int √ (cid:15). (S18)Taylor expanding the square root about (cid:15) = 0 yields h l int i ≈ s πN cells N int (1 + 12 (cid:15) − (cid:15) ) . (S19)Similarly, substituting h a i = 1 + (cid:15) to Eq. S17 andexpanding yields h l sp i = s √ s A box − N cells N sp r − N cells A box − N cells (cid:15) ≈ s √ b (1 − b (cid:15) − b (cid:15) ) , (S20)where b = p ( A box − N cells ) /N sp and b = N cells / ( A box − N cells ). The total energy is now e tot = N cells (cid:15) + γ p πN cells (1 + 12 (cid:15) − (cid:15) )+ 12 k sp N sp s √ b (1 − b (cid:15) − b (cid:15) ) − l ! + k p N cells ( h p i − p ) . (S21)We assume that in this regime, h p i does not depend on h l int i or h l sp i . The cell perimeters are not coupled to thetotal tissue perimeter in the way that cell areas are cou-pled to the total tissue area. Instead, for a given fixed cellarea, the cell perimeter can vary, given the constraints oncell shape.At this point, we simplify things by first taking thelimit of no ECM tension, or, equivalently, k sp →
0. Com-bining powers of (cid:15) and completing the square then gives e tot = ( N cells − γ p πN cells ) (cid:15) + 1 √ πN cells π γ − ! + k p N cells ( h p i − p ) + B, (S22)where B is a collection of constants that simply shifts thetotal energy. Dropping this term, and putting the energyin terms of h a i again, we get e tot = ( N cells − γ p πN cells ) h a i−
1+ 1 √ πN cells π γ − ! + k p N cells ( h p i − p ) , (S23)which we rewrite as e tot = k a ( h a i − a ) + k p N cells ( h p i − p ) (S24)and define k a = N cells − γ p πN cells (S25)and a = 1 − √ πN cells π γ − . (S26)We now divide both sides by a (dimensionless) con-stant, k A a to arrive at˜ e tot = ( ˜ h a i − + ˜ k p ( ˜ h p i − ˜ p ) , (S27)which is the same as Eq. 2 in the main text. The effec-tive model parameters are denoted by tildes and definedin Table SIII. They are each functions of k a and/or a ,which are defined in Eqs. S25 and S26. It is this effectivemodel that we use to produce the results in Fig. 2 (c)and (d) of the main text.If we do not ignore the third term in Eq. S21, we canexpand it, drop terms greater than O ( (cid:15) ), again collectterms that are linear and quadratic in (cid:15) , complete thesquare, and divide by a new a . This yields Eq. S27, butwith new definitions of ˜ p and ˜ k p . We do not detail thealgebra here, but show the results of this reformulationin Fig. S5, which demonstrates the behavior of the meancell area, h a i , as a function of γ , for a set of l valuesand fixed p = 3 .
91. The measured values from simu-lations are shown as solid curves. We assume that thecells reach their preferred areas ( h a i = a ) and therefore Effective Original˜ e total e total k a a ˜ h a i h a i a ˜ h p i h p i √ a ˜ p p √ a ˜ k p k p N cells k a a TABLE SIII. The mapping from the new, effective modelparameters to the original, dimensionless model parameters.FIG. S5. The mean cell area, h a i , as a function of interfa-cial tension, γ , for varying spring rest length, l . Simulationresults are shown in solid lines and the theory is shown indashed lines. plot a ( γ, l , k sp , N cells , N springs , A box ) from our model,on top of the observed values, using dashed curves. Al-though our analytical model ignores the disorder andfluctuations in our simulations, we find good agreementbetween the predicted and measured values for the cellareas. F. Analytical model: predicting the criticalinterfacial tension
In Fig. S5, it is clear that the model described in Sec.E, represented by dashed lines, breaks down below a spe-cific value of γ for each curve. We call this special point γ c , and define it as the value corresponding to the peakof the mean cell area, h a i , vs. γ curve, for a fixed p and l . For γ > γ c , the tissue is compact, and compact tissueregime model works well to describe h a i . For γ < γ c , thetissue is irregular and cavities are present, and specificallyfor γ ≈
0, our irregular tissue regime model works wellto describe the cell shapes (see Sec. G). At γ c , the sys-tem clearly switches sharply from one of these regimes tothe other, and we hypothesize that this sharpness is dueto a sudden change in whether cavities are energetically-favorable. In other words, γ c signifies the onset of aninstability in the tissue-ECM interface.To test this hypothesis, we devise a simple toy model,in which the tissue is circular, with radius, R , and con-tains a cavity, also circular, of radius, r . The energeticcontribution from the boundary is E int = γ (2 πr + 2 πR ) . (S28)We approximate the energetic contribution from thespring network as E sp.net = 12 k sp N sp ( h l sp i − l ) . (S29)The total area of the spring network is related to R as A sp.net = A box − πR , (S30)and therefore, using the same description of the springnetwork as in Sec. E and applying Eq. S15, the meanspring length, h l sp i , is h l sp i = s √ s A box − πR N sp . (S31)We assume the energetic contributions from the ten-sion in the spring network and the interfacial tensiondominate, and therefore deviations in cell areas andperimeters are small. Specifically, we let h a i = 1, andtherefore A tissue = πR − πr = N cells , (S32)which allows us to write R as a function of r : R = r N cells π + r . (S33)The total energy is now E tot ( r, γ, l ) = 2 πγ ( r + R ( r ))+ 12 k sp N sp ( h l sp i ( r ) − l ) . (S34)We minimize this energy with respect to r , for fixed l and γ , and solve for r ∗ – the value of the cavity radius, r ,that minimizes the energy. The results are shown as thesolid curves in Fig. S6. For each fixed l value, we findand plot r ∗ for a range of γ values.For each l value, there is a jump in r ∗ at a particularvalue of γ . For γ above this point, the energetically-favorable cavity radius is zero, and the system preferscircular, compact tissue geometries. Immediately belowthis point, however, the system suddenly prefers a cavityof finite radius, which continues to grow as γ decreasesfurther. To compare the γ value at which we see a jumpto γ c , we also plot the location of the peaks of h a i vs. γ curves as Xs. We find good agreement between the γ val-ues associated with the peaks and the points at which r ∗ becomes non-zero, meaning that our instability hypoth-esis successfully predicts γ c . r * l FIG. S6. The instability model prediction for the criticalcavity radius, r ∗ , as a function of interfacial tension, γ , isshown in solid lines for varying spring rest length, l , values.Each X indicates the critical γ value, γ c , extracted from themaximum of the h a i vs. γ curve for the corresponding l . γ c values and their error bars are found by binning data pointsnear each maximum and taking the largest bin mean and itsstandard deviation. G. Analytical model: irregular tissue regime
When the surrounding spring network tension is high,the minimum energy configurations for the tissue be-come extremely irregular. We observe an increase in theperimeter of the cells with increasing spring network ten-sion and with increasing cavity area between the springnetwork and tissue. To understand this, we employ yetanother toy model. We start with a collection of squarebricks, with side length b , compactly arranged as in Fig.S7, (a). In our tissue-spring network simulations, thesystem is driven to decrease the spring network area anddoes this by creating more tissue-ECM boundary. In ourtoy model system, we can increase the total boundaryby simply arranging the blocks in a line, as in Fig. S7,(b). If we now want to increase the boundary even more,we must start to deform the blocks. To do this withoutchanging the area of each individual block (which wouldincur an energetic cost in the actual system), we simplyexpand each block along the line and shrink it perpen-dicular to the line in such a way that the area does notchange, as shown in Fig. S7, (c). As the sides of theblocks perpendicular to the line shrink and the parallelsides grow, the perimeter of each block becomes propor-tional to the length of the line.If we now imagine wrapping this line around a cavity,its length becomes the circumference of this cavity, andeach block’s perimeter is therefore proportional to thiscircumference, and in turn, the radius, of the cavity. Thismeans that each block’s perimeter is proportional to thesquare root of the area of the cavity. Multiplying thiscavity area by some value gives us the total cavity area, (a)(b)(c) b b b b 2b 0.5b FIG. S7. Illustration of toy model used to explain the increasein mean cell shape in the irregular tissue regime . A cav , and we can write h p i = B p A cav + B . (S35)Since the shape, s , is defined as s = p/ √ a , and we assume a = a = 1, we can replace h p i with h s i and we arrive atEq. 3 of the main text. H. Defining and computing the cell-cell alignmentparameter
To quantify the cell-cell alignment, we define and com-pute Q , as in Ref. [S5]. The value of this quantity indi-cates how aligned a set of cells are, modulated by theircircularity. In other words, unlike a simple nematic or-der parameter, Q weighs the alignment of highly polar-ized shapes more than that of less polarized shapes, evenif the orientations of their axes are identical. This as-sures that regular, symmetric shapes are assigned zeroalignment, regardless of their orientation relative to eachother.To compute Q for a set of cells, one first triangulatesthe system with each vertex shared by three cells. Byconnecting the centers of these three cells, we produce atriangle. Doing this for each vertex in our system pro-duces a set of triangles– the dual of the original networkof vertices. The general idea is now to determine the de-formation of each triangle from some reference triangle,and compare these deformations to see how similarly-elongated the triangles in the set are. The details are asfollows:1. Define a reference triangle. For example, an equi-lateral triangle with area 1. Its side length is givenby d ref = 2 p √ m = d ref × × (cid:20) −√ −√ − (cid:21)
2. For each triangle, list the coordinates of its cor-ners as [ r A , r B , r C ], such that they are in counter-clockwise order.3. Construct a matrix that defines the current trian-gle: m = (cid:20) r x,B − r x,A r x,C − r x,A r y,B − r y,A r y,C − r y,A (cid:21)
4. Compute the shape tensor, S , of the current trian-gle, which is defined as S = m × m
5. Compute the trace part, t , of S : t = 12 Tr[ S ] × (cid:20) (cid:21)
6. Compute the symmetric, traceless part of S : S symm = 12 ( S + S T ) − t
7. Compute the antisymmetric part of S (this is al-ready traceless): S asymm = 12 ( S − S T )8. Compute the rotation angle, θ , of the triangle. Inprogramming languages with the arctan2 function,this can be found as θ = arctan2( S asymm,yx , t xx )9. Compute the “norm” of the symmetric, tracelesspart of S . Note this is not defined in the standardway. | S symm | = ( S symm,xx + S symm,xy ) /
10. Compute the determinant of S , det( S ).11. Finally, compute the elongation tensor of the tri-angle, q : q = 1 | S symm | sinh − ( | S symm | p det( S ) )( S symm × R )where R is the 2D, clockwise rotation matrix: R = (cid:20) cos θ sin θ − sin θ cos θ (cid:21) Q of the system is defined as the weighted averageof all q tensors (for all triangles): Q = P i a i q i P i a i where the sum is over triangles, i , and a i is the areaof triangle i . The area of triangle with its cornerslisted in counter-clockwise order can be computedusing the “shoelace formula” as: a i = 12 ( r x,A ( r y,B − r y,C )+ r x,B ( r y,C − r y,A )+ r x,C ( r y,A − r y,B ))13. Q is a tensor, and we report the norm of it, givenby | Q | = ( Q xx + Q yy ) / This procedure can be done for any group of cells (andcorresponding group of triangles). In the irregular tissueregime , our tissue is being expanded roughly radially bythe tension in the surrounding network. The cells are notall being deformed in the same direction and we thereforedo not expect the system to be globally aligned. Instead,high cell-cell alignment in our system means high local alignment. To measure this, we define, for each cell, a setof triangles using the centers of the cell and its neighbors.We compute Q for this set and then average over all sets.This gives a measure of the mean local alignment for agiven tissue configuration. I. Measuring T1 energy barriers
In addition to computing the effective diffusivity, D eff ,of the cells at finite temperature, the phase of the tis-sue can be measured at zero temperature, by analyzingits response to a non-linear deformation, such as a T1transformation. In other words, if there is a non-zeroenergy barrier to undergoing a T1, the system is solid-like [S12, S13]. Although the trends of increasing shapeand alignment with increasing fluidity remain intact, wefind that in the irregular tissue regime , the quantitativemeasures of mean cell shape and cell-cell alignment aresensitive to temperature. Moreover, D eff is sensitive tofinite-size effects, while the calculation of T1 energy bar-riers is not (see Sec. J). Therefore, in an attempt todecouple the effects of temperature and system size fromtissue fluidity, we compute the energy barriers for cells toundergo T1s at zero temperature, and compare the re-sults to our finite-temperature measure of fluidity, D eff . To measure a T1 barrier, we minimize the total energyof the system using the FIRE algorithm (see sec. D 1),choose a cell edge at random, and then incrementallyshorten the edge, minimizing the system energy againafter each shortening step, and stop when the edge lengthequals 0.006. The difference between the final and initialenergies is taken to be an estimate of the barrier height toexecuting that T1. Averaging over many edges gives usa mean barrier height for a system with given parametervalues, such as p of the tissue and l of the surroundingspring network.We find that for a floppy ECM, the energy barrier isapproximately 1 × − for a tissue with p = 3 .
95, andincreases as p increases, to approximately 1 × − for p = 3 .
71, as seen previously for bulk tissue [S13]. As theECM stiffens, however, the energy barriers for all p val-ues approach 1 × − , indicating a solidification of thetissue. This is shown in Fig. S8 (a). For each l - p pair,we also measure the mean cell shape and cell-cell align-ment parameter, Q , and find that high cell shapes andlow alignment correspond to low energy barriers (fluid-like tissue), whereas high cell shapes and high alignmentcorrespond to higher barriers (more solid-like tissue), asshown in Fig. S8 (b). These results corroborate our find-ings for D eff reported in the main text. J. Finite-size effects
To understand the system-size dependence of our re-sults, we first compute the mean squared difference be-tween the average cell shapes for our spheroid embed-ded in ECM and the average cell shapes of a pure vertexmodel, for floppy network and varying interfacial tension, γ , as a function of total system size. In other words, westudy how the mean squared difference between the redand blue curves in Fig. 2 (d) of the main text dependson system size, for fixed ratio of tissue area to box areaof 20%. The results are shown in Fig. S9. We find thatthe difference between the mean cell shape of a spheroidand that of an all-cell vertex model decreases as approx-imately 1 /N system .Throughout our work, in order to reduce the numberof free parameters in our model, we’ve maintained a con-stant ratio of N cells /N total = 0 .
2, where N total is thetotal number of polygons in the simulation box. How-ever, we expect our results to depend on this choice. Tounderstand the extent of this effect, we repeat the aboveanalysis, now for fixed N cells = 800 and varying ratios of N cells /N total . As shown in Fig. S10, the squared differ-ence of our model from the bulk vertex model does not infact depend significantly on the ratio of tissue size to boxsize. For a fixed tissue size of N cells = 800, the tissue cellbehavior is similar to the bulk vertex model, regardlessof the relative size of the box.We also find that D eff depends on system size. With0 ! " (a) - . (b) FIG. S8. T1 energy barrier as a function of spring rest length,mean cell shape, and cell-cell alignment. (a) The energy bar-rier, E barrier , for a tissue cell edge to undergo a T1 transi-tion, as a function of the rest length, l , of the springs inthe ECM. For l > .
62, the ECM is floppy, and E barrier de-creases from approximately 1 × − for p = 3 .
71 to 1 × − for p = 3 .
95. For approximately l < .
62, the ECM isrigid, and E barrier approaches 1 × − for all p values. (b)A heatmap of E barrier as a function of the measured meanshape, h s i , and cell-cell alignment, Q , of the tissue cells. Thisdiagram is made by binning the results of simulations acrossranges of p and l . N total ( s b u l k s v m ) FIG. S9. The squared difference of the mean cell shapes,between spheroid and all-cell simulations, for varying systemsizes. In each spheroid simulation, the tissue takes up 20% ofthe total system size. The dashed line has a slope of -1 onthis log-log plot, illustrating the roughly N − behavior. N cells / N total ( s b u l k s v m ) FIG. S10. The squared difference of the mean cell shapes, be-tween spheroid and all-cell simulations, for varying box sizes,for fixed N cells = 800.FIG. S11. Effective diffusion coefficient, D eff , vs. system sizefor a tissue of p = 3 . D eff decreases, although all othersystem parameters are fixed. a finite boundary, D eff decreases as the system size isdecreased, when surrounded by floppy network, as showin Fig. S11.As we increase the network tension, cavities are formedand the tissue spreads into “channels” that surroundthese cavities. To avoid conflating the finite-size effecton D eff with the potential rigidifying effect of the sur-rounding network tension, we also quantify the phase ofthe tissue using T1 energy barriers (see Sec. I), which donot suffer from finite-size effects, as show in Fig. S12. Wefind that as l decreases and the tension in the networkincreases, the energy barriers increase, supporting the ar-gument that the simultaneous decrease in D eff is due tothe interaction with the rigid ECM and not a finite-sizeeffect.1 FIG. S12. The energy barrier, E barrier , for a cell edge un-dergoing a T1 vs. system size for a tissue of p = 3 .
9, forboth floppy ( l = 0 .
65) and rigid ( l = 0 .
58) surroundingECM. For a rigid network, E barrier is about two orders ofmagnitude greater than for a floppy network, regardless ofthe system size.FIG. S13. Effective diffusion coefficient, D eff , as a functionof simulation temperature and p , for a “balanced” pair of l and γ values. When γ and l are balanced, temperatureeffects the phase of the tissue, as a function of p , as it doesin the bulk vertex model. K. Temperature dependence
Increasing the simulation temperature, T , increasesthe magnitude of fluctuations, adding forces of increasingstrength and random direction to each vertex (see Sec.D 2). In the bulk vertex model, adding thermal energycan enable cells to overcome T1 energy barriers, even for p < .
81, and induce high cell shapes, as thermal forcesperturb the positions of the cell vertices. In our model,we expect to see a similar trend in the cell shape as afunction of temperature. However, these thermal forcesare now also competing with the forces due to interfacialtension and ECM tension. In order to understand thiscompetition, we first check that temperature is playinga fundamentally similar role in our model. We identifya pair of interfacial tension, γ , and equilibrium spring length, l , values for which the tissue is neither expandedor compressed. In other words, for a chosen l , we iden-tify the γ value for which the pressures generated by thesurface and by the external network are balanced, mean-ing that the tissue is compact and the mean cell shape iswhat we would expect for bulk tissue for the current p .Fig. S13 demonstrates the dependence of the diffusivity, D eff , on T and p , for a balanced pair of γ and l of0.25 and 0.6, respectively. The behavior at T = 1 e − γ and l , as there we see a transition in D eff around p ≈ .
81, matching the bulk vertex model.We find that as we increase T , the transition point in p decreases. This indicates that, like the bulk vertexmodel, thermal fluctuations effectively decrease the en-ergy barriers for cell rearrangements, fluidizing the tissue“sooner” (at lower values of p ).Having confirmed the role of temperature alone on thetissue, we now analyze the competition between increas-ing temperature and increasing ECM tension. At T = 0,as l decreases and ECM tension increases, we see themean cell shape, h s i , and cell-cell alignment, Q , bothincrease. Simultaneously increasing T in this case hastwo competing effects: first, as in bulk tissue and whenthe ECM tension is balanced by interfacial tension, largefluctuations drive tissue cell irregularity, which increases h s i and D eff . Second, though, fluctuations mitigate theeffect of ECM tension, therefore reducing the increase in h s i and Q driven by external tension. These results aredemonstrated in Fig. S14, (b) and (c), for fixed p = 3 . l > l ∗ ≈ .
62, the ECM is floppy, and increasingthe temperature results in an increase in h s i and littlechange in Q . However, for l < l ∗ , increasing T attenu-ates the increase in h s i and decreases Q . In other words,at higher temperatures, higher ECM tensions are neededto observe the same h s i and Q . Despite the mitigating ef-fect of temperature on h s i and Q , the effective diffusivity, D eff , is not strongly effected by increasing temperature,as shown in Fig. S14 (a), at least for the temperaturerange studied.2 (a) (b) (c) FIG. S14. The temperature-dependence of diffusivity, D eff , mean cell shape, h s i , and cell-cell alignment, Q . Here, γ = 0 and l varies from values corresponding to rigid ECM (purple to yellow curves) to floppy network (orange to red curves). For a fixedtemperature, increasing the rigidity (decreasing l ) decreases D eff and increases both h s i and Q . The differences in h s i and Q across l values are less pronounced at higher temperatures, while the effect on D eff is relatively constant. While increasingthermal fluctuations mitigates the influence of external network tension on the tissue, it also promotes higher cell shapes andhigher diffusion. [S1] R. Farhadifar, J.-C. Rper, B. Aigouy, S. Eaton, andF. Jlicher, Curr. Biol. (2007).[S2] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning,Nat. Phys. (2015).[S3] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Phys. Rev. X (2016).[S4] D. M. Sussman and M. Merkel, Soft Matter (2018).[S5] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan, M. L. Manning, and K. E. Kasza, PNAS (2020).[S6] J. L. Shivers, S. Arzash, A. Sharma, and F. C. MacK-intosh, Phys. Rev. Lett. , 188003 (2019).[S7] J. C. Maxwell, The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science (1864).[S8] A. Sharma, A. J. Licup, K. A. Jansen, R. Rens,M. Sheinman, G. H. Koenderink, and F. C. MacK-intosh, Nat. Phys. (2016).[S9] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L.Manning, PNAS (2019).[S10] D. M. Sussman, Comput. Phys. Commun. (2017).[S11] E. Bitzek, P. Koskinen, F. G¨ahler, M. Moseler, andP. Gumbsch, Phys. Rev. Lett. (2006).[S12] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning,Soft Matter (2014).[S13] P. Sahu, D. M. Sussman, M. Rbsam, A. F. Mertz,V. Horsley, E. R. Dufresne, C. M. Niessen, M. C.Marchetti, M. L. Manning, and J. M. Schwarz, SoftMatter16