Motility-driven glass and jamming transitions in biological tissues

Cell motion inside dense tissues governs many biological processes, including embryonic development and cancer metastasis, and recent experiments suggest that these tissues exhibit collective glassy behavior. To make quantitative predictions about glass transitions in tissues, we study a self-propelled Voronoi (SPV) model that simultaneously captures polarized cell motility and multi-body cell-cell interactions in a confluent tissue, where there are no gaps between cells. We demonstrate that the model exhibits a jamming transition from a solid-like state to a fluid-like state that is controlled by three parameters: the single-cell motile speed, the persistence time of single-cell tracks, and a target shape index that characterizes the competition between cell-cell adhesion and cortical tension. In contrast to traditional particulate glasses, we are able to identify an experimentally accessible structural order parameter that specifies the entire jamming surface as a function of model parameters. We demonstrate that a continuum Soft Glassy Rheology model precisely captures this transition in the limit of small persistence times, and explain how it fails in the limit of large persistence times. These results provide a framework for understanding the collective solid-to-liquid transitions that have been observed in embryonic development and cancer progression, which may be associated with Epithelial-to-Mesenchymal transition in these tissues.

Recent experiments have revealed that cells in dense biological tissues exhibit many of the signatures of glassy materials, including caging, dynamical heterogeneities and viscoelastic behavior [2,3,36,43,44]. These dense tissues, where cells are touching one another with minimal spaces in between, are found in diverse biological processes including wound healing, embryonic development, and cancer metastasis.
In many of these processes, tissues undergo an Epithelialto-Mesenchymal Transition (EMT), where cells in a solidlike, well-ordered epithelial layer transition to a mesenchymal, migratory phenotype with less well-ordered cell-cell interactions [52,53], or an inverse process, the Mesenchymalto-Epithelial Transition (MET). Over many decades, detailed cell biology research has uncovered many of the signaling pathways involved in these transitions [17,34], which are important in developing treatments for cancer and congenital disease.
Most previous work on EMT/MET has focused, however, on properties and expression levels in single cells or pairs of cells, leaving open the interesting question of whether there is a collective aspect to these transitions: Are some features of EMT/MET generated by large numbers of interacting cells? Although there is no definitive answer to this question, several recent works have suggested that EMT might coincide with a collective solid-to-liquid jamming transition in biological tissues [18,36,38,42]. Therefore, our goal is to develop a framework for jamming and glass transitions in a minimal model that accounts for both cell shapes and cell motility, in order to make predictions that can quantitatively test this conjecture.
Jamming occurs in non-biological particulate systems (such as granular materials, polymers, colloidal suspensions, and foams) when their packing density is increased above some critical threshold, and glass transitions occur when the fluid is cooled below a critical temperature. Over the past 20 years these phenomena have been unified by "jamming phase diagrams" [29,55].
Building on these successes, researchers have recently used self-propelled particle (SPP) models to describe dense biological tissues [4,15,20,45,48,51]. These models are similar to those for inert particulate matter -cells are represented as disks or spheres that interact with an isotropic soft repulsive potential -but unlike Brownian particles in a thermal bath, self-propelled particles exhibit persistent random walks.
SPP models typically exhibit a glass transition from a diffusive fluid state to an arrested sub-diffusive solid that is controlled by (1) the strength of self-propulsion [15,20,35] and (2) the packing density φ [6,13,14,20,35]. Just like in thermal systems, a jamming transition occurs at a critical packing density φ G , but this critical density is altered by the persistence time of the random walks [6,13,14,20,35].
During many biological processes, however, a tissue remains at confluence (packing fraction equal to unity) while it changes from a liquid-like to a solid-like state or vice-versa. For example, in would healing, cells collectively organize to form a 'moving sheet' without any change in their packing density [26], and during vertebrate embryogenesis mesendoderm tissues are more fluid-like than ectoderm tissues, despite both having packing fraction equal to unity [44].
Recently, Bi and coworkers [7] have demonstrated that the well-studied vertex model for 2-D confluent tissues [8,12,22,31,33,49] exhibits a rigidity transition in the limit of zero cell motility. Specifically, the rigidity of the tissue vanishes at a critical balance between cortical tension and cell-cell adhesion. An important insight is that this transition depends sensitively on cell shapes, which are well-defined in the vertex model. While promising, vertex models are difficult to compare to some aspects of experiments because they do not incorporate cell motility.
In this work, we bridge the gap between the confluent tissue mechanics and cell motility by studying a hybrid between the vertex model and the SPP model, that we name Self-Propelled-Voronoi (SPV) model. A similar model was introduced by Li and Sun [28], and cellular Potts models also bridge this gap [24,50], although glass transitions have not been carefully studied in any of these hybrid systems.

