Hydrodynamic instabilities, waves and turbulence in spreading epithelia

Abstract

We present a hydrodynamic model of spreading epithelial monolayers described as polar viscous fluids, with active contractility and traction on a substrate. The combination of both active forces generates an instability that leads to nonlinear traveling waves, which propagate in the direction of polarity with characteristic time scales that depend on contact forces. Our viscous fluid model provides a comprehensive understanding of a variety of observations on the slow dynamics of epithelial monolayers, remarkably those that seemed to be characteristic of elastic media. The model also makes simple predictions to test the non-elastic nature of the mechanical waves, and provides new insights into collective cell dynamics, explaining plithotaxis as a result of strong flow-polarity coupling, and quantifying the non-locality of force transmission. In addition, we study the nonlinear regime of waves deriving an exact map of the model into the complex Ginzburg–Landau equation, which provides a complete classification of possible nonlinear scenarios. In particular, we predict the transition to different forms of weak turbulence, which in turn could explain the chaotic dynamics often observed in epithelia.

---

1 Introduction

In recent years, in vitro epithelial cell monolayers have become a key model system to investigate mechanical aspects of collective cell migration, a generic situation that is directly relevant to a variety of biological processes in living organisms, including morphogenesis^1–3 or regeneration.^4–6 In particular, much attention has been focused on the collective mechanisms by which cohesive advancing cell sheets are capable of transmitting and building up intracellular stresses over distances of hundreds of microns.⁷ Cells are able to exert actively driven forces to a substrate underneath and migrate towards the maximum principal stress direction, and simultaneously the instantaneous monolayer stress maps may trigger signaling pathways that affect the mechanical state of individual cells.^8,9 This suggests a strong interplay between the physical properties of a tissue and the internal structure of the constituting cells.^10,11 Consequently, it becomes crucial to develop a solid theoretical framework to interpret the force and kinematic maps currently available for spreading epithelial monolayers. From a physical standpoint, a key point is to elucidate to what extent the observed phenomena, even if strongly regulated biologically, can be tackled in purely mechanical terms. In this context, in recent years, there has been an increasing interest in designing experiments in vitro to probe the mechanics of epithelial tissues in controlled situations.^9,12–16

One issue that remains a matter of debate on the theoretical side is whether a continuum description of a tissue must assume a viscous^17–20 or an elastic^21–24 constitutive equation in a given range of time scales of observation, in particular for the long-time regime. This question is nontrivial for living tissues in particular because a given type of cell may respond differently in different environments on the same time scale. For instance, MDCK cells in suspended monolayers under external pulling seem to respond elastically¹⁴ on time scales for which a freely spreading monolayer on a substrate seems to be flowing like a viscous fluid.²⁵ In addition to such phenotypic variations of the mechanical properties of cells in response to the environment, additional confusion may arise when comparing the mechanical response of the tissue to an external force as tested in Harris et al.,¹⁴ to the relationship between the stress and strain variables when the stress is autonomously induced by the tissue, implying that the two observables may be related by some additional constraint, either biological or otherwise, that prevents from establishing a direct causal, stimulus-response relationship between both. A clear example of this point is the observation that both stress and strain in the central region of a spreading monolayer have been shown to grow linearly with time in some initial range. This obviously allows establishment of a linear relationship between stress and strain which in turn yields an effective elastic modulus.²⁶ However, the spatial extension of the material may well be similar to that of a fluid that is accomodating to a moving boundary (the leading edge), while the internal active stress may be growing due to some internal process with no direct causal relation with the advanced speed of the leading edge. Consequently, the fact that an effective elastic modulus can be defined is not necessarily informative of the passive rheology of the material. Indeed, it can be shown that the same data of Vincent et al.²⁶ can be made consistent with a purely viscous constitutive equation for an active fluid, once taken into account that the effective viscosity as the monolayer spreads is generically time-dependent. This point is discussed in Blanch-Mercader et al.²⁵

In this context, our study is directly motivated by the experiments of Serra-Picamal et al.¹² on spreading epithelia where ultra-slow elastic-like waves have been reported, on time scales of several hours, where one could argue that on the basis of the time scales of intracellular processes (∼1 minute) and cell–cell adhesion kinetics (∼10 minutes) one expects viscous behavior. Independent observations of individual displacement of cells and their relative sliding also suggest that the relaxation of stresses is fluid-like. However the same experiments reported a phase lag between stress and strain-rate measurements that is characteristic of elastic waves. The different attempts to model these phenomena so far (for instance, by Serra-Picamal et al.¹² and Banerjee et al.²³) have assumed an elastic constitutive equation for the tissue and different additional hypotheses to account for the emergence of waves.

In this paper we present a continuum model of a epithelial monolayer spreading on a substrate that is based on a viscous constitutive equation for the medium and combines two sources of active stresses: bulk contractility and traction forces at the contacts with the substrate. We will elucidate a hydrodynamical instability that can explain the emergence of elastic-like waves in the range of time and length scales of the experiments. The mechanism and the physical scenario that account for the waves are completely different from those of elastic models, and the waves exhibit distinctive features with no counterpart in those models. In particular we will discuss the observations of the experiments from Serra-Picamal et al.,¹² Vedula et al.,¹³ and Deforet et al.¹⁵ in light of our approach. Our model is based on the phenomenological continuum approach of active gel theory,^27,28 where the medium is treated as polar and the equations are imposed solely by symmetry considerations and linear irreversible thermodynamics. Most of the phenomenological parameters of the theory can be estimated from independent experimental observations and those with no direct evidence will be indirectly inferred with the use of the model. In particular the large values obtained for the flow alignment coefficient will provide interesting insights into the phenomenon of plithotaxis.⁹

The layout of the paper is as follows. In Section 2 we present and discuss the continuum model. In Section 3, we discuss the linear stability analysis of a homogeneously polarized state. Section 4 is devoted to the nonlinear regime, and includes both numerical simulations and the mapping of the problem into a complex Ginzburg–Landau equation. In Section 5 we revisit some experiments on MDCK cell monolayers, and interpret their results in light of our theoretical framework. Finally we sumarize our results in the concluding Section 6.

2 Continuum model

Cells are assumed to have a macroscopic polar order described by the vector field p. At the free-edges of expanding cell sheets, they tend to develop lamelipodium-like structures that require a globally oriented actin cytoskeleton, while cells that are hundreds of microns away from the interface extend basal cryptic lamellipodia underneath the neighboring cells.²⁹ In addition, epithelial cells exhibit other types of complex structures, such as stress fibers connecting them. From a coarse-grained point of view, to mimic the tendency of the polarity field to align with the neighbors and thus avoid large polarity gradients, we introduce an effective free energy for these degrees of freedom of the form of the standard free energy of a polar nematic,³⁰where we use the Einstein summation convention. The polynomial part favors the emergence of a finite polarity vector p_α of modulus p = 1, while the second term penalizes energetically the formation of large gradients. The energy scale of the nematic elasticity is fixed by the parameter ρ > 0, which in our 2d model has dimensions of stress. The balance between the two terms defines a characteristic scale of spatial variation of the polarity, the so-called nematic correlation length . The conjugated field of the polarity, the so-called molecular field, is given by h = −δ[scr F, script letter F]/δp. In other situations where cells are not globallly migrating and adopt in-plane elongated shapes, an apolar nematic order has been reported,^31,32 suggesting that complete characterization of these systems may require the use of two ordering fields with both polar and apolar symmetries.³³

A key assumption of our model is that, at sufficiently long time scales, the medium can be described as a viscous fluid. A more realistic description of the rheology of these materials need to address in detail the inter- and intra-cellular mechanics,¹¹ but this remains beyond the scope of our study and henceforth we omit those effects. Our focus here is rather on counterposing the viscous-like to the elastic-like character of the constitutive equation. A viscous hypothesis is consistent with the direct observation of the fluid-like relative motion of cells on those time scales, and assumes that all the processes that control the elastic properties of the medium relax on much shorter time scales. The kinetics of the cell–cell adhesion, for instance, has a turnover time scale that is in the range of ∼10 min,³⁴ while the time scales of observation for the phenomena studied here are at least one order of magnitude larger. Consequently, the viscous hypothesis seems at least plausible for spreading monolayers, even though there is no direct fully-conclusive evidence. Cell division has also been shown to generate a fluidization mechanism that would support further the viscous scenario, together with the introduction of another source of active stresses.³⁵ However, in the experiments of Serra-Picamal et al.¹² cell proliferation is not very significant at the time scale of observation of 100 min so in our model we ignore it for simplicity, neglecting this source of active stresses.

