Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

C.
Blanch-Mercader
*^{ab} and
J.
Casademunt
^{ac}
^{a}Departament de Física de la Matèria Condensada, Universitat de Barcelona, Barcelona 08028, Spain. E-mail: carles.blanch-mercader@curie.fr
^{b}Laboratoire Physico Chimie Curie, Institut Curie, PSL Research University, CNRS, 26 rue d' Ulm, 75005 Paris, France
^{c}Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain

Received
6th June 2017
, Accepted 7th August 2017

First published on 7th August 2017

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.

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^{14} on time scales for which a freely spreading monolayer on a substrate seems to be flowing like a viscous fluid.^{25} 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.,^{14} 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.^{26} 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.^{26} 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.^{25}

In this context, our study is directly motivated by the experiments of Serra-Picamal et al.^{12} 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.^{12} and Banerjee et al.^{23}) 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.,^{12} Vedula et al.,^{13} and Deforet et al.^{15} 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.^{9}

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.

(1) |

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,^{11} 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,^{34} 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.^{35} However, in the experiments of Serra-Picamal et al.^{12} 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.^{25} 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_{0}p_{α}, | (2) |

(3) |

(4) |

h_{α} = ρ(1 − p^{2})p_{α} + K∇^{2}p_{α}. | (5) |

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.^{30}

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, γ_{1} 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.^{37} The form is limited by symmetry considerations^{38} and can be derived from microscopic models with linker kinetics.^{39} 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.^{43} The active traction must be related to the consumption of ATP, as the active contractility term, but the two parameters T_{0} and ζ originate at different mechanisms and can be considered as independent parameters. T_{0} > 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.^{25} 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

(6) |

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 γ ≡ γ_{1}/ρ). We use the nematic stress scale ρ to define a dimensionless contractility = ζ/ρ and a dimensionless traction _{0} = T_{0}L_{γ}/ρ. In addition to ν, and _{0}, 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.

(7) |

λ = λ_{p} + Re[λ_{ζ}] + iRe[λ_{T0}] | (8) |

(9) |

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^{2}. These read

(10) |

(11) |

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 || > |^{L}_{o}| where the threshold value of contractility for the onset of the longitudinal instability is given by

(12) |

(13) |

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

(14) |

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} ≠ 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.

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

(15) |

(16) |

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_{0}/ξ, we find that

(17) |

(18) |

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} = 0, which reads, with q ≡ q_{y},

(19) |

(20) |

In this case the model equations reduce to

∂_{x}σ_{xx} = ξv_{x} − T_{0}p_{x}, | (21) |

(22) |

(23) |

h_{x} = ρ(1 − p_{x}^{2})p_{x} + K∂_{x}^{2}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,

(25) |

(26) |

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.^{12} and Notbohm et al.^{16} 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} ≠ 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

(27) |

In the particular case of 2ν(1 + q_{c}^{2}/2) ≪ || and ||q_{γ} ≪ _{0}, eqn (27) takes the simple form

(28) |

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_{0}/ξ 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^{6}, ξ = 100, T_{0} = 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. |

δp(x,t) ∼ A(X,T)e^{i(qox−ωot)} | (29) |

∂_{T}A = A + (1 + ib)∂_{X}^{2}A − s(1 + ic)|A|^{2}A. | (30) |

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,^{46} 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 ≡ _{0}L_{γ}/L_{η}, G ≡ L_{c}/L_{η} and H ≡ ν^{2}L_{γ}^{2}/L_{η}^{2}, 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) |

From eqn (12) we thus have

(34) |

(35) |

(36) |

(37) |

In the particular case of T_{0} = 0, the bifurcation is of the stationary periodic type,^{44} with finite q_{o} and ω_{o} = 0. Then the instability will lead to the formation of spatial patterns such as in Bois et al.^{47} and the corresponding amplitude equation will be the so-called Real Ginzburg–Landau Equation, with b = c = 0.^{44} In this case, the dynamics of the amplitude equation is variational, that is, a Lyapunov functional exists such that ∂_{T}A = −δ/δ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,^{51} 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} ≠ 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^{2}. 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^{4}, 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^{6}, ξ = 5, T_{0} = 40, L_{c} = 50, _{1} = 600, ν_{1} = 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.

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.

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.^{52} 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.^{52} 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} = 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.

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 cos2θ = 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.^{9} 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.

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.^{55} 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.^{14} 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.^{15} Recent results under different conditions (Saw et al.^{31}) 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.^{33} 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.

Close to the threshold and expanding the linear growth rate around its maximum at q_{o}, we get

(38) |

(39) |

(40) |

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

(41) |

(42) |

(43) |

(44) |

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}_{1}, v^{m}_{1} and p^{m}_{1}. As an example these coefficients solved as a function of amplitude p^{1}_{1} read

(45) |

∂_{T}p_{11} = p_{11} − s(1 + ic)|p_{11}|^{2}p_{11} + (1 + ib)∂_{X}^{2}p_{11}, | (46) |

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} = 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

(47) |

(48) |

(49) |

(50) |

(51) |

(52) |

(53) |

