Collective chemotaxis and segregation of active bacterial colonies

Still recently, bacterial fluid suspensions have motivated a lot of works, both experimental and theoretical, with the objective to understand their collective dynamics from universal and simple rules. Since some species are active, most of these works concern the strong interactions that these bacteria exert on a forced flow leading to instabilities, chaos and turbulence. Here, we investigate the self-organization of expanding bacterial colonies under chemotaxis, proliferation and eventually active-reaction. We propose a simple model to understand and quantify the physical properties of these living organisms which either give cohesion or on the contrary dispersion to the colony. Taking into account the diffusion and capture of morphogens complicates the model since it induces a bacterial density gradient coupled to bacterial density fluctuations and dynamics. Nevertheless under some specific conditions, it is possible to investigate the pattern formation as a usual viscous fingering instability. This explains the similarity and differences of patterns according to the physical bacterial suspension properties and explain the factors which favor compactness or branching.

which asymptotically gives Eq.(1) whenr → ±∞. Let us consider now the inner expansion in the local coordinate system ( e r , e s ) shown in Fig.(1). In this system, along the interface, using the normal reducedr coordinate and the arclength s, the Laplace operator and the time derivative applied to c give the following relations: where κ is the local curvature of the front (positive when it is convex) and ∂ t r = −V b . In the boundary layer, Eq.(1) is then transformed into: Solving this equation for c order by order reads: c = c 0 + αc 1 + α 2 c 2 + ... and ρ = ρ 0 + αρ 1 + α 2 ρ 2 + ....
and we derive to leading order: where C b is the far-field value of the concentration, eventually a function of the arclength s f. Notice that M 0 is related to the concentration jump at the interface since which is the Dirichlet condition. Considering the next order, we get: which gives forr → ±∞ Eq. (7), which is the Neumann relation, corresponds to a jump in the normal gradient, N 0 playing the role of the latent heat in dendritic growth. Contrary to dendritic growth of alloys [3] where the jumps at the level of the Neumann Eq. (7), and the Dirichlet, Eq.(5), conditions are dictated by thermodynamics, at this stage M 0 and N 0 are postponed. However, the similarity with dendritic growth indicates the physical meaning of N 0 being the auto-chemotaxis rate of production by unit time and volume.
) and the density profile, represented in the inset is ρ = (1 − tanh(r))/2. Notice for Λ0 = 4, the strange behavior of the density function confirming that Λ0 = 4 is the limit of the stability. (1 − ρ)) and the density profile, represented in the inset is found numerically. The same code of colors applies to the inset. A singularity occurs here for Λ0 = 4.

DENSITY BOUNDARY LAYER AND THE PRESSURE EQUATION
In the same way, we consider a continuous version of the velocity equation: For passive bacteria, we represent the velocity field by the gradient of the pressure deduced from Darcy law (here for simplicity, we restrict on Newtonian hydrodynamics) and we add to the pressure a density dependent contribution, different in the bacteria (D b ) and the water (D w ) baths: which gives a more simple scalar equation for the pressure field coupled to the concentration field such that : where the mobility coefficient M = 1+βρ 1+β . According to the Einstein relation, passive bacteria solutions are more viscous than water and correspond to β < 0. Positive values of β indicate that the bacteria are active. Expanding p as c in power of α in the local coordinate frame, we get to zero order: This relation shows that it exists a pressure jump due to chemotaxis. For our mobility choice, the jump in pressure is : The next order reads: At infinity, forr → ±∞, ∂ρ 1 /∂r = 0 as ρ 0 (ρ 0 − 1). Then, the continuity of the normal velocity: v r = V b is checked in both domains as physically required. To conclude, our model for the boundary layer is consistent with the continuity of the normal velocities and with v r = V b , but the pressure field is discontinuous at zero order.

