Mechano-chemical active feedback generates convergence extension in epithelial tissue

Convergence extension, the simultaneous elongation of tissue along one axis while narrowing along a perpendicular axis, occurs during embryonic development. A fundamental process that contributes to shaping the organism, it happens in many different species and tissue types. Here we present a minimal continuum model, that can be directly linked to the controlling microscopic biochemistry, which shows spontaneous convergence extension. It is comprised of a 2D viscoelastic active material with a mechano-chemical active feedback mechanism coupled to a substrate via friction. Robust convergent extension behaviour emerges beyond a critical value of the activity parameter and is controlled by the boundary conditions and the coupling to the substrate. Oscillations and spatial patterns emerge in this model when internal dissipation dominates over friction, as well as in the active elastic limit.


INTRODUCTION
Convergent extension (CE) is a morphogenetic process that occurs during development.It is conserved across many different species, types of tissues and stages of development [1][2][3][4].During convergent extension, a region of sheet-like tissue (an epithelium) elongates in one direction (the long-axis) and contracts perpendicular to the long-axis.Convergent-extension plays a key role in a variety of developmental processes, such as primitive streak formation in chick embryos [5,6] and drosophila germ band extension [7,8].The formation of the primitive streak is an important part of gastrulation, the topological inversion process shared by nearly all multicellular animals and some plants [9] that leads to cells taking up their correct positions within the embryo.
CE in epithelia is driven by cell intercalations [2,10], i.e. local cell rearrangements akin to the well-known topological T1 transitions of two-dimensional (passive) foams [11,12].Such T1s in passive systems relax stresses that build up from external driving and underlie the rheology of foams, which are typically yield stress materials [13].However in epithelia, active T1 transitions [7,8,14] can generate stresses locally even in the absence of external driving and can even develop local stresses that oppose external boundary forces.These are only possible due to motor-driven contractile stress generation, i.e. because the epithelial tissue is active.To obtain macroscopic strain against applied tension instead of a disordered response, the question then becomes how such events coordinate orientations with each other.
Until recently, the accepted answer has been a preexisting morphogenetic gene expression pattern that bias local mechanical properties [10,15,16].However, the actomyosin fibres of the cytoskeleton are themselves susceptible to mechanical feedback [17], and in e.g. the chick embryo there is no evidence for pre-patterning.Therefore, more recent work has begun to include active feedback into models of one or several coupled junctions [18-FIG.1.Schematic of a 2D sheet of tissue.The cell junctions are coloured purple/yellow/orange, to indicate the concentration of ActoMyosin: darker colour means less ActoMyosin.The active stress βM quantifies the concentration of and anisotropy of distribution of ActoMyosin in cells.The tissue is viscoelastic with a viscous relaxation time τv link to the bulk/shear viscosities ηp and ηs, and the bulk and shear moduli B and µ via τv = ηp/B = ηs/µ.20], and the response in tissues without T1s has also been investigated [21,22].
Then it becomes paramount to construct models of active tissue rheology with such feedback.While active models of cell sheets have a long tradition and include active gel theory [23] and active nematic theories [24], the focus there has been on active instabilities and topological defect motion, including in in-vitro experiments [25], but not on the response of the full tissue.Spatial patterns or direction of such feedback can be imposed [26][27][28][29] and the flow quantified, but so far a broader understanding of the emergent active relation between applied stress and strain rate, i.e. the tissue rheology, is missing.
In this letter, we present and analyse a continuum description of an epithelium with an active feedback mechanism where motor-driven contractile stress builds up, rather than relaxes, in response to applied tension as in e.g. a catch-bond [30,31].This model is the continuum counterpart to a microscopic cell junction model [32] driven by stresses generated by molecular motors (myosin-II) and cytoskeletal filaments (F-actin) which is able to generate active T1s and limited C-E flow in a tissue patch.We formulate the model in terms of the distribution of ActoMyosin (myosin-II bound to F-actin) within cells, passive viscoelastic stress and the velocity via momentum balance with a substrate (see Fig. 1).Numerically solving the continuum equations shows that above a critical activity C-E states appear, characterised by flow against externally applied stress which acts like a mechanical signal (Fig. 2).We explain this using a steady-state approximation of the feedback dynamics, where high (low) boundary stresses select a high (low) ActoMyosin fixed point in the interior, and then build up a spatial gradient leading to flow.Separately, we find oscillating and patterned states in this model in the active elastic and low substrate friction limits.To linear order, we are able to show that our equations describe an active nematic coupled to a stress field above a critical activity, but an isotropic material below it.
Model.We write down continuum equations for the epithelium as a 2D viscoelastic material that generates active stresses internally, and is coupled to a substrate via friction.The fundamental quantity of our framework is the anisotropic spatiotemporal distribution of Acto-Myosin within cells, quantified by the 2nd rank tensor M (r, t) (see Fig. 1).It is symmetric but not traceless, i.e. for cells with isotropic ActoMyosin distributions, M is proportional to the identity I. ActoMyosin, which in real tissues is distributed on the apical surface and along cell junctions, generates the stresses needed for cell-cell junction remodelling, which allows for active T1 transitions to occur.The other fields that characterise the material are local velocity v(r, t) and the local passive stress π(r, t).The total stress σ(r, t) in the material is the sum of the passive stress and an active stress proportional to M : σ = π + β(M − m 0 I), where β is the activity parameter and m 0 is the reference concentration for ActoMyosin.We use m 0 = 1/2 throughout.
The dynamics of M (r, t) is based on a model for a single contractile active junction that can remodel itself (see [32] and SI eq.S1-4).In this model, the myosin dissociation constant decreases exponentially with tension, controlled by the susceptibility k 0 .Changing the precise functional form does not affect behaviour qualitatively.In addition, the ActoMyosin tensor is convected and rotated by the flow and we write The over circle represents the corotational derivative is the vorticity tensor.We also include ActoMyosin diffusion with diffusion constant D. The cell-cell junctions within the tissue are viscoelastic [33], as is the tissue as a whole, with a time scale of stress relaxation [34].We use a convected compressible Maxwell model for the passive stress, superimposing separate Maxwell models for compression and shear deformations, where τ v is the viscous relaxation time scale, η p is the bulk viscosity, η s is the shear viscosity, and the strain rate tensor is related to the velocity field via γ = (1/2)(∇v + (∇v) T ).The bulk and shear moduli of the system are related to the relaxation time scale and the viscosities via The tissue is coupled to a substrate with friction coefficient ζ via momentum balance in the over-damped limit, Results.We integrated the equations in time using the forwards Euler method, and approximated spatial derivatives using second order accurate finite difference on a square grid.The unit of time is set by the substrate elastic relaxation time scale τ el = ζ/B, with both ζ = 1 and bulk modulus B = 1, and we use a shear modulus µ = 0.5.The myosin feedback strength is set by k 0 = 8, with the exponential form of eq. ( 1) limiting the components of M to the range 0 − 1; then the active coupling β sets the active stress scale.We use a system size of L = 50 cell units with a linear grid spacing of 0.25 and we fix D = 1.We simulate a patch inside a larger tissue by imposing a constant total stress σ ext on the boundary while material flows through freely.Note that our model is neither incompressible nor density conserving as cells can reshape in the 3rd dimension, and also can divide or be extruded [5].We complement this with the equilibrium value M = (1 + e −k0σ ext ) −1 of the ActoMyosin tensor at the boundary and invert for the passive stress π = σ − β(M − m 0 I).To study C-E, we impose pure shear boundary conditions with simultaneous tension along x and compression along y as shown in Fig. 1, i.e. σ ext xx = −σ ext yy and σ ext xy = 0, resulting in a pure shear stress σ s ≡ σ ext xx − σ ext yy .Figure 2 summarises our findings.We measure tissue response using the spatially averaged pure shear strain rate γs ≡ γxx − γyy of the steady state.We can see in Fig. 2B that at zero activity, γs = σ s /η s i.e. the tissue indeed behaves as a viscous liquid.Below a threshold activity β c = 0.5, the tissue continues to flow in the direction of applied stress, and the effective viscosity stays positive.However, ActoMyosin is built up by active tension feedback and eventually overwhelms pulling when β > β c .Fig. 2C-D show how M xx builds up in the tension direction, while M yy symmetrically drops in the compression direction.Above β c , the tissue then shows convergence extension (Fig. 2A) with the axis of elongation along the direction of compression, i.e. the tissue flows against the applied force.ActoMyosin gradients and hence tissue flow is strongest near the boundary and decays into the bulk.The C-E rheological curves above β c in Fig. 2B are highly unusual: the tissue responds to σ s → 0 with a strongly symmetry broken C-E, show-ing that the applied stress acts like a mechanical signal.When σ s increases, the C-E response diminishes, until at a β-dependent value the tissue flow reverses into the direction of pulling.For pure stretch / compression boundary conditions, the tissue also contracts / expands above β c (see SI Fig. 1-4 for full spatial profiles).
Analysis.We can understand the observed spontaneous CE by approximating the steady-state solutions of eq.(1-3).From setting • M = 0, we can derive the ActoMyosin nullcline equations where α = x, y and the off-diagonal components decay to zero (Fig. 3B).The only fixed point of the viscoelastic passive stress is π αα = 0, resulting in the transcendental equation π αα (M αα ) = 0. Below β c = 0.5, this equation has one stable solution, M αα = m 0 .Above the critical activity, there is a pitchfork bifurcation with two stable branches M αα = m + > m 0 and M αα = m − < m 0 , while the M αα = m 0 branch becomes unstable (Fig. 3A).During C-E, the equal and opposite imposed boundary stresses select a pair of points (stars) on the nullcline that break symmetry, and at the center of the tissue, we find M xx = m + and M yy = m − .The boundary conditions determine the branches in the sense that if we reverse tension and compression directions, we have M xx = m − and M yy = m + instead.Convergence extension flows are generated by the spatial gradients in stress between boundary and centre via eq.( 3).We empirically observe that the values of these stresses interpolate between boundary and centre points along the π αα (M αα ) nullcline (Fig. 3B).
We can derive an approximate solution for the C-E steady state by linearly expanding around the stable π αα (m ± ) = 0 fixed points and write and set the off diagonal components to zero.We thus eliminate the ActoMyosin equation and once we use (3) to write the strain rate, the stress equation to linear order in m x and m y becomes where constants A ± , B ± , C ± are given by SI eq.S9.If we work in the limit t τ v , we can neglect the time derivative resulting in coupled PDEs in x and y for m x (x, y) and m y (x, y).The final solution takes the form of a hyperbolic cosine in x and in y, with the full solution and derivation given in SI eq.S10-S12 and where the prefactors of each term are set by the boundary conditions.The length scale Λ + ∼ (η s + η p )/ζ, and we can derive the precise decay length (eq.SI S14) Λ d ∼ √ τ v , independent of τ m of the M αα profiles.Fig. 3C shows that the analytic approximation closely matches the numerical solution, and Fig. 3E shows that both the value of Λ d and the fact that it is τ m -independent are good prediction.The same ratio of viscosity to substrate friction then determines the penetration length of the gradient and therefore the extent of C-E flow, as can be seen in the numerical velocity profiles in Fig. 3D.We can derive an analytical prediction for the C-E strain rate γs , SI eq.S15, shown as a dashed line together with the numerics in Fig. 3F, showing the same scaling.Broader context.In phase diagram Figure 4A-B, we show that other solutions than C-E emerge in our model when we strongly increase either τ v or τ m .For very small τ v if τ m is large, the solution as expected localises near the boundaries -however the state is expanding as M drops to the m − solution throughout (Fig. S7).If we instead take the limit τ v → ∞ at small τ m , we observe pattern formation and regular oscillations in the system, including for the first time a significant M xy component (Fig. 4C-D).This is the active elastic limit where our system behaves as an elastic solid with moduli B and µ coupled to the substrate with friction ζ.It has recently been shown that active instabilities and odd (offdiagonal) responses are a characteristic feature of active elastic systems with feedback [35,36] and that they also formally arise in viscoelastic systems [37].Here we show that they arise in a model for a biological tissue, rais- ing the intriguing possibility that pattern formation in development could make use of such mechanisms.
In the limit of where τ m , τ v τ el (corresponding to the 'wet' limit where substrate friction can be neglected), we observe a spatial destabilisation of the C-E pattern with slow dynamics that we have yet to fully explore (Fig. S6).
Our governing equations bear some similarities to (and differences from) the equations of active nematic systems.However unlike in those systems, we do not find the generic instabilities usually observed.We observe instead robust and steady CE flows which is clearly very useful for biological functionality and control.While the friction with the substrate and the viscoelasticity act as stabilisers on short times and lengthscales, the key feature that keeps robust control is the interplay between the non-zero stress boundary conditions and the mechanochemical feedback of the ActoMyosin dynamics: To explore this nice feature of the model and to compare with the equations of classical nematodynamics, it is helpful to consider the dynamics of the traceless part of the ActoMyosin tensor, Q, i.e.M = Q + 1  2 Tr(M )I.By expanding the matrix exponential in equation (1) to linear order in Q, one can show that (see SI eq.S16-S27) (7) where π is the traceless part of π, π shows that there is an isotropic to nematic transition at β = 1/2, the critical activity derived from our theory.For β < 1/2, we have a stable isotropic material.For β > 1/2, we have a nematic, however the second term coming from the feedback from the passive stress at leading order resembles an applied field that will depend on the boundary conditions.This field will in general suppresses instabilities.The higher order terms will decorate this base state and can lead to a variety of interesting dynamical states, and consistent with the simulations.In summary, here we have introduced a continuum model of developmental tissues where convergence-extension flows arise wholly from mechanical feedback.We find robust C-E flows where applied tension acts like an external field to determine the flow direction, based on breaking the symmetry of spontaneous ActoMyosin polarisation.C-E then arises from the active stress profile between a central fixed point solution and the imposed boundaries.Our model also shows pattern formation and spontaneous oscillations in the active elastic limit.
Mechanochemical active feedback generates convergent-extension in a continuum model for epithelial tissue: Supplemental Material