- T. Lecuit and L. Le Goff, Nature, 2007, 450, 189–192 CrossRef CAS PubMed.
- D. M. Bryant and K. E. Mostov, Nat. Rev. Mol. Cell Biol., 2008, 9, 887–901 CrossRef CAS PubMed.
- P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol., 2009, 10, 445–457 CrossRef CAS PubMed.
- K. J. Sonnemann and W. M. Bement, Annu. Rev. Cell Dev. Biol., 2011, 27, 237–263 CrossRef CAS PubMed.
- 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 CrossRef CAS PubMed.
- O. Cochet-Escartin, J. Ranft, P. Silberzan and P. Marcq, Biophys. J., 2014, 106, 65–73 CrossRef CAS PubMed.
- D. E. Leckband, Q. le Duc, N. Wang and J. de Rooij, Curr. Opin. Cell Biol., 2011, 23, 523–530 CrossRef CAS PubMed.
- 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 CrossRef CAS PubMed.
- 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 CrossRef CAS.
- G. F. Weber, M. A. Bjerke and D. W. DeSimone, J. Cell Sci., 2011, 124, 1183–1193 CrossRef CAS PubMed.
- 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 CrossRef PubMed.
- 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 CrossRef CAS.
- 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 CrossRef CAS PubMed.
- 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 CrossRef CAS PubMed.
- M. Deforet, V. Hakim, H. Yevick, G. Duclos and P. Silberzan, Nat. Commun., 2014, 5, 3747 CAS.
- 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 CrossRef CAS PubMed.
- J. C. Arciero, Q. Mi, M. F. Branca, D. J. Hackam and D. Swigon, Biophys. J., 2011, 100, 535–543 CrossRef CAS PubMed.
- P. Lee and C. Wolgemuth, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061920 CrossRef PubMed.
- A.-K. Marel, M. Zorn, C. Klingner, R. Wedlich-Söldner, E. Frey and J. O. Rädler, Biophys. J., 2014, 107, 1054–1064 CrossRef CAS PubMed.
- P. Recho, J. Ranft and P. Marcq, Soft Matter, 2016, 12, 2381–2391 RSC.
- S. Banerjee and M. C. Marchetti, Phys. Rev. Lett., 2012, 109, 108101 CrossRef PubMed.
- M. H. Kopf and L. M. Pismen, Soft Matter, 2013, 9, 3727–3734 RSC.
- S. Banerjee, K. J. C. Utuje and M. C. Marchetti, Phys. Rev. Lett., 2015, 114, 228101 CrossRef PubMed.
- S. Tlili, E. Gauquelin, B. Li, O. Cardoso, B. Ladoux, H. Delanoë-Ayari and F. Graner, 2016, arXiv:1610.05420, arXiv preprint.
- C. Blanch-Mercader, R. Vincent, E. Bazellières, X. Serra-Picamal, X. Trepat and J. Casademunt, Soft Matter, 2017, 13, 1235–1243 RSC.
- R. Vincent, E. Bazellières, C. Pérez-González, M. Uroz, X. Serra-Picamal and X. Trepat, Phys. Rev. Lett., 2015, 115, 248103 CrossRef PubMed.
- 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 CrossRef CAS PubMed.
- F. Jülicher, K. Kruse, J. Prost and J.-F. Joanny, Phys. Rep., 2007, 449, 3–28 CrossRef.
- R. Farooqui and G. Fenteany, J. Cell Sci., 2005, 118, 51–63 CrossRef CAS PubMed.
- 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 Search PubMed; E. Olbrich, O. Marinov and D. Davidov, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 2713, 48 Search PubMed.
- 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 CrossRef CAS PubMed.
- K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327–331 CrossRef CAS PubMed.
- 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 CrossRef CAS.
- M. Lambert, O. Thoumine, J. Brevier, D. Choquet, D. Riveline and R.-M. Mége, Exp. Cell Res., 2007, 313, 4025–4040 CrossRef CAS PubMed.
- 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 CrossRef CAS PubMed.
- P. Lee and C. W. Wolgemuth, PLoS Comput. Biol., 2011, 7, e1002007 CAS.
- I. S. Aranson, A. Sokolov, J. O. Kessler and R. E. Goldstein, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 040901 CrossRef PubMed.
- F. Jülicher and J. Prost, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 29, 27–36 CrossRef PubMed.
- D. Oriola, R. Alert and J. Casademunt, Phys. Rev. Lett., 2017, 118, 088002 CrossRef PubMed.
- L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101 CrossRef PubMed.
- S. Sankararaman and S. Ramaswamy, Phys. Rev. Lett., 2009, 102, 118107 CrossRef PubMed.
- P. Marcq, Eur. Phys. J. E: Soft Matter Biol. Phys., 2014, 37, 29 CrossRef PubMed.
- U. S. Schwarz and S. A. Safran, Rev. Mod. Phys., 2013, 85, 1327–1381 CrossRef CAS.
- M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys., 1993, 65, 851 CrossRef CAS.
- I. S. Aranson and L. Kramer, Rev. Mod. Phys., 2002, 74, 99 CrossRef.
- H. Chate, Nonlinearity, 1994, 7, 185 CrossRef.
- J. S. Bois, F. Jülicher and S. W. Grill, Phys. Rev. Lett., 2011, 106, 028103 CrossRef PubMed.
- S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101 CrossRef PubMed.
- L. Giomi, Phys. Rev. X, 2015, 5, 031003 Search PubMed.
- R. Ramaswamy and F. Jülicher, Sci. Rep., 2016, 6, 20838 CrossRef CAS PubMed.
- R. Voituriez, J. F. Joanny and J. Prost, EPL, 2005, 70, 404 CrossRef CAS.
- A. Rossberg, A. Hertrich, L. Kramer and W. Pesch, Phys. Rev. Lett., 1996, 76, 4729 CrossRef CAS PubMed.
- G. F. Weber, M. A. Bjerke and D. W. DeSimone, Dev. Cell, 2012, 22, 104–115 CrossRef CAS PubMed.
- B. Aigouy, R. Farhadifar, D. B. Staple, A. Sagner, J.-C. Röper, F. Jülicher and S. Eaton, Cell, 2010, 142, 773–786 CrossRef CAS PubMed.
- 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.^{16} and Banerjee et al.^{23} is of second order in time and includes viscoelasticity and dispersivity, but is still local in space. |

This journal is © The Royal Society of Chemistry 2017 |