Since the cell monolayer is a quasi-two dimensional system, we will assume an effective 2d description, extending the approach of Blanch-Mercader et al.²⁵ We thus take the hydrodynamic equations describing a 2d (compressible) active polar fluid, with nematic elasticity, that are consistent with symmetries and include active and passive contact forces with the substrate. Our model is completely specified by the set of equations

∂_βσ_αβ = ξv_α − T₀p_α, (2)

h_α = ρ(1 − p²)p_α + K∇²p_α. (5)

where σ_αβ is the traceless stress tensor, and v_αβ and ω_αβ are the traceless symmetric and antisymmetric components of the velocity gradient tensor, respectively.

Eqn (2) expresses the force balance in the absence of inertia, with the total external force on the rhs, including a passive contribution describing friction with the substrate, and an active one describing the traction force.

Eqn (3) is the constitutive equation for the total stress of an active polar liquid. It is assumed that on the time scales of observation, elastic effects can be neglected. The first term on the rhs accounts for viscous stresses with η being the shear viscosity. The second term accounts for active stresses, with the activity parameter ζ < 0 for contractile stress, which is the relevant case for epithelial monolayers. The rest of the terms are the usual ones describing nematic elasticity.³⁰

Eqn (4) describes the dynamics of the polarity field. The lhs is the total co-moving co-rotating derivative, and the rhs describes the rotational relaxation of the polarity, γ₁ being the rotational viscosity. The last term, which couples the polarity and the flow, is the so-called flow-alignment contribution, and by virtue of the Onsager reciprocity relations must be characterized by the same coefficient ν appearing in eqn (3).

Finally, eqn (5) specifies the molecular field h in terms of the polarity consistent with eqn (1).

In addition to the active contractility, we have introduced an additional active term that accounts for traction forces exerted on the substrate. The proposed form of this external force has been used before in the context of cell monolayers,^24,25,36 as well as for bacterial suspensions.³⁷ The form is limited by symmetry considerations³⁸ and can be derived from microscopic models with linker kinetics.³⁹ Unlike previous continuum models of active polar gels,^40–42 the polar term is in the cell–substrate interactions rather than the material constitutive equations. As we demonstrate below, active traction forces introduce important qualitative changes in the dynamics of the system. The form of the friction force, with an effective friction coefficient ξ is standard.⁴³ The active traction must be related to the consumption of ATP, as the active contractility term, but the two parameters T₀ and ζ originate at different mechanisms and can be considered as independent parameters. T₀ > 0 defines the scale of the cellular traction forces.

An effective 2d model such as the one proposed can be derived from a 3d thin layer using the lubrication approximation and averaging over the monolayer thickness, as described in Blanch-Mercader et al.²⁵ for a 1d reduction. The 3d incompressibility of the fluid then allows eliminating the pressure from the description. In the reduced description, however, the fluid is compressible, and for a reduction to 2d, the constitutive equation for the trace of the stress tensor σ must be specified. For simplicity we assume

with the three coefficients [small eta, Greek, tilde], [small zeta, Greek, tilde], [small nu, Greek, tilde] each being zero. Strictly speaking, this choice gives up 3d incompressibility, since these parameters are not free within that condition. However we may still omit pressure gradients in the force balance eqn (2) by assuming that the 2d effective fluid has a large compressibility, which accounts for the fact that an in-plane compression offers no significant resistance because it can be accommodated by an expansion in the third dimension.

The parameters of the problem can be grouped into different useful combinations. In addition to the nematic correlation length , we introduce two friction lengths defined by and (with γ ≡ γ₁/ρ). We use the nematic stress scale ρ to define a dimensionless contractility [small zeta, Greek, macron] = ζ/ρ and a dimensionless traction T̄₀ = T₀L_γ/ρ. In addition to ν, [small zeta, Greek, macron] and T̄₀, the model contains the three length scales L_c, L_η and L_γ and one time scale γ. Using two of them to scale length and time, we are left with a set of five independent dimensionless parameters.

3 Linear instability of a homogeneously polarized state

3.1 Linear stability analysis

The set of eqn (2)–(5) have a trivial homogeneous solution with a uniform polarity field and a uniform velocity V in the direction of the polarity field, with V = T₀/ξ. The direction of these fields has the continuous degeneracy of the rotational invariance of the problem, which is spontaneously broken when the direction of the polarity is fixed. Without lack of generality, we chose the polarity and the velocity to be oriented in the x-direction so that p_x = 1, p_y = 0, v_x = V, v_y = 0. In the particular case T₀ = 0 we recover the uniform state from previous studies.^40,41 In this section we address the linear stability analysis of this base state, that is, we will obtain the growth rate ω(q) of sinusoidal perturbations of the form exp(ω(q)t + iq·r) under the linearized dynamics around the base state. We distinguish two types of modes, transverse and longitudinal, which designate perturbations of the polarity with q parallel to the x-direction that are perpendicular (δp_x = 0) and parallel (δp_y = 0) respectively, or equivalently, that modulate the direction and the modulus, respectively. At the linear level these two types of modes are decoupled, so we find two branches for the linear growth rate ω^L and ω^T, which take the formFrom the anisotropy of the problem, this growth rate depends on the modulus q² ≡ q_x² + q_x² and on terms of the form qⁿ cos nθ where θ is defined by cos θ = q̂·x̂. To make the notation more compact we introduce the complex wave vector Q ≡ q_x + iq_y, such that Re[Qⁿ] = qⁿ cos nθ. Equivalently, we have Re[Q²] = q_x² − q_y², Re[Q³] = q_x(q_x² − 3q_y²), and Re[Q⁴] = q_x⁴ − 6q_x²q_y² + q_y⁴. We can now split the contributions to the dimensionless growth rate ω/Ω according to passive, contractility and traction terms, such that we have

λ = λ_p + Re[λ_ζ] + iRe[λ_T0] (8)

If we adopt the compact notation q_c ≡ L_cq, q_η ≡ L_ηq and q_γ ≡ L_γq, and Q_γ ≡ L_γQ, we have
and

The above expressions contain a wealth of physical information about the dynamics of the system already at the linear level. The spontaneously broken isotropy is captured by the terms containing Q. The prefactor Ω(q) in eqn (7) is a Lorentzian propagator that expresses the nonlocal character of the dynamics (see also discussion of Section 3.3). In addition, the contributions contained in are also nonlocal, since they involve all orders in q. Note also that, because of the symmetries of the problem, the linear growth rate satisfies ω(−q) = ω*(q), where the asterisk denotes complex conjugate, ω(−q_y) = ω(q_y) for any q_x, and, for q_x = 0 we have Im[ω] = 0. For further reference it is interesting to write explicitly the growth rates in the long wavelength limit, as an expansion up to order q². These read

The hydrodynamic limit q → 0 is fundamentally different for the two types of modes. The transverse zero-mode is marginal, i.e. ω^T(q = 0) = 0, because of the rotational invariance of the problem, since this mode accounts for an infinitesimal homogeneous rotation of the base state, which has a continuous degeneracy. This soft mode associated with the rotational symmetry plays an important role in the discussion of possible routes to chaos in this problem in Section 4. By contrast, the longitudinal modes relax in a time of order one in this limit, and hence they are not hydrodynamic. The separation of time scales between the relaxation of the modulus of the polarity vector and its direction is often invoked to justify an adiabatic elimination of the dynamics of the modulus. In our case, however, the region of interest for the experiments is not in the hydrodynamic limit and longitudinal modes are essential to interpret the observations.