I. THE SPV MODEL
While the vertex model describes a confluent tissue as a polygonal tiling of space where the degrees of freedom are the vertices of the polygons, the SPV model identifies each cell only using the center (r r r i ) of Voronoi cells in a Voronoi tessellation of space (Dirichlet domains) [27]. The observation that Voronoi tessellations can describe cellular patterns in epithelial tissues was first proposed by Honda [21]. For a tissue containing N-cells, the inter-cellular interactions are captured by a total energy which is the same as that in the vertex model. Since the tessellation is completely determined by the {r r r i }, the total tissue mechanical energy can be fully expressed as E = E({r r r i }): The term quadratic in cell area A(r i ) results from a combination of cell volume incompressibility and the monolayer's resistance to height fluctuations [22]. The term involving the cell perimeter P(r i ) originates from active contractility of the acto-myosin sub-cellular cortex (quadratic in perimeter) and effective cell membrane tension due to cell-cell adhesion and cortical tension (both linear in perimeter). This gives rise to an effective target shape index that is dimensionless: p 0 = P 0 / √ A 0 . K A and K P are the area and perimeter moduli, respectively. For the remainder of this manuscript we assume p 0 is homogenous across a tissue, although heterogeneous properties are also interesting to consider [59].
In the vertex model [7], a rigidity transition takes place at a critical value of p 0 = p * 0 ≈ 3.81. When p 0 < p * 0 , cortical tension dominates over cell-cell adhesion and the energy barriers for local cell rearrangement and motion are finite. The tissue then behaves as a elastic solid with finite shear modulus. When p 0 > p * 0 , cell-cell adhesion dominates and the energy barriers for local rearrangements vanish, resulting in zero rigidity and fluid-like behavior. While the energy functional for cell-cell interactions is identical in the vertex and SPV models, the two are truly distinct: the local minimum energy states of the vertex model are not guaranteed to be similar to a Voronoi tessellation of cell centers, although we do find them to be very similar in practice. Therefore, we are also interested in whether a rigidity transition in the SPV model coincides with the rigidity transition of the vertex model.
We define the effective mechanical interaction force experienced by cell i as F F F i = −∇ ∇ ∇ i E (see Appendix A for details). In contrast to particle-based models, F F F i is non-local and non-additive: F F F i cannot be expressed as a sum of pairwise force between cells i and its neighboring cells. Nevertheless, one can show that momentum is still precisely conserved by this energy functional in the absence of the additional selfpropulsion forces introduced below.
In addition to F F F i , cells can also move due to self-propelled motility. Just as in SPP models, we assign a polarity vector n n n i = (cos θ i , sin θ i ) to each cell; alongn n n i the cell exerts a selfpropulsion force with constant magnitude v 0 /µ, where µ is the mobility (the inverse of a frictional drag). Together these forces control the over-damped equation of motion of the cell centers r r r i dr r r i dt = µF F F i + v 0n n n i .
The polarity is a minimal representation of the front/rear characterization of a motile cell [50]. While the precise mechanism for polarization in cell motility is an area of intense study, here we model its dynamics as a unit vector that undergoes random rotational diffusion where θ i is the polarity angle that definesn n n i , and η i (t) is a white noise process with zero mean and variance 2D r . The value of angular noise D r determines the memory of stochastic noise in the system, giving rise to a persistence time scale τ = 1/D r for the polarization vectorn. For small D r 1, the dynamics ofn is more persistent than the dynamics of the cell position. At large values of D r , i.e. when 1/D r becomes the shortest timescale in the model, Eq. (2) approaches simple Brownian motion.
The model can be non-dimensionalized by expressing all lengths in units of √ A 0 and time in units of 1/(µK A A 0 ). There are three remaining independent model parameters: the self propulsion speed v 0 , the cell shape index p 0 , and the rotational noise strength D r . We simulate a confluent tissue under periodic boundary conditions with constant number of N = 400 cells (no cell divisions or apoptosis) and assume that the average cell area coincides with the preferred cell area, i.e. A i = A 0 . This approximates a large confluent tissue in the absence of strong confinement. We numerically simulate the model using molecular dynamics by performing 10 5 integration steps at step size ∆t = 10 −1 using Euler's method. A detailed description of the SPV implementation can be found in the Appendix Sec. A.

