Hydrodynamic theory of tissue shear flow

We propose a hydrodynamic theory to describe shear flows in developing epithelial tissues. We introduce hydrodynamic fields corresponding to state properties of constituent cells as well as a contribution to overall tissue shear flow due to rearrangements in cell network topology. We then construct a constitutive equation for the shear rate due to topological rearrangements. We identify a novel rheological behaviour resulting from memory effects in the tissue. We show that anisotropic deformation of tissue and cells can arise from two distinct active cellular processes: generation of active stress in the tissue, and actively driven cellular rearrangements. These two active processes result in distinct cellular and tissue shape changes, depending on boundary conditions applied on the tissue. Our findings have consequences for the understanding of tissue morphogenesis during development.


Introduction
During morphogenesis, epithelial tissues grow and reshape to form different organs in the adult animal. These tissue shape changes result from external stresses acting on the tissue as well as from autonomous force generation by cells [1,2,3]. Cellular forces induce cell deformations and topological rearrangements of the network of bonds joining the cells. Topological rearrangements occurring in tissue morphogenesis include neighbour exchanges through T 1 transitions, cell divisions and cell extrusions. During T 1 transitions, an edge joining two cells shrinks and two neighbors loose their contacts, resulting in a 4-fold vertex (Figure 1 A). A new bond can then form, establishing a contact between two cells which were not neighbors before. During cell divisions, a new bond is formed between the two daughter cells, and during cell extrusion, an entire cell leaves the tissue (Figure 1 A). In passive systems such as foams, topological rearrangements occur as a response to an external force deforming the system [4,5]. In biological tissues, which work out-of-equilibrium, topological rearrangements can however be internally driven by the system to generate deformation. In germ-band elongation of Drosophila embryos for instance, T 1 transitions are preferentially oriented, with cell bonds removed along the dorso-ventral axis of the tissue and added along the antero-posterior axis, and are actively driven by the cells [6,7], leading to tissue-scale reorganisation and flows.
During tissue morphogenesis, tissue deformation can either arise from deformation of individual cells of the tissue, or from topological rearrangements of the tissue: for instance, cellular neighbour exchange can result in shear by rearranging the cell positions, without overall cell deformation (Figure 1 B). The contribution of topological rearrangements to tissue flows has been measured to be a significant component of total tissue deformation in some morphogenetic processes in Drosophila [8,9,10,11]. It is however unclear how cell deformation, cellular rearrangements and the overall tissue flow are physically coupled to each other.
Different models have been proposed to describe the collective behavior of cells in an epithelium from the mechanics of single cells [12,13,14,15]. Vertex-model numerical simulations for instance attempt to capture the mechanics of epithelial tissues by representing the tissue by a network of bonds [15]. The vertices are then subjected to forces derived from an effective mechanical energy, taking into account line tensions acting on the edges of the networks and an area elasticity. Cell-based simulations of tissues however require to make specific hypothesis on cell mechanics. Alternatively, hydrodynamic theories have been proposed [16,17,18] to describe tissue growth and deformation due to cell proliferation.
In this work, we propose a coarse-grained, hydrodynamic theory for the deformations and flows at tissue scale. Hydrodynamic fields in the theory correspond to observables characterising the average cell state in the tissue (Figure 1 C). Topological rearrangements are taken into account by introducing a tensor of shear due to topological rearrangements. We construct generic constitutive equations for a polar tissue, characterising the tissue stress as well as the shear created by internally driven topological rearrangements. We show that topological rearrangements fluidify the epithelium on a timescale which depends on the amplitude of the response of shear rate due to topological rearrangements to cell deformation. Similar to active gels [19,20], an active stress can exist in the system from the forces generated inside cells by the cytoskeleton. In addition, we introduce an additional active term, distinct from an active stress, which describes active anisotropic topological rearrangements internally driven by the system. Furthermore, we discuss the properties of a novel rheological equation for the tissue, supported by experimental observations [9]. The equation introduces memory effects in the tissue on a timescale τ d , corresponding to the timescale of the dominant internal processes responsible for producing topological rearrangements in elongated cells. Finally, we discuss the dynamics of the deformation of a homogeneous  Figure 1. A) Cellular processes contributing to tissue deformation. Cell elongation: Deforming cells reshape tissue proportionally to cell area and elongation change. T1 Transition: During a T1 transition two cells lose a bond they share (red bond) and after passing trought 4-fold vertex configuration (middle) a new bond (green) is created between the other two cells taking part in T1 transition. Cell division: Cell division produces two new daugher cells from a mother cells. Cell extrusion: During a cell extrusion a single cell (red) is removed from the tissue. B) Tissue deformation can arise from changes in cell elongation (top) and from topological rearrangements (bottom). Note that in the figure topological rearrangements were represented by T1 transitions but in general they can include both cell divisions and cell extrusions. C) Cell area, cell elongation and anisotropy of internal cell properties are treated as state variables characterizing the tissue. D) Active anisotropic processes in a tissue can result in anisotropic active stress (top), or drive actively oriented topological rearrangements (bottom). Active stress and shear due to active oriented topological rearrangements produce different cell and tissue behavior (see Section 4.4).
tissue, and discuss qualitative differences between flows driven by the active stress and flows driven by autonomous topological rearrangements.

