{'fraction_non_alphanumeric': 0.04603246851767562, 'fraction_numerical': 0.013867394932483689, 'mean_word_length': 4.461716937354988, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From cells to tissue: A continuum model of epithelial mechanics

A continuum model of epithelial tissue mechanics was formulated using cellular-level mechanical ingredients and cell morphogenetic processes, including cellular shape changes and cellular rearrangements. This model can include finite deformation, and incorporates stress and deformation tensors, which can be compared with experimental data. Using this model, we elucidated dynamical behavior underlying passive relaxation, active contraction-elongation, and tissue shear flow. This study provides an integrated scheme for the understanding of the mechanisms that are involved in orchestrating the morphogenetic processes in individual cells, in order to achieve epithelial tissue morphogenesis.


I. INTRODUCTION
During tissue morphogenesis, tissues acquire their unique shape and size through a series of deformations. Morphogenesis occurs at multiple levels, and molecular, cellular, and tissue level changes are interdependent. At cellular level, tissue deformation is accounted for by changes in cell shape, position, and number ( Fig. 1(a); hereafter named cell morphogenetic processes), which are triggered by biochemical signaling and forces generated by cells [1][2][3]. While the tissue stress can affect cell morphogenetic processes through the changes in molecular activity and localization, cell morphogenetic processes generate stress [4][5][6][7][8][9][10]. However, the mechanisms by which the shape of a tissue emerges from these multi-scale feedback processes remain unclear.
In order to clarify this, a coarse-grained description and modeling of cellular and tissue dynamics at an appropriate length scale is required: while the position and timing of cell morphogenetic processes are stochastic at single cell level, the averaging of values obtained over a larger length scale yields a smooth spatial pattern that is reproducible among different samples. We previously determined the appropriate averaging length scale for describing epithelial tissue dynamics (several tens of cells in a patch), and developed coarse-grained methods for measuring stress and kinematics [4][5][6]11]. A force inference method was used for the quantification of cell junction tensions and cell pressures ( Fig. 1(b), (c)), which can be integrated to calculate a stress tensor [4,5,15,16]. A texture tensor method was used for the measuring of different cell morphogenetic processes (e.g., cell division, cell shape changes) in the same physical dimension, which can be further integrated to obtain tissue scale, spatio-temporal maps of tissue growth and cell morphogenetic processes [6].
Together, these methods provide the information on the amplitude, orientation, and anisotropy of tissue stress, tissue growth, and cell morphogenetic processes, and correlations between them [6].
A modeling scheme capable of accommodating the quantitative data described above is still lacking [17]. Cell-based models, such as the cell vertex model (CVM) [18] and the cellular Potts model (CPM) [19] are often employed for the simulation of epithelial tissue morphogenesis ( Fig. 1(d); [3,21,22]), and have proven useful for including experimental data obtained at cellular level, such as the laser ablation of cell junctions or subcellular distribution of proteins. However, the relationship between cell morphogenetic processes and tissue scale deformation and rheology emerges from numerical simulations without being directly tractable. Continuum models allow the indepth analysis of tissue rheology [1,23], yet in many cases do not include the information of the cellular structure by construction, and thus fail to discriminate between different cell morphogenetic processes. A limited number of studies considered the degrees of freedom that represent cell morphogenetic processes and cell polarity [1,[25][26][27], but do so in the context of macroscopic models, which do not incorporate cell-level mechanical parameters explicitly. The finite-element model introduced in [28] includes at a coarse-grained level the contributions of cellular rheology, shape changes, rearrangements and divisions. A continuum model has been derived from the CVM previously but without considering cell rearrangements [29] (see also [30,31] in 1D).
The main aim of our study was to develop a two-dimensional hydrodynamic model of the epithelial tissue. This model included a field that represents coarse-grained cell shape, which enabled us to treat different types of cell morphogenetic processes distinctively. First, kinematics was identified by decomposing tissue deformation into cell shape changes and cell rearrangements. Following this, by introducing an energy function deduced from CVM/CPM, thermodynamic formalism was employed to determine kinetics. The model we derived here describes tissue deformation through stress and deformation tensors, which can be compared with the data obtained experimentally [4][5][6], and can incorporate active terms smoothly. We solved the model for several conditions typical for deforming planar tissues during development, and demonstrated that the model predicts the relaxation of cellular shape following the tissue stretching, the relation of the direction of cell elongation and tissue flow during active contraction-elongation (CE), and shear-thinning in tissues.
Our approach provides a theoretical framework that enables to assess how cellular level mechanical parameters and cell morphogenetic processes are integrated to realize tissue-scale deformation.

A. Cellular shape tensor
To construct a continuum model, we approximated each cell by an ellipse, characterized by a 2 × 2 symmetric tensor M . Cellular shape can be expressed as ( where r c represents the center of a cell and superscript T denotes the transpose (Fig. 1(e)). Since the eigenvalues of M are the square lengths of ellipse semi-axes, cell area can be presented as where |M | is the determinant of M . By coarse-graining over a representative surface element comprising a sufficient number of cells [6,11], we obtained a spatially smooth tensor field M ( r). Similar to a previously described texture tensor [6], the symmetric tensor M ( r) represents measure of tissue scale deformation in our model, with the physical dimension of square length, which can be experimentally quantified from segmented images of two-dimensional (2D) epithelia.

B. Kinematics
Total tissue deformation rate can be represented by the tensor ∇ v, in which v is the velocity field.
Here, we used (∇ v) ij = ∂ j v i , where indices i and j denote cartesian coordinates. The deformation rate ∇ v represented the sum of its symmetric part D = (∇ v + [∇ v] T )/2 and its antisymmetric part We decomposed tissue deformation rate into the sum of contributions, due to the cellular shape alterations and other cell morphogenetic processes, and here, we considered cell rearrangement, division, and death: The quantity Ω + D s represents the tissue deformation rate stemming from cellular shape alterations, while D r denotes the deformation rates that involve topological changes in a network of cell junctions, i.e. cell rearrangement, division, and death. We assumed that these processes are irrotational, so that D r is symmetric. In practice, these tensors can be experimentally determined by cellular shape tracking [6]. In cellular materials, D r is kinematically associated with cell shape whereṀ ≡ ∂ t M + v∇M represents the Lagrange derivative of M ( Fig. 1(f)). The kinematic relationship (B3) was derived for the non-affine deformation in the rheology of polymer melts [32][33][34] and foams [35,36], where D r can be interpreted as due to the slippage between polymer molecules and between bubbles, respectively. When D r = 0, the left-hand side of (B3) becomes the co-deformational upper-convected derivative [37]. In SI Section 3, we demonstrate that the conservation of cell number density ρ = 1/π|M | 1/2 can be presented as where Tr denotes the trace, and Tr D r coincides with the variation rate of ρ. The effects of cell division and cell death should be investigated in future studies, while here we considered a situation in which individual cells only deform elastically and/or intercalate, and the tissue area is invariant and the deformation rate is traceless (Tr D r = 0).

C. Energy function and elastic stress
In CVM and CPM, cell geometry can be determined by minimizing energy function [2,3] where A i and L i represent the area and perimeter, respectively, of cell i, and ℓ ij is the length of the interface between cells i and j. The first term represents cell-area elasticity with elastic modulus K and reference cell area A 0 . The second and third terms represent cell junction tensions, where the tension is given by γ 0 + κ 0 (L i + L j ).
Here, by using the cell shape tensor M ( r), we considered an energy function for the continuum model, which is comparable with that of the CVM/CPM, (3). For arbitrary semi-axes a, b, the perimeter of an ellipse can be given by Euler's formula (see Suppl. Fig. S1) [8]: where F is a hypergeometric function. Using the Cayley-Hamilton equation, we derived the iden- Upon coarse-graining (3), we obtained energy density per unit area where L(M ) is defined for arbitrarily large cellular shape anisotropy by The total elastic energy was obtained by integration F = F (M ( r)) d r.
We focused, for simplicity, on conditions close to the isotropic case (see SI Sec. 1.2 for higherorder approximations). Expanding F close to a = b, or (Tr M ) 2 = 4|M |, the zero order energy can be written as from which the elastic stress can be further derived as where I is the unit tensor (SI Sec. 1, general case). The first and second terms represent the isotropic pressure −P ce I due to the area-elasticity of cells, and the cellular shape-dependent stress σ T due to cell junction tensions, respectively. (A4) is consistent with the expression of the Batchelor stress tensor [4,5,40] relating tissue-scale stress to cell pressures and cell junction tensions. Note that σ e and M commute, and therefore, share the same eigendirections, which is consistent with our previous observation showing that cells are elongated along the inferred maximal stress direction in Drosophila epithelia [4,6,6].
We performed numerical simulation of CVM under the isotropic and anisotropic stress environments, and compared the coarse-grained stress with the true one. Coarse-grained cell shape M was evaluated by averaging the second moment of cell shape, and Euler expansion up to the second order was considered (SI Sec. 1.2). The results obtained here confirmed that the coarse-grained stress values agree with those obtained for the true one ( Fig. 2(a,b), Suppl. Fig. S2).
Factorizing the tensor M as M = M 0 e cΘ can be useful [41], and here, M 0 and c are scalar fields and Θ is a symmetric, trace-less tensor field parameterized by the angle θ as Since existing cell-based models, including CVM and CPM, use ad hoc prescriptions for kinetics, we considered the thermodynamic formalism [43] in order to derive generic hydrodynamic equations. The total stress tensor σ can be given by the sum of elastic and dissipative stresses, as The entropy production rate of an isothermal process was calculated as [33,34]: where s is the entropy density, and T is the temperature. Here, a:b ≡ Tr [ab T ] denotes the scalar product of two arbitrary tensors a and b, and a ′ ≡ a−(Tr a/2)I is the deviatoric part of a. Because D r is traceless, we can replace σ e by σ e ′ in (11). By identifying conjugate flux-force pairs as σ p -D and D r -σ e ′ , the fluxes (σ p , D r ) can be given by using the linear combinations of the forces (D, σ e ′ ): The coefficients χ ss , χ sr , χ rs , and χ rr are fourth-order tensors that satisfy Onsager's reciprocity (e.g., χ sr ijkl = χ rs klij ). Note that Maxwell's model was obtained for χ ss = χ sr = 0 (σ = σ e , σ p = 0), and that Kelvin-Voigt's model was obtained for χ rs = χ rr = 0 (D = D s , D r = 0) [33,34]. The term χ ss D characterizes dissipative stress due to the tissue strain rate, and reduces it to the usual bulk and shear viscous terms for isotropic material. According to (13), cell rearrangements may be driven both by the tissue strain rate and by its elastic stress [5,6,44]. The presence of the (12) is a non-trivial prediction of the model, which cannot be represented by a rheology diagram in terms of a combination of springs and dashpots.
In the absence of external forces, the above equations can predict the relaxation to the steady state, and are not sufficient to address the active phenomena, such as the sustained epithelial flow [9] or self-organized spatio-temporal patterning [46]. Therefore, we included the formalism of active gels [47,48], and added the term r∆µ to the entropy production rate (11): where ∆µ represents the change in chemical potential associated with a chemical reaction that supplies energy to the system, and r is the reaction rate. By identifying r-∆µ as an additional force-flux pair, additional terms σ a and D a , which we refer to as active stress and active cell rearrangement, respectively, appeared in (12) and (13), both of which are proportional to ∆µ.
The coupling coefficients can depend on M ′ . Including lowest-order non-linearities, with the condition that D r is traceless, the generic form of the force-flux relationships can be written as: where the coupling coefficients are scalar in an isotropic system. In (B9), η and η ′ denote the tissue shear and bulk viscosity, respectively, and µ, µ ′ , and µ ′′ denote their anisotropic correction depending on the cellular shape. The term ν 1 D ′ of (B10) plays a role similar to that of the Gordon-Schowalter process in the rheology of polymer melts, which describes the relaxation of polymer deformation, and, consequently, stress, by slippage [37]. The cross terms, including ν 2 , and ν 3 , are possible in general. We introduced the active stress tensor as σ a = − (ζ 1 I + ζ 2 M ′ ) ∆µ. Using the terminology of active nematic liquid crystals, a negative and positive ζ 2 values correspond to a contractile and extensile, respectively, material [48]. These activities are often attributed to myosin contractility, for which ATP is consumed. The terms coupling D r and σ e ′ in (B10) underlie the Maxwellian dynamics of the system, and include the positive coefficient η 1 with the dimension of viscosity. In (B10), active cell rearrangements contribute to the constitutive equation for D r as Both σ a and D a are symmetric second-order tensors. Our treatment of active stresses and active cell rearrangements was similar to that suggested previously [27,44], since both approaches are inspired by the active gel models [47,48].
Finally, the force balance equation was used to close the system: where f ex represents the external force field, supplemented with the appropriate boundary conditions. Collectively, the constitutive equations (A4,10,B9,B10) with the kinematic relationships (1,B3) and the force balance equation (17) determine hydrodynamic equations of a tissue (SI Sec. 2).

III. APPLICATIONS
We investigated three simple examples of dynamical behavior predicted by our model, including the passive response following the axial stretch induced by an external force, the deformation of a tissue due to the active internal forces, and the generation of shear flow. Two assumptions were used for simplicity to obtain the following analytical solutions: 2D incompressibility of a tissue (M 0 is constant and Tr D = div v = 0; SI Sec. 2) and spatial homogeneity of all relevant fields.
A. Passive relaxation following the axial stretching In Drosophila pupal wing, an external force from the proximal part of the body is responsible for the stretching of the wing along the proximal-distal (PD) axis. Upon the tissue stretching, wing cells elongate along the PD axis, while the tissue relaxes during several hours when cells intercalate and adopt a less elongated shape [5,6,6,44]. Below, we demonstrated that our model can recapitulate this process, and determined the characteristic relaxation time in terms of cell mechanical parameters.
We considered a tissue with an initial size L x × L y = L 0 × L 0 and in an initial isotropic, uniform state where the cell shape tensor is M = M 0 I. From time t = 0, the tissue elongates along the x-axis at the constant rateL x /L x = ∂ x v x = λ, and consequently, contracts along the y-axis at the rateL y /L y = ∂ y v y = −λ ( Fig. 3(a)). When the tissue size reaches aL 0 × a −1 L 0 , the stretching . We attempted to identify a uniform solution to the problem, so that (17) is automatically verified when f ex = 0. All tensor variables are diagonal (M = M 0 e cΘ with θ = 0 ). The M -dependent terms in Eqs (B9-B10) are either canceled, or are isotropic tensors that may be absorbed into the pressure. Since we consider a passive process here and ignore active terms, we set ∆µ = 0. Using (B3), the time evolution of c can be written as: where Γ(c) = γ 0 /2 √ M 0 cosh c + 2πκ 0 represents the strength of cell junction tensions as a function of cell mechanical parameters γ 0 and κ 0 . The normal stress components were where the pressure p was determined in order to satisfy the incompressibility condition.
Setting a = 5, the time course of cell shape anisotropy c(t) is presented for several values of λ in Figure 3(b). When |c| ≪ 1, the temporal evolution becomes Maxwellian, CE denotes the simultaneous shrinkage and expansion of a tissue along two orthogonal axes [52], often controlled by the anisotropic localization/activity of signaling/driving molecules, such as molecular motors [1,3,11,53,54]. During the CE, cells are often elongated along the axis perpendicular to the direction of the tissue flow ( Fig. 4(a)) [52][53][54][55] which may occur in order to facilitate the force transmission along the axis of tissue contraction, through the formation of multicellular myosin cables through the mechanosensing of neighboring cells [54]. However, the mechanisms whereby tissue deformation due to the cellular shape alterations counteract that due to the cell rearrangements remain unclear.
Here, we investigated the CE by extending our model to include active stress and rearrangements provided by signaling molecules oriented along a fixed direction, as represented by a traceless tensor Q = n ⊗ n − Tr ( n ⊗ n)/2, where n = (cos φ, sin φ) T represents a unit vector field pointing to the direction φ. Possible feedbacks on the signal activity were ignored. We considered the lowest-order active contributions where the parameters ζ Q and β Q , respectively, quantify the strength of the active stress and of the active rearrangements. Negative and positive ζ Q values correspond to the contractile and extensile, respectively, stress along the direction n. For positive β Q , D a drives cell rearrangements where a cell junction parallel to n shrinks and is remodeled to form a new cell junction perpendicular to n.
As above, −ζ 1 ∆µ I is absorbed into the pressure term when the tissue is incompressible.
We considered a uniform and fixed signal n = (0, 1) T . We set ζ 2 = β 2 = 0 to focus on the activity induced by Q (SI Sec. 3.1 for full calculation). Assuming as above that θ = 0, Eqs. (18)(19)(20) becomeċ In isotropic stress conditions (σ xx = σ yy ), both the cellular shape anisotropy at steady state (ċ = 0), determined by: and the tissue deformation rate ∂ remain non-zero, indicating that the tissue anisotropy and flow are maintained through the active stresses and cell rearrangements. In Eqs. (27) and (28), cellular shape anisotropy c and velocity gradient ∂ x v x may adopt either an identical or an opposite sign depending on the numerical values of the active coefficients β Q and ζ Q . Therefore, cell elongation occurs either parallel or perpendicular to the direction of flow, depending on the parameter values (phase diagram in Fig. 4(b), with ν 1 < 1). Considering the contractile effect of myosin motors (upper right quadrant of Fig. 4(b), −ζ Q > 0 and β Q > 0), cell elongation occurs mostly perpendicular to the tissue flow, except below the green line of the slope (1 − ν 1 )/η. An earlier study using CPM suggested that the differential cell adhesion accounts for CE with cell elongation orthogonal to tissue flow, in which only the outer tissue boundary contributes to the driving of tissue deformation [56,57]. Our model provides an alternative mechanism in which activities play an essential role, in agreement with recent observations of elevated myosin activity in the elongated cell junctions orthogonal to the tissue flow [53,54].

C. Shear flow
A fundamental geometry for investigating rheology [58], shear flow is commonly found in many developmental tissues [6,9,54]. Here, we considered the simple geometry given in Fig. 5(a), inspired by the plane Couette flow [58]. The flow, with shear rateγ = ∂ y v x , is driven by an external shear stress σ b acting in the opposite directions on the boundaries. The effective shear viscosity η eff can be calculated as: where the cell shape anisotropy c and orientation θ depend on the driving stress σ b (SI Sec. 3.2). In the presence of a coupling between the cellular rearrangement rate and the elastic stress (η 1 = 0), η eff depends on σ b , which makes the tissue non-Newtonian (Fig. 5). The shear rate is an increasing function of the external stress σ b , whereas the effective shear viscosity decreases withγ, indicating that the model predicts shear-thinning ( Fig. 5(b-c)). Cellular shape anisotropy c increases with σ b orγ, to converge to a finite value for large driving. Cells turn in the direction of the applied stress as they elongate (Fig. 5(d-f)). Shear-thinning was reported in vitro, using cellular spheroids [59], and was shown to be related to stress-dependent barriers that may control cell rearrangements (see [48] for a mechanism leading to shear-thinning in active materials that does not involve topological effects). To the best of our knowledge, there is no experimental evidence for the shear-thinning in epithelial tissues in vivo, which is a non-trivial prediction of our model, obtained assuming only linear force-flux couplings.

IV. DISCUSSION
We formulated a continuum model of epithelial mechanics based on the definition of a tissue-scale field M representing cellular shape. The advantages of using this approach are as follows. Most importantly, our model is designed to connect cellular level mechanical ingredients (e.g., cell area elasticity and cell junction tension) and cell morphogenetic processes (e.g., cell rearrangements), in order to drive tissue mechanics and deformation. This was achieved by defining the energy function and kinematic relationship in terms of the cell shape field M , which distinguishes our work from the previous continuum models [1,27,28]. The model describes time-dependent flows, and allows the evaluation of time scales as a function of material parameters. Large and nonaffine deformations can be treated. In addition, the model can also incorporate a signal field, for instance, the axial tensor Q, which, here, orients active stresses and cellular rearrangements. Many relevant fields can be experimentally determined, including tissue stress, tissue deformation, cell morphogenetic processes, and chemical signaling fields, such as the concentrations of cell polarity molecules or the orientation field describing the spatial distribution of myosin molecules. Once their dynamics are quantitatively characterized by the relevant scalar, vector, or tensor variables [4-6, 6, 11, 44, 60, 61], comparison between the model and experiments is feasible.
The data presented here demonstrate that using our model, we can predict the dynamical behaviors underlying epithelial tissue morphogenesis. In the future, the quantitative comparison of the model predictions with experimental data should help us evaluate material parameters and validate constitutive equations. For example, since D, D r [6] and σ e [4,5] are measurable quantities, the validity of (13) can be tested by comparing it with the experimental data. Furthermore, the relevance of the non-linear terms in (B10) for epithelial rheology may be tested.
The current approach can be extended in several ways. Other cell morphogenetic processes, such as cell division and cell death, should be incorporated to the current modeling scheme [62].
Plastic behavior [59] may be considered as well, either within the dissipation function formalism [1], or by considering non-linear constitutive equations to effectively incorporate a yield stress.
In analogy to the recent adaptations of the CPM [63] or the CVM [64], a cell polarity field may be included to describe collective cell migration [62]. Let each material point at the position r move to r ′ = r + u( r), thus defining the displacement field u. The center of a cell changes as r ′ c = r c + u( r c ) and the cell shape changes as Upon coarse-graining, this indicates that M changes as whereby, at order O(|∇ u|), This equation represents the relationship between the change in the cell shape and tissue displacement field, in the absence of cell rearrangements. It has been derived rigorously for a cellular material in [1].
By the virtual displacement u, the total energy F(M ) = F (M )d x changes as follows: where I represents the unit tensor and δM = M ′ − M . The first term vanishes at the boundary of the system, and the elastic stress σ e is given as: Since M is symmetric, this expression further simplifies to Using the relations and given Eq. (7) for the expression of the energy, Eq. (A3) leads to

Comparison of macroscopic and microscopic stress
To check the validity of coarsening, we conducted numerical simulation of CVM with an energy function given by Eq. (3) in the main text [2,3]: and compared two expressions of the stress tensor.
The first one is the 'microscopic' expression directly calculated from CVM [4][5][6] σ CVM = 1 where A i is the area of cell i, and ℓ ij is the length of the interface between cells i and j. P i is the pressure of i-th cell originating from cell elasticity, and T ij is the tension of cell interface [ij]. P i and T ij are determined from the above energy function as The second one is the 'macroscopic' expression of stress obtained by coarse-grained cell shape M , Eq. (A4), estimated from the given geometry of the cell in a CVM simulation. In practice, we calculated the centroid and the second moment µ i 2 of each polygon (cell) i [7], and then averaged it over N cells: The factor 4 is needed since the second moment of an ellipse with major and minor radii a and b has eigenvalues a 2 /4 and b 2 /4. With the estimated M , we calculated the cell area as π|M | 1/2 , and the cell perimeter by using Euler's formula for the ellipse perimeter to the second order [8]: Taking into accout the factor c h = 2 √ 3/π ∼ 1.05, the ratio between the perimeters of a circle and an hexagon with the same area, slightly improves the fitting. The precision of Euler's expansion for the ellipse perimeter is illustrated in Fig. S1.
CVM simulations were conducted by minimising the energy (A5). An external stress σ ex was applied on the boundary of the system, for which we took σ ex xy = σ ex yx = σ ex yy = 0, while σ ex xx was controlled to stretch the system along the x-axis. After the system relaxed, we confirmed that the force was balanced and that the stress σ CVM converged to coincide with the external stress σ ex .
In Fig. S2, we distinguish the stress that comes from cell elasticity (−P ce I, where P ce denotes the pressure) and from cell junction tension (σ T ), respectively (i.e., the first and the second terms of Eqs. (A6) and (A4)). In the simulations, parameters are set as K = 10.0, A 0 = 1.0 and κ 0 /KA 0 = 0.02, 0.04. The results are summarized in Fig. 2(a,b) in the main text and detailed in The pressure P is obtained by differentiating f with respect to A = πM 0 (equivalently P = −Tr σ e /2) Note that cosh c ≥ 1, where the equality holds at c = 0. The cell area A depends on parameters and boundary conditions. For a given cell area, the described condition does not hold for large γ 0 values, indicating that the homogenous cell size state becomes unstable. Taking higher order in approximating ellipse perimeter does not change the condition.
b. Cell shape instability f is an even function with respect to c, and is expanded as f (c) = f takes its minimal value at c = 0 as long as f ′′ (0) > 0 i.e., If γ 0 is smaller than the threshold value −4κ 0 πM 1/2 0 , the circular shape is no longer stable, and cells preferentially take an elongated shape. Using a non-dimensionalized parameter κ 0 = κ 0 /KA, the above condition can be written as This condition is unchanged when higher orderer correction of ellipse perimeter is taken into account.

Appendix B: Hydrodynamic equations of epithelial mechanics
For convenience, we list the equations describing tissue hydrodynamics derived in the main text.
a. Force balance equation reads: where f ex represents the external force.
b. Kinematic equations are given as: Here, Ω = (∇ v − [∇ v] T )/2 represents the asymmetric part of the deformation rate tensor, while D = D s + D r represents its symmetric part. Without cell division and death, D r represents the deformation rate caused by cell rearrangement, the trace of which is zero.
c. Cell number balance equation Given Eq. (B3), we calculate Here, d|M |/dt is Lagrange derivative, i.e., d|M |/dt ≡ ∂|M |/∂t + v∇|M |. The cell number density field ρ is defined by ρ = 1/π|M | 1/2 , and its evolution equation reads where Tr D r indicates the cell number variation rate. Without cell division and death, cell number density does not change through cell rearrangement and Tr D r = 0 (B4).
d. Constitutive equations are given as: where σ e and σ p are the elastic and the dissipative stress respectively, and η and η ′ are the tissue shear and bulk viscosity.
To obtain Eqs. (B9-B10) (Eqs. (15) and (16) in the main text), we set the fourth-order tensors χ sr and χ rr as follows: Here δ ij is the Kronecker tensor. χ rs ijkl is determined by Onsager's reciprocity χ rs ijkl = χ sr klij : The constitutive equations are replaced by where p represents the tissue pressure.
In the main text, we considered the case ζ 2 = β 2 = 0. Here we take into account these terms and consider their role. For active contraction-elongation with signal along n = (0, 1), the governing equations becomeċ for which Γ(c) = γ 0 /2 √ M 0 cosh c + 2πκ 0 and the incompressibility condition ∂ x v x + ∂ y v y = 0 is taken into account. With isotropic boundary condition σ xx = σ yy , ∂ x v x is given by At steady state (ċ = 0), cell shape anisotropy c is determined by The tissue exhibits steady flow, with ∂ x v x given by

Shear flow
We will consider shear flow for which three kinds of driving are taken into account. The first is a shear stress acting on the boundary, the other two are the cell-intrinsic active stresss and rearrangements. Cell vertex model simulations have shown that directed cell rearrangements may produced self-driven shear flow [9]. In addition, the properties predicted by the following analysis will give opportunities to test the model in the future.

a. Assumptions
We look for a solution with steady and uniform shear velocity gradient in the form of In the incompressible case, σ ′ e , σ ′ p , and D r are given as follows Here, for simplicity, we omit possible dependences of the coefficients χ ss , χ sr , χ rs , and χ rr on M ′ .
With an orientation along n = (cos φ, sin φ) T , the external signal reads Writing D r as d r and δ r are given as follows which determine cell shape (c, θ) for a given shear rateγ.

c. Stress
The total stress σ = σ ′ e + σ ′ p − pI reads The stress boundary condition at y = +L (top) is given as σ xy = σ b , where σ b is the force per unit length applied at the boundary to drive the shear flow. This condition reads The active stress ζ Q ∆µ sin 2φ/2 plays a role equivalent to the external driving stress σ b in the sense that it shifts σ b by a constant as σ ′ b = σ b − ζ Q ∆µ sin 2φ/2. d. Shear rate From Eqs. (C20), (C21) and (C23), we can evaluate how the shear rateγ depends on the driving stress σ b .
For the right hand sides of these equations to remain finite in the limitγ → ∞, c and θ converge to c → c ∞ and θ ∼γ −1 → 0, respectively, where c ∞ is a solution of the following equation: Considering small deviations c = c ∞ − ∆c and θ = ∆θ, Eq. (C26) reads thus ∆c is of the order of ∆c ∼γ −2 . The difference η eff − η for largeγ is evaluated from Eq. (C24), as which is of the order ofγ −2 .   [52,53]. Due to the cell rearrangement, the tissue simultaneously shrinks along the medio-lateral, y axis and elongates along the anterior-posterior, x-axis, whereas the cells adopt an elongated shape along the y-axis.