II. CHARACTERIZING GLASSY BEHAVIOR
We first characterize the dynamics of cell motion within the tissue by analyzing the mean-squared displacement (MSD) of cell trajectories. In Fig. 1(a), we plot the MSD as function of time, for tissues at various values of p 0 and fixed v 0 = 0.1 and D r = 1. The MSD exhibits ballistic motion ( slope close to 2 on a log-log plot) at short times, and plateaus at intermediate timescales. The plateau is an indication that cells are becoming caged by their neighbors. For large values of p 0 , the MSD eventually becomes diffusive (slope =1), but as p 0 is decreased, the plateau persists for increasingly longer times. This indicates dynamical arrest due to caging effects and broken ergodicity, which is a characteristic signature of glassy dynamics.
Another standard method for quantifying glassy dynamics is the self-intermediate scattering function [56]: F s (k,t) = e i k·∆ r(t) . Glassy systems possess a broad range of relaxation timescales, which show up as a long plateau in F s (t) when it is analyzed at a lengthscale q comparable to the nearest neighbor distance. Fig 1 (b) illustrates precisely this behavior in the SPV model, when | k| = π/r 0 , where r 0 is the position of the first peak in the pair correlation function. The average ... is taken temporally as well as over angles of k. F s (t) also clearly indicates that there is a glass transition as a function of p 0 : at high p 0 values F s approaches zero at long times, indicating that the structure is changing and the tissue behaves as a viscoelastic liquid. At lower values of p 0 , F s remains large at all timescales, indicating that the structure is arrested and the tissue is a glassy solid. Fig 1 (d) demonstrates that at the structural relaxation time, the cell displacements show collective behavior across large lengthscales suggesting strong dynamical heterogeneity. This is strongly reminiscent of the 'swirl' like collective motion seen in experiment in epithelial monolayers [2,3,15,39,40].

A. A dynamical order parameter for the glass transition
Although the phase space for this model is three dimensional, we now study the model at a fixed value of D r = 1.
We then search for a dynamical order parameter that distinguishes between the glassy and fluid states as a function of the two remaining model parameters,(v 0 , p 0 ). A candidate order parameter is the self-diffusivity D s : D s = lim t→∞ ∆r(t) 2 /(4t). For practicality, we calculate D s using simulation runs of 10 5 time steps, chosen to be much longer than the typical caging timescale in the fluid regime. We present the self-diffusivity in units of D 0 = v 2 0 /(2D r ), which is the free diffusion constant of an isolated cell. D e f f = D s /D 0 then serves as an accurate dynamical order parameter that distinguishes a fluid state from a solid (glassy) state in the space of (v 0 , p 0 ), matching the regimes identified using the MSD and F q . In Fig. 2, the fluid region is characterized by a finite value of D e f f and D e f f drops below a noise floor of ∼ 10 −3 as the glass transition is approached. In practice, we label materials with D e f f > 10 −3 as fluids indicated by an orange dot, and those with D e f f ≤ 10 −3 as solids indicated by blue squares. Importantly, we find that the SPV model in the limit of zero cell motility shares a rigidity transition with the vertex model [7] at p 0 ≈ 3.81, and that this rigidity transition controls a line of glass transitions at finite cell motilities. Typical cell tracks (Fig. 2) clearly show caging behavior in the glassy solid phase.

B. Cell shape is a structural order parameter for the glass transition
In glassy systems it can be difficult to experimentally distinguish between a truly dynamically arrested state and a state with relaxation times longer than the experimental time window. Similarly, in tissues it is experimentally challenging to quantify a glass transition through the measurement of a dynamical quantity such as the diffusivity D s . Identifying a static quantity that directly probes the mechanical properties of a tissue would therefore be a powerful tool for experiments. Puliafito et al. have suggested that shape changes accompany dynamical arrest in proliferating tissues [41]. Similarly, a structural signature based on cell shapes -the shape index q = p/ √ a -was previously shown to be an excellent order parameter for the confluent tissue rigidity transition in the vertex model [38]. In a model where cells were not motile (v 0 = 0) we found that when p 0 < 3.813, q is constant ∼ 3.81 and when p 0 > 3.81 q grows linearly with p 0 . Quite surprisingly, we found that the prediction of q = 3.813 works perfectly in identifying a jamming transition in in-vitro experiments involving primary human tissues, where cells are clearly motile (v 0 = 0) [38]. At that time, we did not understand why the v 0 = 0 theory worked so well for these tissues.
The prediction of a solid-liquid transition in the SPV model presented here provides an explanation for this observation.We find that q (which can be easily calculated in experiments or simulations from a snapshot) can be used as a structural order parameter for the glass transition for all values of v 0 , not just at v 0 = 0. Specifically, the boundary defined by q = 3.813, shown by the blue dashed line in Fig. 2(A) co-incides extremely well with the glass transition line obtained using the dynamical order parameter, shown by the round and square data points. The insets to Fig. 2 also illustrate typical cell shapes: cells are isotropic on average in the solid phase and anisotropic in the fluid phase. This highlights the fact that q can be used as a structural order parameter for the glass transition at all cell motilities, providing a powerful new tool for analyzing tissue mechanics.