3.2 Bifurcation into travelling waves

We find that, in general, both λ and D are complex numbers, so the growth rates ω^L,T contain both a real and an imaginary part. The condition Re[ω] = 0 defines an instability boundary. If Im[ω] ≠ 0 when Re[ω] = 0, then the instability is oscillatory (i.e. a Hopf bifurcation), allowing for travelling waves.

From the explicit expression of the growth rates, we see that the passive contribution to Re[ω] is always negative, which makes the base state stable in the absence of active forces. However, if the system is sufficiently active, there may exist modes with Re[ω] > 0 and thus an instability sets in. Hereinafter we will restrict ourselves to the case ν > 0 which is the one relevant for the discussion of spreading epithelial monolayers. Then, for longitudinal modes, the instability can only happen for contractile active stress (i.e. ζ < 0) and for a band of modes with finite q around the critical wave number q_o for which Re[ω^L(q_o)] = 0, excluding the neighborhood of the mode q = 0, which is always stable. An unstable band thus exists for |[small zeta, Greek, macron]| > |[small zeta, Greek, macron]^L_o| where the threshold value of contractility for the onset of the longitudinal instability is given by

In this case, the onset of instability occurs at a finite q_o ≠ 0. This situation is usually referred to as an oscillatory periodic instability (see Cross and Hohenberg⁴⁴). For the 1d case (i.e. q_y = 0) the critical mode where the instability sets in isSince Re[ω^L(q)] is arbitrarily small close to the threshold, while the frequency Im[ω^L(q)] remains finite, slightly above the threshold a localized perturbation is advected faster than it grows, so the perturbed fields eventually relax to the unperturbed values, a situation that is referred to as convective instability. Only over a finite distance from the threshold is the system said to be absolutely unstable, that is a localized perturbation at a given location does grow in amplitude at that location. Note that for the longitudinal modes, the instability requires a finite value of the flow alignment coefficient ν. Large values of ν, and small values of friction favor the instability.

For the transverse modes, the instability is controlled by both active parameters. The condition for the instability then reads

Remarkably, in this case, the instability may be driven solely by traction forces and may occur even for extensile activity (i.e. ζ > 0). The transverse instability is a long-wave length one, that is with the critical wave vector q_o = 0. At a finite value above the threshold, Re[ω^T] is peaked at finite q but remains marginal at q = 0, as imposed by the rotational invariance of the problem. Note also that the transverse instability sets in before the longitudinal one as we increase |ζ|, since condition (14) is satisfied for [small zeta, Greek, macron]^L_o < [small zeta, Greek, macron], even if T̄₀ = 0.

We remark that, for any q_x ≠ 0, both longitudinal and transverse modes have a non-zero imaginary part of ω^L,T provided that T₀ ≠ 0. Therefore, the presence of traction forces is directly associated with oscillatory behavior. In addition, the sign of Im[ω] is always opposite to that of q_x, implying that waves propagate only in the direction of the polarity p (from negative to positive x in our case), reflecting the anisotropy of the base state. Similarly, the fact that for q_x = 0 we have Im[ω] = 0 implies that waves travelling along the y-axis are not possible.

The linear instability criterion does not exclude that finite-amplitude nonlinear waves or other types of solutions can exist even in the linearly stable regions, in particular near the linear stability boundary. This nonlinear instability of the base state, is characteristic of subcritical bifurcations and implies that sufficiently large finite-amplitude perturbations of the base state may grow even if smaller amplitudes do decay. In the analysis of the following sections we will determine exactly the boundaries where the bifurcation in our model is subcritical.

In Fig. 1 we plot the real and imaginary parts of both branches ω^L,T, for parameter values relevant to MDCK cells.^12,13,15 We observe that the transverse instability appears first and the growth rate is peaked at finite q_y and q_x = 0. By contrast, the longitudinal instability peaks at finite q_x and q_y = 0.

Fig. 1 Linear growth rate of a small perturbation of a homogeneously polarized state. The upper (lower) row displays the positive real (imaginary) part of the growth rate in the (q_x,q_y) plane. The left (right) column represents the growth rate of the longitudinal (transverse) modes, with color-coded amplitude. The values of parameters are η = 10⁶, ξ = 100, T₀ = 10, L_c = 50, γ = 600, ν = 10, ρ = 10 and ζ = −2000, in units of Pa, μm and s.

The explicit expressions for the growth rate on the x, y-axes take a simpler form. For q_y = 0 and writing q ≡ q_x and ω_x(q) ≡ ω(q_x), we have

According to these results, the characteristic frequency scale of the oscillations w_o = |Im[ω(q_o)]| is given by ω_o ∼ T₀q_o/ξ. Therefore, the origin of the slow time scale of oscillations is not intrinsic of the material properties but depends on the interaction with the substrate, which fixes the friction and traction forces. In Section 4 we will see that the observed waves in experiments can indeed be interpreted as those resulting from this instability.

The phase velocity of a wave with wave-vector q_x for the two types of modes is given by v^L,T = |Im[ω^L,T_x]|/q_x. Using that the flow velocity of the base state is V = T₀/ξ, we find that

This result allows us to infer a value of the parameter ν, which is usually unknown for tissues, just comparing the wave velocity and the front velocity in the experiments of Serra-Picamal et al.¹² (see Section 4).

In the y-axis, q_x = 0, we have Im[ω] = 0. The presence of traction quantitatively modifies the growth rate but this is qualitatively similar to that for the case T₀ = 0, which reads, with q ≡ q_y,

3.3 Physical origin of waves and phase lag

In most of the subsequent analysis, we will pursue the study of the 1d-case corresponding to longitudinal modes with q_y = 0 (eqn (15)). This is the simplest case and at the same time the most interesting to gain insights into the physical mechanism behind the waves, to analyze in depth the nonlinear dynamics of the problem, and to compare with experiments.

In this case the model equations reduce to

∂_xσ_xx = ξv_x − T₀p_x, (21)

h_x = ρ(1 − p_x²)p_x + K∂_x²p_x. (24)

We remark that the emergence of elastic-like waves appears naturally in our model for a purely viscous constitutive equation as long as active traction is present, without invoking any additional time scale related to an extra coupling to internal variables, as is usually required in models based on elastic constitutive equations.^12,16,23 For the sake of discussion, let us consider the simple case of ν = 0 and ζ = 0. At the linear level, the evolution of a small perturbation δp is coupled to traction through the advective term,

with V = T₀/ξ, implying that a perturbation of the polarity proportional to exp(ω(q)t + iq·r) will decay but at the same time be advected with the velocity V = T₀/ξ. We thus have a dispersion relation with Im[ω] = −Vq. This is reminiscent of that of elastic waves (i.e. Im[ω] = ±V|q|), but with the fundamental difference that we get only one propagation direction, due to the fact that rotational symmetry (or parity in the 1d case) is spontaneously broken in the base state. The key observation is that we can close an equation for δp that is of first order in time derivatives. This is in contrast to the case where the medium is assumed elastic. As pointed out in Notbohm et al.¹⁶ and Banerjee et al.,²³ in that case an effective inertia must be invoked to obtain elastic waves. In Banerjee et al.,²³ this is achieved by introducing an additional coupling to a slow variable, such that, at the linear level, a wave equation (i.e. with second order time derivatives) is obtained for the strain field. Note that the resulting elastic waves are apolar, that is, insensitive to the sign of p. In a purely viscous medium, however, we do obtain propagating waves through the advective coupling, as long as the force balance equation includes a finite traction T₀, and these waves are polar. The other observables can be directly related to δp obtained from eqn (25). Note that the combination of the force balance equation and the viscous constitutive equation implies that the relationship between δp and both δv and δσ is nonlocal withand δσ(x) = η∂_xδv(x).

