Tissue fluidization by cell-shape-controlled active stresses

Biological cells can actively tune their intracellular architecture according to their overall shape. Here we explore the rheological implication of such coupling in a minimal model of a dense cellular material where each cell exerts an active mechanical stress along its axis of elongation. Increasing the active stress amplitude leads to several transitions. An initially hexagonal crystal motif is first destabilized into a solid with anisotropic cells. Increasing activity further, we find a re-entrant transition to a regime with finite hexatic order and finite shear modulus, in which cells arrange according to a rhombile pattern with periodically arranged rosette structures. The shear modulus vanishes again at a third threshold beyond which spontaneous tissue flows arise. In this last regime, we observe the emergence of cell shape patterns called topological defects, with flow and stress fields around defects agreeing with those observed in epithelial tissue experiments. We further provide a testable prediction of cell-cell rearrangement hotspots near topological defects. Overall, our work connects seemingly distinct features - e.g. rosettes and topological defects - observed across various types of epithelial tissues.

Biological cells can actively tune their intracellular architecture according to their overall shape. Here we explore the rheological implication of such coupling in a minimal model of a dense cellular material where each cell exerts an active mechanical stress along its axis of elongation. Increasing the active stress amplitude leads to several transitions. An initially hexagonal crystal motif is first destabilized into a solid with anisotropic cells. Increasing activity further, we find a re-entrant transition to a regime with finite hexatic order and finite shear modulus, in which cells arrange according to a rhombile pattern with periodically arranged rosette structures. The shear modulus vanishes again at a third threshold beyond which spontaneous tissue flows arise. In this last regime, we observe the emergence of cell shape patterns called topological defects, with flow and stress fields around defects agreeing with those observed in epithelial tissue experiments. We further provide a testable prediction of cell-cell rearrangement hotspots near topological defects. Overall, our work connects seemingly distinct features -e.g. rosettes and topological defects -observed across various types of epithelial tissues.
Connecting the single-cell behavior to large-scale mechanical properties of biological tissues is key to understand development, regeneration, and disease [1]. Growing experimental evidence supports the idea that biological cells actively tune their intracellular architecture according to their overall shape. For instance, cortical actomyosin [2][3][4] and microtubules (whether in Drosophila [5,6] or in plants [4,7]) tend to align along the direction of cell shape elongation. Such oriented fibers are known to generate anisotropic stress [8].
Here, we explore the consequences of such cell-shape feedback on the tissue-scale behavior. We study a minimal model where a bulk cellular active stress σ (act) is created by filaments, which in turn align with cell shape, represented by a tensor Q. To lowest order, which we incorporate within a computational model for dense epithelilal tissue. For β > 0 (resp. β < 0), cells actively push (resp. pull) on their neighbours along their direction of elongation [9,10]. Relations like Eq. (1) have been considered in previous models. For instance, in active hydrodynamic theories, where both σ (act) and Q are defined by averages over several cells, Eq. (1) gives rise to a classical flow instability [11][12][13]. In the absence of confining boundaries, this instability occurs at arbitrarily small activities β > β c ≡ 0. An active stress as in Eq. (1) is also included in cellbased phase field simulations [14], where Q is defined at the single-cell level. In these simulations the transition to spontaneous flows occurs for a critical value of β c , which (i) is finite, β c > 0, and appears to be independent of the system size, contrasting with active hydrodynamic theory results, and (ii) scales with a cell surface tension. So far it is unclear why the critical activity β c is different in these two models, and why it scales with the surface tension in the phase field model Ref. [14].
Here we address these questions by combining analytical arguments and vertex model simulations. Vertex models describe epithelial tissues as networks of polygons [15]. Forces on the polygon vertices are defined by a mechanical energy: E = 1/2 n K A (A n − A 0 ) 2 + K P (P n − P 0 ) 2 where the sum is over all cells n of the tissue, A n and P n are cell area and perimeter, respectively. The parameters A 0 and P 0 are preferred cell area and perimeter with the associated rigidities K A and K P , respectively. A transition occurs at P 0 = P * 0 , with the tissue behaving as a yield stress solid for P 0 < P * 0 and as a fluid for P 0 ≥ P * 0 [16,17]. The numerical value of P * 0 is in the range 3.72 . . . 3.94, which a value that depends on the disorder in the cellular packing [16][17][18][19][20] (SM [21], Sec. I). Cell-based active polar forces at the individual cell level were also shown to drive a solid-to-fluid transition through an intermediate hexatic regime [22,23].
We introduce Eq. (1) in the vertex model framework, and we find that three transitions occur with increasing β. For β < β 1 the tissue is solid with isotropic cells. At β = β 1 , cells lose their isotropic shape. At a second transition β = β 2 , a regime with rhombile cell shapes, many-fold vertices, long-range crystalline order and finite shear modulus emerges. Finally, at β = β 3 , the tissue turns into an active fluid and starts to display persistent small-scale chaotic flow. We discuss the relation of these transitions to the P 0 -dependent vertex model transition between yield stress solid and fluid. This also allows us to demonstrate the appearance of a finite activity threshold β c in cell-based models such as Ref. [14]; the latter behave like yield stress solids, requiring finite β c for active flows to appear. Meanwhile, active hydrodynamic theories as in [11,12] describe tissues as fluids, such that active flows can appear without threshold. We then discuss potential applications to the modelling of biological tissues.   1. (a, b) Morphology of (a) a four-cell system and (b) a cell sheet at different activities β, with an initially hexagon (top) or Voronoi pattern (bottom). In (b), we mark vertices with coordination number Zi > 3.
(c) Cell elongation parameter q, Eq. (3), and T1 transition rate kT1 versus cell activity β, for an initially hexagonal (H) or Voronoi (V) pattern. The solid black line refers to an analytical approximation of q. (d) Hexatic order parameter ψ6 and average vertex coordination number Z versus cell activity β. (e) Long-time shear modulus Gxy and yield stress σ yield versus cell activity β.
Inset: red shaded area indicates the mean ± standard deviation of the shear modulus for the initially hexagonal pattern, estimated from n = 5 simulations. Default parameter values with P0 = 1.0.
Method We implement Eq. (1) through the following friction-based dynamics for the vertex positions r i Here, γ is a friction coefficient, F (svm) i = −∂E/∂r i are the standard gradient-based vertex model forces, and F (act) i are the active forces induced by the active anisotropic bulk stresses σ (act) defined for each cell according to Eq. (1). There are different ways to translate the cellular bulk stresses σ (act) into the vertex forces F (act) i [24][25][26]. Here we use the approach proposed by Tlili et al. [24,26], which relies on Cauchy's stress definition. For the cell shape anisotropy tensor Q in Eq. (1), we use the symmetric, traceless tensor where the sum is over all sides k of the cell, while l k and t k denote length and unit tangent vector of side k, respectively. The eigenvalues and principal directions of Q provide metrics for the cell shape and cell orientation. We measure the cell shape anisotropy by q = 2tr(Q 2 ) ∈ [0, 1), with q = 0 for round cells and q → 1 for increasingly elongated ones. We initialize the system with cells arranged according to either (i) a regular hexagonal pattern with small random deviations in the vertex positions, or (ii) random Voronoi tessellations (SM [21], Sec. I). We use periodic boundary conditions with fixed system size. We set K A = 1, A 0 = 1, K P = 0.02, γ = 1, N = 10 3 cells and P 0 = 1 if not otherwise stated.
Results Increasing β, we observe several rheological and structural transitions (Fig. 1a,b; Movies S1-S3). In Fig. 2 we show the dependence of these transitions on both β and P 0 ; however, in the following, we focus on the case P 0 = 1 (Fig. 1). For small β < β 1 ≈ 0.20 the vertex model tissue is solid (Fig. 1e) with isotropic cell shapes (Fig. 1c), where the average cell elongation is q = 0 (resp. q ≈ 0.16) when using a hexagonal (resp. Voronoi) initial state. When β increases beyond β 1 , cell shapes become anisotropic, as indicated by an increase in q (Fig. 1c). This is accompanied by a decrease in the hexatic order parameter ψ 6 (Fig. 1d), defined as ψ 6 = | Ψ j /N |, where Ψ j = k∈neighbors exp (i6θ jk )/N j and θ jk = arg (r k − r j ) [27,28]. While in this regime the shear modulus vanishes for the hexagonal initial state, the tissue remains solid, as verified through the examination of the yield stress (Fig. 1e, SM [21], Sec. I). The transition point β 1 decreases with increasing P 0 up until the critical point P * 0 (Fig. 2). The cell shape transition at β = β 1 occurs because the β term in Eq. (1) effectively corresponds to a negative shear modulus. As a consequence, when β > β 1 , the total cellular shear modulus decreases to zero, destabilizing the isotropic cell shape. To see this, we start from the Batchelor stress of a vertex model cell with perimeter P and area A [26,[29][30][31], whose anisotropic partσ is (SM [21], Sec. I):σ To obtain the global tissue shear modulus G aff in an analytical mean-field picture, we apply an affine pure shear strain to an isotropic cell, which creates the cell shape anisotropy q = 3 /2 to lowest order in (SM [21], Sec. II). Comparing Eq. (4) toσ = 2G aff , we obtain: Here, we have used A ≈ A 0 = 1, which corresponds to the limit of incompressible cells. Testing Eq. (5) numerically, where we also include all non-affinities, we find the same result, except for a prefactor: G non−aff ≈ 2G aff /3.
The behavior of cell-shape elongation q is different in the regime P 0 > P 0 , where we observe a discontinuous increase in q and coordination number Z as soon as β is increased above zero (Fig. 2a,b, SM [21], Sec. VI). The discontinuity in q can be understood by noting that in the regime P 0 > P 0 the standard vertex model is floppy with vanishing energy E(q) = 0 until some critical q value q crit ∼ P 0 − P * 0 , beyond which cells and the vertex model tissue start to attain a finite shear modulus [20]. As a consequence, as soon as β is set to a positive value, the energy below q crit becomes E(q) ∼ −βq 2 and so any state q < q crit becomes unstable. Thus, for small positive β, cell elongation will make a jump from zero to a value close to q crit .
For P 0 < P * 0 , the prediction for q in Eq. (7) breaks down beyond some β > β 2 , where β 2 ≈ 0.24 for P 0 = 1 (Fig. 1c). This corresponds to a regime marked by an increasing hexatic order, which, for tissues initiated in the near hexagonal pattern, appears to be even long-ranged (SM [21], Sec. I). (Fig. 1b-d). We call this the rhombile regime because domains appear where cells attain a rhombic shape, which arrange into a periodic arrangement of six-fold vertices (Fig. 1b). Correspondingly, the rhombile regime is marked by the increase in the average vertex coordination number Z (Fig. 1d). We observe such behavior within an intermediate range of activities β for all P 0 < P * 0 , as revealed through an increase in Z and ψ 6 (Fig. 2b,e).
The emergence of the rhombile pattern and manyfold vertices can be understood from single-and four-cell systems. Indeed, at a critical value of β (1) 2 ≈ 0.27 for the single-cell system (β (4) 2 ≈ 0.25 for the four-cell system), the two shortest cell edges shrink to length zero, resulting in the observed rhombic cell shapes ( Fig. 1a; Movie S4-S5; SM [21] Sec. II).
We find that the shear modulus in the rhombile regime is finite and peaks when the rhombile domain extension is maximal at around β * ≈ 0.3, as reflected by the local maximum in the hexatic order ψ 6 and coordination number Z (Fig. 1d). We understand this by considering a single cell: at β * = 0.3, the cell reaches a regular diamond shape with two π/3 and two 2π/3 angles. This shape (also called calisson [34]) is the building block of the crystal rhombile pattern (SM [21], Fig. S10b). When β < β * = 0.3 (resp. β > β * = 0.3), acute angles above (resp. below) π/3 form in the single cell, leading to frustration in the crystal rhombile pattern, which can destabilize the rhombile domains.
Despite a finite shear modulus, the rhombile regime exhibits a finite steady-state T1 transition rate ( Fig. 1d; SM [21] Sec. III; Movie S2). This may appear paradoxical at first sight, since T1 transitions are expected to relax the applied stress [35], leading to a long-time fluid material response. While such argument is true in a passive material at equilibrium, it may not hold for an active system at steady state, as observed here. Further, we notice that these T1 transitions are generated at the interfaces, rather than in the bulk, of the rhombile crystal domains (Movie S2).
The shear modulus eventually vanishes at β 3 ≈ 0.37. We find that this value matches the one at which the shear modulus of the rhombile crystal vanishes (SM [21], Sec. III). This suggests that the finite shear modulus of the rhombile regime is created by the rhombile crystal domains.
For β > β 3 , the system flows continuously, where the steady-state T1 transition rate increases with activity ( Fig. 1e, 2d and Movie S2). Such flow already appears in a four-cell system for an intermediate activity range (Movie S6). This flowing regime exhibits features of an active nematic material. For instance, coarse-graining of the cell orientation field reveals the presence of ±1/2 topological defects in the fluid regime (SM [21], Sec. V) [36,37]. The stress and velocity patterns around ±1/2 topological defects are consistent with those predicted in an incompressible material with extensile nematic activity [13,14,38] (SM [21], Sec. V). Our vertex model simulations also allow us to make new predictions at the We are now in the position to explain the apparent contradiction between the previous cell-based model [14], showing the onset of active flows only beyond a finite critical activity β c > 0, and active continuum models [11,12], where active flows appear at any finite activity, i.e. β c = 0. In our vertex model, such a difference arises as a function of P 0 . For P 0 < P 0 , where the vertex model behaves as yield stress solid, we showed that β acts as a negative shear modulus, destabilizing the tissue beyond a critical value that increases with the distance to P 0 . In contrast, we observe that for the fluid vertex model regime, P 0 > P 0 , any positive β destabilizes the tissue, inducing a sharp increase in flow (Fig. 2d, SM [21], Sec. VI). Comparing both cases, we conclude that while any positive activity leads to flows in a fluid material, a yield stress creates a finite activity threshold to flow. This is likely the case in Ref. [14], because the employed phase field model essentially describes a foam that likely displays a yield stress.
Finally, we tested how much our results depend on the definition of the Q tensor used in Eq. (1) (SM [21], Sec. IV). We find that the transitions in cell shape and rigidity are generic, but not the appearance of the rhombile crystal domains. We further observe that for two different Q definitions chevron patterns of local smectic order for the cell orientation emerge within the solid regime with anisotropic cells, β 1 < β < β 2 .
Conclusion Here, we studied how a feedback of cell shape on the cellular active stress generation affects collective tissue dynamics. We show that such a feedback can fluidify the vertex model tissue through a series of intermediate steps displaying hexatic order, finite shear modulus and spontaneous T1 transitions.
Perspectives We anticipate several applications to the understanding of cellular materials. First, in most instances, active cellular materials have been studied in the vertex model through the introduction of traction forces against a resting substrate [39,40]. In contrast, we introduce activity through momentum-conserving forces F (act) i [26]. This kind of activity could be useful to describe tissue flows in conditions of low environmental friction e.g. free-standing tissues, floating mouse embryos or intestinal organoids [41]. Second, our work also connects through a single active stress parameter seemingly distinct features encountered in epithelial tissues, namely rosettes [20,42] and topological defects [10]. We look forward to test our prediction of T1 transition hotspots near topological defects from experimental data of epithelial tissue. The cell dynamics is dictated by the motions of vertices, which obey the force balance equation (Eq. (2) in the main text). In addition, we implemented T1 transitions once the length of a cell edge decreases below a critical value ∆ T1 [1][2][3].

Supplemental Material Tissue fluidization by cell-shape-controlled active stresses
Simulations were started from a regular hexagonal pattern (or random Voronoi pattern if stated) with small Gaussian perturbations (magnitude = 0.01) applied to the vertices' positions at t = 0, and run for enough steps (typically N simulation = 100, 000) to arrive dynamic steady state, which can be checked by examining the evolutions of the cell elongation parameter q, the hexatic order parameter ψ 6 , the average vertex coordination number Z, see Fig. S1 for example. We systematically consider system prepared with zero residual stress (see Sec. I B). In our simulations, the governing equations and corresponding parameters are rescaled by the length scale = √ A 0 , the time scale τ = γ/(K A A 0 ) , and the stress scale σ = K A A 0 . The following non-dimensional parameter values are adopted: K P = 0.02, ∆ T1 = 0.01, and ∆t = 0.01. If not stated otherwise, we consider a system of N cell = 1, 000 and an initially hexagonal cell pattern. In Sec. III B, we verified that our results were consistent for a wide range of system sizes.

B. Generation of a global stress-free cell sheet
The stress-free configuration is defined as a configuration where the global tissue stress is zero. We here give the details of how to generate the stress-free configuration for a cell sheet without active stress.
Once we set up a hexagonal cell pattern or random Voronoi cell pattern, we compute the global stress in the cell sheet. It is possible to compute cell-based stresses from the forces applied on the vertices using the Batchelor formula [4][5][6]. In the absence of active stress, the cell stress can be expressed as: where the first term is the isotropic stress and the second term is the anisotropic stress. The global stress in a cell sheet can be calculated by, The initial tissue state correspond to a global stress σ (tissue) that is non-zero. Here we propose a two-step protocol to achieve a near zero-stress state.
1. As a first step, we apply a uniform affine transformation to the cell sheet, {r i } → {λr i }, where λ > 0 here is a scaling factor to be determined in the following. In the new configuration, {λr i }, the global stress in the tissue can be expressed as 2. Although Eq. (S3) does not strictly imply that σ (tissue) = 0, such condition yields an initial tissue configuration that is relatively close to a zero stress state. To achieve such zero stress state, we let the system evolve through the classical vertex model simulation. We find that this is efficient to decrease the residual global tissue stress. At the end of the simulation, we evaluate whether the stress is small enough; more precisely, we check whether the maximal absolute component of the stress, denoted |σ (tissue) |, is smaller than a threshold value set to σ 0 = 0.01, i.e. |σ (tissue) | < σ 0 . If the condition is met, we consider that a stress-free state is obtained; otherwise, we repeat the above steps, until the condition |σ (tissue) | < σ 0 is met.
Our protocol was always found to converge.

C. Simple shear simulations
Our simple shear simulations are performed using the Lees-Edwards periodic boundary conditions ( Fig. S2 and [2,7]). Note that before shearing, the cell sheet is relaxed to a steady state (see Fig. 1(b) in the main text for examples). In the following we detail our implementation protocol.

Lees-Edwards periodic boundary conditions
In the sheared configuration, the calculations of the relative position and distance of vertex (cell) pairs should obey the Lees-Edwards periodic boundary conditions, as shown in Fig. S2. For example, the position vector of the neighboring vertex pair A and B, r A→B , cannot be simply measured as r B − r A , but instead should be measured as r B − r A with A being the mirror of vertex A. Specifically, denote r can be calculated as where R [·] represents the function round to nearest integer and γ xy is the simple shear strain.

Tissue under a constant simple shear strain
Definition To examine the shear modulus of a cell sheet, we perform simple shear simulations as below: 1. First, we apply a small simple shear deformation γ xy = 0.01 1 to the cell sheet by mapping all the vertices to new positions, 2. Second, the simple shear strain is maintained for a sufficiently long time for the sheared configuration to relax. We keep track the global shear stress σ (tissue) xy during the relaxation and check that the evolution of the global shear falls below a threshold value. The shear modulus of the cell sheet then reads where σ (tissue) xy B and σ (tissue) xy A are the average tissue shear stress before and after simple shear strain, respectively.
Verification We checked that, within the solid regime, the shear modulus is not significantly changed for a range of values γ xy = (0.001, 0.10) for all considered β.

Tissue under quasi-statically increasing simple shear strain
Definition We also perform simple shear simulations by a quasi-static increasing of the amplitude in the simple shear strain γ xy (t). The procedure is the following: 1. First, the simple shear strain is increased in a step-by-step fashion where the incremental shear strain ∆γ xy is small and set to be ∆γ xy = 0.01 in our simulations.
2. Second, after updating the simple shear strain, vertices are mapped to new positions as where x ref = x − γ xy y and y ref = y are the reference coordinates in the reference configuration without shear.
3. Third, the new simple shear strain is kept and the cell sheet is relaxed for a duration τ relax = 100. After relaxation, the global shear stress is computed to construct the stress-strain relation σ (tissue) xy (γ xy ) during such continuous shearing process.
We then repeat the above three steps up to a maximal strain set at γ xy = 3.0. The tissue-scale shear modulus can be estimated as the derivative at zero shear strain of the shear stress curve with respect to the initial shear strain.
Verification We checked that the relaxation time τ relax = 100 after each step increase is sufficient long to provide stable estimates of the shear modulus and yield stress σ yield ), see Sec. I D 4.