III. A THREE-DIMENSIONAL JAMMING PHASE DIAGRAM FOR TISSUES
Having studied the glass transition as function of v 0 and p 0 at a large value of D r , we next investigate the full threedimensional phase diagram by characterizing the effect of D r on tissue mechanics and structure. D r controls the persistence time τ = 1/D r and persistence length or Péclet number Pe ∼ v 0 /D r of cell trajectories; smaller values of D r correspond to more persistent motion.
In Fig. 3(A), we show several 2D slices of the three dimensional jamming boundary. Solid lines illustrate the phase transition line identified by the structural order parameter q = 3.813 as function of v 0 and p 0 for a large range of D r values (from 10 −2 to 10 3 ). (In Appendix B 2 we demonstrate that the structural transition line q = 3.813 matches the dynamical transition line for all studied values of D r .) In contrast to results for particulate matter [6], this figure illustrates that the glass transition lines meet at a single point (p 0 = 3.81) in the limit of vanishing cell motility, regardless of persistence. Fig. 3(B) shows an orthogonal set of slices of the jamming These slices can be combined to generate a threedimensional jamming phase diagram for confluent biological tissues, shown in Fig. 3(C). This diagram provides a concrete, quantifiable prediction for how macroscopic tissue mechanics depends on single-cell properties such as motile force, persistence, and the interfacial tension generated by adhesion and cortical tension.
We note that Fig. 3(C) is significantly different from the jamming phase diagram conjectured by Sadati et al [42], which was informed by results from adhesive particulate matter [55]. For example, in particulate matter adhesion enhances solidification, while in confluent models adhesion increases cell perimeters/surface area and enhances fluidization. In addition, we identify "persistence" as a new axis with a potentially significant impact on cell migration rates in dense tissues.
To better understand why persistence is so important in dense tissues, we first have to characterize the transitions between different cellular structures. In the limit of zero cell motility, the system can be described by a potential energy landscape where each allowable arrangement of cell neighbors corresponds to a metastable minimum in the landscape. There are many possible pathways out of each metastable state: some of correspond to localized cell rearrangements, while others correspond to large-scale collective modes. The maximum energy required to transition out of a metastable state along each pathway is called an energy barrier [8].
We observe that tissue fluidity can increase drastically with increasing D r at finite cell speeds. This suggests that different pathways (with lower energy barriers) must become dynamically accessible at higher values of D r .
One hint about these pathways comes from the instantaneous cell displacements, shown for different values of D r in Fig. 4. At high values of D r , (p 0 = 3.78, v 0 = 0.1) the instantaneous displacement field is essentially random and largely uncorrelated, as shown in Fig. 4, and the material is solidlike. There is no collective behavior among cells, and each cell 'rattles' independently near its equilibrium position.
However, as D r is lowered, the instantaneous displacement field becomes much more collective (Fig. 4) and the tissue begins to flow, presumably because these collective displacement fields correspond to pathways with lower transition energies.
Two obvious questions remain: How does a lower value of D r generate more collective instantaneous displacements? Why should collective instantaneous displacements generi-cally have lower energy barriers? The first question can be answered by extending ideas first proposed by Henkes, Fily and Marchetti [20] to explain why motion in self-propelled particle models seems to follow the 'soft modes' of a solid. This argument is based on a simple, yet powerful observation: in the limit of zero motility (v 0 = 0), a solid-like state will have a well-defined set of normal modes of vibration (with frequencies {ω ν }), and a corresponding set of eigenvectors ({ê ν }) that forms a complete basis. At higher motilities (v 0 > 0) near the glass transition, the motion of particles in the system can be expanded in terms of the eigenvectors. As discussed in Appendix B 1, one can use this observation to show that in the limit of D r → 0, motion along the lowest frequency eigenmodes is amplified -the amplitude along each mode is proportional to 1/ω 2 ν . These low-frequency normal modes are precisely the collective displacements observed for low D r .
The second question is more difficult to answer because it is impossible to enumerate all of the possible transition pathways and energy barriers in a disordered material. However, a partial answer comes from recent work in disordered particulate matter showing that low-frequency normal modes do have significantly lower energy barriers [30,61] than higher frequency normal modes. A similar analysis could potentially be performed in vertex or SPV models.