Cell density balance
We consider a tissue consisting of cells with local cell number density n. The cell number can change trough cellular events such as division and extrusions. Denoting the rates of cell divisions and extrusion per cell k d and k e respectively, the balance of cell number can be expressed as where we have introduced the cell velocity field v i .

Velocity gradient
Deformations of the tissue result from spatial inhomogeneities of cell velocity v i , described by the velocity gradient tensor The velocity gradient matrix can be uniquely decomposed into a sum of isotropic, traceless symmetric and antisymmetric terms where d is the number of dimensions. In the following we will discuss two-dimensional tissues, but extension to the d = 3 case is straightforward. The trace of the velocity gradient v kk describes local changes in tissue area, the traceless symmetric partṽ ij = (v ij +v ji )/2−δ ij v kk /d is the pure shear rate and the antisymmetric part ω ij = (v ij −v ji )/2 corresponds to local rotation of the tissue.

Cell properties
We propose a description of tissue flows at the scale much larger than the typical cell size, but we retain information on the properties of cell shapes. The cell shape is characterised by the average cell area a and the cell elongation nematic tensor Q ij (Fig. 1 B). The average cell area of a cellular patch can be defined as a = 1/n, where n is the cell number density. The cell elongation nematic has a magnitude Q characterising the strength of cell elongation, and an angle ϕ characterising the orientation of the elongation axis (see Appendix A). Different definitions to calculate cell elongation based on cell outlines have been proposed before [9,11,21,22]. In the framework of the hydrodynamic theory we propose, we expect these various definitions to slightly modify phenomenological coefficients but to leave the hydrodynamic equations unchanged. However, we impose that the tensor Q ij is defined such that homogeneous tissue deformations in the absence of topological changes give rise to the cell elongation change where D/Dt is a corotational convected derivative (see Appendix B.1).
In addition, intracellular components can be distributed anisotropically inside the cell. Proteins of the planar cell polarity pathways, for instance, are known to be distribuded across opposite cell edges [23]. Here we take into account the cell polarisation by introducing a nematic tensor field q ij that describes the orientation of anisotropic structures in the cell (Fig. 1 B). If the average cell polarity is characterised by a vector field p, the corresponding nematic tensor is obtained from

Cellular contributions to tissue flows
Isotropic and shear tissue flows can be decomposed in contributions reflecting changes of cellular properties. First, we note that Eq. (1) can be rewritten as an equation for the isotropic shear in terms of the average cell area a and cell division and extrusion rates where d/dt is the convected derivative (see Appendix B.1). Therefore, the relative change in tissue area is equal to the relative cell area change plus the relative change in cell number. Anisotropic tissue deformation stems from two sources, (1) cellular deformations which can be captured by a change of cell elongation Q ij , and (2) topological rearrangements which include T 1 transitions, cell divisions and cell death (see Fig. 1 B). Therefore, the tissue shear flow rateṽ ij can be decomposed according to the following equationṽ where R ij is a shear rate due to topological rearrangements. Note that collective correlated movements of cells can also contribute to R ij [9].

Force balance
During embryonic morphogenesis, viscous forces dominate over inertial contributions. We therefore write force balance in the low Reynolds number limit Here, σ ij is the tissue stress, and f ext is an external force density acting on the tissue.

Constitutive relations
In Eqs. 7 and 8 we have introduced the tensorial quantities R ij and σ ij , which are determined by tissue properties. We now propose constitutive relations for these quantities, using general assumptions about physical processes acting in the tissue. In addition to passive terms, we account for active processes in cells, which give rise to additional contributions in constitutive relations. In particular, they can lead to actively generated stresses and influence topological rearrangements. If cells are polarised, cell force generation can be anisotropic, giving rise to anisotropic active stress and anisotropic topological rearrangements. Both effects are taken into account in our constitutive relations through the cell nematic polarity tensor q ij introduced in Section 2.3.

