Chemotaxis migration and morphogenesis of living colonies

Development of forms in living organisms is complex and fascinating. Morphogenetic theories that investigate these shapes range from discrete to continuous models, from the variational elasticity to time-dependent fluid approach. Here a mixture model is chosen to describe the mass transport in a morphogenetic gradient: it gives a mathematical description of a mixture involving several constituents in mechanical interactions. This model, which is highly flexible can incorporate many biological processes but also complex interactions between cells as well as between cells and their environment. We use this model to derive a free-boundary problem easier to handle analytically. We solve it in the simplest geometry: an infinite linear front advancing with a constant velocity. In all the cases investigated here as the 3 D diffusion, the increase of mitotic activity at the border, nonlinear laws for the uptake of morphogens or for the mobility coefficient, a planar front exists above a critical threshold for the mobility coefficient but it becomes unstable just above the threshold at long wavelengths due to the existence of a Goldstone mode. This explains why sparsely bacteria exhibit dendritic patterns experimentally in opposition to other colonies such as biofilms and epithelia which are more compact. In the most unstable situation, where all the laws: diffusion, chemotaxis driving and chemoattractant uptake are linear, we show also that the system can recover a dynamic stability. A second threshold for the mobility exists which has a lower value as the ratio between diffusion coefficients decreases. Within the framework of this model where the biomass is treated mainly as a viscous and diffusive fluid, we show that the multiplicity of independent parameters in real biologic experimental set-up may explain varieties of observed patterns.


Introduction
Chemotaxis plays an important role in many biological processes like wound healing, cancer metastasis and fertilization. In colonies or tissues, chemotaxis leads to collective motions regulated by biological processes like cell-cell interactions, adhesion to the substrate, secretion of growth factor or inhibitors. At the origin of this process morphogens are long-range molecules that induce concentration cellular responses, activating different transcriptional pathways [1]. In embryogenesis, diffusion of these species is at the origin of the "French flag model" where cell fates are simply explained by their concentration gradients. Here, we aim to investigate the role of these chemoattractants in two-dimensional migration of bacteria or cells in an epithelium. Our scope is to give a unified picture of these two Contribution to the Topical Issue "Physical constraints of morphogenesis and evolution", edited by Vincent Fleury, Paul François and Marie Christine Ho Ba Tho. a e-mail: benamar@lps.ens.fr different experimental in vitro situations and to present a continuous formalism of the shape of the colony during migration.
Morphogens diffuse while up-taken at the surface of the cells by complex chemo-receptors [2]. Chemotatic pathways for bacteria allow motile bacteria to navigate according to higher concentration of attractants. Colonies exhibit various patterns going from rather compact and cohesive to dendritic (even fractal in extreme cases) structures. In a pioneering review [3] Ben Jacob and collaborators have underlined the close similarity between bacteria colonies (in a Hele-Shaw cell) and patterns in physics experiments like dendritic crystallization [4][5][6], viscous fingering [7,8] and electro-chemical deposition [9,10]. However, the bacterial pattern differs depending on the species, the nutrient conditions and also the substrate, exhibiting diverse morphologies including Eden-like growth, fractal and dense branching morphologies (see the review by Ben-Jacob and Levine [11]). Epithelia are more compact, cells spread as a thin layer over a flat surface in a Petri dish. Fronts are denser compared to bacteria [12] even if some fingering instabilities may also appear.
All these patterns are reminiscent of diffusive instabilities but the relation between solidification or viscous fingering experiments is not obvious a priori. Here we present a continuous model of chemotaxis migration in the framework of the mixture models [13,14] with two phase constituents. In principle, this model can incorporate an arbitrary number of constituents described by a density and a velocity field defined on an intermediate scale between the consituent and the experimental set-up size. In this work, one phase is water while the other phase is a mixing of water and cells: bacteria or epithelium cells. We assume the horizontal size of the Petri dish or the Hele-Shaw cell infinite as the domain of the colony and we examine the evolution of this colony under a chemotactic gradient. We limit ourselves to the case of chemoattractants (the cells swim up gradients of morphogens), not the case of chemorepellants [15]. The mixture model is suitable for numerical investigations but here we use it to transform the front dynamics into a free-boundary problem. The boundary of the colony evolves according to a diffusion equation with an uptake term while at the unknown border boundary conditions are applied. Indeed when the transition zone is sharp enough, it appears as a zone of discontinuities between two domains. Boundary conditions depend on the biological processes involved in the model such as an increase of mitosis at the front as observed in recent experiments [16].
The derived sharp interface model is then a freeboundary problem as in dendritic growth or viscous fingering formulation. It has the advantage to be simple enough to allow analytical treatments for the existence and stability of fronts. It can appear as a new variant of the Keller-Segel [17] model of chemotaxis, originally written to explain aggregation of amoebae and which has been widely studied in mathematical biology. It is known to generate finite time singularities showing its intrinsic instability [18][19][20]. Once transformed into a free-boundary problem, our model indicates that planar fronts remain unstable for chemoattractants due to the existence of a Goldstone mode and a behavior in k 2 for the growth rate of weak perturbations at long wavelengths. This is mainly true for weak chemotaxis and hence, weak velocities. Although capillary effects stabilize the front at short wavelengths, it is not enough to circumvent the intrinsic instability of the Keller-Segel [17] model. However when chemotaxis increases above a critical threshold for the mobility coefficient, surprisingly the front recovers the stability. We examine also other physical effects present in practice in experiments like the diffusion in 3 dimensions or nonlinearities which have been extensively introduced in mathematical biology.
Nonlinearities occurring in the chemoattractant diffusion itself, the uptake or the mobility law introduce complexity in the analytical predictions and we consider only the last two cases. For the uptake process of the chemoattractants by the cells, the Michaelis Menten [21] law is very commonly used in bacteria simulations [22]. This nonlinearity does not suppress the Goldstone mode but may change the behavior at long wavelengths. Surprisingly a much more stabilizing effect is induced by a nonlinearity in the chemotactic mobility as a function of the chemoattractant concentration. In this case we demonstrate that planar traveling-wave stable solutions exist. The mixture model can be extended to incorporate short-range or elastic interactions like it has been done for tumor growth with adhesion to a soft substrate [23]. Our work will open new possibilities for growth of epithelia once introduced more solid interactions in the growth process. It might be a way to explain the relative stability of epithelium growth in vitro [24,25].
The paper is organized as follows: Section 2 presents the continuous model for chemotaxis, including diffusion, interaction between cells, the morphogen uptake and mitosis. Section 3 gives the asymptotic analysis allowing the transformation of coupled partial differential equations into a free-boundary problem which is easier for analysis. Section 4 gives the results of the model: the existence of traveling wave solutions, the stability of these solutions according to a restricted number of dimensionless parameters. The end of sect. 4 concentrates on the possibility of suppressing the Goldstone mode by introducing various realistic effects like mitosis at the border and nonlinearities. Section 5 compares the results to experiments. Finally in conclusion we sum up the main theoretical results.