IV. A CONTINUUM MODEL FOR GLASS TRANSITIONS IN TISSUES
Although continuum hydrodynamic equations of motion have been developed by coarse graining SPP models in the dilute limit, there is no existing continuum model for a dense active matter system near a glass transition. Here we propose that a simple trap [32] or Soft Glassy Rheology (SGR) [47] model provides an excellent continuum approximation for the phase behavior in the large D r Brownian regime, but fails in the small D r limit.
For large D r it is known that particle behave like Brownian particles with an effective temperature T e f f = v 2 0 /2µD r [14]. This mapping becomes exact when D r → ∞ at fixed "effective inertia" (µD r ) −1 [13]. In other words, like in granular systems [1,9], the effective temperature in SPP is dominated by kinetic effects. Guided by this result we conjecture that in our model the temperature also scale quadratically with the velocity, Physically, this effective temperature gives the amount of energy available for individual cells to vibrate within their cage or 'trap'. The next important question is how to characterize the 'trap depths', or energy barriers between metastable states. In the Brownian regime (large D r ) there is no dynamical mechanism for the cells to organize collectively, and therefore a reasonable assumption is that the rearrangements are small and localized.
In [8], some of us explicitly calculated the statistics of energy barriers for localized rearrangements in the equilibrium Comparison between SPV glass transition and an analytic prediction based on a Soft Glass Rheology (SGR) continuum model. The dashed line corresponds to an SGR prediction with no fit parameter based on previously measured vertex model trap depths [8]. Data points correspond to SPV simulations with D r = 10 −3 and where we have defined T e f f = cv 2 0 with c = 0.1 as the best-fit normalization parameter. Blue points correspond to simulations which are solid-like, with D e f f < 10 −3 , and the boundary of these points define the observed SPV glass transition line. (Inset) L 2 difference between SPV glass transition line (at best-fit value of c) and the predicted SGR transition line at various values of D r . The SGR prediction based on localized T1 trap depths works well in the high D r limit, but not in the low D r limit. vertex model. In the 2D vertex model, one can show that localized rearrangements must occur via so-called T1 transitions [58]. Using a trap model [32] or Soft Glassy Rheology(SGR) [47] framework, we were able to use these statistics to generate an analytic prediction, with no fit parameters, for the glass transition temperature T g as function of p 0 .
To see if the SGR prediction for the glass transition holds for the SPV model in the large D r limit, we simply overlay the data points corresponding to glassy states from the SPV model with the glass transition T g line predicted in [8]. There is one fitting parameter c that characterizes the proportionality constant in Eq. 4. Fig 5, shows that the SPV data for D r = 10 3 is in excellent agreement with our previous SGR prediction.
Because T e f f ∼ v 2 0 , and the glass transition line scales as T g ∼ p * 0 − p 0 for p 0 p * 0 and D r → 0, the glass transition line scales as v 0 ∼ (p * 0 − p 0 ) 0.5 in those limits. The reason the effective temperature SGR model works here is that, like in SPP models of spherical active Brownian colloids, the angular dynamics of each cell evolves independently of cell-cell interactions and of the angular dynamics of other cells. An additional alignment interaction that couples the angular and translational dynamics may therefore modify this behavior.
To our knowledge, this is the first time that a SGR / trap model prediction has been precisely tested in any glassy system. This is because, unlike most glass models, we can enumerate all of the trap depths for localized transition paths in the vertex model.
However, for small values of D r , we have shown that cell displacements are dominated by collective normal modes, and therefore the energy barriers for localized T1 transitions are probably irrelevant in this regime. The inset to Fig 5 shows the deviation (L 2 -norm) between glass transition lines in the SPV model and T1-based SGR prediction as a function of D r . We see that the SGR prediction fails in the small D r limit, as expected. A better understanding of the energy barriers associated with collective modes will be required to modify the theory at small D r .