Coordination number Z
To quantify the coordination number Z represented in Figure 1, we first (1) find clusters of vertices, denoted C i , that are closer to one another than a threshold distance 0.01 and (2) estimate the mean coordination number Z = Ci N (C i ) where N (C i ) is the number of vertices in the cluster C i .

T1 transition rate
Definition To quantify the dynamics of cell rearrangements, we examine the rate of T1 transitions after the system arrives at a dynamic steady state. We define the T1 transition rate as, where N T1 is the total number of T1 transitions during the observation time τ . Localization of T1 transitions We also examine the spatial distribution of T1 transitions, as shown in Fig. S3. It demonstrates that at the rhombile regime, T1 transitions are highly localized (with hotspots of frequent T1 transitions corresponding to cell-cell rearrangement oscillations) and sub-domains of rhombile patterns without any T1 transitions (see Movie S2). In the fluid regime β = 0.5, T1 transitions are more uniformly distributed, indicating a completely disordered fluid regime.
Method for counting T1 transitions We exclude repeated T1 transitions from the count of the T1 transition rate in the main text

Bond-orientational correlation function
Definition As in the main text, the tissue hexatic order parameter is defined as where Ψ j = k∈neighbors exp (i6θ jk )/N j , with θ jk = arg (r k − r j ) [8,9]. To quantify the spatial correlation of tissue architecture, we compute the bond-orientational correlation function, where denotes the conjugate complex number and r = |r i − r j |. Results Figure S4(a) shows the bond-orientational correlation function g 6 (r) for different values of β. For the solid regime with isotropic cells, g 6 (r) = 1; for both the solid regime with anisotropic cells and the rhombile regime, g 6 (r) decreases as r is increased and approaches a finite positive value for large r, suggesting long-range orientational order when the tissue is initialized according to a near perfect hexagonal pattern. We confirm such result for the three different definitions of Q. In contrast, when the tissue is initialized according to a Voronoi pattern, g 6 (r) decreases exponentially to zero. In the fluid regime β > β 3 (see main text), g 6 (r) decreases exponentially to zero at large r.