SINGLE JUNCTION MODEL
The mechano-chemical feedback in the ActoMyosin equation is inspired by a model for a contractile active junction wirtten down by SH and colleagues.Let us consider a single junction of length l and rest length l 0 .The junction is contractile and has myosin concentration m, and is subject to an external pulling force σ ext .Assuming the passive mechanical and active component of the junction are acting in parallel, the equation of motion is where β is the activity, ζ is the friction between the junction and the surrounding medium, k is the stiffness of the junction, and m 0 is the reference value for m.The rest length l 0 of the junction undergoes viscous relaxation where λ is the viscous relaxation time of the junction.Finally, we include active feedback by making the unbinding rate of myosin decrease with tension where τ −1 m is the myosin binding rate, f (σ)/τ m is the myosin unbinding rate, and one-dimensional stress σ = k(l − l 0 ) + β(m − m 0 ).The myosin unbinding rate decreases exponentially with tension as (S4)

APPROXIMATE SOLUTION FOR CONVERGENCE EXTENSION
We have derived an approximate solution for the convergence extension steady state by expanding around the fixed points π αα (M αα ) = 0 of the mean field to linear order.We write and expand π yy (M yy (r, t)) = π yy (m − )m y (r, t), where we have used the fact that π αα (m ± ) = 0. From equation ( 3), the velocity is given by We set the off diagonal components of active and passive stress equal to zero.We look at times long compared to the myosin unbinding time τ m .We write equation ( 2) to linear order in m x and m y , which means that we can ignore the term v • ∇π αα because its lowest order terms are proportional to m 2 x , m 2 y , m x m y .The other terms in the corotational derivative vanish because the off diagonal components of the passive stress are zero.This leaves us with where we have used equation ( 3) to write the strain rate, and where If we work in the limit t λ, then we can neglect the time derivative since λ∂ t m x , λ∂ t m y → 0. This gives coupled partial differential equations in x and y for m x (x, y) and m y (x, y).Equation (S8c) suggests that we looks for solutions of the form m x (x, y) = X 1 (x) + Y 1 (y), m x (x, y) = X 2 (x) + Y 2 (y).Plugging this into (S8a) and (S8b), it is straightforward to show that X 1 (x) ∼ e ±x/Λ+ , Y 1 (y) ∼ e ±y/Λ± , X 2 (x) ∼ e ±x/Λ∓ , Y 2 (y) ∼ e ±y/Λ− , where We approximate π xx (m + ) and π yy (m − ) by differentiating equation ( 4) with respect to M αα and evaluating at m + and m − .This solution is not compatible with constant ActoMyosin on the boundary, however there is good agreement with simulations (figure 3C and 3D compared to 2C and 2D) if we apply the boundary conditions like so: yy , and M yy (0, ±L/2) = M B yy , where M B xx and M B yy are the isotropic equilibrium components of of the ActoMyosin tensor consistent with our chosen value σ ext at the boundary: We find that where .
The gradients of the passive stress are obtained by differentiating (4) This solution depends on the viscous time scale via η p = λB and η s = λµ, where B and µ are the bulk and shear moduli of the material.The velocities and strain rates are given by where σ is the total stress.We can get the dynamics for Q by taking the traceless part of (1) where ω = (∇u − (∇u) T )/2 is the anti-symmetric vorticity tensor.Writing this in terms of Q, we have From here, we will use the fact that Rearranging terms slightly, and using equation (S22), we have e −k0σ = I − k 0 (π + β(Q + m0 I)) The quantity e −k0σ • Q is also useful:

FIG. 2 .
FIG. 2. A: Steady state convergence extension velocity field at β = 0.7, σs = 0.08, and τm = λ = 20.0.B: Pure shear strain as a function of pure shear stress at the boundary for various values of activity.C: xx component and D: yy component of the ActoMyosin tensor in the same convergence extension steady state as panel A.

8 FIG. 3 .
FIG. 3. A: Pitchfork bifurcation of the mean field ActoMyosin concentration obtained via παα(Mαα) = 0 as function of activity.B: Mean field nullcline of Ṁ , παα(Mαα) (red) and corresponding simulation data, with blue dots for πxx(Mxx) along y = 0, green dots for πyy(Myy) along x = 0.The central (boundary) points are marked with a triangle (stars), and the three παα(Mαα) = 0 solutions with circles.C: xx and yy components of the simulated (solid) and mean-field (dashed) M tensor for the same parameters as fig.2C.D: Simulated vx velocity profiles in the C-E state as a function of viscous time scale and for ratio R = τv/τm ∈ [0.4 − 0.5], showing decay length.E: Decay length and F: pure strain rate in the C-E phase as a function of λ together with mean field prediction (dashed).

FIG. 4 .
FIG. 4. Phases observed for other time scale ratios.A: Observed phases as a function of viscous relaxation time scale λ and myosin time scale τm, in units of elastic substrate relaxation time scale τ el = ζ/B = 1.B: Characteristic velocity fields of convergent-extension (C-E), instability, oscillating and localised expanding states.C-D: ActoMyosin wave pattern excited in the oscillating state for the Mxx (top) and off-diagonal Mxy (bottom) components.
FIG. S1.Convergence extension phase with parameters β = 0.7, σ ext xx = 0.04, σ ext yy = −0.04,τm = 20.0,τv = 20.0.Tension is applied to the patch of the tissue parallel to the x axis, and compression parallel to the y axis.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S2.Convergence extension phase with parameters β = 0.7, σ ext xx = −0.04,σ ext yy = 0.04, τm = 20.0,τv = 20.0.Compression is applied to the patch of the tissue parallel to the x axis, and tension parallel to the y axis.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S3.Convergence extension phase with parameters β = 0.7, σ ext xx = 0.04, σ ext yy = 0.04, τm = 20.0,τv = 20.0.The patch of tissue is being pulled uniformly.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S4.Convergence extension phase with parameters β = 0.7, σ ext xx = −0.04,σ ext yy = −0.04,τm = 20.0,τv = 20.0.The patch of tissue is being compressed uniformly.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S5.Oscillating phase with parameters β = 0.7, σ ext xx = 0.04, σ ext yy = −0.04,τm = 5.0, τv = 40.0.Tension is applied to the patch of the tissue parallel to the x axis, and compression parallel to the y axis.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S6.Complex spatial pattern phase with parameters β = 0.7, σ ext xx = 0.04, σ ext yy = −0.04,τm = 200.0,τv = 500.0.Tension is applied to the patch of the tissue parallel to the x axis, and compression parallel to the y axis.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.
FIG. S7.Local expansion phase with parameters β = 0.7, σ ext xx = 0.04, σ ext yy = −0.04,τm = 20.0,τv = 20.0.Tension is applied to the patch of the tissue parallel to the x axis, and compression parallel to the y axis.The top row has the components of the ActoMyosin tensor: Mxx, Myy, Mxy, the middle row has the components of the passive stress tensor: πxx, πyy, πxy, and the bottom row has components of velocity and the xx component of the strain rate tensor: vx, vy, γxx.