V. DISCUSSION AND CONCLUSIONS
We have shown that a minimal model for confluent tissues with cell motility exhibits glassy dynamics at constant density. This model allows us to make a quantitative prediction for how the fluid-to-solid/jamming transition in biological tissues depends on parameters such as the cell motile speed, the persistence time associated with directed cell motion, and the mechanical properties of the cell (governed by adhesion and cortical tension). We define a simple, experimentally accessible structural order parameter -the cell shape index -that robustly identifies the jamming transition, and we show that a simple analytic model based on localized T1 rearrangements precisely predicts the jamming transition in the large D r limit. We also show that this prediction fails in the small D r limit, because the instantaneous particle displacements are dominated by collective normal modes.
This model makes several experimentally-verifiable predictions for cell shape and tissue mechanics: • The order parameter q = 3.81 is a structural signature for the glass transition, even in tissues with significant cell motility or dynamics. This prediction has already been tested in epithelial lung tissue [38], but it should be much more broadly applicable. We have performed a rudimentary shape analysis of a small number of images from other systems that have been previously published, including proliferating MCDK monolayers [10] and convergent extension in fruit fly development [11] and found that the shapes are consistent with this prediction. A much more careful analysis with full data sets should be performed to further validate this prediction or understand where it breaks down.
• In the limit of vanishing cell motility, shape and pressure fluctuations should vanish when the jamming transition is approached from the solid side, and remain zero in the fluid. A finite motility v 0 will induce such fluctuations in the fluid phase, as confirmed by preliminary calculation of cellular stresses and pressure in the SPV model [62]. This could be studied by combining measurement of cell shape fluctuations with traction force microscopy (TFM) in wound healing assays. After locating the glass transition by imaging cell shape changes, it may be possible to extract information on cell motility v 0 from cellular stresses and pressure inferred from TFM in the fluid phase near the glass tran-sition. This suggests that one may estimate cell motility by examining the changes in cellular stresses and pressure in the cell monolayer near the unjamming transition and assuming that the local velocity of the monolayer is very small just above the transition. The latter assumption can also be verified independently via particle image velocimetry (PIV).
• Cell proliferation, so far neglected in our model, causes an increase in cell number density in confluent tissues. Often this is accompanied by a reduction in individual cell motility v 0 , via contact inhibition of locomotion. In cases where this is the dominant effect and changes to the ratio between A 0 and P 0 are negligible, our work predicts that proliferation would drive the system towards jamming. This is consistent with existing reports in the literature [10], although more work is required to test the prediction carefully. In tissues where v 0 remains low at all times [41], our model predicts that proliferation can either cause jamming or unjamming, depending on whether cell divisions are oriented in such a way to decrease or increase cell shape anisotropy.
• Spatial correlations and fluctuations of the cell displacement field, such as swirl sizes [2,3,15,39,40] , should grow as a tissue approaches the glass transition from the fluid side. Very recently, a similar prediction for displacements and correlation lengths based on a particle-based model has been verified in one cell type [15]. The SPV model, which makes predictions for cell shapes in addition to displacements and correlation lengths, could be tested simply by compiling detailed statistics about cell shapes and cell motion in epithelial monolayers.
Although all of the work presented here focuses on the SPV model, which tracks cell centers and therefore has only two degrees of freedom per cell, we found that in the limit of zero cell motility it exhibits the same rigidity transition as the vertex model which has two degrees of freedom per vertex. We have also checked that an "active vertex model", where active motile forces are added to the vertex model vertices, also exhibits a robust glass transition characterized by the shape order parameter q. The fact that two models with ostensibly different degrees of freedom share the same transition suggests that there is a deeper universality, perhaps generated by isostaticity, that remains to be understood.
Another result of this work is the surprising and unexpected differences between confluent models (such as the vertex and SPV models) and particle-based models (such as Lennard-Jones glasses and SPP models). For example, works by Berthier [6], Fily and Marchetti [14] in SPP models suggest that the location of the zero motility glass transition packing density φ G (defined as the density at which dynamics cease in the limit of v 0 → 0) depends the value of noise, D r . This is also related to the observation that the jamming and glass transition are not controlled by the same critical point in nonactive systems [23,37]. We find this is not the case in the SPV model. Fig. 3(a & b) show that while the glass transition point p * 0 shifts with D r at finite values of v 0 , in the limit of vanishing motility, all glass transition lines merge on to a single point in the limit v 0 → 0, namely p * 0 = 3.81. Given these differences, it is important to ask which type of model is appropriate for a given system. We argue that SPV models are maybe more appropriate for many biological tissues. Whereas SPP models interact with two-body interactions that only depend on particle center positions, both SPV and vertex models naturally incorporate contractility as a key property of living cells and capture the inherently multibody nature of intercellular forces due to shape deformations. Unlike equilibrium vertex models, SPV models account for cell motility, and they are also much easier to simulate in 3D (which is nearly impossible in practice for the vertex model.) Recent work by Li and Sun [28] also models a confluent cell as a Voronoi tessellation of the plane. An important difference between this work and ours is that in Ref. [28] cellcell adhesion is captured via a potential that is quadratic in the distance between cell centers, just as in particle models. We might guess that stronger cell-cell adhesion in their model will result in stiffening of the tissue, which is common for particle based models, although that remains to be tested in active systems. In contrast, adhesion enters our model through the coupling of the shape energy to the cell perimeter. Increasing cell-cell adhesion (or decreasing cortical tension) yields a larger value of p 0 , which leads to the tissue becoming softer.
We expect that other shape-based models of confluent tissue dynamics will also yield the glass transition described here. For example, it has been reported in recent works based on the Cellular Potts model [24,50] and in a modified SPP model [15] when the cell motile force is decreased beyond a certain threshold, the motion of cells transitions from diffusive to sub-diffusive. This is similar to crossing the glass transition line in the SPV model by decreasing the value of v 0 .
In this work and in previous work based on the vertex model, the cell volume is generally assumed to be fixed. While this is a good assumption in developmental systems such as drosopholia [16,25] and zebrafish [31], epithelial tissues cells can show significant volume fluctuations, as reported recently [63,64]. Therefore, it will be important to incorporate volume fluctuations in future iterations of the vertex model or the SPV model, as they introduce another source of active shape fluctuation and could therefore lead to jamming or unjamming of the tissue locally and potentially shift the location of the rigidity and glass transitions.
In our version of the SPV model, we have assumed that cell polarity is controlled by simple rotational white noise. It is also possible to include more complex mechanisms. For example, external chemical or mechanical cues could be modeled by coupling v 0 andn n n i to chemoattractant or mechanical gradients, allowing waves or other pattern formation mechanisms to interact with the jamming transition. Similarly, simple alignment rules (such as those in the Viscek model [57]) could lead to collective flocking modes that also affect glassy dynamics.
Another interesting extension of the SPV model would be to study the role of cell-cell friction, which has already been shown to be important in controlling collective dynamics in particle-based tissue models [15]. Our current model includes viscous frictional coupling of cell to the 2D substrate and cellcell adhesion enters as a negative line tension on interfaces. However, it would be possible to add a frictional force between cells proportional to the length of the edge shared between two cells, and we know from previous work on particulate glasses that these localized frictions can change the location of jamming/glass transition and the nature of spatial correlations in a glass [19,46].
It is also tempting to speculate about the relationship between the unjamming transition captured by our model and the epithelial-mesenchimal transition (EMT) that precedes cell escape from a solid tumor mass. The EMT involves significant changes in cell-cell adhesion and cytoskeletal composition, with associated changes in cell shape and motility. This suggests that escape from the tumor mass is controlled not just by the chemical breakdown of the basement membrane, but also by specific changes in mechanical properties of both individual cells and the surrounding tissue [60]. One could then hypothesize that the collective unjamming described here may provide the first necessary step towards the mechanical changes needed for cell escape from primary tumors.
In particular, recent work suggests that cancer tumors are mechanically heterogeneous, with mixtures of stiff and soft cells that have varying degrees of active contractility [59]. Our jamming phase diagram suggests that the soft cells, which often exhibit mesenchymal markers and presumably correspond to higher values of p 0 , might unjam and move towards the boundary of a primary tumor more easily than their stiff counterparts. Examining the effects of tissue heterogeneity on tissue rigidity and patterns of cell motility is therefore a very promising avenue for developing predictive theories for tumor invasiveness and metastasis.

