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

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 k0. 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.

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 e for the morphogens inside the cellular domain Ω e and in the water bath Ω w . In Ω e , the morphogens diffuse but are also captured Figure 1. Human melanocytic cells (primary tumour) were plated on 33 cm 2 Petri dish with at the center a disk (32 mm 2 ) 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 12 with permission.
by the cells and disappear with a characteristic uptake τ 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: (1) t e e e e e e t w e w w Let us define dimensionless physical quantities as capital letters and choose τ e as time unit, τ D e 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: t e e e t w w At the front of normal N , oriented towards Ω w , the jump conditions for both the concentration field and its normal gradient read: e w e w i 0 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: where K e is the mobility coefficient.  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][33][34] . In the laboratory frame, we transform Darcy's law into = −∇ V 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: where the chemotactic migration constant is introduced: Λ = λ 0 c i /D e , suppress the cell diffusion coefficient: D and the coefficient of proliferation: α 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), ρ 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: Δ (P + Λ C e ) − α 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 In grey, the cell sheet (Ω e ), in colours, the bath of water and nutrients (Ω w ). From left to right, a steady front with constant velocity, a wavy front at threshold of stability, then the spatio-dynamical pattern.
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 subtile equilibrium by imposing: If the size of the transition zone remains small, the mechanical equilibrium of the front fixes the boundary condition for the pressure: where 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 κ is chosen positive for convex interface and σ represents the dimensionless capillary parameter σ τ = S K D /( ) T e e e 3/2 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 V s and a convection term ∇ V C s must be added into the left-hand side of equations (2). Also a modification of V i into − V V i 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 V into + V V e s y and transforming ρ and C into ρ(X, Y − V s t, t) and C(X, Y − 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: ( 1) and C 0 = r e U(1 + N 0 ). The pressure results from proliferation and chemotaxis: i y e y ry e (0) 0 2 2 2 e which gives the interface velocity through an implicit relationship, function of the chemotactic forcing Λ and the proliferation rate α: 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 Λ (1 + N 0 ) and αP i . From equation (8), we recover the two limiting cases: A front driven only by chemotaxis without proliferation, having velocity: ) and a proliferative front with velocity: U = − αP i . Let us consider now the stability of such fronts. The stability analysis of such  (8)).
Scientific RepoRts | 6:33849 | DOI: 10.1038/srep33849 fronts is performed by introducing a sinusoidal perturbation induced by a wavy frontier of the cellular domain: ζ = Ut + εe ikx e Ωt (in analogy with the Fig. 2, middle panel) and expanding all fields according to: e e e P P p y e e ; ; At linear order in ε, we have: where the growth rate Ω reads: 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 χ 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 = α = 1 and σ = 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 17 ). 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Ω/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 Λ 0 . Fixing the front velocity and the parameters σ, N 0 , α, 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 1,2,12 . 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. 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 Λ 0 . 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 36,37 , it is less true for free boundary problems, except for the side-branching instability in dendritic growth [38][39][40] . 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 Λ 0 and wavenumber k 0 . Indeed, when Λ steps across the threshold Λ 0 , one can approximate locally the dispersion relation, equation (11), by: where = t k t 0 6 . How to select the nonlinearity in equation (13)? Above threshold, the first nonlinearity is − β 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 u/ eliminates the unknown constant β from the equation for  u and we recover equation (13). However a more fancy way consists in applying the Galilean invariance symmetry (1/2∂ u 2 /∂ x): x → x − vt and u → 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 |μ| of order ε 2 , from the dispersion relation, we deduce that the lengths scale as 1/ε while, from the right-hand-side of equation (13), the velocity u scales also as 1/ε. The left-hand-side gives us the time-scale as 1/ε 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 η 1 , indicating a noisy threshold for bifurcation and an additive noise η 2 arising from all other sources. Both sources of noise are chosen Gaussian in space, white in time and are not correlated so 〈 η i (x, t)η j (x′ , t′ )〉 = C j (x)δ ij δ(t − t′ )δ(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 η 1 and η 2 . It reads: is the differential operator, right-hand-side of equation (13). The correspondence between the biochemical noise and η 1 and η 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 , (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 μ = 0. The effect of the multiplicative noise can be approximated by using Novikov's theorem 37 leading to a new contribution proportional to  C S k ( ) 1 which modifies the threshold μ becoming μ + C 1 . It contributes to a possible jump into the unstable domain even if the parameter μ is negative. Finally the quadratic nonlinearities do not contribute, in opposition to the Swift-Hohenberg equation 41 . 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 μ values, the velocity u results from the competition between both modes at =  k 0 and at =  k 1, respectively, which are coupled by the nonlinear terms. This competition suggests the following decomposition ε ε = + u x t Re A X T e B X T ( , ) 2 [6 ( , ) ( , )] ix , B representing a small change in the front velocity, X and T being rescaled variables: X = εx and ε = T k t 2 0 6 . A classical multi-scale analysis 41,42 gives us the coupled amplitude equations for both modes A and B: The last terms of the first equation (15) represent the advection of the pattern, =  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 43,44 and such coupling has been called Hopf-Turing bifurcation 42  , 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 The nonlinear analysis proves that there exists a band of stationary spatial oscillations in the neighbourhood of =  k 1, with amplitude − q 1 4 2 and wavenumber µ + k q (1 ) 0 , where q h < q < 0. At the frontier, q = q h , the frequency of the time dependent oscillations becomes µ ω k p h 0 6 . Also the averaged front velocity oscillates with the same frequency µ ω k p h 0 6 . 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 , µ 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 ε − − ~s 410 5 1 , which identifies to . As shown in ref. 45, τ e varies between 45 minutes to two hours while D e = 0.3(μm) 2 /s. Taking the lowest estimation for τ e , we find p ~ 0.15 which gives a period T h = 400s for k 0 ~ 7.2. From Fig. 4, on right, such k 0 value indicates a velocity U ~ 10, giving a front velocity 10(D e τ e ) 1/2 ~ 0.1μm/s while the wavelength is about 25μm. These estimations are in good agreement with the experimental results 16 . A phase-diagram (Λ 0 , U) is presented in Fig. 5 for fixed values of N 0 = α = 1 and σ = 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 Ω in the dispersion relation 46 . Although we cannot quantify the level of chemical signaling, any diffusive front is a strong indication of auto-chemotaxis 26,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 e (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 η 1 and η 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 where, the noise correlation function, for a white noise in time but coloured in space η 2 is : , which is positive below the threshold of instability (μ being negative). The velocity correlation function, , where u * denotes the complex conjugate of u, is evaluated explicitly in the following. Indeed, from equation (16), we derive: The long time behaviour, t 0 → − ∞ , cancels the first term in the right-hand side of equation (17). The function ′  S k t t ( , , ) satisfies In order to evaluate the quantity 〈 η k (t)u −k (t′ )〉 , from equation (17), we derive: where θ(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: If , 0 and , then ( ) At long times, the structure function reads: Close to the threshold, given by μ = 0, for negative values of μ, the only root of γ  k ( ) is =  k 1, which increases the velocity structure function for this value. Therefore, at linear order, fingering with =  k 1 is selected by the noise below the hydrodynamic threshold. However, at threshold, the structure function diverges for =  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 Λ and the two boundary conditions in equation (3). This leads to three stochastic variables that we consider uncorrelated.