THE MASS FLUX EQUATION
Taking into account Eq.(8), the mass flux equation is reduced to To leading order, it reads: where Λ 0 = (ΛM 0 )/D. We require that ρ 0 (r) behaves like a kink with or without chemotaxis. Obviously it will depend on the function F s (ρ 0 ) which must cancel for ρ 0 = 0 and ρ 0 = 1. This equation can be solved exactly for arbitrary F s function by defining dρ/dr as a function of ρ transforming the nonlinear second order O.D.E, Eq.(14), into 2 linear O.D.E of first order and we get: Eq.(15) shows that the bacteria repartition in the boundary layer is entirely determined by the surface proliferation rate F s (ρ). This function scaled by γ s ∼ 1/α 2 must vary sharply at the front to provide a controlled size of the front.To be a kink solution requires that, in Eq.(15), the integral vanishes for ρ 0 = 1. Imposing that the kink exists with and without chemotaxis gives 1) A kink solution pre-exits without chemotaxis and F s (ρ) vanishes for ρ = 0 or 1. The bacteria do not change their proliferative strategy but the density ρ is modified. Assuming F s is fixed one can suspect a singularity for Λ 0 ≥ 4.
2) Due to chemotaxis, F s is modified. It means that chemotaxis affects the proliferation rate.
3) The chemotaxis destroys the front making it more and more diffuse as time goes on.
A standard function for F s giving a kink solution [7] and G(0) = 0. It will give a kink solution whatever the value of Λ 0 , provided that the integral in the right-hand-side of Eq.(15) remains always negative. Considering hypothesis 1, and making the traditional and simple choice as F s = −4W (1 − 2ρ), we get in absence of chemotaxis ρ = (1 − tanh(r))/2 which corresponds to the solution in usual phase-field models, derived from the Landau approach for phase-transition. With chemotaxis, it becomes difficult to perform a direct integration of Eq.(15) and a numerical analysis is required.
For hypothesis 2, the transformation of the proliferative rate into Fig.(2), where F s is represented) . The parameter Λ 0 is fixed by the proliferation rate and is not an unknown. It is not the only choice. For example, we can choose F s = −4W (1 − 2ρ)/(1 − Λ 0 W ) and in this case we get : All values of Λ 0 are possible, but these values so M 0 are fixed by the proliferation rate at the front separation (see F s and the associated profile functions in Fig.(2)) .
Let us consider the next order in α. We first introduce the linear second-order differential operator: Using Eq.(12) we derive, after simplification, L(ρ 1 ) = 0 accepts as a solution ρ 1 = dρ 0 /dr, which can be checked by simple differentiation of Eq.(14) and whose adjoint is A necessary condition for the convergence of the expansion in α is given by the adjoint theorem and reads: Integrating the first integral by part, noticing that: and introducing the following integrals: we get σ being the capillary parameter as defined usually in phase-field models and 2J 2 = 1 − Λ 0 /6. This solvability condition gives a relation between the jump in concentration gradient N 0 , the jump in concentration M 0 , the jump in pressure T (s) due to the proliferation rate represented by γ when the front keeps its kink shape. When the surface proliferation rate defined in Eq.(2) of the main manuscript allows such kink formation it is possible to treat Eq.(2) as a free-boundary problem provided that suitable boundary conditions are applied for the continuous field. In absence of chemotaxis Λ 0 = 0, the coupling between morphogens and bacteria disappear but a diffuse front may exist. In this case we have: and the pressure jump at the interface, equivalent to Laplace law, becomes: The sign of κ has been chosen positive here for convex interface in agreement with the Laplace law, the pressure inside being larger than outside according to this law, in absence of growth. However, there is also a constant jump of pressure due to the volumetric growth.

PARTICULAR CHOICES FOR THE PROLIFERATION RATE
According to hypothesis (1), we fix the proliferation rate F s to its value in absence of chemotaxis. Then the border proliferation rate is given by 4ρ 0 (1 − ρ 0 )(2ρ 0 − 1), and we can calculate only J 0 exactly. It does not depend on ρ 0 (difficult to evaluate), contrary to J 1 and σ. In absence of chemotaxis having ρ 0 = (1 − tanh(r))/2, all integrations can be done and we get: A linear expansion for weak Λ 0 value can be achieved for J 1 and σ. For hypothesis (2) with F s = 4w(2ρ 0 − 1)(1 − 2Λ 0 w), being w = ρ 0 (1 − ρ 0 ) and ρ 0 = (1 − tanh(r))/2, all integrals can be achieved exactly giving: Adding chemotaxis modifies the Laplace relation giving the so-called kinetic effects of dendritric growth and a modification of the surface tension: γv γs is a constant related to the growth. As expected, the sign of the kinetic effect λ is negative since J 1 is negative behaving as a stabilizing effect as the surface tension.