The scenario for the emergence of waves in our model is thus fundamentally different from that provided by an elastic medium with an effective inertia. Our scenario is not essentially modified if we reintroduce non-zero values of ν and ζ. Then the closed equation for δp becomes nonlocal but still of first order in time derivatives.† The physical picture is essentially the same with two important additional features. First, the presence of contractility can reverse the damping and generate the growth of the wave amplitude, thus allowing for sustained (nonlinear) waves. Second, the presence of either one of the two parameters is sufficient to introduce an elastic phase shift between the stress and strain rate. The presence of such a phase lag in the experiments of Serra-Picamal et al.¹² and Notbohm et al.¹⁶ has been usually interpreted as a signature of an elastic constitutive equation of the medium. Here we show that this inference is not really justified, since it is also possible to have the same phase lag in a purely viscous medium, provided that active traction is present (T₀ ≠ 0), together with at least one of the two parameters ν or ζ being nonzero.

Solving the linearized equations eqn (22)–(24) in Fourier space, the relative phase between the stress and strain rate can be found exactly and reads

The rhs of eqn (27) being real would be characteristic of a purely viscous medium. The presence of traction forces T₀, however, introduces an imaginary part that produces a phase shift, which would be typical of elastic behavior. Note that the presence of nematic elasticity in the constitutive equation (i.e. the term proportional to ν in eqn (22)) would not be sufficient to introduce the elastic-like phase shift in the absence of traction.

In the particular case of 2ν(1 + q_c²/2) ≪ |[small zeta, Greek, macron]| and |[small zeta, Greek, macron]|q_γ ≪ T̄₀, eqn (27) takes the simple form

The phase lag in the spatial dependence will produce a phase lag in the time oscillations of the two observables at any given location, provided that travelling waves are sustained, that is, for |ζ| > |ζ^L_o|, and consequently q ∼ q_o. We thus obtain that, if |ζ^L_o| ≫ ω_oη = T_oq_oη/ξ the phase lag for the observed waves will typically be that of an elastic medium, even though the passive rheology of the system is that of a purely viscous material.

4 Nonlinear waves

In the previous section we have analyzed the linearized dynamics around the homogeneously polarized state, and have identified broad regions of parameters where this state is unstable. The amplitude of unstable modes will grow until saturation by nonlinearities. As we will show, the observed waves in a variety of experiments can be identified with such nonlinear waves. In this section we pursue the numerical and analytical study of these nonlinear waves. We will focus mostly on the nonlinear analysis in the case of 1d longitudinal waves, and address only briefly more general situations.

4.1 Numerical analysis of longitudinal waves

The numerical exploration of the five-dimensional parameter space of our problem is obviously prohibitive, so in a first tentative study of the highly nonlinear regime we will construct a phase diagram ζ − ρ, relating the main parameter driving the instability and the stress scale of the nematic elasticity, a phenomenological parameter which is rather elusive for cell tissues. With the rest of the parameters fixed, we have that ρ ∝ L_γ².

We have integrated numerically the full nonlinear dynamics of the 1d model for longitudinal modes, eqn (21)–(24), using a semi-implicit algorithm. For parameter values in the range relevant to experiments, we have identified different types of solutions resulting in the phase-diagram plotted in Fig. 2. Beyond the threshold of the longitudinal oscillatory instability (i.e. |ζ| > |ζ^L_o|) we find three classes of nonlinear solutions. In domain A, we observe that the polarity field, after a random perturbation, develops a transient array of localized pulses, spaced by the typical length set by the most unstable mode, which after some time coalesce giving rise to the formation of an isolated pulse of polarization propagating through a nonpolarized medium, as shown in Fig. 3. The maximal value of the polarity field is of the order of and the transition from the homogeneously polarized state is discontinuous, implying that the propagating solutions appear at a finite amplitude right on the threshold (subcritical bifurcation). In domains B and C, the physical observables exhibit a nonlinear travelling periodic pattern. In both cases, the travelling wave speed is of the order of T₀/ξ and the spatial periodicity is of the order of 1/q_o. The transition from the homogeneous polarized state is continuous (supercritical) for domain B and discontinuous (subcritical) for domain C. The rest of the diagram, plotted in yellow, corresponds to the linearly stable region.

Fig. 2 Phase-diagram of different nonlinear dynamics in the ζ − ρ plane. The yellow domain corresponds to the linearly stable region while the other domains correspond to nonlinear oscillatory solutions displayed in Fig. 3. The boundary of the yellow domain is given by the longitudinal instability threshold ζ^L_o. The other borders are determined by the complex Ginzburg–Landau eqn (30). The parameters are η = 10⁶, ξ = 100, T₀ = 10, L_c = 50, γ = 600 and ν = 20, in units of Pa, μm and s.

Fig. 3 Representative steady oscillatory nonlinear profiles for the different domains of the phase diagram of Fig. 2, with the same set of parameters. In the first column (A) ρ = 0.1 and ζ = −1013, in the second column (B) 10 and −1273 and in the third column (C) 100 and −3398, both coefficients in units of Pa. The solid red and dashed blue curves are spaced by a time lapse of 6 min. Periodic boundary conditions are assumed.

4.2 Complex Ginzburg–Landau equation

The above solutions correspond to regions arbitrarily far from the instability threshold and are illustrative of typical solutions that could be identified with experimental observations, where both wave trains and solitary fronts have been reported by Serra-Picamal et al.,¹² Vedula et al.,¹³ and Deforet et al.¹⁵ However, it is obviously unpractical to explore numerically the different qualitative scenarios of nonlinear dynamics in a five-dimensional parameter space. Alternatively, one may exploit the universality of the dynamics of any nonlinear system close to an instability, to develop a reduced description of the nonlinear dynamics of our model near threshold. Such a center-manifold projection does capture the essential nonlinear features of the problem.⁴⁴ The resulting description depends only on symmetries and the nature of the bifurcation, and thus allows classifying the nonlinear dynamics regardless of the physical mechanisms responsible for the instability. This is a powerful approach that can be carried out analytically in the framework of a formal expansion on a small parameter ε defined by the normalized distance to the threshold, and exploits the separation of length and time scales of the spatiotemporal variations of the bifurcating modes with respect to the original scales of the system. In this framework, the so-called amplitude equation appears as a solvability condition at the lowest nontrivial order within that expansion. In our case ε ≡ (ζ − ζ^L_o)/ζ^L_o, which is positive (above threshold) when ζ < ζ^L_o< 0. We define the complex envelope A that describes the modulus and phase modulations of the bifurcating mode in the form

δp(x,t) ∼ A(X,T)e^{i(q_ox−ω_ot)} (29)

where q_o is the critical wave vector and ω_o its corresponding frequency and the slow variables denoted by capital letters are defined as T = εt and X = ε^1/2(x + V_gt), where V_g is the group velocity of the envelope. The modulations of the other physical fields, for instance the cell velocity δv, are proportional to δp, hence they exhibit the same qualitative dynamics. The amplitude of the wave modulation of the base state is δp(x,t) ∼ ε^1/2. Since in our case, the symmetry x → −x is broken and waves travel only in one direction, we will have only the amplitude equation for the right-traveling wave. Accordingly, we can formulate the problem in a reference frame travelling with the group velocity. The normal form is then that of the so-called uniform oscillatory instability (i.e. ω_o ≠ 0, q_o = 0),⁴⁴ even though we have q_o ≠ 0. After proper rescaling, the normal form for the range above the threshold is the so-called complex Ginzburg–Landau equation CGLE,^44,45

∂_TA = A + (1 + ib)∂_X²A − s(1 + ic)|A|²A. (30)

The parameter in front of the cubic nonlinearity is s = ±1. If s = 1, the equation corresponds to the case of the bifurcation being supercritical (continuous). For s = −1 the bifurcation is subcritical (discontinuous) and the equation must be supplemented by a quintic term.

It is worth remarking that, once an additional parameter is eliminated by the condition of being at the instability threshold, the original five-parameter problem gets reduced to a two-parameter problem. Accordingly, since the phase diagram of eqn (30) is well-known,⁴⁶ to classify all possible nonlinear scenarios it suffices to obtain the map of the physical parameters of the original problem into the parameters of the normal form c and b. This mapping always exists but it may be difficult to obtain in practice. We have performed this analysis following standard methods and with the help of the Mathematica software. The extent to which the results obtained through this weakly nonlinear analysis apply to situations more deeply in the nonlinear regime is not known a priori, but it is plausible to expect that the qualitative behavior will be similar relatively far from the theshold.