The continuous model for chemotaxis
We consider a two-dimensional cell population of density ρ s immersed in a morphogenetic gradient. This population can be bacteria in a Hele-Shaw cell, bacteria in a thin film put on a solid substrate or an epithelium grow- ing on a substrate. Initially located in a domain Ω, with a front separation δΩ (see figs. 1 and 2), the colony moves in the direction of the highest chemoattractant concentration. When they are also nutrients, the colony also grows by cell division. Both processes are considered but we will focus mostly on migration.

Description of the cellular population
We adopt here a continuous model for the cells immersed in a liquid bath (see figs. 1 and 2). For growing epithelia on a substrate, like in wound healing mimetic experiments [25], the fact that the tissue remains a monolayer ensures the 2D characteristic of the migration even if the experiment takes place in a 3D Petri dish. For bacteria, on the contrary, the 2D hallmark is obtained by confinement in a Hele-Shaw cell or when bacteria grow making a thin film. We will not consider this last possibility here although adaptation of our model can be done for biofilms. In both cases (epithelia or bacteria), the colony is assumed to live in a flat land saturated by nutrients. We assume that the cellular population is located in an infinite domain Ω represented by y < 0 in Cartesian coordinates. Our sample being a mixture of two species with the same mass density, the density of water, we call ρ l and ρ s the concentration in each phase with ρ l + ρ s = 1. Far away from δΩ, for y positive ρ l = 1 and ρ s = 0. For y → −∞ ρ l and ρ s reach constant values: ρ s = ρ ∞ and ρ l = 1 − ρ ∞ . For epithelia only, when the tissue is at confluence, ρ s = ρ ∞ = 1 and the layer thickness h 0 is then constant. At the front it decreases and this variation [16] can be taken into account in our model by an increase of ρ l . Each phase has locally a velocity (v l or v s ) and the average velocity is then given by V = ρ l v l + ρ s v s . We consider now the diffusion of the nutrients with concentrationc.