STABILITY OF THE BOUNDARY LAYER VERSUS TIME DEPENDENT PERTURBATIONS
We perform now a stability analysis of the boundary layer in the spirit of Kasyap and Koch [4,5] treatment for passive and active bacteria. We focus first on the passive case, the analysis being rather technical in our case since, first, the layer is curved and, second, has no fixed boundaries. We consider a dynamical perturbation of amplitude , being smaller than 1 and we solve the dynamics up to order α 2 . The relevant physical quantities are then: R(t, s) = R b + e Ωt e iks C = c (st) (r) + c(r)e Ωt e iks P = p (st) (r) + p(r)e Ωt e iks ρ = ρ (st) (r) + ρ(r)e Ωt e iks (27) Here, a steady state variable is mentioned by the label exponent (st) contrary to a dynamical one represented by a capital letter (except ρ). Dynamical quantities will be expended in modes up to linear order in . These modes are due to the waviness of the center line of the boundary layer (see Fig.(1)). Contrary to usual stability analysis and to the work by Kasyap and Koch [4], we only know the steady state ρ (st) , P (st) , C (st) up to order α which adds a degree of difficulty to the following treatment. From Eq.(13) we derive: Making the expansion in and taking the average according to the definition < f >= ∞ −∞ f (y)dy simplifies Eq.(28) by eliminating the first 2 terms of the left-hand-side combined with the first 3 terms of the right-hand side and givingρ = dρ (st) /dr, the analysis being similar to the determination of ρ 1 given in the previous section. So we get Considering modes having a wavelength smaller than the radius of curvature of the boundary eliminates the second term of the left-hand side of Eq.(31). Let us evaluateΠ s knowing that to leading order we have: v s = −M ∂p ∂s = −Λ 0 Dρ ∂ρ ∂s , we obtain: Using again the hypothesis of large k value compared to relative variation of ρ (st) we finally get: gives The stability of the boundary layer requires J 2 = 1 − Λ 0 /6 > 0 which limits the value of Λ 0 . For Λ 0 < 6, the boundary layer is stable versus perturbations of short wavelengths. However, it does not mean that the boundary layer is stable for arbitrary wavelength. Indeed it depends on the proliferation rate itself which can be dependent of the arclength and of the local curvature. In this case, an instability mode may exist at finite k.