Appendix A: Simulation algorithm for the SPV model
To create an initial configuration for the simulation, we first generate a seed point pattern using random sequential addition (RSA) [54] and anneal it by integrating Eq. 2 with v 0 = 0 for 100 MD steps. The resulting structure then serves as an initially state for all simulations runs. The use of (RSA) only serves to speed up the initial seed generation as using a Poisson random point pattern does not change the results presented in this paper.
At each time step of the simulation, a Voronoi tesselation is created based on the cell centers. The intercellular forces are then calculated based on shapes and topologies of the Voronoi cells (see discussion below). We employ Euler's method to carry out the numerical integration of Eq. 2, i.e., at each time step of the simulation the intercellular forces is calculated based on the cell center positions in the previous time step.
In a Delaunay triangulation, a trio of neighboring Voronoi centers define a vertex of a Voronoi polygon. For example in Fig. 6, ( r i , r j , r k ) define the vertex h 3 , which is given by where the coefficients are given by In the vertex model, the total mechanical energy of a tissue depends only on the areas and perimeters of cells: In a Voronoi tessellation, the area and perimeter of a cell i can be calculated in terms of the vertex positions where z i is the number of vertices for cell i (also number of neighboring cells) and m indexes the vertices. We use the convention h z i = h 0 . With these definitions, the total force on cell-i can be calculated using Eq. A3 here µ denotes the cartesian coordinates (x, y). The first term on the r.h.s. of Eq. A5 sums over all nearest neighbors of cell i. It is the force on cell i due to changes in neighboring cell shapes. The second term is the force on cell i due to shape changes brought on by its own motion. It maybe tempting to treat ∂E j ∂r iµ as the force between cells-i and j, but since the interaction is inherently multi-cellular in nature and interactions between i and j also depend on k and l (see Fig. 6).
For the typical configuration shown in Fig. 6, the first term in Eq. A5 can be expanded using the chain rule and calculated using Eq. A1 In Eq. A7, only terms involving h 2 and h 3 are kept since E j does not depend on other vertices of cell i. ν is a cartesian coordinate index. The energy derivative in Eq. A7 can be calculated in a straightforward way, by using Eqs. A3 and A4 Similarly, the second term on the r.h.s. of Eq. A5 can be calculated in a similar way.
Appendix B: Cell displacements and structural order parameter as a function of D r

Expanding cell displacements in an eigenbasis associated with the underlying dynamical matrix
In the absence of activity (v 0 = 0), the tissue is a solid for p 0 < p * 0 = 3.81. As v 0 is increased, the solid behavior persists up to v 0 = v * 0 (p 0 ), which is given by the glass transition line in Fig. 2. In order for the tissue to flow, sufficient energy input is needed to overcome energy barriers in the potential energy landscape, which are a property of the underlaying solid state at v 0 = 0. In this limit, the instantaneous cell center positions { r i (t)} can be thought of as a small displacement { d i (t)} from the nearest solid reference state { r 0i } [20] where d i (t) = r i − r 0i . The r 0i correspond to positions of cell in a solid, which has a well-defined linear response regime [7]. The linear response is most conveniently expressed as the eigen-spectrum of the dynamical matrix D i jαβ . Since the eigenvectors {ê i,ν } of D i jαβ form a complete orthonormal basis, the cell center displacement can then be expressed as a linear combination of For simplicity, we will adopt the Bra-ket notation and express the eigenbasis simply as |ν and Eq. B1 becomes whereD and ω 2 ν are the eigenvalues of the dynamical matrix. The polarization vectorn i can also be expressed as a linear combination of eigenvectors Since the polarization vector and eigenvector are both unit vectors, it follows that b ν (t) = n|ν = cos(θ ν − ψ). Where ψ is the angle of the polarization and θ ν the angle of the eigenvector.
(B11) In the limit of D r → ∞ and Eq. B11 becomes a ν (t) = a ν (0)e −kt . (B12) This suggests that while normal modes control the rate of decay, they do no affect the long-time behavior. However as D r → 0, Eq. B11 becomes a ν (t) = a ν (0)e −µω 2 ν t + v 0 µω 2 ν cos (θ ν − ψ(0)) 1 − e −µω 2 ν t (B13) The second term in this equation scales as ∼ 1/ω 2 ν . Therefore, at short times (corresponding to instantaneous response), the mode amplitude a ν is much larger for modes at lower frequencies. Since the reference state is an elastic solid with Debye scaling D(ω) ∼ ω as ω → 0 [7], this suggests that the displacement will be heavily dominated by the lowest frequency modes that are spatially more collective in nature.