Morphogen diffusion in a Hele-Shaw cell and definition of the physical units
We begin by the bacteria (Hele-Shaw cell). The nutrients may play different roles and we will try to involve all of them in the same formalism. Being chemoattractants, they are the motor of the migration but they will be absorbed by the cells for proliferation or simply because of surface adhesion. Choosing the uptake time τ c as time unit and calling D l the morphogen diffusion coefficient in the liquid phase, we define a length scale by l 0 such that l 2 0 = D l τ c /ρ ∞ and a velocity unit given by The validity of the model requires l 0 to be much bigger than the cell size. All the equations presented hereafter are dimensionless. The chemoattractant concentration c =c/c ∞ satisfies a 2D diffusion equation given by ∇ being the horizontal dimensionless gradient ∇ = (∂ x , ∂ y ) and D(ρ s ) = (β − 1)(ρ s /ρ ∞ ) + 1. Clearly the morphogens are consumed only if members of the colony exist, i.e. for nonvanishing ρ s . c and ρ s are space and time dependent. Let us examine now the epithelium migration in a Petri dish according to figs. 2 and 3.

Morphogen diffusion in a Petri dish
The same equation in 3D for nutrient diffusion in a Petri dish with a localized absorption at z = 0 shows an additive term and we get with c 3d = c 3d (x, y, z, t) being now a function of all coordinates while ρ s remains a function of the horizontal position only. In this case the saturation conditionc 3d = c ∞ = 1 has to be applied for z = H, the position of the upper bath surface and we assume that H is a large quantity (much larger than our length unit and the thickness of the epithelium h 0 ). This hypothesis is checked in the experiment [25]. In this context, considering a distance d intermediate h < d < H, one can expand the density c as Above the distance d we assume that the saturation is provided by a continuous vertical flux, a valid assumption if H is much smaller than the lateral size of the Petri dish. This flux being amplified when the absorption process increases, we have ∂ z c 3d = α(1 − c 0 )d. Averaging the concentration gives then in the layer d eq. (4) is transformed into where a correction of the length unit has been introduced to take into account the transverse flux in 3D: (l 0 = l 0 / (1 + α)). The average concentration satisfies a 2D diffusion equation with an additive contribution due to the transverse flux. The depletion of the nutrients has to be decomposed into two contributions: one occurring at the border due to possible mitosis activity (not yet introduced in eq. (6)), the second present in eq. (6) acting on the scale of the colony. Both effects may have also different amplitudes. A more detailed description will be given hereafter. To solve one of these diffusion equations (1), (6) one needs to know the colony density ρ s and the velocity V . This is the aim of the next section.

Cellular density repartition and dynamics
We present here the mixture model mainly used in the context of tumor growth which constitutes a rather intuitive formalism when dynamics, growth, interactions between species and dissipation occur. We describe this model in the following and applied it to chemotaxis migration. This model [13,14] concerns a mixture of species in interaction which move according to mechanical laws. Being continuous it is valid for processes occurring at a scale larger than the size of the cells. For our purpose we limit ourselves to a two-phase model, mainly the cells and water, but the model can be extended easily to an arbitrary number of species. The species are then described by an average density ρ i and average velocity v i , the average being made at an intermediate scale between the cellular and the experimental sizes. These quantities are then space and time dependent. When the mixture contains much smaller elements like morphogens or nutrients, they are not considered as part of the mixture. However, they will modify either the dynamics or the species repartition. We represent chemotaxis by a simple linear relation between the cell velocity and the gradient relative to the concentration at the origin of a forcing velocity: More complex choices with nonlinearities can be found in the abundant literature on the subject [18,19] and will be discussed in sect. 4. The mass balance equation written in dimensionless units reads γ being the cell proliferation rate. For the liquid phase an equivalent equation with ρ l = 1 − ρ s reads Since the mass density of cells and water is very similar, we have γ l = −γ. The average velocity is given by with the incompressibility condition transformed into We assume a free energy F s for representing the attraction between cells where a penalty energy has been added for sharp density gradients. Dissipation in this system is due to the friction with the substrate or friction between phases. Assuming a viscous drag between phases, the variation of energy per unit times called Q is then The pressure p ensures the incompressibility of the mixture, the coefficient M refers to the friction between both phases, M l and M s represent either the fluid friction with a stiff static substrate or the internal viscous dissipation. Q has been called Rayleighian in [26] and appears as the equivalent of the Lagrangian for dissipative systems. As a consequence, our formalism is purely mechanical and does not introduce any thermodynamic concept as phase-field models for example [27,28]. The first term of eq. (12) exists only if both phases exist at an arbitrary point which explains the factor ρ l and ρ s , the same argument being applied to the second and third term of eq. (12). This energy rate Q must be an extremum which gives by elementary variation analysis with respect to both velocities with Introducing the incompressibility condition eq. (10) into eq. (14) allows to calculate the pressure p as usually done in hydrodynamics. Once p is introduced into eq. (13) we deduce v s as a function of ρ s and finally eq. (7) gives the density ρ s for a given interaction potential Ψ (ρ s ) and the growth rate γ. At this stage, the presentation of the model is complete and ready for numerical simulations. It is an alternative to phase-field models [27,28]. Another possibility is to derive the asymptotic sharp interface limit and transform this set of coupled partial differential equations into a free-boundary problem. Free-boundary problems allow to predict analytically the front morphology in simple initial geometries: planar and circular fronts, while continuous models are more easy to implement in numerical front simulations.