The crucial piece of information is thus the explicit map that relates c and b to the physical parameters, and the region of parameters for which s = 1. In terms of the conveniently redefined set of dimensionless parameters F ≡ T̄₀L_γ/L_η, G ≡ L_c/L_η and H ≡ ν²L_γ²/L_η², in Appendix A we derive the explicit maps

b = b(ν,F,G,H), (31)

c = c(ν,F,G,H), (32)

s = s(ν,F,G,H). (33)

Note that the dimensionless contractility [small zeta, Greek, macron] does not appear, because it has been eliminated by the additional constraint of being at the instability threshold [small zeta, Greek, macron] = [small zeta, Greek, macron]^L_o(ν,G,H).

From eqn (12) we thus have

We then findwhere f_i(ν,F,G,H) are functions given explicitly in Appendix A. The parameter b has a much simpler expression because it depends solely on the linear part of the dynamics. The genuinely nontrivial part of the dynamics is contained in the parameters c and s. The regions where the bifurcation is subcritical (s = −1) must be analyzed separately and will not be addressed here.

In the particular case of T₀ = 0, the bifurcation is of the stationary periodic type,⁴⁴ with finite q_o and ω_o = 0. Then the instability will lead to the formation of spatial patterns such as in Bois et al.⁴⁷ and the corresponding amplitude equation will be the so-called Real Ginzburg–Landau Equation, with b = c = 0.⁴⁴ In this case, the dynamics of the amplitude equation is variational, that is, a Lyapunov functional [script L] exists such that ∂_TA = −δ[script L]/δA*. Interestingly, this is no longer true if at least one of b and c is nonzero. In this case, the dynamics is said to be persistent and does not relax asymptotically to a specific pattern. This opens the way to different forms of spatio-temporal chaos. The two-dimensional phase diagram of the 1d CGLE has been established in detail,^45,46 and is indeed extremely rich. Different complex dynamical regimes were identified such as the so-called phase turbulence, amplitude turbulence, spatio-temporal intermittency, and bistable chaos, in addition to regimes with more regular behavior. We may refer generically to the above classes of irregular persistent dynamics as weak turbulence. This type of turbulence, arising from secondary instabilities of nonlinear waves, is fundamentally different from the active turbulence of nematic fluids,^48–50 which originates at the spontaneous flow instability of active fluids,⁵¹ and is usually mediated by the dynamics of topological defects. In our case, defects may or may not be present, but the most salient distinction is the crucial role of active traction forces (T₀ ≠ 0).

The boundaries of the different dynamical regimes are usually determined numerically. However, there is an exact boundary that is relevant to our analysis, which locates the so-called Benjamin–Feir (BF) instability. This is a long wavelength instability of the nonlinear travelling-wave solutions of the CGLE of the form , with Ω_Q = c + (b − c)Q². Such waves are unstable when the Newell criterion 1 + bc < 0 is satisfied. Beyond the BF line, one possibility is to have phase turbulence, where the wave phase changes in an irregular manner but preserving the winding number. In that regime, close to the BF line, the phase dynamics can be approximated by the Kuramoto–Sivashinsky equation.^44,46 In other regions, crossing the BF line may lead to amplitude turbulence where the wave amplitude can reach zero values giving rise to non-conservation of the winding number. On the BF-stable side, however, one may also find regions with spatio-temporal intermittency, where patches of regular and chaotic behavior coexist. All these qualitative behaviors are contained in our original model provided that the corresponding values of c and b can be reached by changing the model parameters. In Fig. 4 we show an example of phase turbulence for a situation in the BF-unstable region, obtained from numerical simulations of the original model in the appropriate parameter region.

Fig. 4 Phase turbulence. Numerical simulation of waves in the 1d model given by eqn (21)–(24) for a region of the parameter space beyond the Benjamin–Feir instability line. The total size of the system is 8 × 10⁴, hence the kymograph represents only ∼10% of it. The x-axis is the position in a reference frame moving at the group velocity V_g = 1.41 and the y-axis is time. V_g has been estimated from the average of the phase velocity of polarity peaks. The colour bar labels the modulus of the director field p_x. The parameters are η = 10⁶, ξ = 5, T₀ = 40, L_c = 50, [small gamma, Greek, macron]₁ = 600, ν₁ = 20, ζ = −662.3 and ρ = 10, in units of Pa, μm and s.

To illustrate the complexity of the phase diagram in the four-dimensional space of the model projected in the region close to the instability threshold, we plot it in Fig. 5 for some ranges of the physical parameters. For simplicity we distinguish only three regions, namely the subcritical one, the supercritical BF-stable region and the supercritical BF-unstable one. In the latter, one may find phase turbulence, amplitude turbulence and bistable chaos. In the BF-stable region, traveling-waves are linearly stable to long wavelength modulations, but spatio-temporal intermittency can also be found.

Fig. 5 Exact subcritical–supercritical and Benjamin–Feir instability boundaries in a 4-dimensional plot of the model parameters, at the onset of the primary instability of the problem. In the yellow region the primary instability is subcritical (discontinuous). In the red region the primary instability is supercritical and waves are BF-stable (1 + bc > 0). In the blue region, the primary instability is supercritical and waves are BF-unstable (1 + bc < 0). Each figure represents a cross-section of the 4-dimensional parameter space in the plane given by T₀γ/ξL_η and L_c/L_η, which are varied 4 decades each. The dimensionless parameters γρ/η and ν vary as indicated.

We remark that the explicit knowledge of the exact maps eqn (35)–(37) in a problem with so many parameters and with such a rich variety of complex nonlinear behaviors, is extremely valuable. Indeed, for any set of physical parameters of the model, one can immediately find the expected nonlinear behavior by checking the corresponding point in the known b–c phase diagram, which is universal, and determined once for all.

4.3 Transverse modes and soft-mode turbulence

The previous weakly nonlinear analysis is suitable for the case where q_o ≠ 0, implying that near the threshold, there is a narrow band of nearly marginal modes that excludes q = 0. For transverse modes, however, q_o = 0, and above the threshold the band of unstable modes extends all the way to q = 0, which remains marginal because of rotational invariance. A 1d amplitude equation for the transverse mode along the y axis can also be derived, now for a real amplitude field. We do not address this case here but, as discussed in Cross and Hohenberg,⁴⁴ we remark that this scenario also includes the possibility of phase turbulence, in the form of the Kuramoto–Sivashinsky equation. A similar form of turbulence is reported in Marcq.⁴²

In general, for a 2d case the longitudinal and transverse modes will be coupled at the nonlinear level. A combined weakly nonlinear analysis of this case is beyond the scope of this paper. However, there is a particular case where this coupling may be worked out more easily. This is the coupling of the longitudinal modes to the q = 0 transverse mode. The coupling of a finite q_o mode with a Goldstone mode has been discussed in the literature of liquid crystal electroconvection with homeotropic alignment.⁵² This case was shown to be a remakable case of a direct transition to spatio-temporal chaos at the onset, due to the nonlinear coupling between the Goldstone mode and the bifurcating mode.⁵² The scenario was called soft-mode turbulence, and was demonstrated for a stationary periodic instability (i.e. w_o = 0, q_o ≠ 0), so it would correspond to our case with T₀ = 0, for the nonlinear coupling of the soft transverse mode and the longitudinal q_o mode. We are not aware of any study of soft-mode turbulence for a periodic oscillatory instability, but it is again plausible to expect that the dynamics will be no less chaotic. Consequently, we have every indication that the behavior of the system in sufficiently extended 2d regions will generically contain different forms of weak turbulence, possibly in parameter regions where the 1d modes are not yet turbulent.

5 Application to experiments on epithelial monolayers

We will now discuss the applicability of the present approach to interpret existing data and also to extract quantitative estimates of parameters, mostly from three series of experiments on MDCK cells under different confinement conditions, but all exhibiting some kind of slow oscillatory dynamics, with characteristic periods of several hours and wavelengths in the range of hundreds of micrometers.^12,13,15

5.1 Dynamics of freely spreading epithelia

The present work was directly motivated by the in vitro experiments on freely spreading epithelia described in Serra-Picamal et al.,¹² where apparently elastic ultraslow waves were reported. Those experiments study a wound-healing assay, where a wide planar front of the monolayers spreads at an approximately constant speed. The reported data were averaged over the transversal y-direction, so the possible structure of the fields along the y coordinate, if present, was averaged out. More specifically, the transverse modes with finite q_y and q_x = 0 will be present but averaged out in the data, while the transverse modes with finite q_x and q_y = 0 will generically be present if the instability is sufficiently above the threshold for them to be unstable. At distances of the order of L_c from the leading edge, the systems are manifestly polarized. Note, however, that since the transversal dimension is much larger than L_c, the monolayer could well be polarized in regions further apart from the leading edge, but such polarization would be averaged out in the 1d projection. The waves reported from those experiments can be interpreted as the nonlinear waves reported here, Fig. 6. The origin of the long period of the waves and their elastic-like phase lag were indeed puzzling and several possible explanations have been proposed so far, such as the nonlinear viscoelastic spring model described by Serra-Picamal et al.¹² or the model by Banerjee et al.,²³ where a feedback between local strain, polarization, and contractility is postulated to endow the elastic medium with an effective inertia. Both cases assume that the medium is constitutively elastic, at least partially, and that the origin of the phase-lag between the stress and strain-rate is due to the dominance of the elastic relaxation. By contrast, our result shows that the elastic-like phase lag can be entirely associated with the presence of active traction forces, and thus be observed in a medium with a purely viscous constitutive equation. This remarkable result is consistent with similar considerations discussed in Blanch-Mercader et al.,²⁵ where a purely viscous continuum model for spreading epithelia was shown to explain other apparently elastic behaviors, such as the emergence of an effective elastic modulus.²⁶ In contrast to the above two approaches, in our model the wave frequency turns out to be extrinsic to the material properties, since it is fixed essentially by the parameters characterizing the contact forces with the substrate, friction and traction. This is an interesting feature that could be used to discriminate between different theories, since these parameters can be in principle changed by modifying the properties of the substrate or the molecular complexes that interact with it. Finally, in these experiments the waves seem to propagate in the direction of the polarization and not backwards, consistently with our prediction of polar waves, Fig. 6.

Fig. 6 Time evolution of mechanical quantitates for spreading MDCK monolayers. Each panel is an x − t kymograph of (a) velocity v_x, (b) strain rate ∂_xv_x, (c) traction forces T_x and (d) total stress σ_xx (normalised), with periodic boundary conditions.

5.2 Flow alignment and plithotaxis

In order to fit the data from Serra-Picamal et al.¹² with our model, we can obtain estimates from experimental data of all parameters of our model except for the flow alignment coefficient ν and ρ. The values we obtained are listed in Table 1. The parameter ν which couples the flow and the polarity is difficult to measure in living tissues, and is usually not known. We are only aware of values inferred from data for the epithelium of the Drosophila wing, which are negative with |ν| in the range 3–10.⁵⁴ The use of our model to fit the experimental data gives in our case ν ≈ 10 or larger. Values of |ν| for liquid crystals are only slightly larger than 1.³⁰ The stronger coupling in the case of tissues may effectively entail an active response of the cells to the stress environment. The estimated value of ν follows from eqn (17) which relates the underlying tissue velocity V to the phase velocity v_L of the waves. Although the presence of transverse modes also travelling along x cannot be ruled out in the experiments, the presence of the longitudinal modes is clear because the phase-lag of the stress vs. strain rate, would not be present for purely transverse modes. In any case, for large values of ν the phase velocities predicted for both types of waves are very similar.

Table 1 Order of magnitude of the model parameters for spreading epithelial monolayers. The parameters T₀, L_c, η and ξ can be extracted for instance from Blanch-Mercader et al.²⁵ The coefficients ζ and γ are estimated from Roure et al.,⁸ Weber et al.⁵³ The coefficient ρ is extracted from Lee and Wolgemuth.³⁶ The parameter ν is estimated to make our model consistent with the rest of the parameters and the experimental observations

T 0 [Pa min]	L c [μm]	η [Pa min]	ξ [Pa min μm^−2]
10	10	10^5–10^6	10

---

ζ [Pa]	γ [min]	ρ [Pa]	ν [—]
10^3	10–100	10	10

---

The fact that ν is large has an interesting interpretation. It is known from the hydrodynamics of nematic liquid crystals that this parameter sets the orientation angle θ of the polarity (director) field with respect to a pure shear flow, such that cos 2θ = 1/ν. This simple relation implies that, for |ν| ≫ 1 the polarity of cells orients with an angle of θ ≈ π/4 with the shear. Taking into account that, for a pure viscous shear, the principal stress directions are precisely at this angle, we conclude that cells tend to reorient along the directions of maximal principal stress, along which the shear vanishes. This tendency has been observed in different situations and has been named plithotaxis.⁹ Regardless of whether this response of cells to intercellular stress is an active, regulated process, we find that it can be naturally encoded in the parameter ν. Accordingly, this phenomenon, usually seen as an emergent collective property, can be effectively described as the result of a passive local hydrodynamic coupling between the flow and the polarity.

5.3 Spreading with lateral confinement

In Vedula et al.¹³ epithelial monolayers migrate along adhesive strips with a controllable width. For the most narrow channels, the problem is the closest experimental situation to our specific case of 1d, with purely longitudinal waves. The modes q_y that are present may be limited to a rather large modulus, and therefore the transverse instability may be suppressed (not averaged out as in Serra-Picamal et al.¹²). Similarly, the transverse modes with finite q_x and q_y = 0 will also be suppressed by the boundary conditions on the lateral walls that enforce a fixed orientation of the polarity along them. Consequently, in very narrow channels only longitudinal modes with q_y = 0 are expected to be relevant. For the narrow channels, we associate the contraction–elongation caterpillar-like motion as a signature of our longitudinal waves. In addition, they seem to propagate only in one direction, as predicted by our model. When the channel width is gradually increased, unstable transverse modes are expected to appear and yield the progressively more complex dynamic scenario. The change of behavior for increasing channel width is thus qualitatively and quantitatively consistent with the prediction of our linear analysis. As in all the other cases, the suppressing effect of the treatment with blebbistatin is consistent with the instability mechanism that we propose for the phenomena. Finally, the complexity of the flow patterns observed for wide channels seems to be qualitatively consistent with the scenarios of weak turbulence predicted by our model.

5.4 Oscillations in totally confined monolayers

Finally, our model can be used to revisit the experiments of Deforet et al.¹⁵ where the monolayers are totally confined in circular islands, but nevertheless exhibit oscillatory collective modes. However the fact that the time and length scales are the same as in the other experiments is suggestive that the mechanism behind such a collective mode could be a similar instability adapted to the confined geometry. Whether our model yields oscillatory modes in these situations is an open question that we do not address here. However, we can show that the reported linear dependence of the oscillation period with the tissue radius R is consistent with our linear dispersion relation analysis. In fact, the oscillation frequency is ω_o ∼ qT₀/ξ, then assuming that the range of wave numbers allowed by the geometry is such that q > q_o, increasing the radius sets the most unstable mode available as the minimum q_min ∼ R⁻¹. Consequently, the period will grow linearly with R as reported in Deforet et al.¹⁵

5.5 Collective modes and turbulence in epithelia