Tissue stress
The total two-dimensional tissue stress can be decomposed in a pressure P and shear stressσ ij The anisotropic stress in the tissue depends on the cell elongation Q ij and cell polarity q ij . We consider a linear elastic response of the tissue stress to tissue elongation described by Q ij . In addition we introduce an active anisotropic stress capturing anisotropic force generation in the cell where the memory kernels φ K and φ ζ have units of a two-dimensional elastic modulus divided by a time. They characterise the response of anisotropic tissue stress to cell elongation Q ij and cell nematic polarity q ij ( Fig. 1 D -top), respectively. In general, the memory kernels are fourth order tensors, however, here we consider the case when all anisotropies have been accounted for by Q ij and q ij and thus the memory kernels are isotropic.
For completeness, we express the isotropic stress in the tissue as a linear response to the natural strain ln (a/a 0 ) of cell area a where K is the isotropic elastic modulus and a 0 is the cell area in a pressure-free tissue.
In what follows, we focus on the role of anisotropic stress and cell elongation.

Shear rate due to topological rearrangements
Tissue rheology is governed by cell rearrangements described by the tensor R ij . In the spirit of linear response theory, we express R ij in terms of other relevant nematic quantities, the average cell elongation Q ij and the internal cell anisotropy q ij . Taking into account memory effects, we write to linear order The memory kernels φ τ (t) and φ λ (t) have units of inverse squared time and characterise the response of shear produced by topological rearrangements to cell elongation and active cellular processes, respectively ( Fig. 1 D -bottom). Eq. 12 corresponds to the underlying assumption that oriented, anisotropic cellular rearrangements depend on the average cell elongation and the orientation of planar cell polarity. Note that similarly to the constitutive equation for the tissue stress Eqs. 10 and 11, the constitutive equation 12 and the memory kernels φ τ and φ λ depend on the material or tissue considered. We expect that passive, apolar foams can be characterised by the function φ τ , which describes how bubbles in the foam rearrange as a response to shear stress applied to the foam. In a biological process driven by anisotropic force generation on cell bonds, such as germ band elongation [6,7], we expect the function φ λ to contain information characterising how internal tissue anisotropy results in oriented topological rearrangements. The functions φ τ and φ λ therefore carry key information about cellular processes in the tissue, and can be measured experimentally in different morphogenetic tissues, in the same way that a tissue elastic modulus or viscosity can be measured in a rheological experiment.

Shear flows
Eqs. 6 to 12 constitute a system of equations which can be solved for the velocity field of the tissue v, the average cell area a and the average cell elongation field Q ij , once the rates of cell division and extrusion k d and k e , the cell nematic polarity tensor q ij , the phenomenological coefficients and memory kernels, and boundary conditions are specified. In the next section, we discuss simple limits of the hydrodynamic theory of tissue flows we propose here.

Tissue shear rheology
We first discuss the rheology of a homogeneous passive tissue subjected to external forces driving its deformation in the absence of active anisotropic cellular processes q ij = 0. Moreover, we consider the case where the memory in the response of topological rearrangements to cell deformation arises from one dominating underlying relaxation process: in that case the long time behavior of memory kernel φ τ , introduced in Eq. 12, which can also be written in a differential form This equation, together with the shear flow decomposition Eq. 7, describes the coupled dynamics of cell elongation and topological rearrangements, schematically represented in Fig. 2 A. Here, τ d is the delay timescale over which cells integrate changes in cell elongation and the changes in R ij are delayed by τ d to the changes in cell elongation Q ij . τ is a characteristic timescale of topological rearrangements, corresponding to the sensitivity of topological rearrangements to cell elongation (Figure 2 B). This form of response function was found to describe cellular rearrangements during the morphogenesis of the Drosophila pupal wing [9]. In the limit τ d → 0, the tensor of cellular rearrangements is simply proportional to the tensor of cell elongation R ij = Q ij /τ . Similarily, we discuss here a particular choice of shear stress constitutive relation in Eq. 10 where memory effects relax on a timescale much shorter than other relevant timescales. Thereforeσ where K is the anisotropic elastic modulus.
Assuming that cellular rearrangements are described by Eq. 14, a relationship between shear stress and shear flow can be derived by combining the shear flow decomposition Eq. 7, the constitutive relation describing shear due to cellular rearrangements Eq. 13 and the shear stress constitutive relation Eq. 15 For simplicity we ignore here corotational terms. Eq. 16 can be rewritten in the frequency domain asσ where η = Kτ is a viscosity andṽ ij (ω) andσ ij (ω) are Fourier transforms ofṽ ij (t) and σ ij (t), respectively (see Appendix B.2). In Eq. 17 we have introduced the frequency dependent mechanical response function χ(ω), which characterises the rheology of the tissue. The function χ(ω) is plotted on Fig. 2 D. The poles of χ(ω) are in the upper half of complex plane, as required by causality. We now discuss the form of the response function.
For zero delay timescale, τ d = 0, the response function reduces to χ(ω) = 2η/(1 + iωτ ), corresponding to a viscoelastic Maxwell material with relaxation time τ and long-time viscosity η. The viscoelastic behaviour can be understood as follows: on short time scale, tissue deformation results in cell elongation, and the emergence of an elastic stress in the tissue. On a time scale larger than τ , cell elongation is relaxed by topological rearrangements, resulting in the relaxation of elastic stress and a fluid behaviour of the tissue.
For τ d < τ /4, the poles have vanishing real parts and the system response is an exponential relaxation, similar to the case with zero delay timescale.
For large enough delay time, τ d > τ /4, an original, oscillatory rheological behaviour arises. In that limit, the poles of the response function χ have non-zero real parts and the system exhibits damped oscillations. As a result, the coupled dynamics between cell elongation and topological rearrangements (Fig. 2 A) results in an oscillatory response of the tissue. To demonstrate this, we consider the stress response σ ij (t) to a step function in imposed shear v ij (t) =ṽ 0 ij Θ(t), with Θ(t) the Heaviside step function. The resulting stress responseσ ij (t) reads for t > 0 where we have introduced β = 4τ d /τ − 1. For β 2 < 0 the stress relaxes exponentially, while for β 2 > 0 it exhibits damped oscillations (Fig. 2 C).