Morphogenesis free-boundary problem
We consider now the asymptotic limit of our system of eqs. (1), (6), (14) far away from the interface. It corresponds also to the equations we have to solve in the sharp interface limit.

Bulk equation
In Ω, the cell density ρ s has reached its constant value ρ ∞ . Cell divisions if they occur are localized at the interface border, and hence γ vanishes far away from this border and constant values are reached for ρ l and ρ s . Taking V from eq. (14) gives where C ∞ is the asymptotic limit of C with ρ s = ρ ∞ and ρ l = 1 − ρ ∞ . Note that for epithelia, C ∞ is equal to 1/M s . Hereafter we will use capital quantities for the asymptotic values or for the sharp interface model. Equation (16) means simply a Darcy flow both for a dilute bacteria bath but also for tissues when the hydrostatic part dominates the shear stresses. Far away from the front but in water now, ρ s = 0 and ρ l = 1. We have also a Darcy law with a constant permeability M l . Defining P = C ∞ · p the hydrodynamic flow equation in Ω becomes while in water it is When the ratio of permeability C ∞ M l is weak we can neglect the flow in the liquid and have P = 0. For the concentration field, we need simply to consider eq. (1) or eq. (6) with ρ s = 0 in the liquid and ρ s = ρ ∞ in the colony. We examine now the boundary conditions.

Boundary conditions
The correct boundary conditions to apply in the sharp interface limit come from the inner domain around ∂Ω when → 0. The asymptotic analysis is not difficult but a little technical. We give here the outline of the proof. First, as for phase-field models [27,28], we define a dimensionless order parameter φ which will vary between −1 and 1 across the tiny boundary layer, the interface itself being defined by φ = 0. So we get and ρ l = 1−ρ s = (1+φ)ρ ∞ /2. We define a local coordinate system on the line δΩ (see fig. 1), the unit vectors being the tangent and normal to this line. We solve eqs. (6), (14) in this frame of coordinates assuming the following scaling for φ, p and c as a function of R = r/ and S Here we make the assumption that the concentration field c(R, S) is continuous at leading order and may present a sharp variation only at first order.

Pressure jump at the interface
A sharp transition zone is obtained for typical energy interactions with two wells such as in the Cahn-Hilliard model [29,30]: . This choice has the advantage to give analytical results [31]. Combining eq. (10) with eq. (14), taking into account the scaling defined by eqs. (20), the leading contribution is then: which gives The next order gives with κ the local curvature, positive for the convex part of the interface and L 1 (φ 1 ) denotes Equation (23) is an equation for φ 1 once p 0 and c 0 are known. However one needs to ensure convergence of the expansion given by eq. (20), so one needs that φ 1 and all its derivatives cancel at both ±∞. Since ∂φ 0 /∂R is a solution of L 1 = 0, the theorem of the adjoint imposes Noticing that the integrand of the left-hand side is well behaved at both infinities, we integrate by part to extract the main contribution of the left-hand-side and we derive: where {p 0 } means the pressure jump between its value at the interface (φ = 0) and its value at R → −∞ (φ = −1), the same definition applies to {c 0 }. We transform eq. (26) into Equation (26) gives for our "mixture" the equivalent of a Laplace contribution with an additive correction when the concentration is discontinuous at the interface. Coming back to P we get the pressure jump as which is the boundary condition for the pressure.

Conditions on the normal velocities at the moving front
The moving front is given by δΩ defined by φ = 0. Examining eq. (7) written in the local moving frame of δΩ, we get where γ scales as 1/ so has a noticeable amplitude at the front but cancels in the liquid and also when the tissue is at confluence ρ s → ρ ∞ . γ depends also on the availability of the nutrients. We make the simplest linear choice and assume that at the front C 0 (S) the concentration of nutrients depends only on the arc-length (smooth variation of the concentration field). Taking into account eq. (19), the right-hand side can be transformed into Integrating between the interface (φ = 0) and −∞ (with φ = −1) we get the interfacial boundary condition for the cell velocities in the normal direction of the front

