# Onset of nonlinearity in a stochastic model for auto-chemotactic advancing epithelia

Martine Ben Amar, Carlo Bianca

<figure><img src="image_1.png" /></figure>

## To cite this version:

Martine Ben Amar, Carlo Bianca. Onset of nonlinearity in a stochastic model for auto-chemotactic advancing epithelia. Scientific Reports, 2016, 6, pp.33849. 10.1038/srep33849. hal-01382784

HAL Id: hal-01382784

https://hal.sorbonne-universite.fr/hal-01382784

Submitted on 17 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Onset of nonlinearity in a stochastic model for auto-chemotactic advancing epithelia

Martine Ben Amar<sup>1,2,*</sup> & Carlo Bianca<sup>1,2,*</sup>

Received: 12 July 2016

Accepted: 01 September 2016

Published: 27 September 2016

We investigate the role of auto-chemotaxis in the growth and motility of an epithelium advancing on a solid substrate. In this process, cells create their own chemoattractant allowing communications among neighbours, thus leading to a signaling pathway. As known, chemotaxis provokes the onset of cellular density gradients and spatial inhomogeneities mostly at the front, a phenomenon able to predict some features revealed in *in vitro* experiments. A continuous model is proposed where the coupling between the cellular proliferation, the friction on the substrate and chemotaxis is investigated. According to our results, the friction and proliferation stabilize the front whereas auto-chemotaxis is a factor of destabilization. This antagonist role induces a fingering pattern with a selected wavenumber *k*₀. However, in the planar front case, the translational invariance of the experimental set-up gives also a mode at *k* = 0 and the coupling between these two modes in the nonlinear regime is responsible for the onset of a Hopf-bifurcation. The time-dependent oscillations of patterns observed experimentally can be predicted simply in this continuous non-linear approach. Finally the effects of noise are also investigated below the instability threshold.

Recently moving epithelium on a substrate has motivated many experimental<sup>1,2</sup> and numerical works<sup>2,3</sup> with the main aim to understand the collective migration of cell sheets occurring during development, inflammation, wound-healing<sup>4-7</sup> or tumour metastasis<sup>8,9</sup>. Consideration of the movement of several thousand cells in a tissue of constant thickness can explain the mechanism of cell communication and migration. Intercellular interactions and substrate adhesion<sup>10</sup>, i.e mechanical forces, together with sensitivity to biochemical signals<sup>11</sup> are responsible for these phenomena. One may argue that the importance of these morphogenetic signals depends on the experimental context, some cells being more chemotactic than others. However, in physiological conditions, chemotaxis cannot be ignored, e.g. as in the wound-healing process: Just after injury<sup>4</sup>, a battery of morphogens originating from the immune system organizes the migration of several connected cell layers of the basal epithelium. In particular, it has been recognized that cells can produce their own chemoattractant, as melanocytic cells. For these cancerous cells, it is a way to escape the tumour and undergo metastasis<sup>9</sup>. As shown in Fig. 1, in the circular geometry<sup>12</sup>, *in vitro* advancing epithelium of melanocytic cells exhibits simultaneously a noisy front and a diffuse interface, which strongly suggests the local production of chemo-attractants.

This paper is concerned with the so-called auto-chemotaxis process where cells themselves produce their own chemo-attractants. Auto-chemotaxis can be considered as a mechanism for neighbour cell communication, which by transforming biochemical signaling into mechanical stresses induce collective migration<sup>13</sup>. This process has also been investigated in the dynamical aggregation of dense micro-swimmer colonies<sup>14</sup>. It is worth stressing that, even in extremely well controlled experiments, it appears difficult to distinguish between chemotaxis and mechanotaxis<sup>10,13,15</sup> (including all forces and stresses that a cell can feel in a moving epithelium). In recent theoretical approaches for *in vitro* advancing epithelia, chemotaxis has not been considered and discrete cellular models with short-range interactions, treated with simulations eventually including noise, are preferred<sup>2,16</sup>. However, here the choice is made of a continuous model to explain experimental features, such as the origin of waves occurring during tissue expansion<sup>16</sup>.

<footnote n="1"><sup>1</sup>Laboratoire de Physique Statistique, Ecole Normale Supérieure, PSL Research University; Sorbonne Universités UPMC Univ Paris 06; CNRS; , 24 rue Lhomond, 75005 Paris, France. <sup>2</sup>Institut Universitaire de Cancérologie, Faculté de médecine, Université Pierre et Marie Curie-Paris 6, 91 Bd de l'Hôpital, 75013 Paris, France. <sup>*</sup>These authors contributed equally to this work. Correspondence and requests for materials should be addressed to M.B.A. (email: benamar@lps.ens.fr)</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Human melanocytic cells (primary tumour) were plated on 33 cm² Petri dish with at the center a disk (32 mm²) at confluence. The day after, the disk was removed and cells were allowed to keep migrating and dividing, before they were fixed with 3.7% paraformaldeide and stained with hematoxillin/eosin solution, at different times. On top, typical experiment, at two different times, on bottom, details from the same images with X5 magnification, compared to the top (Leica MZFLIII mounted with a camera Leica DFC320). Observe the diffuse front and an increase of noise after 4 days. From publication¹² with permission.</figcaption></figure>

Pattern formation has a rather long history in physics, in diffusion of growing crystals¹⁷,¹⁸, viscous fingering in fluids¹⁹,²⁰, electrostatics²¹,²² and soft matter²³. Theoretical methods have been developed to understand the front dynamics originating from the interplay between interfacial and bulk processes. If continuous models have been proposed in the past for advancing epithelia³,²⁴,²⁵, simulations are preferred in most cases, although nonlinear analysis of pattern formation can help in the understanding of the collective motion of cells when the experimental size is much larger than the cell size. This paper aims at explaining the pattern formation observed in recent experiments¹²,¹⁵,²⁶, by means of the definition of a continuous model of auto-chemotactic migration, including also cell mechanics and cell proliferation. We show that the nonlinearity coupled to the inherent stochasticity of any biochemical signaling induces pattern formation below the instability threshold and we evaluate the velocity structure function which is a measurable parameter.

The paper is organized as follows: after a presentation of the system of study, a continuous theoretical model is presented and an advancing steady front is theoretically predicted. A stability study restricts the range of physical parameters allowing periodic patterns that are determined at the linear approximation. Then the stochasticity and the nonlinearities are introduced to predict the wavy dynamics observed in the quasi-infinite geometry.

**The model**

**Objectives and Hypotheses.** A theoretical model is proposed to investigate the dynamics of a cellular sheet, controlled by auto-chemotactic gradients and cellular proliferation. A front separates the experimental cells into two domains Ωₑ (the epithelium) and Ω<sub>w</sub> (the water and nutrients), as shown in Fig. 2. Physical quantities characterizing the front, such as velocity, stress repartition, stochastic and spatio-temporal patterns, will be predicted with a minimal set of parameters according to the following strategy: steady leading edge dynamics, stochastic and wavy edge, spatio-temporal front, (see the three schemes of Fig. 2). Among the three scenarios, from left to right, there is an increase of the chemotactic gradients which controls the evolution: first an advancing linear epithelium edge with constant velocity, second a wavy edge, also advancing with constant velocity, then above the stability limit, a spatio-temporal pattern.

The biological description is based on simple assumptions:

* The cellular proliferation occurs as soon as the level of tension inside the cells is higher than the homeostatic one, according to a mechano-sensitive growth process⁷,¹²,²⁷,

* On the substrate the motility law is dominated by friction,

* The cells communicate via chemical signals giving rise to auto-chemotaxis.