In the experiments previously addressed, when the tissue is strongly active, highly disordered flow patterns are observed, often described as noisy.¹⁵ Noisy data of local measurements are often seen as reflecting inherent strong fluctuations of the physical variables. In particular in experiments with very large monolayers, simple visual inspection shows an apparently turbulent behavior.⁵⁵ Whether this apparent chaos is the manifestation of intrinsic noise in the system, or some form of collective modes in a turbulent regime is an interesting open question. In this paper we have seen that secondary instabilities after a Hopf bifurcation do generically lead to different forms of spatio-temporal chaos or weak turbulence, in particular in large systems, so the scenario of chaotic collective modes seems plausible. The distinction between the two possibilities is not only of theoretical interest, but may have some practical relevance. Indeed, if a turbulent state results from an instability, it can be suppressed or triggered at convenience by tuning a single parameter across the appropriate boundary in parameter space. In contrast, if the complex dynamics reflects intrinsic noise, this is virtually impossible to control or regulate.

6 Conclusion

In this paper we present a general framework to account for the mechanics of epithelial monolayers. The model is built on the idea that at sufficiently long length and time scales, a continuum hydrodynamic approach can capture a large variety of mechanical aspects of such monolayers, encoding their complex biological regulation in a set of (possibly time-dependent) physical parameters. Our model includes a polarity field that is not locally aligned with the velocity but coupled to the flow as in nematic hydrodynamics. The contact forces with the substrate contain two contributions, a passive friction force aligned with the velocity, and an active traction force aligned with the polarity. The material also exerts active contractile stresses. A key feature of our model is that the constitutive equation of the medium is taken as that of a viscous (polar) fluid for the slow dynamics. This seems to be at odds with the apparently elastic behavior reported in the literature. In the spirit of Blanch-Mercader et al.,²⁵ where it was first emphasized that such behavior was not inconsistent with a purely viscous constitutive equation provided that the fluid was active, here we show that the presence of waves with an elastic-like phase lag between the stress and strain rate is not necessarily a signature of elasticity, but can also occur in viscous fluids with active traction. Even though the support for the fluid-like scenario based on direct observation of relative cell movements and on arguments based on time scales associated with short-time elasticity is not fully conclusive, our results imply that the assumption of long-time elasticity based on the observation of waves is not justified. In any case, the issue of the passive rheology of epithelial monolayers in different situations still requires further scrutiny.

In addition, our framework provides insights into the physics of collective cell migration. The idea that interacting cells collectively set the local stress environment and the motion of individual cells is incorporated into our physical picture through the nonlocal character of the interaction, which establishes that the flow velocity at a given point is determined by an integration over a region of the size of the friction length. Similarly, the tendency of cells to align with the principal stresses (plithotaxis) appears as a consequence of the large values of the flow alignment coefficient ν, which are obtained here independently to fit observations on the propagating waves. This parameter is crucial to explain the instability leading to longitudinal waves, which can be identified in the experiments by the elastic-like phase lag.

The test of the quantitative predictions of this type of continuum model is not simple because of the difficulty in determining the model parameters and also because these may be changing with time due to the ongoing biological regulation. Nevertheless, we remark that our model is predictive also in qualitative aspects. For instance, we predict that stress waves must be polar, in the sense that they only propagate along the polarization of the medium and not backwards, as opposed to elastic waves. We also predict that the wave frequency, as well as the phase lag are essentially determined by contact forces (friction and traction), which are more amenable to experimental control than material parameters.

In addition, we have pursued our study of waves into the nonlinear regime, and have shown that different forms of weak turbulence are generically present in the nonlinear waves that emerge in our model. We speculate that this chaotic dynamics of the waves may be at the root of the noisy dynamics of tissues. In particular, experiments in very large spreading monolayers exhibit what could be loosely described as turbulence.⁵⁵ Within this picture, the dynamic disorder could well be an expression of highly unstable collective modes, and not a signature of intrinsic noise. This idea could also be tested qualitatively by observing sudden changes between regular and irregular collective behavior through the change of a single parameter, implying the transition from a non-chaotic to a chatic region in parameter space.

The model is not expected to apply to epithelial tissues that are not moving on a substrate.¹⁴ It is unclear to what extent the model can be adapted to situations where there is no global flow, such as in the fully confined experiments of Deforet et al.¹⁵ Recent results under different conditions (Saw et al.³¹) have also reported that the same type of cells can exhibit a nematic (apolar) ordering, based on cell in-plane elongation rather than traction forces, and even extensile stresses rather than contractile. Our approach cannot apply to such a situation, implying that a more complete description of such tissues could possibly require the introduction of two ordering fields with both polar and apolar symmetries, as in other contexts of active matter.³³ However, for the slow spreading dynamics of epithelia on substrates, the physical scenario unveiled here is expected to be generic. Indeed, even though the explicit exact calculation presented here for both linear and nonlinear dynamics refers to a specific choice of the terms included in the model, which has been kept as simple as possible, our central results are robust and do not depend on the model details. The model could be enriched with more parameters, and new physical ingredients, such as effects of cell division, short time elasticity or more complex rheology. However, from the generic nature of the linear and nonlinear analysis discussed here, which relies to a large extent on symmetries, it is expected that two basic results are robust to changes in model details, namely, the mechanism that controls the traction-driven instability of an active viscous polar fluid, leading to polar nonlinear travelling waves, and the nonlinear instabilities that lead generically to weak turbulence scenarios. The need of at least two sources of activity, and a coupling between polarity and flow seem also well established. We expect that further experimental investigations will eventually test the ideas developed here and clarify the appropriate mechanical framework for a continuum description for collective cell migration in tissues.

Conflicts of interest

There are no conflicts to declare.

Appendix A

Here we provide details of the weakly nonlinear analysis leading to the CGLE in our physical model for the case of 1d longitudinal modes. We use the corresponding physical model given by eqn (21)–(24) and assume a system of units such that γ = η = ξ = 1.

Close to the threshold and expanding the linear growth rate around its maximum at q_o, we get

withThe weakly nonlinear analysis is a formal expansion on the small parameter ε ≡ (ζ − ζ^L_o)/ζ^L_o that measures the normalized distance to the threshold. We will refer to ε > 0 as the system being (slightly) above the threshold. Then, a narrow band of modes with size are unstable but nearly marginal, since for them Re[ω] ∼ ε. In contrast, modes with |q − q_o| ∼ ε⁰ will relax much faster, with Re[ω] ∼ ε⁰. Accordingly, long wavelength spatial modulations of order ε^−1/2 with slow relaxation times of order ε⁻¹ are expected to dominate the dynamics, and slave the rest of the (fast) modes. This separation of time and length scales is at the root of the universal description in terms of an amplitude equation for the envelope of the bifurcating mode.

Consequently, in general the perturbations about the reference ordered state can be expressed as a superposition of plane waves with wavenumber multiples of q_o and phase velocity v_o = −Im[ω(q_o)]/q_o plus an envelope wave with slow spatio-temporal dynamics. In our particular case the most general solution reads

where the envelope waves of the corresponding physical fields p_nm, v_nm and σ_nm are functions of the long spatial variable X ≡ ε^1/2(x + V_gt) and the slow temporal variable T ≡ εt. V_g is the travelling speed of the wave envelope (the group velocity) and in general is a power series in ε^1/2, whose coefficients are treated as unknowns.

The physical solution valid near the vicinity of the transition (41–44) is inserted into the PDE's (21–24) and the different terms are sorted in powers of ε^1/2. Note that the zeroth order leads to the ordered uniform solution. The first order turns into an undetermined set of linear equations for σ^m₁, v^m₁ and p^m₁. As an example these coefficients solved as a function of amplitude p¹₁ read

Consequently, the different mechanical fields will exhibit the same qualitative dynamics. Similarly the second order can also be arranged into a linear set of equations for the coefficients σ_2m, v_2m, p_2m and V₁. Through a solvability condition these second order coefficients are connected to p_1m. These conditions are often used to construct self-consistent solutions through perturbative analysis. Lastly, the third order solvability condition is analogous to the complex Ginzburg–Landau equation. Thus after rescaling the variables conveniently, it can be expressed in the form

∂_Tp₁₁ = p₁₁ − s(1 + ic)|p₁₁|²p₁₁ + (1 + ib)∂_X²p₁₁, (46)

s, c and b being parameters that in general depend on the details of the mechanical properties of the system. The coefficient s is either ±1 and it controls whether the transition is continuous or discontinuous. With respect to eqn (30), we have replaced the variable p₁₁ by A.