Boundary conditions for the chemoattractant concentration
Mitotic activity localized at the front induces a supplementary absorption of nutrients at this front. Our work hypothesis is that this additional depletion will not perturb too much the concentration field which will remain continuous even in the sharp interface limit. It means that the additive term we must add to eq. (1) or eq. (6) will have the same dependence as in eq. (29) but with a different constant Γ 2 . In the local coordinate frame both diffusion equations reduce to with Expanding the concentration field according to eq. (20), we get that the main order is a slowly varying function of the orthogonal coordinate r: C 0 (r, S). Next order in gives which gives for R → −∞ (or r → 0 − ) and for R → ∞ (or r → 0 + ) where capital letters mean functions defined out of the boundary layer, f (S) being an arbitrary function of the arc-length. So at zero order we get a continuous concentration field C 0 (r, S), with a discontinuity in the normal gradient, the discontinuity in the concentration appearing only at first order in . It will be discarded. Restricting to the zero order, we finally derive the boundary condition when a localized mitosis production occurs at the front and Finally we present a summary of the set of equations which constitutes the sharp-interface limit.

Summary for the sharp-interface model
We take the convention to use capital letters for physical quantities in the sharp-interface limit. We detail now these equations in the liquid and in the colony phase. Neglecting the fluid convection, the diffusion equation in the liquid is given by with δ = α/(ρ ∞ + α). In the colony we have Boundary conditions to apply at δΩ is the continuity of C and the discontinuity of the normal gradient N · ∇C C l = C s and N · ∇C l = βN · ∇C s + Γ 2 C s on δΩ, C l = 1 when y → ∞, with N the normal to the interface. The front velocity is given by the displacement of the cells We assume the viscosity of water negligible. Boundary conditions for hydrodynamics are then If we neglect mitosis, and focus on front migration, there exits 2 independent parameters which are Λ and σ in the 2D case, 3 in the 3D case with δ. The two first are related to physical and biological quantities that we recall here: α or δ are parameters much more difficult to evaluate. Indeed it is the result of viscous and chemical boundary layers. It changes the effective diffusion coefficient which becomes D l → D L /(1 + α), decreasing the length and the velocity units. Let us study the role of these different dimensionless parameters for the existence of planar solutions.

Steady planar front 4.1 Conditions for the existence of steady fronts
Examining more carefully the set of equations, one notices that our growth problem couples a diffusion equation with consumption and a Laplace equation typical of viscous fingering. Being a more complex front dynamics than dendritic or viscous fingering front, we examine the existence of steady solutions moving along the y-axis at constant velocity U and we analyze their stability. In the moving frame, going towards the highest concentration of the chemical signal, the diffusion equation setting C ∞ = 1 for y → +∞, the continuity of C and the discontinuity of the chemical flux represented by eq. (39) at the free surface gives with C 0 (0) = (r l + δr s )/(2Γ 2 + r l + r s ), A solution for the pressure P given by solving eq. (40) is P 0 (Y ) = −Λ(C 0 (Y )−C 0 (0)) and the front velocity is then For arbitrary values of the physical parameters, this relation gives an implicit equation for U and must be solved numerically. Only in the case where δ = 0, one can find an explicit solution with showing that contrary to viscous fingering or dendritic growth, the velocity of the steady front is fixed by the complete set of equations and boundary conditions. Moreover a necessary condition for advancing front is that the chemotatic effect must overcome the diffusion inside the cellular domain but also the depletion of nutrients due to mitosis represented by Γ 2 . Mitosis plays the role of a motor (proliferation represented by Γ 1 ) but is also a brake since it consumes nutrients. To study the 3D effect of diffusion, let us make an expansion for small values of δ when Γ 1 = Γ 2 = 0 for simplification. We get indicating a decrease of the eigenvalue U with a transverse flux, for U 0 not too small. Moreover one must keep in mind the decrease of the length scale of our problem l 0 with decreasing α and δ, so the decrease of the unit for velocities. The true physical velocities depend on both variations. For very weak velocity, much smaller than √ β and √ δ, one gets If δ = 0 we refer to eq. (45).