Yield shear stress
Definition In Sec. I C, we proposed a protocol to implement a quasi-static increase of a simple shear strain. We represent typical stress-strain curves for different values of β for tissues either initialized in a near hexagonal state ( Fig. S5(a)) or according to a Voronoi pattern (Fig. S5(c)). We define the yield stress σ yield as the maximum of the shear stress σ xy curve. The yield stress is represented in the main text Figure 1(e).
Discussion When the initial pattern is near hexagonal (Fig. S5(a)), the shear stress sharply drops beyond a critical yield strain. Such a sharp drop of σ xy is due to a large number of simultaneous T1 transitions. The curve is more smooth in the case for a tissue initialized according to a Voronoi pattern.
Verification We check that the measured shear modulus and yield stress are insensitive to the relaxation time duration of simple shear strain increment in a broad range 10 < τ relax < 1000, Fig. S5(b).

E. Simulation check for the passive case
Here we set β = 0 (passive case) and check that our simulation framework agrees with previous studies [2,7,[10][11][12]. We first check it for an initially hexagonal cell pattern. We represent the cell elongation q and shear modulus G xy as function of the target shape index P 0 in Fig. S6, which clearly indicates a shape and rigidity transition at P * 0 ≈ 3.723. We also check it for an initially Voronoi cell pattern. Figure S6 shows that such a cell pattern undergoes a rigidity transition at P * 0 ≈ 3.846.

A. Shape stability analysis of a single cell
Here we consider a regular hexagonal cell with shape-dependent active stress σ (act) = −βQ and ask whether the cell will elongate. To theoretically address such question, we consider a uniform deformation of the cell, as shown in Fig. S7. The cell deformation can be described by the stretches {λ 1 , λ 2 }: {x i , y i } → {λ 1 x i , λ 2 y i }, where i is the index of cell vertices and λ 1 and λ 2 are cell stretches along x-axis and y-axis, respectively.

Equilibrium radius of a hexagonal cell without active stress
We first examine the equilibrium state of a hexagonal cell without active stress. We express the area and perimeter of the cell as A = 3 √ 3 2 R 2 and P = 6R in terms of a quantity R, that we call the cell radius. We define the dimensionless equilibrium cell radiusR h = R h / √ A 0 as the solution of the mechanical equilibrium condition σ = 0, which reads In the incompressible limit case, i.e.,K P 1, the equilibrium cell radius can be easily computed asR h 2 √ 3/3 ≈ 0.6204.

Equilibrium configuration of a hexagonal cell with active stress
We now consider the effect of an active stress σ (act) = −βQ, in which case the total cell stress reads Under an affine deformation: the cell area, perimeter and anisotropy tensor of an initially hexagonal cell read The equilibrium state {λ 1 , λ 2 } of a cell corresponds to σ = 0. Substituting the expressions of A, P and Q into Eq. (S13), the condition σ = 0 leads to: (S19) whereβ = β/(K A A 0 ) corresponds to the dimensionless strength of the active stress. Equations (S18) and (S19) shows that {λ 1 , λ 2 } = {1, 1} is always one of the equilibrium state of the single cell system for arbitrary β. Yet increasing β, the system undergoes a pitchfork bifurcation and the cell start to elongate, as shown in Fig. S8.

Stability analysis of a hexagonal cell
Now we turn to examine the stability of the equilibrium state {λ where is the Jacobian matrix at the equilibrium state {λ 1 , λ 2 } = {1, 1}; we further defined the quantities ε 1 = λ 1 − 1 and ε 2 = λ 2 − 1. The stability of an equilibrium hexagonal cell is equivalent to the positive definiteness of J σ , which results in the set of conditions The first inequation represents the stability condition for the single cell system without active stress. Therefore, under the shape-dependent active stress σ = −βQ, the stability condition of a hexagonal cell at rest radius R h is whereR h depends on (K P ,P 0 ), as given by Eq. (S12). When β > β 1 , the regular hexagonal cell is unstable: under small perturbations, the cell will elongate and cannot recover to a regular hexagonal shape. Our numerical calculations suggest that such cell shape transition belongs to the pitchfork bifurcation class (see Fig. S8). Further, in the incompressible limit case, i.e.,K P 1, the force balance equation (S12) givesR h 2 √ 3/3 ≈ 0.6204. Then the critical cell activity can be further simplified as: This results is confirmed by our numerical simulations (see main text Fig. 2).

B. Estimation of the cell elongation beyond the cell-shape instability critical point
Here we get an analytical prediction for the cell elongation parameter q for β β 1 . For simplicity, we consider the incompressible limiting case, i.e.,K P 1. In that case, A A 0 , λ 2 1/λ 1 (as defined in Eq. (S14)). Based on the expression of Q (Eq. (S17)), we obtain q = 2tr (Q 2 ) = 2λ 1 Further assuming λ 1 1, we get λ 1 = . We expand λ 1 and λ 2 to second order terms in q: Substituting Eqs. (S28) and (S29) into Eq. (S16), we obtain Consequently, the anisotropic stress can be expanded to higher order terms of q as, Combining it withσ = ∂E/∂ε /(2A) with ε = λ 1 − 1 ≈ 2q/3, we obtain an effective potential energy of the cell, E eff (q), which reads to fourth order terms in q as: which corresponds to the expression derived in the main text. Minimizing E eff (q) with respect to q yields an approximation of the cell elongation q min : In Figure S9(a), we represent several E eff (q) curves for several β values; this suggest that an instability upon increasing β beyond β 1 ≈ 0.2. Figure S9(b) shows the comparison of the theoretical predicted cell elongation q in Eq. (S33) to the exact cell elongation q obtained by solving Eqs. (S18) and (S19) and numerical simulations of a single cell (see Sec. II C). We can see that all match very well for β close to the transition point β 1 ≈ 0.2.

C. Simulations of single-cell system and four cell system
To gain physical insight into the tissue fluidization driven by cell-shape-controlled active stress, we also perform simulations in a single-cell system (Fig. S7) and a four-cell system (see main text Fig. 1(a), β = 0.150). In both systems, initial configurations are regular hexagonal patterns with small perturbations applied to each vertex.
We find that both the single-cell system and the four-cell system exhibit an elongation transition at around β ≈ 0.20 for P 0 = 1 (Fig. S10(a,c)), as in larger scale simulations (main text Fig. 1(c)). These two systems also undergo a second transition where some edges vanish. Specifically, the four-cell system undergo such a transition at β ≈ 0.25 (Fig. S10(c)), which corresponds to the end of the solid anisotropic cell regime identified in our larger scale simulations (see also Movie S4 and S5 to visualize these transitions upon increasing β).
Further, we measure the minimum cell angle θ M in the single-cell system. Interestingly, we find that the single cell system exhibits a perfect rhombile shape with the minimum cell angle θ M ≈ 60 o at β ≈ 0.3. This can explain why we observe the maximum fraction of rhombile regime near β ≈ 0.30.
In addition, we find that at the higher activity β = 0.4, there are many metastable states in the four-cell system (Fig. S10(d)).
Finally, we compare the cell elongation parameter q for simulations of the system-cell system, the four-cell system and the cell sheet system in Fig. S10(e); all agree very well for β close to the first transition point.

A. The rhombile regime exhibits a finite shear modulus
We examine more closely the rheological properties of the rhombile regime defined in the main text. In Fig. S11(a), we represent a typical evolution curve of the shear stress before and after the application of a small simple shear deformation for β = 0.3 (P 0 = 1). The sustained increase in the measured shear modulus is consistent with the interpretation that this tissue displays a finite shear modulus. We further perform a statistics in the measured shear modulus for different initial tissue configurations, see Sec. I A; we find that the mean shear modulus is positive, Fig.  S11(b), which is also consistent with the interpretation that the tissue displays, on average, a finite shear modulus.
To confirm that such finite shear modulus is not due to numerical effects, we further investigated the roles of • the T1 transition length threshold ∆ T1 . We performed simulations with varied T1 transition length threshold ∆ T1 and we find that the shear modulus is not sensitive to ∆ T1 ; specifically for β = 0.3, setting ∆ T1 = 10 −3 , 10 −2 , 10 −1 yielded consistent shear modulus values G xy /K P = 0.3244, 0.3000, 0.3321, respectively.
We conclude that the tissue is solid within the identified rhombile regime. To explore the effect of system size on the active stress driven tissue fluidity transition in cell sheets, here we vary the system size from a small one (N = 100) to a large one (N = 10, 000). We find similar tissue fluidity transition regimes: as the cell activity β is increased from small values to large values, the tissue transitions from the solid regime (isotropic cells), the solid regime (anisotropic cells), the rhombile regime, to finally the fluid regime, as shown in Fig.  S12(a,b). Moreover, the critical β values between these regimes for a large system (N = 10, 000) reads: β 1 ≈ 0.20, β 2 ≈ 0.24, and β 3 ≈ 0.37. These critical β values are close to those presented in the main text. Finally, we present the rhombile regime in a large system (N = 10, 000) in Fig. S12(c).

C. Stability of a perfect rhombile pattern
Here we examine the stability of a global rhombile tiling pattern, as shown in Fig. S13(a), where cells are of identical diamond shape and display both 3-way and 6-way vertices. Here, the 6-fold vertices, being considered as a single vertex, cannot split into two separate 3-way vertices.
Our simulations show that the rhombile tiling pattern becomes unstable as the activity parameter β is increased beyond the threshold value β cr,3 ≈ 0.38, see Fig. S13(a,b). The value β cr,3 ≈ 0.38 matches the value mentioned in the main text as delineating the rhombile regime from the disordered fluid regime.
We further estimated the shear modulus of the rhombile tiling pattern through numerical simulation, see Fig. S13(c) and find that the rhombile tiling pattern is of finite shear modulus for β < β cr, 3 . The presence of nucleated rhombile crystals may explain why a residual shear modulus can be observed within the rhombile regime.

A. Vertex-based inertia tensor
Definition Here we first consider an alternative cell shape description to Q defined in terms of a vertex-based inertia tensor where we have introduced the relative vertex position vectors ρ i = r i − r C with respect to the cell center r C = k r k /N . By analogy to Q defined in main text, which is a deviatoric tensor, we define the following vertex-based inertia Q-tensor: Results We find that the hexagonal pattern is destabilized at a critical coupling value β c ≈ 0.2, see Fig. S14(a,b) in favor of a solid regime with anisotropic cells in the range β ∈ (0.2, 0.35). Within the latter range of activity, we observe buckling along tissue-scale stripes, reminiscent of the chevron pattern within smectic A liquid crystal [13]. Such chevron pattern is of finite shear modulus, Fig. S14 Discussion Overall, the results for the vertex-based inertia cell shape tensor are similar to the one presented in the main text. Following Ref. [6], one can show that the inertia tensor M v displays the same principal orientations as the cell nematic shape tensor Q defined in the main text; yet the associated eigenvalues for these two tensors are, in general, different.

B. Area-based inertia tensor
Definition We next consider an area-based inertia tensor M a , defined as, where r is the position vector relative to the cell centroid. Similarly, we define the following Q-tensor: Results -similarities with other schemes We find that the hexagonal pattern is destabilized at a critical coupling value β c ≈ 0.2, see Fig. S15(a,b) in favor of a solid regime solid with anisotropic cells. Within the range β ∈ (0.2, 0.27), we also observe buckling along tissue-scale stripes (Fig. S15(a)) that are similar to those obtained with the vertexbased inertia tensor (Eq. (S35)); the corresponding tissue state display a finite shear modulus (Fig. S15(d)). At a second critical value β 2 ≈ 0.27, the vertex coordinate number Z increases, see Fig. S15(c), which marks the end of the anisotropic solid regime. For larger value of the activity parameter β > β 3 = 0.32, spontaneous flows emerge.
Results -differences with other schemes In the intermediate activity (β 2 , β 3 ), we observe some higher order vertices, but less than in the other schemes for Q. We do not observe well-defined rhombile crystal domains, Fig. S15(a,c).
Discussion With the area-based inertia definition, the rhombile crystal pattern is not observed as β is increased. We interpret that through a single cell shape analysis, which leads to less of shrinkage at cell edges in the case of the area-based inertia tensor. At the tissue level, this is reflected by a much smaller T1 transition rate k T1 (Fig. S15(b)) than for the standard inertia tensor Q defined on cell edges and Fig. 1(c)) or the vertex-based inertia tensor (Eq. (S35) and Fig. S14(b)). In presence of cell active stress, for β > β 1 , a cell sheet exhibits elongated cell shape and coordinated orientations in local areas that are conflicting at singular points, called topological defects (see Fig. S16(a)).
To detect the locations and orientations of such topological defects, we first compute the cell geometric centers r J = i∈neighbor r i /N J with N J being the number of vertices belonging to the J-th cell, and the nematic tensors Q J of cells. Based on r J and Q J , we next construct a coarse-grained nematic field Q(r) on a regular grid, where w(r − r J ) is the weight function and r cut−off = 3σ is a cut-off length. Here we set the weight function as a Gaussian function, where we set the kernel size σ at σ = 0.75 cell length. We set the cut-off length at r cut−off = 3σ = 2.25 cell length.
After obtaining the coarse-grained nematic field Q(r), we detect the location of topological defects based on the calculation of winding number [14,15]. The orientation of topological defects is extracted following the method presented in [16]; for a topological charge k (k = ±1/2), the orientation is given by where · denotes an average along the shortest available loop enclosing the defect core; sgn(·) is the sign function; arctan 2(y, x) is the 2-argument arctangent function.
To check the reliability of such detection scheme, we superimpose the detected topological defects on the original cell sheet morphology and coarse-grained nematic field, as shown in Fig. S16. We can see that the detected topological defects fit well with both the original cell sheet morphology and the coarse-grained nematic field. We mainly observe two kinds of topological defects, i.e., comet-like +1/2 topological defects and trefoil-like −1/2 topological defects.

B. Stress and flow fields near topological defects
We examine the average stress and flow fields near ±1/2 topological defects in our vertex model implementation. The resulting fields are consistent with the predictions of an incompressible active nematic continuum theory [14].
Method -We first examine the average isotropic stress (σ iso = (σ xx +σ yy )/2) field. We first define a local coordinate system of a topological defect: its origin locates at the defect core, with x-axis along the defect orientation. The average isotropic stressσ iso at position x (local coordinate system) is estimated through the following expression σ iso (x) = |x−r J,α |<r cut-off w (x −r J,α ) σ (J) iso |x−r J,α |<r cut-off w (x −r J,α ) .
Here σ (J) iso is the isotropic stress of the J-th cell; w(·) is a Gaussian weight function (Eq. (S39));r J,α represents the position vector of the geometric centre of the J-th cell relative to the α-th topological defect, and can be calculated according to the following expression:r where θ α refers to the orientation of the α-th topological defect; R[r; θ] represents a rotation of the vector r by an angle θ. Similarly, the average flow velocityv(x) at position x (local coordinate system) is estimated aŝ v (x) = |x−r J,α |<r cut-off where v J = dr J /dt is the motion velocity of the J-th cell. Result - Figure S17 shows the coarse-grained average isotropic stress and flow fields near ±1/2 topological defects, for β = 0.4 > 0 (extensile active nematic case). We observe that both the stress and flow fields are consistent with previous experimental measurements (including MDCK cell monolayer system [14,17], human bronchial epithelial cell monolayer system [18] and neural progenitor cell monolayer system [15]) and predictions of a continuum active gel theory [14,19].

C. Cell-cell rearrangement field near topological defects
In the main text Figure 3, we examine the cell-cell rearrangement rate field near topological defects. We define the cell rearrangement rate as the number of T1 transitions per unit time per unit area, estimated according to the following formula χ (x) = |x−r T ,α |<r cut-off δ (x −r T,α ) wherer T,α = R [(r T − x α ) ; −θ α ] represents the position vector of the T -th T1 transition relative to the α-th topological defect; the observation period spans a duration τ preceeding the end of simulation (here τ = 200).

VI. THRESHOLDLESS FLOWS
For P 0 > P * 0 ≈ 3.72 (hexagonal cell pattern), we observe spontaneous tissue flows for very small cell activity, suggesting that the transition to flow is thresholdless, i.e. at exactly β = 0 in the limit of a very large system N → ∞.
For illustration, here we set P 0 = 4.2 and vary β. Figure S18(a) shows the morphology of a cell sheet at zero cell activity (β = 0) and small cell activity (β = 0.001), respectively. We observe strikingly different cell morphologies at β = 0 and β = 0.001: in the former, cells are concave polygons, with no rosettes; in the latter, cells are elongated β = 0.001, with many involved in rosettes. In Fig. S18(b) we represent the T1 transition rate k T1 as a function of cell activity β for different system sizes. Notably, when the system size is increased, the activity threshold for spontaneous T1 transitions shifts towards zero, which suggests thresholdless tissue flows. In addition, we observe a discontinuous (first-order like) transition at β = 0 in both the cell elongation parameter q and the average vertex coordination number Z upon increasing cell activity β, as shown in Fig. S18(c). These results suggest that even a very small cell activity leads to a drastic change in the cell shape pattern (including the onset of multicellular rosettes).