Representation by a simple mechanical network
We now discuss whether the rheological properties of the material described by the response function χ(ω) can be mapped to a simple rheological behaviour. We first note that the response function in Eq. 18 cannot be realised by any finite network of parallel and serial connections of springs and dashpots. This can be shown by inspecting the real and imaginary parts of χ(ω) Re Im We note that the imaginary part of the response function can become positive when τ d > τ . As we now explain, a network of springs and dashposts in series and in parallel can only have a positive real part and a negative imaginary part of the response function.
Let us consider two elements in series or in parallel in a rheological network, and assume that these two elements have the same sign of real and imaginary part of the response function at given frequency. One can verify then that the response function of the  Figure 2.
A) Schematics of the relation between cell elongation and topological rearrangements (Eqs. 7, 12). Shear due to topological rearrangements inhibits itself and cell elongation, which in turn induces shear due to topological rearrangements, forming a feedback loop with damping. Tissue shear flow can also modify cell elongation. B) Schematics of the tissue response to a sudden change in cell elongation. The onset of cell elongation triggers topological rearrangements after a delay timescale τ d , which relax cell elongation over a timescale τ . C) Shear stress response of a material subjected to a step in shear rateṽ ij (Eq. 19). Blue line: for β 2 = −0.2 < 0 the shear stress relaxes exponentially to a steady state value. Green line: for β 2 = 3 > 0, the shear stress exhibits damped oscillations relaxing to the steady state value. Red line: Case τ d = 0, β 2 = −1, corresponding to a viscoelastic Maxwell material. D) Mechanical response function χ(ω) for the three cases in C with corresponding line colors. E) A network containing a spring in series with a parallel connection of an inerter and a dashpot has equivalent rheology to a tissue with delayed topological rearrangements.
combined elements has the same sign of real and imaginary parts than the individual elements (Appendix C.1). The response function of a spring with elastic constant k is χ spring = −ik/ω, and the response function of a dashpot with viscosity η is χ dashpot = η. As a result, any combination of spring and dashpot elements in series and in parallel has a positive real part and a negative imaginary part.
However, it is easy to verify that the Laplace transform of the response function χ(t) is a positive real function. A positive real function is a rational complex function which is real for real values of s and has a positive real part for Re s > 0. The Bott and Duffin synthesis theorem [24] for electrical circuits guarantees that a positive real response function can be reproduced by a network of resistors, capacitors and inductors. By drawing a mechanical analogy to electrical networks, one can verify that similarly, any rheological network with a positive real Laplace transform of the response function can be represented by a network of spring, dashpots, and inerters. An inerter is an additional mechanical element corresponding to a capacitor in electrical networks, in the analogy where electrical current corresponds to stress and electric potential to shear rate [25]. The response function of an inerter is χ inerter = iωm where m is called inertance of the inerter and has units of mass for a two-dimensional system. Interestingly, inerters are generally omitted from rheological schemes aimed at describing tissue rheology. Indeed, because biological tissues operate at low Reynolds numbers, inertial terms associated to the mass density of the tissue can be ignored when compared to viscous forces. We find here however that a delay in topological rearrangements introduces an effective inertial term, which does not come from physical masses in the system. We find that the response function of a tissue with delayed topological rearrangements is equivalent to a circuit made of a spring connected in series with a circuit of an inerter and dashpot connected in parallel (Fig. 2 E, Appendix C.2). The inertance of the inerter is given by m = 2Kτ τ d , the effective elastic constant of the spring is k = 2K and the dashpot viscosity η = 2η = 2Kτ . In Ref. [9], we found that Eq. 17 accounted for experimental observations with τ 2 h and τ d 4 h. Assuming a typical tissue three-dimensional elastic modulus K 10 Pa and characteristic cellular length-scale l 10µm, this corresponds to an inertance of m 2 · 10 4 kg. Of course, this very large inertance, arising from memory effects in the system, is not related to the actual mass of the system.