General perturbation analysis in 2D and 3D
The steady front can be observed experimentally only if it is stable. If it is unstable the dispersion relation gives us some insight on possible nonlinear patterns as reported in [4]. Let us study first the case without mitosis at the border. So we imagine a small perturbation of the front δΩ given by ζ = 2 e ikx e Ωt inducing a perturbation of order 2 for the pressure and concentration fields as follows: where P 0 (Y ), C 0 (Y ) are the solutions calculated above eq. (45). Due to the weakness of 2 , our system of equations, once linearized, can be solved and a calculation tedious but without difficulty gives the dispersion relation Ω as a function of the wave number k. First, the diffusion equation in the moving front gives for Y > 0 and while the pressure p − reads given automatically the dispersion relation if Γ 1 = Γ 2 = 0. Ω is a function of r + , r − and c − , themselves function of Ω. So eq. (52) is an implicit equation for Ω which must be solved numerically. However, we cannot predict the characteristic of the eigenvalue spectrum. In particular, the existence and unicity of solutions are not assured without numerical investigation.

Results for 2D diffusion
First, we present the numerical results obtained for Ω as a function of k in true 2D diffusion. According to eq. (4.2.1) and eq. (52), Ω is a function of k given in an implicit way. As a consequence, the number of Ω-eigenvalues is not obvious to guess a priori and we focus on real and positive values of Ω. Asymptotic analysis is doable for k → 0 and for k infinite. Near k = 0, one can show that a Goldstone mode Ω = 0 for k = 0 exits and the limit k close to zero is crucial for global stability of the front. For this limit we get which proves that the process is unstable for small velocities and recovers stability at "large" velocities given by Above this value for Λ the planar front is stable since Ω is a decreasing function of k. Below Λ c or for small values of U , whatever the value of β the existence of a domain of unstable modes is usually the signal of strong instability suggesting dendritic patterns [4][5][6]. For k → ∞, Ω has an imaginary part but a real dominant part given by −σk 3 . Intermediate values are reached numerically and solutions are displayed in figs. 4, 5, 6 where we vary one parameter among (Λ, β, σ), the others remaining fixed. The spectrum is always the same showing a dendritic dispersion relation: strong instability at long wavelengths, stabilization by surface tension. Surface tension defines a cut-off that means a length scale l c above which patterns are unstable. As an example, radial patterns can be obtained if the radius of the colony is smaller than the cutoff given by

Results for diffusion in 3 dimensions
When we take into account the 3D effects of diffusion via a transverse flux, there is an increase of the velocities for small U 0 values which is responsible of an increase of the instability behavior as shown in fig. 7. However stability is recovered for large values of δ that is 0.8, being 1 the limiting value. The fact that an experiment is done in a Petri dish does not suppress the diffusive instability as shown by fig. 7 which is compared to the purely 2D case (with δ = 0). The transverse flux may give an astonishing strength to the instability for a weak value of its amplitude δ. In fact due to the various independent parameters of our model, the variation is not always monotonous and cannot be simply predicted.
The dispersion diagrams established here suggest that in an infinite geometry destabilization of the front gives patterns which can have infinite amplitude such as in viscous fingering or crystal growth [4][5][6]: there exists a continuum of unstable wavelengths in the vicinity of k = 0. Dendritic patterns are observed for sparsely distributed bacteria solutions [3] but are more questionable for epithelia grown on substrate [16,25]. For epithelia, the ratio between the diffusion coefficient in the tissue and in the liquid β is probably weak and there exists a transverse flux represented by δ in eqs. (44), (45). However, the saturation of the solution via the transverse flux that maintains a source of nutrients everywhere is not always sufficient to suppress the 2D diffusive instability, imposed by the migration of the colony. One can argue that perhaps we are in the regime where the 2D process is stable having Λ > Λ c (54). Clearly the instability of the model imposed by the behavior at k = 0 is robust. In the next section, we examine other physical possibilities which may be at the origin of the relative stability of epithelia.