Auto-chemotaxis and proliferation are two factors of migration that occur simultaneously. In order to limit the number of independent parameters in our analysis, we assume the same diffusion coefficient Dₑ for the morphogens inside the cellular domain Ωₑ and in the water bath Ω<sub>w</sub>. In Ωₑ, the morphogens diffuse but are also captured
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The three steps of leading edge evolution as the auto-chemotactif gradients increase. In grey, the cell sheet ($\Omega_e$), in colours, the bath of water and nutrients ($\Omega_w$). From left to right, a steady front with constant velocity, a wavy front at threshold of stability, then the spatio-dynamical pattern.</figcaption></figure>

by the cells and disappear with a characteristic uptake $\tau_e$. The cells of the interface produce more morphogens than inside the epithelium where these chemicals are absorbed, so only a small number of rows are responsible for the morphogen production. Notice however that cells away from the front also produce morphogens but in less quantity that the ones at the frontier. In addition, it is known that epithelia have a constant thickness, except at the boundary, where a thinning can be observed as in Fig. 1 and reported$^{12,28}$. Based on these biological informations and hypotheses, the main equations follow.

**The main equations.** The morphogen concentration satisfies the diffusion-consumption equation in the tissue and the diffusion equation in water:

$$ \left\{ \begin{array}{l} \partial_t c_e = D_e \Delta c_e - c_e / \tau_e \quad \text{if } \mathbf{r} \in \Omega_e \\ \partial_t c_w = D_e \Delta c_w \quad \text{if } \mathbf{r} \in \Omega_w \end{array} \right. \qquad (1) $$

Let us define dimensionless physical quantities as capital letters and choose $\tau_e$ as time unit, $\sqrt{D_e \tau_e}$ as length unit and as concentration unit, the averaged value at the front, $c_i = (c_e + c_w)/2$. Then, the equation (1) are transformed into:

$$ \partial_t C_e = \Delta C_e - C_e \quad \text{and} \quad \partial_t C_w = \Delta C_w \qquad (2) $$

At the front of normal $\vec{N}$, oriented towards $\Omega_w$, the jump conditions for both the concentration field and its normal gradient read:

$$ C_e - C_w = 1 \quad \text{and} \quad \vec{N} \cdot (\vec{\nabla} C_e - \vec{\nabla} C_w) = N_0(\vec{N} \cdot \vec{V}_i) \qquad (3) $$

The first relation of equation (3) gives the normalized jump in the morphogen concentration while $N_0$ indicates the rate of morphogen production at the border, responsible for the flux discontinuity. These relations indicate a discontinuity due to the localization of the morphogen sources, attached to the front and being transported with it. These relationships are similar to the boundary conditions of the so-called one sided-model of solidification for binary mixtures$^{17}$ and they have been mathematically established$^{29}$ and derived in chemotaxis$^{30}$. In both domains, far from the interface, $C_e = C_w = 0$. When the friction on the substrate$^{31}$ dominates the sheet motility, the cell velocity follows a Darcy's law: $\vec{v}_e - \vec{v}_s = -K_e \nabla \rho$ where $K_e$ is the mobility coefficient. $\vec{v}_s$ is usually neglected since most of the experimental set-ups are at rest in the laboratory frame. However, this formulation makes more explicit the Galilean invariance of the Darcy's law and underlines the fact that friction results from velocity differences$^{32-34}$. In the laboratory frame, we transform Darcy's law into $\vec{V} = -\vec{\nabla} P$, once the unit of pressure is chosen as $P_0 = D_e/K_e$. This law can eventually be modified into a power law model used in ref. 30 where the mobility coefficient $K_e$ is itself a function of the pressure gradient. The cells proliferate with a rate $K_p$ as soon as the pressure is below the homeostatic pressure $P_h^{12,27}$, (chosen as our pressure of reference) which leads to the following mass-flux equation:

$$ \partial_t \rho + \vec{\nabla}(\rho(\vec{V} - \Lambda \vec{\nabla} C_e) - D \nabla \rho) = - \alpha^2 \rho P \qquad (4) $$

where the chemotactic migration constant is introduced: $\Lambda = \lambda_0 c_i / D_e$, suppress the cell diffusion coefficient: $D$ and the coefficient of proliferation: $\alpha^2 = K_p D_e / K_e$. Equation (4) superposes the convective flux to the diffusive and chemotactic fluxes$^{35}$. For an epithelium of constant thickness (not true at the interface), $\rho$ is a constant and the model assumes a priori cohesiveness and excluded volume interactions. Contrary to discrete simulations, these two constraints are automatically verified. Moreover, far from the front, we get for the pressure the simplified equation: $\Delta(P + \Lambda C_e) - \alpha^2 P = 0$. However chemotaxis does not allow a constant thickness at the front$^{30}$ and depletion of density is experimentally observed$^{12,28}$ (see also Fig. 1). This density depletion is necessary to ensure the mechanical equilibrium of the interface as it has been demonstrated$^{30}$ in details. Indeed, an abrupt front under chemotaxis does not allow to maintain simultaneously two constraints: on the one hand, the equality of cellular
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left panel (A): Density profile of the morphogens as a function of the distance from the front position $y=0$, in the experimental cell. Velocities are increasing from 3 to 7 while setting $N_0 = 1$. Right panel (B): The front velocity as a function of combined chemotactic parameters, when varying the interface pressure and setting $\alpha = 1$ (see equation (8)).</figcaption></figure>

and interface normal velocities, and on the other hand, the cancellation of chemical flux. The cellular density variation and diffusion, inside the boundary layer, give the necessary degrees of freedom allowing this subtle equilibrium by imposing: $-\Lambda\rho\vec{N} \cdot \vec{\nabla} C_e - D\vec{N} \cdot \vec{\nabla} \rho = 0$. If the size of the transition zone remains small, the mechanical equilibrium of the front fixes the boundary condition for the pressure:

$$ \vec{N} \cdot \vec{V} = -\vec{N} \cdot \vec{\nabla} P = \vec{N} \cdot \vec{V}_i \quad \text{and} \quad P = P_i + \sigma\kappa \qquad (5) $$

where $\vec{V}_i$ is the interface velocity. The interface pressure $P_i$ is the pressure in the water bath, relative to the homeostatic pressure, the curvature $\kappa$ is chosen positive for convex interface and $\sigma$ represents the dimensionless capillary parameter $\sigma = S_T K_e / (D_e^{3/2} \tau_e^{1/2})$, proportional to the surface tension $S_T$. Before going further, it is worth noticing that, for an epithelium advancing in one direction in an infinite experimental set-up, the Galilean invariance is verified for this set of equations and in particular for equation (4). Indeed a change of Galilean frame induces a substrate velocity $\vec{V}_s$ and a convection term $\vec{V}_s \cdot \vec{\nabla} C$ must be added into the left-hand side of equations (2). Also a modification of $\vec{V}_i$ into $\vec{V}_i - \vec{V}_s$ is required into the boundary condition (3). Because it makes more complex the writing of equations, the choice of the laboratory frame is often preferred as in equation (5), destroying apparently (but only apparently) this fundamental symmetry. Changing the velocity $\vec{V}$ into $\vec{V} + V_s \vec{e}_y$ and transforming $\rho$ and $C$ into $\rho(X, Y - \vec{V}_s t, t)$ and $C(X, Y - \vec{V}_s t, t)$ restore the initial equations (2) and (4) by a simple change of coordinates, confirming the Galilean invariance of the model.

**Free boundary problem and stability analysis.** The usual step for the analysis of such free-boundary problem consists in the search of a simple solution moving with a constant velocity and the study of its stability. For a steady traveling front, Fig. 2, left panel, moving with the velocity $U$ along the $y$ axis, both concentrations $C$ vary exponentially in the co-moving frame of velocity $U$. Elementary algebra gives for the epithelium:

$$ \mathcal{C}_e^{(0)} = C_0 e^{r_e y} \quad \text{and} \quad \mathcal{C}_w^{(0)} = (C_0 - 1) e^{-Uy} \qquad (6) $$

where $r_e = \frac{1}{2}(-U + \sqrt{(U^2 + 4)})$ and $C_0 = r_e U(1 + N_0)$. The pressure results from proliferation and chemotaxis:

$$ P^{(0)} = P_i e^{\alpha y} - \Lambda C_0 r_e^2 \frac{e^{\alpha y} - e^{r_e y}}{\alpha^2 - r_e^2} \qquad (7) $$

which gives the interface velocity through an implicit relationship, function of the chemotactic forcing $\Lambda$ and the proliferation rate $\alpha$:

$$ U\{1 - \Lambda r_e^3(1 + N_0)/(r_e + \alpha)\} = -\alpha P_i \qquad (8) $$

Figure 3 summarizes the results for a particular choice of the parameters. The density (Fig. 3A) and pressure fields vary exponentially with the distance from the front, only the morphogen concentration is given, having its maximum at the front. Notice that fast fronts maintain the morphogen concentration on large distances. As $U$ decreases, the morphogen concentration is more localized around the leading edge but for higher velocities (so higher cellular proliferation), one may observe a spreading of the morphogen concentration. In the water, the decrease is abrupt, but this event does not affect the cells. Figure 3B shows the front velocity as a function of the combined parameters $\Lambda(1 + N_0)$ and $\alpha P_i$. From equation (8), we recover the two limiting cases: A front driven only by chemotaxis without proliferation, having velocity: $U = (\Lambda(1 + N_0) - 1)/\sqrt{\Lambda(1 + N_0)}$ and a proliferative front with velocity: $U = -\alpha P_i$. Let us consider now the stability of such fronts. The stability analysis of such
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (A) The growth rate Ω versus the wavenumber $k_0$ when U and Λ vary at fixed values: α = N₀ = 1 and σ = 10⁻³. Note that the chemotactic flux (proportional to Λ) is destabilizing contrary to the advective one, proportional to U. (B) Selected wavenumber $k_0$ versus the front velocity U for σ = 10⁻³, N₀ = 1 and α is varied from 3 to 0.1. (C): α is fixed to 0.5, N₀ varies from 3 to 0.1.</figcaption></figure>

fronts is performed by introducing a sinusoidal perturbation induced by a wavy frontier of the cellular domain: ζ = Ut + εeikx eΩt (in analogy with the Fig. 2, middle panel) and expanding all fields according to:

$$C_e = C_e^{(0)} + \epsilon\chi_e e^{ikx}e^{l\delta}e^{\Omega t}; \quad C_w = C_w^{(0)} + \epsilon\chi_w e^{ikx}e^{-l\omega}e^{\Omega t}; \quad P = P^{(0)} + \epsilon p(y)e^{ikx}e^{\Omega t}. \quad (9)$$

At linear order in ε, we have:

$$l_w = (U + \sqrt{U^2 + 4(k^2 + \Omega)})/2; \quad l_e = (-U + \sqrt{U^2 + 4(k^2 + \Omega + 1)})/2 \\ \chi_e = \frac{-U(N_0 + 1)r_e + N_0(\Omega - U(l_w - U))}{l_e + l_w} \quad (10)$$

where the growth rate Ω reads:

$$\Omega = F(\Omega, k) = -k_t k^2 \sigma - U(k_t - \alpha) + \Lambda \left( \frac{r_e^4 U(N_0 + 1)}{\alpha + r_e} + \frac{\chi_e (l_e^2 - k^2)}{k_t + l_e} \right) \quad (11)$$

with $k_t = \sqrt{k^2 + \alpha^2}$. The dispersion relation, equation (11), which predicts stable steady fronts for Ω < 0, is not easy to analyze owing to the implicit relation for the growth rate Ω, $l_w$, $l_e$ and $\chi_e$ also being functions of Ω. However, one easily notices that, proliferation alone (without chemotaxis) makes the front stable, Ω being negative, and the increase of the parameter Λ may induce an instability. Setting $N_0 = \alpha = 1$ and $\sigma = 10^{-3}$, numerical investigations show the existence of stable patterns when Ω remains negative for all k values (purple curve of Fig. 4A). They also detect fully unstable patterns, out of reach of the present analysis (green and red curves where Ω > 0 for a finite interval of k values). For these two cases, we go directly from a stable planar front to a developed instability without periodic wavy patterns (a situation commonly encountered in diffusive instabilities such as dendritic growth¹⁷). Finally, from Fig. 4A, we can predict steady periodic patterns of finite amplitude (blue curve) since the spectrum exhibits the simultaneous cancellation of Ω and $d\Omega/dk$. This requirement of double cancellation is mandatory for observation of steady undulating patterns. From this analysis, we deduce that a steady wavy-front with prescribed wavelength $l_0$ and wavenumber $k_0$ can be observed, travelling with a constant velocity U above a critical value of Λ called $\Lambda_0$. Fixing the front velocity and the parameters $\sigma$, $N_0$, $\alpha$, as shown in Fig. 4B, or Fig. 4C, a pair of solutions for $k_0$ can be numerically found so a pair of Λ values. This pattern at $k_0$ can explain the fingering pattern observed in experiments¹,²,¹². As for the mode observed at $k = 0$, it results from symmetries as the translational and Galilean invariances (called Goldstone mode). It will play an important role for the nonlinear analysis.

This linear analysis can predict a fingering instability of the leading edge. However, at this stage of linear order, it is not sufficient because these two modes $k = 0$ and $k = k_0$ compete, in the nonlinear regime, thus requiring a careful nonlinear study. Moreover, the relative disorder observed in moving epithelia, leading to fingering and noisy fronts, also suggests that a stochastic study may be helpful to understand the observed patterns.
<figure><img src="image_8.png" /><figcaption>Figure 5. Phase-diagram Λ versus the averaged front velocity U for the parameter values corresponding to Fig. 4, on left: α = N₀ = 1, σ = 10⁻³. The fully stable domain is in white (without noise), in green the estimated weakly nonlinear domains with steady advancing periodic patterns (represented below) and Hopf bifurcation, in yellow, the domain of spatio-temporal pattern, not treated here as the fully unstable domain represented in red. In the inset, schematic representation of the steady pattern, the pressure at the interface is fixed to Pᵢ = -6.76, giving U = 7 and k₀ = 5.</figcaption></figure>

**Nonlinearity and Stochasticity.** The previous analysis was restricted to weak perturbations of the planar front represented in Fig. 2, left panel leading to oscillations of tiny amplitude represented in Fig. 2, middle panel when the parameter Λ is close to Λ₀. We consider now the disorder which always exits in an experiment with cells, especially when they are cancerous. If many theoretical works have been devoted to patterns with noise in experiments of hydrodynamic such as the Rayleigh-Benard³⁶,³⁷, it is less true for free boundary problems, except for the side-branching instability in dendritic growth³⁸⁻⁴⁰. The noise originating from the biochemical signals will modify equation (3) and equation (5) as explained in Section **Stochastic fronts**. However, the analysis of stochasticity will be made easier if we can establish an effective weakly nonlinear equation, close to the threshold Λ₀ and wavenumber k₀. Indeed, when Λ steps across the threshold Λ₀, one can approximate locally the dispersion relation, equation (11), by:

$$ \Omega \sim -\tilde{k}^2 k_0^6 \{-\mu + (\tilde{k}^2 - 1)\}^2 \quad \text{with} \quad \tilde{k} = k/k_0 \qquad (12) $$

where μ represents the gap from threshold: μ = (Λ - Λ₀)/(Λ₀k₀⁴). This suggests a nonlinear equation for the interface velocity deviation u (such as V = U₀ + u) which recovers the dispersion relation at threshold (equation (12)):

$$ \frac{\partial u}{\partial \tilde{t}} = \frac{\partial}{\partial x} \left( \frac{\partial u}{\partial x} (-\mu + 1) + 2 \frac{\partial^3 u}{\partial x^3} + \frac{\partial^5 u}{\partial x^5} - \frac{1}{2} u^2 \right) \qquad (13) $$

where $\tilde{t} = k_0^6 t$. How to select the nonlinearity in equation (13)? Above threshold, the first nonlinearity is $-\beta^2 u^2$ where β is a numerical value and the negative sign indicates a saturation of the velocity amplitude. There is no symmetry $u- > -u$, since the velocity correction u cannot behave symmetrically in the front direction (upstream) and in the opposite direction (downstream). Defining $u = \tilde{u}/\beta$ eliminates the unknown constant β from the equation for $\tilde{u}$ and we recover equation (13). However a more fancy way consists in applying the Galilean invariance symmetry $(1/2\partial u^2/\partial x): x \to x - vt$ and $u \to u + v$. Such consideration allows to avoid the very tedious nonlinear analysis from the free-boundary equations. This quadratic nonlinearity also imposes the scaling of the parameters: Choosing $|\mu|$ of order $\epsilon^2$, from the dispersion relation, we deduce that the lengths scale as $1/\epsilon$ while, from the right-hand-side of equation (13), the velocity u scales also as $1/\epsilon$. The left-hand-side gives us the time-scale as $1/\epsilon^2$. As shown in Fig. 1, chemotaxis, which is responsible for diffuse fronts, increases stochasticity and the disorder. Restricting to the biochemical noise in our model, stochasticity enters at the level of the boundary conditions: The jump condition for both concentration and gradient of concentration and the chemotactic migration constant Λ. Stochasticity from the pressure and the friction of the cells at the substrate are neglected. We are then face simultaneously to a multiplicative noise η₁, indicating a noisy threshold for bifurcation and an additive noise η₂ arising from all other sources. Both sources of noise are chosen Gaussian in space, white in time and are not correlated so $\langle \eta_i(x, t) \eta_j(x', t') \rangle = C_i(x) \delta_{ij} \delta(t-t') \delta(x-x')$. Restricting on the vicinity of the interface, close but below the bifurcation threshold μ < 0, (which corresponds to the white zone of parameters in Fig. 5), we transform equation (13) into a stochastic equation adding η₁ and η₂. It reads:
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left panel (A): Domains of spatial wavelengths $q = (k - k_0)/k_0$ in the nonlinear regime defined by $\mu = (\Lambda - \Lambda_0)/\Lambda_0$. In grey, unstable modes, in green stable modes. The left-hand frontier between grey and green domains corresponds to the Hopf bifurcation. Right panel (B): Schematic representation of the interface profile $\zeta$ and the average front velocity, for $\mu = 0.1$. The interface oscillates in time between the blue and red curves, the curve in purple indicates the mean velocity.</figcaption></figure>

$$ \frac{\partial u}{\partial \tilde{t}} = \mathcal{L}(\mu, u) + \frac{\partial}{\partial x}(\eta_1 u + \eta_2) \qquad (14) $$

where $\mathcal{L}(\mu, u)$ is the differential operator, right-hand-side of equation (13). The correspondence between the biochemical noise and $\eta_1$ and $\eta_2$ is established in Section **Stochastic fronts**. Below the instability threshold, the velocity $u$ has only a stochastic component and neglecting the non-linearity and the multiplicative noise contributions, the velocity structure function reads $S(\tilde{k}) = C_2(\tilde{k})/(2(\tilde{k}^2 - 1) - \mu)$, (see Section **Stochastic fronts**), a result similar to the linearized Swift-Hohenberg equation derived in refs 36 and 41, indicating that periodic patterns can be observed below the threshold $\mu = 0$. The effect of the multiplicative noise can be approximated by using Novikov's theorem³⁷ leading to a new contribution proportional to $C_1S(\tilde{k})$ which modifies the threshold $\mu$ becoming $\mu + C_1$. It contributes to a possible jump into the unstable domain even if the parameter $\mu$ is negative. Finally the quadratic nonlinearities do not contribute, in opposition to the Swift-Hohenberg equation⁴¹. This analysis allows to conclude that below the threshold of instability of the planar front, it is possible to observe a wavy pattern induced by a small level of noise. However, above the threshold, noise may play a role only when it competes with the nonlinear effects so very close to the onset. If its level is too low, nonlinearities dominate.

**Nonlinearity and the Turing-Hopf bifurcation.** For positive $\mu$ values, the velocity $u$ results from the competition between both modes at $\tilde{k} = 0$ and at $\tilde{k} = 1$, respectively, which are coupled by the nonlinear terms. This competition suggests the following decomposition $u(x, \tilde{t}) = 2\varepsilon Re[6A(X, T)e^{iX} + \varepsilon B(X, T)]$, B representing a small change in the front velocity, $X$ and $T$ being rescaled variables: $X = \varepsilon x$ and $T = \varepsilon^2 k_0^6 t$. A classical multi-scale analysis⁴¹,⁴² gives us the coupled amplitude equations for both modes A and B:

$$ \left\{ \begin{aligned} \partial_T A &= (1 + 4\partial_X^2 - |A|^2)A - i(AB) - \partial_X(AB) \\ \partial_T B &= \partial_X^2 B - 36\partial_X |A|^2 - B\partial_X B \end{aligned} \right. \qquad (15) $$

The last terms of the first equation (15) represent the advection of the pattern, $\tilde{k} = 1$, by the change in the front velocity $B$. We also notice that $B$ has only a nonzero value for inhomogeneous values of the amplitudes $A$ and $B$ (see Section **Hopf-Bifurcation**). Equations (15) are similar to those established previously⁴³,⁴⁴ and such coupling has been called Hopf-Turing bifurcation⁴². Exact solutions for $A$ and $B$ in equation (15) can be easily found: $A = A_0 e^{iq^X}$ with amplitude: $A_0 = \sqrt{1 - 4q^2}$ and $B = 0$ (see the domain in green in Fig. 6A). A stability analysis of these solutions is necessary knowing that their domain of existence is constrained by the inequality $q^2 < 1/4$. Once perturbed according to $A = (A_0 + r)e^{iq^X+\phi}$ and $B = \psi$, where $r$, $\phi$ and $\psi$ are small real functions of wavenumber $p$ and growth rate $\eta$, an algebraic equation is established (see Section **Hopf-Bifurcation**) for the growth rate $\eta$ giving modes growing exponentially in time for $q > 0$, a range of stable modes for negative $q$ values, bounded by a Hopf bifurcation for $q = q_h$ (see the schematic representation of the analysis in Fig. 6A where unstable modes are in the grey domain, stable modes are in the green domain and the Hopf-Turing modes are at the frontier on left). For $q_h = 18/13(1 - \sqrt{1 + 533/(36^2)}) \sim -0.260$, waves appear and superpose to the static spatial mode $q_h$ affecting simultaneously the spatially periodic oscillation A and the average velocity $B$ of the front. The interface oscillates between the blue and red curves of Fig. 6B, but also the average velocity of the front (represented in purple). For a number $p$ of the wave, the frequency for the time-dependence is equal to $\omega_h p$ with $\omega_h = 12\sqrt{2}|q_h|(1 + q_h) \sim 7.446$.
The nonlinear analysis proves that there exists a band of stationary spatial oscillations in the neighbourhood of $\tilde{k} = 1$, with amplitude $\sqrt{1 - 4q^2}$ and wavenumber $k_0(1 + \sqrt{\mu}q)$, where $q_h < q < 0$. At the frontier, $q = q_h$, the frequency of the time dependent oscillations becomes $\mu k_0^6 \omega_h p$. Also the averaged front velocity oscillates with the same frequency $\mu k_0^6 \omega_h p$. Let us evaluate an order of magnitude for the oscillation period $T_h$ and compare our results with measurements performed in ref. 16. From the evaluation of the ratio between the amplitude of the fingering and the spatial wavelength$^{16}$, $\sqrt{\mu}$ is estimated as 0.1/12 (the factor 12 comes from the definition of A). The wave-number $p$ can be deduced from the measured strain-rate $\dot{\epsilon} \sim 410^{-5}s^{-1}$, which identifies to $\dot{\epsilon} = \sqrt{\mu}p^2\omega_h/(2\tau_c)$. As shown in ref. 45, $\tau_c$ varies between 45 minutes to two hours while $D_c = 0.3(\mu m)^2/s$. Taking the lowest estimation for $\tau_c$, we find $p \sim 0.15$ which gives a period $T_h = 400s$ for $k_0 \sim 7.2$. From Fig. 4, on right, such $k_0$ value indicates a velocity $U \sim 10$, giving a front velocity $10(D_c\tau_c)^{1/2} \sim 0.1\mu m/s$ while the wavelength is about 25µm. These estimations are in good agreement with the experimental results$^{16}$. A phase-diagram ($\Lambda_0, U$) is presented in Fig. 5 for fixed values of $N_0 = \alpha = 1$ and $\sigma = 10^{-3}$ which shows the domain of a planar front (with velocities oriented towards the leading edge), the homogeneous oscillatory patterns combined eventually, with a Hopf bifurcation (in yellow). We show also the domain of full instabilities with spatio-temporal dynamics (in green) and the fully developed instabilities (in red) which are not studied here.

The main conclusion of this analysis is that temporal global modes of oscillations is the result of the Galilean invariance of the experimental set-up, which occurs in nonlinear systems as soon as their dispersion relation presents a similar behaviour as Fig. 4.

## Discussion

Epithelial monolayer expansion is mostly regarded as the interplay between mechanical forces combined with stochasticity. The samples involved in many experiments$^{1,2}$ have sizes of order 100µm and perhaps cannot be compared with the example shown in Fig. 1 which is of order centimeters and lasts for several days. It is not an easy task to quantify precisely the magnitude of auto-chemotaxis and its contribution to motility. In addition, it probably depends on the nature of the cells. However it deserves to be studied as a possible mechanism for *in vitro* tissue motility. Our macroscopic model averages the mechanical forces acting at the cellular level which is justified by the coordination existing among cells in the bulk, far from the border. The model suggests a different scenario via the exchange of chemical factors made by the cells themselves and gives front undulations and front temporal oscillations. These waves result from the combined effects of the translational invariance (responsible for the mode $k=0$), the existence of a spatial undulation (mode $k=k_0$) detected by the linear stability analysis, equation (11), and the nonlinearities deduced from the Galilean invariance. It is different from other instabilities where waves appear at linear order via a complex growth rate $\Omega$ in the dispersion relation$^{46}$. Although we cannot quantify the level of chemical signaling, any diffusive front is a strong indication of auto-chemotaxis$^{2,6,30}$. The experiment shown in Fig. 1 does not respect the translational invariance required by the model but it shows a diffuse front. Once the initial geometry is modified to correspond to the free planar growth, some aspects of the theory can be verified since the model gives analytical predictions. Different parameters can vary such as the mobility coefficient $K_*$ (by changing the substrate) or the proliferation rate (by adding a biocompatible polymer, dextran, to the solution). It will allow to explore the phase diagram shown in Fig. 6. It has been shown that dextran induces a controlled osmotic pressure$^{12,47}$ so allows to vary the interface pressure $P_i$. Among quantitative measurements, the velocity structure function, predicted here equation (21), is an accessible quantity in most of the experimental set-ups$^{1,2}$. A more systematic study of the waves, detected above threshold in a specific domain of parameters, will be highly valuable.

## Methods

**Stochastic Fronts.** Here, we rely the additive and multiplicative noise term $\eta_1$ and $\eta_2$ of equation (14) to the biochemical and physical sources of noise. Stochasticity has been introduced in nonlinear systems, near threshold of bifurcations via effective equations called stochastic normal forms. The most studied equation is probably the stochastic Swift-Hohenberg equation (for a review, see ref. 37). As for our front nonlinear equation (13), the main advantage is the simplicity of these approaches compared to the full hydrodynamic equations. Neglecting the nonlinearities, the Fourier transform of equation (14) gives

$$ \frac{du_k}{dt} = \tilde{k}^2 (\mu - (\tilde{k}^2 - 1)^2) u_k(t) + i\tilde{k}\eta_k(t) \quad (16) $$

where, the noise correlation function, for a white noise in time but coloured in space $\eta_2$ is: $\langle \eta_k(t) \eta_k(t') \rangle = \delta(\tilde{k} + \tilde{k}') C(k) \delta(t' - t)$. We define $\gamma(\tilde{k}) = (-\mu + (\tilde{k}^2 - 1)^2)$, which is positive below the threshold of instability ($\mu$ being negative). The velocity correlation function, $S(\tilde{k}, t, t') = \langle u_k(t) u_{-k}^*(t') \rangle$, where $u^*$ denotes the complex conjugate of $u$, is evaluated explicitly in the following. Indeed, from equation (16), we derive:

$$ u_k(t) = u_0 e^{-\tilde{k}^2 \gamma(\tilde{k}) (t-t_0)} + i \tilde{k} \int_{t_0}^{t} \eta_k(\tau) e^{-\tilde{k}^2 \gamma(\tilde{k}) (t-\tau)} d\tau \quad (17) $$

The long time behaviour, $t_0 \to -\infty$, cancels the first term in the right-hand side of equation (17). The function $S(\tilde{k}, t, t')$ satisfies

$$ \frac{dS(\tilde{k}, t, t')}{dt} = -\tilde{k}^2 \gamma(\tilde{k}) S(\tilde{k}, t, t') + i\tilde{k} \langle \eta_k(t) u_{-k}(t') \rangle \quad (18) $$
In order to evaluate the quantity $\langle \eta_k(t)u_{-k}(t') \rangle$, from equation (17), we derive:

$$ 
\begin{aligned}
\mathcal{I} &= \langle \eta_k(t) u_{-k}(t') \rangle = -i\tilde{k} \int_{-\infty}^{t'} \langle \eta_k(t) \eta_{-k}(\tau) \rangle e^{-\tilde{k}^2 \gamma(\tilde{k})(t'-\tau)} d\tau \\
&= -i\tilde{k} \int_{-\infty}^{\infty} \langle \eta_k(t) \eta_{-k}(\tau) \rangle \theta(t'-\tau) e^{-\tilde{k}^2 \gamma(\tilde{k})(t'-\tau)} d\tau
\end{aligned} 
\quad (19) $$

where $\theta(X)$ is the Heaviside step function, equal to 1 for $X>0$, to 0 for $X<0$, and to 1/2 for $X=0$. At this step any coloured noise in time and space can be introduced as example a Gaussian distributed noise with zero mean. However, to simplify, we choose a noise term with zero mean, white in time. Accordingly, it reads:

$$ 
\left\{
\begin{array}{ll}
\text{If } t > t', & \mathcal{I} = 0 \quad \text{and} \quad t < t', \quad \text{then} \quad \mathcal{I} = -i\tilde{k}C_a(\tilde{k})e^{-\tilde{k}^2\gamma(\tilde{k})(t'-t)} \\
\text{If } t = t', & \text{we have } \mathcal{I} = -i\tilde{k}C(\tilde{k})/2.
\end{array}
\right. 
\quad (20) $$

At long times, the structure function reads:

$$ S(\tilde{k}, t' - t) = \langle u_k(t) u_{-k}(t') \rangle = \frac{C(\tilde{k}) e^{-\tilde{k}^2 \gamma(\tilde{k})(t'-t)}}{2\gamma(\tilde{k})} \quad (21) $$

Close to the threshold, given by $\mu=0$, for negative values of $\mu$, the only root of $\gamma(\tilde{k})$ is $\tilde{k}=1$, which increases the velocity structure function for this value. Therefore, at linear order, fingering with $\tilde{k}=1$ is selected by the noise below the hydrodynamic threshold. However, at threshold, the structure function diverges for $\tilde{k}=1$ and nonlinearity cannot be discarded. Let us identify the origin of the multiplicative and additive noises for our specific free boundary problem.

**Stochasticity and free boundary problems.** The correspondence between the biochemical or biophysical noise and the stochastic variables introduced into equation (14) requires to come back to the free boundary problem, where the nature of the noise is well identified. The interplay between nonlinearities and stochasticity is a challenging theoretical problem for fully developed instabilities$^{38-40}$ but becomes easier for patterns of finite amplitude like the one investigated here. However the origin of noise and its nature (additive or multiplicative) have to be elucidated. Patterns result from the interplay between fields, solutions of partial differential equations and initial boundary conditions. These problems give rise to strong nonlinear instabilities and noise is not always a pertinent parameter, except close the threshold of instability. Focussing on chemotaxis, the chemical production of morphogenetic gradient is surely a source of noise, however, the cellular activity such as mitosis is also a noisy mechanism. Here we restrict to noise arising from the morphogen production. The morphogens, occurring directly at the border, may affect more the fingering instabilities. Mathematically, morphogens enter into the chemotatic migration constant $\Lambda$ and the two boundary conditions in equation (3). This leads to three stochastic variables that we consider uncorrelated.

**The multiplicative noise.** The multiplicative noise is easily identified. The phase diagram has demonstrated that a possible control parameter is the average front velocity $U$, equation (8), which is itself related to the chemotactic flux. In equation (8), stochasticity appears through the parameter $N_0$ transformed into $N_0+\eta^{(4)}$, and the chemotactic migration constant $\Lambda$ which becomes $\Lambda+\eta^{(5)}$. A straightforward calculation shows:

$$ \eta_1 = - \frac{\Lambda \eta^{(4)} + (1 + N_0) \eta^{(5)}}{\Lambda(1 + N_0)} Q(U) \quad \text{with} \quad Q(U) = \frac{(U-1)(r_e + \alpha)\sqrt{U^2 + 4}}{\sqrt{U^2 + 4}(\alpha + r_e) + U(U-1)} \quad (22) $$

**Numerical estimation.** For our choice of parameters, $\alpha=N_0=1$, $\Lambda=7.65$ and $U=7$, $\eta_1$ can be estimated and gives $\eta_1 = -(0.5\eta^{(4)} + 0.26\eta^{(5)})$ with $Q(U) \sim 0.99$.

**The additive noise noise.** Each boundary condition (equation (3) and equation (5)) is modified by a stochastic contribution $\eta$:

$$ C_e - C_w = 1 + \eta^{(3)} \quad \text{and} \quad \vec{N} \cdot (\vec{\nabla} C_e - \vec{\nabla} C_w) = (N_0 + \eta^{(4)}) \vec{N} \cdot \vec{V}_i \quad (23) $$

We will solve the free-boundary problem linearly, by using the Fourier transform for the three noise sources $\eta^{(3)}$, $\eta^{(4)}$ and $\eta^{(5)}$ as follows:

$$ \eta(x, t) = \int_{-\infty}^{\infty} d\omega \int_{-\infty}^{\infty} d\lambda \, \eta(\lambda, \omega) e^{i\omega t} e^{i\lambda x} \quad (24) $$

and for the interface:

$$ \zeta(x, t) = \int_{-\infty}^{\infty} d\omega \int_{-\infty}^{\infty} d\lambda \, \zeta(\lambda, \omega) e^{i\omega t} e^{i\lambda x} \quad (25) $$

The calculation is made simpler if we do not consider the time variation of the noisy terms in both diffusion equations and we derive $\eta$
$$ \zeta(\lambda, \omega) = \frac{\Lambda(\lambda^2 - r_{ep}^2)(\eta^{(3)}(\lambda, \omega)r_{wp} + \eta^{(4)}(\lambda, \omega)U)}{(k_t + r_{ep})(r_{ep} + r_{wp})D_{\eta}} + \frac{U(1 + N_0)(\lambda^2 - r_e^2)r_e\eta^{(5)}(\lambda, \omega)}{(k_t + re)D_{\eta}} \quad (26) $$

with $D_{\eta}$

$$ D_{\eta} = -i\omega \left( 1 + \frac{\Lambda N_0 (\lambda^2 - r_{ep}^2)}{(k_t + r_{ep})(r_{ep} + r_{wp})} \right) - \gamma \lambda^2 k_t - (k_t - \alpha) U + \Lambda \frac{(1 + N_0) r_0^4 U}{\alpha + r_0} + \Lambda \frac{(r_{ep}^2 - \lambda^2) X_0}{k_t + r_{ep}} \quad (27) $$

where $r_{ep} = (-U + \sqrt{U^2 + 4(\lambda^2 + 1)})/2$, $r_{wp} = (U + \sqrt{U^2 + 4\lambda^2})/2$ and $X_0 = (-U(N_0 + 1)r_0 - N_0 U(r_{wp} - U)/(r_{ep} + r_{wp})$. Using equation (11), we transform $D_{\eta}$ as follows:

$$ D_{\eta} = -i\omega \left( 1 + \frac{\Lambda N_0 (\lambda^2 - r_{ep}^2)}{(k_t + r_{ep})(r_{ep} + r_{wp})} \right) + \tilde{\Omega}(\lambda) \quad (28) $$

and we introduce

$$ \Gamma(\lambda) = \frac{-\tilde{\Omega}}{1 + \Lambda N_0 (\lambda^2 - r_{ep}^2) / \{(k_t + r_{ep})(r_{ep} + r_{wp})\}} \quad (29) $$

Remember that $\tilde{\Omega}$ remains negative below threshold. Then, the noisy interface equation (26), now reads:

$$ \zeta(\lambda, \omega) = \frac{E_3(\lambda)\eta^{(3)} + E_4(\lambda)\eta^{(4)} + E_5(\lambda)\eta^{(5)}}{-i\omega + \Gamma(\lambda)} \quad (30) $$

where $E_3$, $E_4$ and $E_5$ are given by:

$$ E_3(\lambda) = \frac{\Lambda(r_{ep}^2 - \lambda^2)r_{wp}}{\Lambda N_0(\lambda^2 - r_{ep}^2) + (k_t + r_{ep})(r_{ep} + r_{wp})} \quad (31) $$

$$ E_4(\lambda) = \frac{\Lambda(r_{ep}^2 - \lambda^2)U}{\Lambda N_0(\lambda^2 - r_{ep}^2) + (k_t + r_{ep})(r_{ep} + r_{wp})} \quad (32) $$

$$ E_5(\lambda) = \frac{(k_t + r_{ep})(r_{ep} + r_{wp})}{k_t + r_e} \frac{U(1 + N_0)(\lambda^2 - r_e^2)r_e}{\Lambda N_0(\lambda^2 - r_{ep}^2) + (k_t + r_{ep})(r_{ep} + r_{wp})} \quad (33) $$

The Langevin equation for the front position. It is more easy to analyze the free-boundary problem with noise by writing the equivalent Langevin equation for equation (30):

$$ \frac{\partial}{\partial t} \zeta(\lambda, t) = -\Gamma(\lambda)\zeta(\lambda, t) + E(\lambda)\eta(\lambda, t) \quad (34) $$

Since the three sources of noise in equation (30) are independent, we can treat each of them separately. Applying the time Fourier transform to equation (34), we recover the same structure as equation (30), which facilitates the identification of the noise sources occurring in the free boundary problem and in equation (14). Equation (34) can be integrated in time giving:

$$ \zeta(\lambda, t) = \zeta_0 e^{-\Gamma(\lambda)t} + E(\lambda) \int_0^t \eta(\lambda, t') e^{-\Gamma(\lambda)(t-t')} dt' \quad (35) $$

Since $\langle \eta(\lambda, t) \rangle = 0$, then the averaged value of the front position in the moving frame vanishes with time: $\langle \zeta \rangle = \zeta_0 e^{-\Gamma(\lambda)t}$ whether the noise appears at time $t=0$. Let us evaluate the auto-correlation function:

$$ \sigma_{\zeta}(\tau) = \langle \delta\zeta(\lambda, t) \delta\sigma(\lambda, t+\tau) \rangle = E(\lambda)^2 e^{-\Gamma(\lambda)(2t+\tau)} \int_0^t dt_1 \int_0^{t+\tau} dt_2 \langle \eta(\lambda, t_1) \eta(\lambda, t_2) \rangle e^{\Gamma(\lambda)(t_1+t_2)} \quad (36) $$

If we take a white noise in time but coloured in position, we have $C(\tau) = C_0(\lambda)\delta(\tau)$ and then:

$$ \sigma_{\zeta}(\tau) = \frac{C(\lambda)}{\Gamma(\lambda)} E(\lambda)^2 e^{-\Gamma(\lambda)\tau} (1 - e^{-2\Gamma(\lambda)t}) \quad (37) $$

Clearly $\sigma_{\zeta}(\tau)$ increases in the neighborhood of $\lambda = k_0$ and the noise spatial frequency is selected even below the threshold. In the neighbourhood of $\lambda = k_0$ a regular pattern induced by noise is expected. Not surprisingly, near $\lambda \sim k_0$, the auto-correlation function of the position becomes linear in $t$. In addition the structure function, which
corresponds to $S(\lambda) = \sigma_c(0)$ and $t \to \infty$, can be evaluated in the vicinity of $\lambda \sim k_0$. Indeed, once $\Gamma(\lambda)$ is approximated by $\Gamma(\lambda) = \Gamma_0(\lambda^2 - k_0^2)^2 - \mu$, it reads:

$$S(\lambda) = \frac{C_3(k_0)E_3(k_0)^2 + C_4(k_0)E_4(k_0)^2 + C_5(k_0)E_5(k_0)^2}{2\Gamma_0 k_0^4 (\lambda^2/k_0^2 - 1)^2 - \mu} - \frac{C_2(\tilde{k})}{2((\tilde{k} - 1)^2 - \mu)} \quad (38)$$

where $\mu$ is the negative gap from the hydrodynamic threshold of instability. The relationship between the noise introduced in the nonlinear equation (14) and the biochemical noise is now established. Obviously $\eta_2$ turns out to be a nontrivial combination of the biochemical sources produced at the frontier.

**Hopf-Bifurcation.** The non-linear front velocity equation (13) cannot be solved naively with a Fourier series expansion $\sum_m a_m e^{imx}$, where $m$ is integer. Indeed, because of the nonlinear term $-\frac{1}{2}\mu^2$, the related series will include secular terms which produce its divergence$^{48}$. In order to eliminate these terms, one choice consists in considering the following decomposition$^{42}$:

$$u(x, t) = 2\varepsilon Re[6A(X, T)e^{ix} + \varepsilon B(X, T)] \quad (39)$$

where $X = \varepsilon x$, $T = \varepsilon^2 k_0^6 t$, $\mu = \varepsilon^2$, and $\mu = \varepsilon = 0$ at the instability threshold, which transforms the nonlinear equation (13) into the two following amplitude equations:

$$\frac{\partial A}{\partial T} = (1 - |A|^2)A + 4 \frac{\partial^2 A}{\partial X^2} - i(AB) - \frac{\partial(AB)}{\partial X} \quad (40)$$

$$\frac{\partial B}{\partial T} = \frac{\partial^2 B}{\partial X^2} - 36 \frac{\partial |A|^2}{\partial X} - \frac{1}{2} \frac{\partial B^2}{\partial X} \quad (41)$$

The above equations accept trivial solutions: $A = \sqrt{1 - 4q^2} e^{iqX}$ and $B = 0$, which indicates a band of possible wavenumbers $q$ in the close vicinity of $\tilde{k} = 1$. In order to study the dynamical stability of these solutions, we define $A = (A_0 + r)e^{i(qX+\phi)}$ and $B = \psi$ where $r, \phi, \psi$ are small real quantities of $X$ and $T$ compared to $A_0$. The linearization of equations 40 and 41 gives:

$$\left\{
\begin{aligned}
\frac{\partial r}{\partial T} &= -2A_0^2 r + 4 \frac{\partial^2 r}{\partial X^2} - 8A_0 q \frac{\partial \phi}{\partial X} - A_0 \frac{\partial \psi}{\partial X} \\
\frac{\partial \phi}{\partial T} &= 4 \frac{\partial^2 \phi}{\partial X^2} + \frac{8q}{A_0} \frac{\partial r}{\partial X} - (1+q)\psi \\
\frac{\partial \psi}{\partial T} &= \frac{\partial^2 \psi}{\partial X^2} - 72A_0 \frac{\partial r}{\partial X}
\end{aligned}
\right.
\quad (42)$$

In the long wavelength limit, expanding $r, \phi, \psi$ as $e^{i\eta t+ipx}$ ($p$ being small) and solving the above system, one finally obtains:

$$\eta^2 + c_1 \eta p^2 + c_2 p^4 = 0 \quad (43)$$

with $c_1 = (41 + 144q - 52q^2)/(1 - 4q^2)$ and $c_2 = -288q(q+1)$. The two roots for the growth rate $\eta$ are then

$$\eta = \frac{1}{2}(-c_1 p^2 \pm p \sqrt{c_1^2 p^2 - 4c_2}) \quad (44)$$

When $c_2 < 0$, $q > 0$, equation (43) accepts two real solutions of opposite signs, giving a mode with exponential growth. Then, this domain $q > 0$ is forbidden. If $q < 0$ and $c_2 > 0$, two complex roots exit and the stability of the front is thus given by the sign of $c_1$ which has two roots: $q_h = 18/13(1 \pm \sqrt{1 + (533/(36^2))})$, one being negative $q_h \sim -0.260$. In the interval of negative $q$ between $q_h$ and 0, $c_1$ is positive, so static oscillatory modes exist around $\tilde{k} = 1$. This stability stops at $q_h$ but for this particular value, one reaches a Hopf bifurcation with a time dependent oscillation with frequency $\omega_h = 12p\sqrt{2|q_h|(1+q_h)} \sim 7.4463p$.

## References

1. Trepat, X. et al. Physical forces during collective cell migration. *Nat. Phys.* **5**, 426–430 (2009).

2. Sepulveda, N. et al. Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model. *PLoS Comput. Biol.* **9**(3), e1002944 (2013).

3. Koepf, M. H. & Pismen, L. M. A continuum model of epithelial spreading. *Soft Matter* **9**, 3727–3734 (2013).

4. Martin, P. Wound healing-Aiming for perfect skin regeneration. *Science* **276**, 75–81 (1997).

5. Ben Amar, M. & Wu, M. Re-epithelialization: advancing epithelium frontier during wound healing. *Roy. Soc. Interface* **11**, 20131038, 1–7 (2014).

6. Wu, M. & Ben Amar, M. Growth and remodelling for profound circular wounds in skin. *Biomech Model Mechanobiol.* **14**, 357–370 (2015).

7. Bowden, L. G. et al. A morphoelastic model for dermal wound closure. *Biomech Model Mechanobiol.* **15**, 3663–3681 (2016).

8. Bockhorn, M., Jain, R. K. & Munn, L. J. Active versus passive mechanisms in metastasis: do cancer cells crawl into vessels, or are they pushed? *Lancet* **8**, 444–448 (2007).

9. Muinonen-Martin, A. J. et al. Melanoma Cells Break Down LPA to Establish Local Gradients That Drive Chemotactic Dispersal. *PLoS Biol* **12**, e1001966 (2014).
10. Roca-Cusachs, P., Sunyer, R. & Trepat, X. Mechanical guidance of cell migration: lessons from chemotaxis. *Curr Opin Cell Biol.* **25**, 543–549 (2013).

11. Harris, H. Role of chemotaxis in inflammation. *Physiol. Rev.* **34**, 529–562 (1954).

12. La Porta, C. A. M. et al. Osmotic stress affects functional properties of human melanoma cell lines. *EPJ Plus* **130**, 64–79 (2015).

13. Etienne-Manneville, S. Neighborly relations during collective migration. *Curr Opin Cell Biol.* **30**, 51–59 (2014).

14. Lushi, E., Goldstein, R. E. & Shelley, M. J. Auto-chemotactic micro-swimmer suspensions: modeling, analysis and simulations. *arXiv:1310.7614 (2013).*

15. Insall, R. The interaction between pseudopods and extracellular signalling during chemotaxis and directed migration. *Curr Opin Cell Biol.* **25**, 526–531 (2013).

16. Serra-Picamal, X. et al. Mechanical waves during tissue expansion. *Nat. Physics* **8**, 628–634 (2012).

17. Langer, J. S. Instabilities and pattern formation in crystal growth. *Rev. Mod. Phys.* **52**, 1–30 (1980).

18. Cummings, L. M., Hohlov, Y. E., Howison, S. D. & Kornev, K. Two-dimensional solidification and melting in potential flows. *Journ. Fluid Mech.* **378**, 1–18 (1999).

19. Lindner, A. Coussot, P. & Bonn, D. Viscous Fingering in a Yield Stress Fluid. *Phys. Rev. Lett* **85**, 314–318 (2000).

20. Al-Housseiny, T. T., Tsai, P. A. & Stone, H. A. Control of interfacial instabilities using flow geometry. *Nat Phys.* **8**, 747–750 (2012).

21. Ben Amar, M. Void electromigration as a moving free-boundary value problem. *Physica D* **134**, 275–286 (1999).

22. Laurençot, P. & Walker, C. A stationary free boundary problem modeling electrostatic MEMS. *Archive for Rational Mechanics and Analysis* **207**, 139–158 (2013).

23. Debiais, A., Harich, R., Bonn, D., Colin, A. & Kellay, H. Spreading of an Oil-in-Water Emulsion on a Glass Plate: Phase Inversion and Pattern Formation. *Journ. Fluid Mech.* **469**, 237–256 (2002).

24. Gov, N. S. Traction forces during collective cell motion. *HFSP J.* **3**, 223–227 (2009).

25. Mark, S. et al. Physical Model of the Dynamic Instability in an Expanding Cell Culture. *Biophys J.* **98**, 361–370 (2010).

26. Tweedy, L., Knecht, D. A., Mackay, G. M. & Insall, R. H. Self-Generated Chemoattractant Gradients: Attractant Depletion Extends the Range and Robustness of Chemotaxis. *PLoS Biol.* **14**(3), e1002404 (2016).

27. Taber, L. A. Biomechanics of Growth, Remodeling, and Morphogenesis. *Appl. Mech. Rev* **48**, 487–545 (1995).

28. Puliafito, A. et al. Collective and single cell behavior in epithelial contact inhibition. *Proc. Nat. Acad. Sci. USA* **109**, 739–744 (2012).

29. Caginalp, G. & Xie, W. An analysis of phase-field alloys and transition layers. *Arch. Rational Mech. Anal.* **142**, 293–329 (1998).

30. Ben Amar, M. Collective chemotaxis and segregation of active bacterial colonies. *Sci. rep.* **6**, 21269, doi: 10.1038/srep21269 (2016).

31. Callan-Jones, A., Joanny, J. F. & Prost, J. Viscous-Fingering-Like Instability of Cell Fragments. *Phys. Rev. Lett.* **100**, 258106–258110 (2008).

32. Lowengrub, J. S. et al. Nonlinear modelling of cancer: Bridging the gap between cells and tumors. *Nonlinearity* **23**, 1–91 (2010).

33. Preziosi, L. & Vitale, G. A Multiphase model of tumor and tissue growth including cell adhesion and plastic reorganization. *Math. Mod. Meth. Appl. S.* **21**, 1901–1932 (2011).

34. Balois, T. & Ben Amar, M. Morphology of melanocytic lesions in *situ*. *Sci. rep.* **4**, 3622, doi: 10.1038/srep03622 (2014).

35. Kasyap, T. V. & Koch, D. L. Chemotaxis driven Instability of a confined bacterial Suspension. *Phys. Rev. Lett.* **108**, 03101–104 (2012).

36. Carillo *et al*. Spatial coherence resonance near pattern-forming instabilities. *Europhys. Lett.* **65**(4), 452–458 (2004).

37. Sagués, F., Sancho, J. M. & J. Garcia-Ojalvo, J. Spatiotemporal order out of noise. *Rev. Mod. Phys.* **79**, 829–882 (2007).

38. Langer, J. S. Dendritic side branching in the three-dimensional symmetric model in the presence of noise. *Phys. Rev. A* **36**, 3350–3358 (1987).

39. Karma, A. & Rappel, W. J. Phase-field model of dendritic growth side branching with thermal noise. *Phys. Rev. E* **60**, 3614–3625 (1999).

40. Brener, E. & Temkin, D. Pattern formation in three-dimensional dendritic growth *Physica A* **263**, 338–344 (1999).

41. Cross, M. C. & Hohenberg, P. C. Pattern formation outside of equilibrium. *Rev. Mod. Phys.* **65**, 851–1124 (1993).

42. Meron, E. Hopf-Turing instability in *Nonlinear Physics of Ecosystems*, 85–94 (CRC Press, 2015).

43. Coullet, P. & Fauve, S. Propagative Phase dynamics for systems with Galilean invariance. *Phys. Rev. Lett* **55**, 2857–2861 (1985).

44. Coullet, P. & Spiegel, E. A. Amplitude Equations for Systems with Competing Instabilities. *SIAM J. Appl. Math.* **43**, 775–821 (1983).

45. Wartlick, O., Kicheva, A. & González-Gaitán, M. Morphogen Gradient Formation. *Cold Spring Harb Perspect Biol.* **1(3)** (2009).

46. Plapp, M. & Karma, A. Eutectic colony formation: A stability analysis. *Phys. Rev. E* **60**, 6865–6890 (1999).

47. Montel, L. et al. Stress clamp experiments on multicellular tumor spheroids. *Phys. Rev. Lett.* **107**, 188102–188106 (2011).

48. Bender, C. M. & Orszag, S. A. Singular perturbation analysis in *Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory* 549–554 (Springer, 1999).

## Acknowledgements

This research was supported in part by Institut Universitaire de France and by the National Science Foundation under Grant No. PHY11-25915. It is a pleasure to acknowledge very enlightening discussions with S. Fauve and V. Hakim.

## Author Contributions

M.B.A. and C.B. wrote the main manuscript text and M.B.A. prepared all figures. All authors reviewed the manuscript.

## Additional Information

**Competing financial interests:** The authors declare no competing financial interests.

**How to cite this article:** Ben Amar, M. and Bianca, C. Onset of nonlinearity in a stochastic model for auto-chemotactic advancing epithelia. *Sci. Rep.* **6**, 33849; doi: 10.1038/srep33849 (2016).

<figure><img src="image_10.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/