Viscoelastic nematic gel close to equilibrium
In general, a tissue operates far away from equilibrium due to adenosine triphosphate (ATP) consumption. The constitutive equation 14 we propose here is generic and describes tissues as well as active viscoelastic gels. In a gel, the cell elongation tensor Q ij corresponds to the local elastic shear strain. Viscoelastic nematic gels can be close to equilibrium, where Onsager symmetry relations impose relations between the phenomenological coefficients relating fields and their conjugated thermodynamic forces [26].
By deriving generic constitutive equations close to equilibrium (Appendix D), we find following constitutive equation containing an additional term involving a time derivative of elastic shear strain Q ij appears in Eq. 14 In this work, we consider the simple case where α = 0. We discuss in Appendix D parameter regimes of the viscoelastic nematic gel theory where α is negligible.

Activity-induced shear flows
We now discuss the effects of active terms in specific examples of tissue deformation subjected to constant anisotropic active processes. We start from the constitutive equation for the tensor of topological rearrangements, Eq. 12, and assume that the two memory kernels φ and φ λ are exponential with the same relaxation time scale τ d . This leads to a differential equation for R ij where λ characterises the magnitude of shear flow due to active oriented topological rearrangements. Moreover, we assume that the response of tissue shear stress to tissue polarity is instantaneous and has magnitude ζ Let us first discuss a simple but instructive case of a convergent extension process in a stress free, homogeneous, polarised tissue. Two distinct active processes can drive tissue and cell deformation: active cellular rearrangements, characterised by the coefficient λ in Eq. 24, and active stress generation, characterised by the coefficient ζ in Eq. 25. At steady state, using Eqs. 7, 24 and 25, one finds that as a result of these active processes in the tissue, the tissue deforms with the shear ratẽ From Eq. 26, one finds that at long time scales, the tissue is constantly deforming, with a steady-state shear flow controlled by the active shear coefficient γ = −ζ/(2η) + λ. Therefore for γ = 0, non-zero flows are present in the steady state. Such steady-state flows can not arise in a passive system; therefore the constitutive eqs. (25) and (12) describe a passive system only when γ = 0, and in that case the magnitude of the coefficients ζ and λ are related by ζ = 2Kτ λ. We now discuss the behaviour of our hydrodynamic model in a tissue subjected to different combinations of active processes and boundary conditions. We consider a rectangular shaped, homogeneous tissue, oriented such that polarity axis is along the x axis: q xx = 1, q xy = 0 (see Fig. 3A). The tissue is constrained in space by an elastic material connected to a solid frame. The elastic material provides resistance to changes in tissue length and height. We define the natural strain variables L = ln (l/l 0 ) and H = ln (h/h 0 ), where l, h are tissue length and tissue width. Introducing the elastic moduli of the surrounding material k x and k y , the external stress acting on the tissue can be written We assume that for t < 0, the tissue is at rest, not subjected to active anisotropic processes, that cells are isotropic and that there is no stress in the elastic material surrounding the tissue. We then assume that active anisotropic processes are turned on at t = 0 and are constant for t > 0. We then solve for the tissue shape and cell