Compact colonies and the sharp interface model
In the previous example, we conclude that at long times, diffusion will give dendritic patterns, a prediction compatible with bacteria colonies but not with growth of epithelium. Epithelia growing on substrate are more compact, even when they exhibit instabilities. Contrary to physical intuition, these fully nonlinear instabilities are not suppressed by 3D diffusion as we have seen previously. In this paragraph, we want to investigate the modification to limit the strength of the instability. Our scope is to define the conditions for observation of these compact patterns with perhaps weak instabilities as observed in in vitro wound healing experiments. To do so the behavior of the dispersion relation Ω(k) near k = 0 is crucial, at large k, the capillary effect behaving as k 3 will saturate the instability and will allow a control of the size of the instabilities. But the small-k behavior is an indicator of a possible threshold for growth instability.
First of all we want to examine the possibility to have a stable front so to avoid a Goldstone mode in this front problem. To simplify, we restrict to 2D diffusion but we examine the possibility to have a nonlinear diffusion equation for the nutrients in the colony. Indeed nonlinearities in the diffusion equation may have 3 origins: nonlinearities of the diffusion coefficient as a function of nutrient concentration or of the concentration gradient, nonlinearities that can arise from the chemotactic model with a more complex mobility law but also from the process of absorption of nutrients. In biochemistry enzymatic transformations lead to the Michaelis-Menten law [21] very commonly used in studies on bacteria colonies [22]. Then the nutrient concentration in the cellular phase N (Y ) (to avoid confusion with the linear case, we change the notation from C s to N ) is solution of a complicated nonlinear equation, while in water we keep our ordinary linear law of diffusion. There is little chance to find the analytic solution of N (Y ). We simply assume that it exists and that the requirement of a finite value at N (Y ) at −∞ selects a unique solution. Rewriting more formally the dispersion relation Ω(k) (52) than before, we obtain then n is solution of The zero order prescribes the logarithmic derivative of the nonlinear solution at the front. At next order we derive If Ω and k vanish, r + = 2U . Then knowing that n = ν 0 N , we get ν 0 = −1 except for very specific values of our problem so irrelevant. This value of ν 0 means Ω = 0 and nonlinearities cannot suppress the Goldstone mode. However, more attention has to be paid on the behavior near k = 0.
If Ω decreases negatively with k 2 we get a possibility for observation of a quasi-steady front. If mitosis is included, the limit U small is singular but looking for finite values of U we get which gives the same answer as eq. (53), which corresponds to the linear case. Again the front is expected to be highly unstable at low velocities. The effect of nonlinearities which have been extensively studied in mathematics at the origin of all the variants on the Keller-Segel model [17] are perhaps more surprising. Analytical predictions are really difficult to get: for example for the Michaelis-Menten law [21], we can prove only the existence of the Goldstone mode. But we cannot do more without knowing the base solution N 0 , which results from the solution of a nonlinear ordinary differential equation (O.D.E). Recent works [32,33] are concerned by exact derivation of this kind of (O.D.E.) called Fisher-like equation but we are not aware of a solution for the Michaelis-Menten law. We do not consider the case of nonlinear diffusion but the case of nonlinear mobility for the chemoattractants is more easy to handle. Nonlinear mobility has been introduced because cells do not measure only the gradient of N but they can be sensitive to more complex functions. For the chemotaxis, the model is extended to V = ∇χ(N ) the only restriction being Λ(N ) = χ (N ) positive for chemoattractants. In this case, focussing directly on the behavior for k = 0 we get Λ 1 being the first derivative of Λ with respect to N . U is given by a nonlinear equation One can show that if Λ is a slightly increasing function of N so that Λ 1 /Λ 1, then we get so for k = 0, Ω is negative and a planar stable front can be observed. This calculation is explicit only if Ω is weak (see also fig. 8). However for arbitrary values the spectrum may be more complex as shown by fig. 9 which exhibits domains where Ω has an imaginary part. Considering fig. 8, for specific values of the parameters, varying Λ 1 , it is possible to observe planar fronts (purple curve in fig. 8) and also a stable, steady periodic pattern with a prescribed wavelength corresponding to the vicinity of the red curve in fig. 8. One can argue that the nonlinearity which is introduced to stabilize the front is "unphysical" since it has a tendency to increase the linear behavior and usually nonlinearities are responsible of saturation. However bacteria use chemotaxis as natural self-engineering adaptation [11], especially at low nutrient concentration and such behavior is typical of small nutrient concentration.