The longitudinal mechanical transition in our system is oscillatory and periodic, meaning that the critical wavenumber q_o is finite and also the travelling speed of the perturbations v_o. As a result, the coefficients of our amplitude eqn (46) are complex. Only in the particular case of null active traction forces (i.e. T₀ = 0) do both coefficients b and c vanish, reducing eqn (46) to the real Ginzburg–Landau equation for s = 1. This equation is purely relaxational, that is, there exists a Lyapunov function that is maximised over time. Except for some particular cases no such Lyapunov functional can be constructed for the CGLE, giving rise to a richer phenomenology of dynamical states: from travelling coherent states to different forms of spatiotemporal chaotic states.

For completeness, we present the analytical form of the coefficients s, c and b as

where the function Sign returns the sign of its argument. The functions f₁, f₂, f₃ and f₄ readand depends only on the dimensionless physical parameters F = T₀γ/ξL_η, H = ργν²/η and G = L_c/L_η.

Acknowledgements

We thank R. Alert, F. Graner, V. Hakim, J.-F. Joanny, P. Marcq, J. Prost, X. Trepat, R. Vincent and S. Yabunaka for useful discussion. We acknowledge financial support from MINECO under projects FIS2013-41144-P and FIS2016-78507-C2-2-P from Generalitat de Catalunya under project 2014-SGR-878. CBM also acknowledges a FPU fellowship from the Spanish Government.

References

1. T. Lecuit and L. Le Goff, Nature, 2007, 450, 189–192.
2. D. M. Bryant and K. E. Mostov, Nat. Rev. Mol. Cell Biol., 2008, 9, 887–901.
3. P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol., 2009, 10, 445–457.
4. K. J. Sonnemann and W. M. Bement, Annu. Rev. Cell Dev. Biol., 2011, 27, 237–263.
5. A. Brugues, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Munoz, G. W. Brodland, B. Ladoux and X. Trepat, Nat. Phys., 2014, 10, 683–690.
6. O. Cochet-Escartin, J. Ranft, P. Silberzan and P. Marcq, Biophys. J., 2014, 106, 65–73.
7. D. E. Leckband, Q. le Duc, N. Wang and J. de Rooij, Curr. Opin. Cell Biol., 2011, 23, 523–530.
8. O. du Roure, A. Saez, A. Buguin, R. H. Austin, P. Chavrier, P. Siberzan and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 2390–2395.
9. X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, D. A. Weitz, J. P. Butler and J. J. Fredberg, Nat. Phys., 2009, 5, 426–430.
10. G. F. Weber, M. A. Bjerke and D. W. DeSimone, J. Cell Sci., 2011, 124, 1183–1193.
11. S. Tlili, C. Gay, F. Graner, P. Marcq, F. Molino and P. Saramito, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 1–31.
12. X. Serra-Picamal, V. Conte, R. Vincent, E. Anon, D. T. Tambe, E. Bazellieres, J. P. Butler, J. J. Fredberg and X. Trepat, Nat. Phys., 2012, 8, 628–634.
13. S. R. K. Vedula, M. C. Leong, T. L. Lai, P. Hersen, A. J. Kabla, C. T. Lim and B. Ladoux, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 12974–12979.
14. A. R. Harris, L. Peter, J. Bellis, B. Baum, A. J. Kabla and G. T. Charras, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16449–16454.
15. M. Deforet, V. Hakim, H. Yevick, G. Duclos and P. Silberzan, Nat. Commun., 2014, 5, 3747.
16. J. Notbohm, S. Banerjee, K. J. Utuje, B. Gweon, H. Jang, Y. Park, J. Shin, J. P. Butler, J. J. Fredberg and M. C. Marchetti, Biophys. J., 2016, 110, 2729–2738.
17. J. C. Arciero, Q. Mi, M. F. Branca, D. J. Hackam and D. Swigon, Biophys. J., 2011, 100, 535–543.
18. P. Lee and C. Wolgemuth, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061920.
19. A.-K. Marel, M. Zorn, C. Klingner, R. Wedlich-Söldner, E. Frey and J. O. Rädler, Biophys. J., 2014, 107, 1054–1064.
20. P. Recho, J. Ranft and P. Marcq, Soft Matter, 2016, 12, 2381–2391.
21. S. Banerjee and M. C. Marchetti, Phys. Rev. Lett., 2012, 109, 108101.
22. M. H. Kopf and L. M. Pismen, Soft Matter, 2013, 9, 3727–3734.
23. S. Banerjee, K. J. C. Utuje and M. C. Marchetti, Phys. Rev. Lett., 2015, 114, 228101.
24. S. Tlili, E. Gauquelin, B. Li, O. Cardoso, B. Ladoux, H. Delanoë-Ayari and F. Graner, 2016, arXiv:1610.05420, arXiv preprint.
25. C. Blanch-Mercader, R. Vincent, E. Bazellières, X. Serra-Picamal, X. Trepat and J. Casademunt, Soft Matter, 2017, 13, 1235–1243.
26. R. Vincent, E. Bazellières, C. Pérez-González, M. Uroz, X. Serra-Picamal and X. Trepat, Phys. Rev. Lett., 2015, 115, 248103.
27. K. Kruse, J. F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2005, 16, 5–16.
28. F. Jülicher, K. Kruse, J. Prost and J.-F. Joanny, Phys. Rep., 2007, 449, 3–28.
29. R. Farooqui and G. Fenteany, J. Cell Sci., 2005, 118, 51–63.
30. P. G. De Gennes and J. Prost, The physics of liquid crystals, International series of monographs on physics, Oxford University Press, London, 2nd edn, 1997; E. Olbrich, O. Marinov and D. Davidov, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 2713, 48.
31. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
32. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327–331.
33. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
34. M. Lambert, O. Thoumine, J. Brevier, D. Choquet, D. Riveline and R.-M. Mége, Exp. Cell Res., 2007, 313, 4025–4040.
35. J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. Jülicher, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 20863–20868.
36. P. Lee and C. W. Wolgemuth, PLoS Comput. Biol., 2011, 7, e1002007.
37. I. S. Aranson, A. Sokolov, J. O. Kessler and R. E. Goldstein, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 040901.
38. F. Jülicher and J. Prost, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 29, 27–36.
39. D. Oriola, R. Alert and J. Casademunt, Phys. Rev. Lett., 2017, 118, 088002.
40. L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101.
41. S. Sankararaman and S. Ramaswamy, Phys. Rev. Lett., 2009, 102, 118107.
42. P. Marcq, Eur. Phys. J. E: Soft Matter Biol. Phys., 2014, 37, 29.
43. U. S. Schwarz and S. A. Safran, Rev. Mod. Phys., 2013, 85, 1327–1381.
44. M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys., 1993, 65, 851.
45. I. S. Aranson and L. Kramer, Rev. Mod. Phys., 2002, 74, 99.
46. H. Chate, Nonlinearity, 1994, 7, 185.
47. J. S. Bois, F. Jülicher and S. W. Grill, Phys. Rev. Lett., 2011, 106, 028103.
48. S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101.
49. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
50. R. Ramaswamy and F. Jülicher, Sci. Rep., 2016, 6, 20838.
51. R. Voituriez, J. F. Joanny and J. Prost, EPL, 2005, 70, 404.
52. A. Rossberg, A. Hertrich, L. Kramer and W. Pesch, Phys. Rev. Lett., 1996, 76, 4729.
53. G. F. Weber, M. A. Bjerke and D. W. DeSimone, Dev. Cell, 2012, 22, 104–115.
54. B. Aigouy, R. Farhadifar, D. B. Staple, A. Sagner, J.-C. Röper, F. Jülicher and S. Eaton, Cell, 2010, 142, 773–786.
55. R. Vincent and X. Trepat, private communication.

---

Footnote

† Note that the linear equation for the fluctuations of the strain field in the model of Notbohm et al.¹⁶ and Banerjee et al.²³ is of second order in time and includes viscoelasticity and dispersivity, but is still local in space.

---