Rigid boundary conditions
Active shear stress Active oriented topological rearrangements A B Figure 3. A) A rectangular, homogeneous tissue, is attached to an external solid frame by springs. Active processes drive internal anisotropic tension and oriented topological rearrangements in the tissue. B) The polarity coupling parameters ζ and λ, corresponding to active stress and active topological rearrangements, affect the tissue differently. The deformation of a tissue with λ = 0, ζ = 0 (active topological rearrangements, no active stress) is compared to the deformation of a tissue with ζ = 0, λ = 0 (active stress, no active topological rearrangements), while keeping γτ = −ζ/(2K) + λτ = −0.2 in all cases. Other parameters: τ d /τ = 2 and in the third column, the spring constants are k x /K = k y /K = 1. Note that a Dirac delta peak at t = 0 in v xx (t) is not plotted. elongation using Eq. (24). For simplicity, we discuss here only the case k x = k y = k. In this case, the tissue area A = A 0 exp (H(t) + L(t)) is conserved. The general solution is given in Appendix E. We find for t 0; where and the corresponding dynamics of the shear rate v xx and cell elongation Q xx are plotted on Fig. 3 B for a few parameter values. We distinguish an oscillatory solution when s 2 > 0 corresponding to high values of τ d and non-oscillatory dynamics when s 2 < 0 corresponding to low values of τ d . Note that for s 2 < 0, Eq. 28 can be written in terms of hyperbolic functions (see Appendix E). When boundary conditions are stress free, ν/µ = 0, s 2 = −1 and the solution cannot be oscillatory. When imposing rigid boundary conditions, ν/µ = 1 and the parameter s is equal to the parameter β in Eq. 19. In general these two coefficients are related by (s 2 + 1)µ/ν = β 2 + 1. The factor ν/µ is related to boundary conditions; by imposing firmer boundary conditions, this factor is increased. Therefore, the oscillatory behavior disappears for decreased rigidity of the boundary springs. However, oscillatory flow does not necessarily appear at high external rigidity, due to the constraint ν/µ < 1.
Using the solution in Eqs. 28 and 31, we now discuss the cases of free, rigid, and intermediate boundary conditions, and consider the difference between flows and cell deformation induced by the anisotropic active stress ζq ij and anisotropic topological rearrangements λq ij (Fig. 3 A). We take here λ < 0 and ζ > 0, such that the active stress is larger along the x direction and horizontal bonds are actively removed. Tissue flow and cell deformation depend crucially on boundary conditions (Fig. 3B).
For rigid boundary conditions, the tissue can not deform and the active anisotropic stress has no effect. Active topological rearrangements however drive cell elongation along the x direction. The process reaches steady-state when topological rearrangements driven by cell elongation are balanced by active topological rearrangements, and a final cell elongation Q xx = −λτ is reached.
For free boundary conditions, active anisotropic stresses result in cell elongation along the y direction and tissue flow. In the steady state limit, the cell elongation reaches the value Q xx = −ζ/(2K). Active topological rearrangements do not result in cell deformation, but generate tissue shear.
For intermediate values of boundary spring elasticity, the tissue flows until the boundary springs deform sufficiently to balance stresses in the tissue. When the tissue shape reaches steady state, flows vanish and the average cell elongation reaches the same value as in the case of rigid boundary conditions, Q xx = −λτ .
In the cases described above, when boundary conditions are such that the tissue eventually stops flowing, a steady state is reached only when the shear created by topological rearrangements also vanishes. This occurs when topological rearrangements induced by cell elongation and spontaneous, active topological rearrangements balance each other. This process selects a value of the cell elongation tensor, Q ij = −λτ q ij (see Eq. 24). Therefore, by controlling active topological rearrangements, a tissue should be able to establish a fixed value of cell elongation at steady-state.
In addition, we note that although both active shear stress and active topological rearrangements can induce tissue flows (Fig. 3 B), these two processes affect cell shape differently. This suggests that the relative contribution in tissue morphogenesis of these two active processes can be distinguished by observing both cell and tissue shape changes.

Discussion
We have presented a hydrodynamic theory of tissue shear flows which explicitly accounts for topological rearrangements at the cellular scale. We have introduced active terms influencing the shear stress as well as oriented topological rearrangements of cells, coupled to the nematic field q ij . The theory applies to both effectively two-dimensional epithelia, which have been most widely studied [1,2,3], and to three-dimensional tissues, which can also undergo cellular rearrangements [27,28]. We have introduced phenomenological parameters characterising the response of shear stress to cell shape anisotropy and cell polarisation (φ K and φ ζ in Eq. 10) and shear due to topological rearrangements (φ τ and φ λ in Eq. 12). These parameters can be experimentally measured, similarly to an elastic modulus or a viscosity, and we expect that future analysis of morphogenetic processes will involve the quantification of their values.
We find that the dependency of topological rearrangements on cell shape generically leads to tissue fluidification. Without delay between cell elongation and topological rearrangement, the long-time scale tissue viscosity depends on cellular elasticity and on a characteristic time scale of topological rearrangements (Eq. 18) .
We have also introduced a delay in the response of topological rearrangements to cell shape change (Eq. 13) trough an exponential memory kernel, motivated by experimental observations in the Drosophila wing epithelium during pupal morphogenesis [9]. For delay time scale sufficiently smaller than the characteristic time scale of topological rearrangements, τ d < τ /4, the qualitative behavior of the tissue is not different from the one without delay in the response. However, sufficiently large delay time scale τ d > τ /4 leads to a novel rheological behavior, with an oscillatory response of the tissue to the imposed shear. The response of the tissue can be described by a simple rheological scheme, involving a spring, a dashpot, and an effective inertial element (Fig. 2 E). Such an effective inertial response is not generally taken into account as tissues operate at low Reynolds number, but we show that complex tissue time-dependent behaviour can indeed show an effective inertial behaviour, unrelated to the tissue's real mass.
The theory we propose here also accounts for autonomously produced stress and oriented topological rearrangements that are a consequence of active processes in the tissue. They drive tissue flows and cell elongation changes on long time scales even in the absence of external forces. In general, active stress and active topological rearrangements may occur together in tissue morphogenesis, and both drive tissue and cell shape changes. Quantification methods must be developed to characterise the effects of these active processes [9].
We have focused here on anisotropic flows. We expect that future work combining a description of isotropic and anisotropic flows in tissues, incorporating the effect of cell division and cell extrusion, will allow to obtain a full picture of tissue morphogenesis.

Appendix B.1. Convected and corotational derivatives
The convected derivative of a scalar field S is defined as The convected corotational derivative of a vector field U i is defined as and for a tensor field V ij which can be written in terms of real and imaginary parts Re χ serial = Re χ 1 Re χ 2 (Re χ 1 + Re χ 2 ) + Re χ 1 [Im χ 2 ] 2 + Re χ 2 [Im χ 1 ] 2 (Re χ 1 + Re χ 2 ) 2 + (Im χ 1 + Im χ 2 ) 2 (C.5) Im χ serial = Im χ 1 Im χ 2 (Im χ 1 + Im χ 2 ) + Im χ 1 [Re χ 2 ] 2 + Im χ 2 [Re χ 1 ] 2 (Re χ 1 + Re χ 2 ) 2 + (Im χ 1 + Im χ 2 ) 2 (C.6) By inspecting Eqs. C.2 and C.5 we conclude that if the real parts of the response functions of two elements in a serial or parallel connection have the same sign, this sign is preserved in the real part of the response function of the connection. The same is true for the imaginary parts (Eqs. C.3 and C.6).
Appendix C.2. Response function of a rheological network with a spring in series with a parallel connection of an inerter and a dashpot We calculate here the response function of the mechanical network shown in Fig. 2E.
Using Eq. C.1 we find that the response function of a parallel connection of an inerter with inertance m and dashpot with viscosity η reads The full response function is found using Eq. C.4 for a serial connection of a spring with elastic constant k and parallel connection of an inerter and dashpot Comparing with the Eq. 18 we can identify η = 2Kτ , k = 2K and m = 2Kτ τ d .

Appendix D. Viscoelastic nematic gel close to equilibrium
Here we derive hydrodynamic equations for a viscoelastic nematic gel close to equilibrium and we discuss parameter regimes which reproduce Eq. 14. Our derivation is similar to the one given in the Ref. [20], but for a one-component, nematic gel. We start by writing conservation equations for density and momentum where g i = ρv i and σ tot ij is the total two-dimensional stress. We treat the local elastic shear strain Q ij as a thermodynamic state variable. We also include in our description a nematic order parameter q ij describing other internal anisotropies, motivated by the cell polarity in tissues. Here we will discuss the case when q ij is not spontaneously generated in equilibrium. Assuming that the dynamics of elastic shear strain Q ij and nematic order parameter q ij are sufficiently slow, we introduce the free energy density of the gel as f (ρ, Q ij , q ij ).
We define the elastic shear stress as being the thermodynamic conjugate quantity to the elastic shear strain Q ij We also allow for ATP driven active processes described by the chemical reaction rate r and the chemical potential difference between ATP molecule and hydrolysis products ∆µ.
Using the free energy density f (ρ, Q ij , q ij ), we find the expression for the density of entropy production rateσ is the symmetric part of velocity gradient tensor and H ij is the nematic molecular field conjugate to q ij defined as For simplicity, we ignored here terms associated to the gradients of chemical potential and elastic shear strain, and we discuss traceless symmetric components. Identifying the thermodynamic fluxesσ ij , DQ ij /Dt, Dq ij /Dt and r, and the corresponding forcesṽ ij , −σ el ij , H ij and ∆µ, we write the phenomenological equations where Onsager symmetry relations have been taken into account. We have included ATP consumption terms coupling to fluxes of different tensorial order trough Onsager coefficients proportional to the elastic shear strain Q ij and nematic order parameter q ij . Note that since elastic shear strain is conjugate to the elastic shear stress, the Onsager coefficient relatingṽ ij to the change of shear strain Q ij can be set to 1 without loss of generality, see [20]. One can already note that Eq. D.7 is in the form of Eq. 7 if we identify R ij = Γσ el ij − H ij /β 2 − ψ∆µq ij − ψ∆µQ ij . We now show that this general description can result in the delayed response of R ij discussed in the main text, produced by the relaxation dynamics of q ij . In the remaining part of this appendix, we consider the following simple choice for the elastic stress and molecular field σ el ij = 2KQ ij (D.10) Let us first consider the case of a passive gel which would not consume ATP. Using Eqs. D.7 and D.8, the tensor R ij can be expressed as a function of the shear strain Q ij We find that the choice τ d = 1/[γκ(1 − β 1 /(γβ 2 ))] and τ = [1 − β 1 /(γβ 2 )]/[2K(Γ − 1/(β 2 2 γ))] allows to identify Eq. D.12 with Eq. 14. The last term on the right-hand side of Eq. D.12, involving the derivative of Q ij , is not present in Eq. 14. The dimensionless prefactor in front of this term, α = [2KΓ + κβ 1 /β 2 ]/[γκ(1 − β 1 /(γβ 2 ))], can be set to 0 if β 1 = −2KΓβ 2 /κ. This choice implies the following inequality where we have used γ, κ, K, Γ ≥ 0 and Γ − 1/(γβ 2 2 ) ≥ 0. Therefore, if the parameter β 1 is chosen such that α is negligible, the gel can never be in the oscillatory regime τ d > τ /4 described in Section 4.1.
A finite value of β 1 implies a dependence of the shear stress on the molecular field H ij . If we assume that the main contribution to the shear stress comes from the viscous and elastic stresses and other contributions are negligible, one would set β 1 = 0. In this case it is no longer possible for α to be arbirtrarily small. Indeed, due to the positive semi-definiteness of Onsager coefficient matrix Γ − 1/(γβ 2 ) ≥ 0, the coefficient Γ cannot be arbitrarily small. For Γ < 1/(γβ 2 2 ) the timescale τ would become negative and the thermodynamic equilibrium state would become unstable. Moreover, a lower limit to the prefactor α is given by (D.14) Other possible choices of parameters which would allow neglecting the term involving the shear strain are K → 0 and γκ → ∞. However, both choices inevitably remove a term in Eq. (D.12), in such a way that the form of Eq. (14) can not be reproduced. Therefore, a passive gel with only viscous and elastic stresses can reproduce Eq. (14) only with an additional term involving time derivative of the tensor Q ij .
Since the shape of the tissue is constrained to be a rectangle, the velocity gradient tensor has only diagonal components: Combining these relations with the constitutive relation for the isotropic stress in Eq. 11, we obtain a relation between L and H from the isotropic shear component where we have used the fact that in the absence of cell division and extrusion, (da/dt)/a = (dA/dt)/A with A being the area of the tissue. Using the constitutive relation for the shear stress Eq. 25, the constitutive equation for the shear due to topological rearrangements Eq. 24, and the shear decomposition Eq. 7, we obtain a second relation between L and H: where γ = −ζ/(2η) + λ. Integrating Eq. E.3 from an arbitrary lower bound t < 0 yields We now express H in terms of L in Eq. E.4 and we obtain the second order equation where L(0) and ∂ t L(0) are initial conditions, and The parameter s 2 is positive only for high enough values of the memory timescale τ d . When s 2 becomes negative, the solution is equivalent in form to Eq. E.9, with trigonometric functions replaced by their hyperbolic counterparts and s 2 replaced with −s 2 . We now determine the initial conditions for an experimental setting in which the tissue is initially at rest, cells in the tissue are not elongated, external elastic connections are not under tension, and there is no tissue polarity. At t = 0 the polarity is activated on a timescale much shorter than τ and τ d , so that we can treat the activation as instantaneous. First, considering Eq. 24 we can conclude that R xx (0 + ) = 0. Then, using Eqs. E.1, E.2, E.3, 7 and 25 we can show that for t < 0 + (E.11) Then, using ∂ t q xx = δ(t) and integrating over a small finite time interval around t = 0, we obtain the initial conditions L(0 + ) = − ζ µK ∂ t L(0 + ) = 0 .

(E.12)
When the elastic connections around the tissue are isotropic (k x = k y ), one obtains from Eqs. E.5, E.7 and E.8: Since L(t) + H(t) = 0, the tissue area is constant. Therefore, area changes arise only when the surrounding material is anisotropic. Using these relations, the initial conditions discussed above and Eq. E.9, we obtain Eq. 28 in the main text. Finally, we note that in the limit τ d → 0, Eq. E.9 describes a simple exponential relaxation of the tissue length: (E. 16) with a characteristic relaxation timescale equal to (kτ )/(2K + k) for isotropic external springs. It differs from the Maxwell timescale of the tissue τ by a factor describing a competition of internal and external elasticities.