Discussion
Recent published experimental papers on migration and growth of an epithelium [16,12,34,35] do not mention the possibility of chemotaxis as a possible driving force. So we use the published data simply to check that it may be a possible scenario. Let us give an evaluation of the dimensionless parameters entering into the model. We refer to eq. (42) which makes the link between dimensionless parameters and physical quantities. Values of the parameter Λ, motor of our migration process both for the case of epithelium or bacteria are difficult to find. We discuss first an epithelium made of MDCK cells experiments [12,34]. Taking into account the fact that the average velocity of MDCK cells on solid substrate is 40 μm per hour [12], this value is consistent with our velocity unit estimated to be D l /(τ c (1 + α)). Indeed according to [1] the uptake time is of order half an hour, D l in water is 10 −10 m 2 /s (that we decrease to take into account α) so our effective D l ∼ 10 −11 m 2 /s which gives approximatively 0.7 μm/s so 260 μm per hour. The dimensionless velocity U = (Λ − β)/ √ Λ (eq. (46)) is then U ∼ 40/260 ∼ 0.15. β being the ratio between the diffusivity in the tissue and the water is about β = 1/100, then we deduce that Λ ∼ 0.04 giving a typical wavelength for instability of order 16σ in dimensionless units. Let us evaluate now the dimensionless capillary number. In [34], the tension at the border of a MDCK epithelium is estimated to beσ ∼ 0.07 nN μm −1 while the friction on a rigid substrate is about ζ = 1 nN s/(μm) 3 which allows to determine the permeability K = l 0 /ζ. K being given by √ Dτ c /ζ = 1.410 5 , SI finally gives σ ∼ 3.510 −3 . Knowing that the length unit is 140 μm we get the typical length scale for the instability about 8 μm which is the cell size for this particular experiment [12]. This is coherent with experimental observations but at the limit of validity of the model.
In a different experiment involving a migrating cell sheet of Madin-Darby canine kidney, Trepat et al. [35] also observe an undulation of the advancing front with a wavelength of order 100 μm. They measure the forces and conclude that there is a traction force perpendicular to the front, a result consistent with our pressure model P = −Λ(C s − C sf ), C sf being the front concentration. Estimating the pressure in physical units, we get P ∼ ΛD l ζ/l 0 which gives the correct order of magnitude that is 1000 Pa. Moreover the stress grows quite exponentially with the distance from the front also in agreement with our model.
For bacterial colonies, our model may give a new way to explain the quasi-exponential behavior of the mass transfer due to chemotaxy according to eq. (44). Indeed experimental results aiming to measureΛ are not obtained directly but based on a model for the concentration. When mitosis is subdominant, our prediction for the length decay of chemoattractant concentration is √ Λ/(Λ−1) √ D l τ c , our Λ being itself related to its physical value by eq. (42). It varies linearly with the bacterium and the chemoattractant densities. In other experimental geometry, the model for regular fronts can be adapted without difficulty [36].

Conclusion
We have established an analytical model for studying the stability properties of a class of coupled P.D.E. governing front evolution due to chemotaxis. Our formulation based on the mixture model where the interface is diffuse is suitable for numerical investigation of living colonies moving in a chemoattractant field. We consider only chemoattractants which are not a by-product of the colony itself and which are consumed by the cells of the colony. Our model allows to derive the correct boundary conditions when the front between the colony and the environment is abrupt. Then the migration of the colony is transformed into a free-boundary problem which can be addressed analytically in simple geometries. Here we have investigated planar fronts and their stability which gives some insights on the nonlinear behavior of the front. Chemotaxis is known to be essentially an unstable process which is confirmed here when the chemotaxis attraction is weak. Surprisingly, when it increases, the stability of diffusion occurs as for dendritic growth and the front recovers stability. Also surface tension stabilizes short wavelength perturbations fixing a capillary length, the limit of stability of circular patterns. We investigate also other possible stabilization processes. These effects are physically relevant as the saturation of the solution in chemoattractants or a strong but localized mitosis at the front.
However, if we are successful in explaining the dendritic patterns for sparsely distributed bacteria (having β ∼ 1) deriving a dispersion relation typical of dendritic growth [4], perhaps the relative compactness of epithelia as observed in [12,16] may involve more complex interactions. We examine also the possibilities of nonlinearities of the chemotaxis process concerning the law of attraction and the uptake. Our conclusion concerning the Michaelis-Menten law for the uptake, in the absence of known exact solutions for the planar front is that the Goldstone mode is maintained and probably this nonlinearity is not enough. More success concerns nonlinearity in the chemoattraction itself which does not affect the diffusion equation and the boundary conditions. Here we can show that typical mobility laws give stable fronts or weekly unstable fronts. The main conclusion however is that if we want to interpret epithelia or a priori biofilms [37], elastic effects must be introduced which is possible in the mixture model framework [23] as it has been done for tumor growth. Ultimately, the comprehension of the shape instability and the pattern formation in multiphase models is of utmost importance for a wide range of problems relating to tissue growth, remodeling and morphogenesis under chemotaxis.
This work was supported in part by "AAP Physique Cancer 2012". We would like to thank Pascal Silberzan for sending an unpublished figure of epithelium growing on a substrate. We are very indebted to Pasquale Ciarletta and Julien Dervaux for helpful discussions.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.