STABILITY ANALYSIS OF THE BOUNDARY LAYER WITH ACTIVE BACTERIA
A continuous model of active bacteria flow has been proposed by Kasyap and Koch [4] who show that the random walk of these bacteria is perturbed by the chemotacting forcing imposed transversally. It results an anisotropic and active stress in the flow which differs according they are pushers like E.Coli or B. Subtilis or pullers. Assuming a constant chemotacting forcing, represented by a velocity U 0 which induces a variation of bacteria density in the film, they prove the existence of a bio-convection process, recovering the experimental results of Sokolov and Aranson [6]. For weak concentration of chemicals and for dilute bacterial solutions, it is possible to average the bacterial motion on time scale larger than the time separation between 2 tumbling events T ∼ 1s and on a length scale larger than V p T (V p being the individual bacterium velocity) and the Stokes equation (in their case) turns out to be modified by an active stress, average of the force-dipole interaction.
The active stress per bacterium is S = −C a /16µV p L 2 ζ(e r · e r − 1/2I), with ζ = αΛV p ∂ r c a coupling parameter with the chemotactic gradient, C a is a geometrical factor, positive for pullers, negative otherwise, estimated to be |C a | ∼ 0.57, L means the effective size of the bacteria (taking into account the bundle). Individual velocity of bacteria V p is larger than the flow velocity given by v r and V p ζ ∼ 6v r . Defining a new parameter δ as the ratio between the average active stress and the pressure at zero order in the boundary layer leads to δ = (3C a /8)µv r L 2 α < n 0 > /(12µv r αR b /b 2 ), with R b the radius of curvature of the pattern, b the thickness of the film or the Hele-Shaw cell, < n 0 > the average number of bacteria per unit volume. So we get .
One can estimate this parameter taking L ∼ 12µm, b ∼ 10 −3 m, R b ∼ 10 −2 m, α ∼ 10 −1 which gives : δ ∼ 0.510 −14 < n 0 >. In Sokolov et al. experiments < n 0 > is of order 10 16 with a film thickness of order the µ. So δ is a parameter of order 1 or even larger. Because of the anisotropic active stress, the Darcy law is transformed into : which does not affect our treatment if we define an effective pressure Q = p(r) − δ 2 ΛM 0 ρ 2 dρ dr , playing exactly the same role as p. Then we get for the tangential velocity: with a change of sign in front of δ because of the anisotropy of the active stress. It modifies Eq.(32) into: With our choice for the kink solution we get for : which finally gives for the kink solution stability, the following growth rate for perturbations: which increases the probability to get an instability of the boundary layer due to bio-convection in the case of pushers. The dispersion law, Eq.(40), have an effective diffusion coefficient which becomes negative for finite values of Λ 0 and δ. There is no prediction of the wavelength in this range of parameters. As shown in [7,8], such a diffusive process may induce a spinodal decomposition in a full domain with domains of density 1 and others of density 0.
Up to now we focus our attention on the boundary layer which separates a quasi-homogeneous domain of bacteria and the bath containing the chemo-attractant. We show that chemotaxis is responsible of a diffusive layer and prohibits sharp front of separation. If the proliferation rate F s at the level of the front satisfies some properties, the bacteria remain confined in a growing domain and the colony expands regularly. However this adjustment may fail in practice on long times. In addition, for active bacteria of type pushers and morphogen attractant, a wrinkling instability of the boundary layer itself may occur, mostly controlled by the bacteria density < n 0 >. In the next section, we will assume that it does not occur, the parameter Λ 0 being far from 6 and δ being small. This allows to transform the continuous set of equations into a free-boundary problem with the suitable boundary conditions for the morphogen concentration that we have found Eqs.(5,7).

STABILITY ANALYSIS OF THE FREE-BOUNDARY PROBLEM
The aim of this section is to study bulk instabilities, different from the one originated from the boundary layer (compare Fig.(4) with Fig.(1)). We consider the free-boundary problem in radial geometry, taking into account the correct boundary conditions derived previously. In radial geometry the zero-order is: The linear expansion requires the linear order of the morphogen concentration to calculate v 1 (r) as p 1 (r). Then, the continuity of the velocity and the discontinuity of the pressure at the perturbed interface: R = R b + e Ωt cos(mθ) and normal velocity at the interface V = V b + Ωe Ωt cos(mθ) will allow to determine completely the flow field. But first let us evaluate the pressure field to linear order for a p-Laplacian according to Eq. Ψ is an arbitrary function that we choose as Ψ = ψ(r)e Ωt sin mθ to expand Eq.(42) to linear order in . Then we get: with ∂ r P being given by the zero-order solution calculated previously: Elimination of ψ is trivial by solving the first equation and replacing ψ into the second. Then we get the second order O.D.E for p 1 : Although linear this equation cannot be solved exactly except for η = 0 which corresponds to the standard Darcy's equation. In this case we get where r m Cos(mθ) are the free modes which satisfy the Laplace equation. Another asymptotic solution can be found in the W.K.B limit for m 2 /(1 − η) having a large value. Then we get to leading order: with 2a(m, η) = (η + 4m 2 (1 − η) + η 2 )/(1 − η). Defining Ω as we can then separate the stabilizing effect giving a negative contribution to the growth rate Ω n /Ω d from the destabilizing effect with Ω p /Ω d . Of course it depends on the sign of Ω d which is controlled by auto-chemotaxis.
where as in the main manuscript, I j means the ratio between 2 successive Bessel functions evaluated at the front: I j = I j−1 (R b )/I j (R b ).
We consider now the case of an expanding drop or a moving epithelium sliding on a substrate. It means that our parameter β → −1 to cancel the viscosity, outside the living matter domain. Then, the dispersion relation simplifies and reads: