CONGESTION-DRIVEN DENDRITIC GROWTH

In order to observe growth phenomena in biology where dendritic shapes appear, we propose a simple model where a given population evolves feeded by a diffusing nutriment, but is subject to a density constraint. The particles (e.g., cells) of the population spontaneously stay passive at rest, and only move in order to satisfy the constraint ρ ≤ 1, by choosing the minimal correction velocity so as to prevent overcongestion. We treat this constraint by means of projections in the space of densities endowed with the Wasserstein distance W2, defined through optimal transport. This allows to provide an existence result and suggests some numerical computations, in the same spirit of what the authors did for crowd motion (but with extra difficulties, essentially due to the fact that the total mass may increase). The numerical simulations show, according to the values of the parameter and in particular of the diffusion coefficient of the nutriment, the formation of dendritic patterns in the space occupied by cells.


Introduction
Dendritic growth is very common in nature: solidified alloys, snowflake formation, corals, bacterial colonies, and some viscous fluid flows, present similar branching patterns. Those patterns reflect a wide variety of mechanisms. In most cases, branches result from an external supply of a diffusing substance. The simplest model for such a phenomenon is the following (Diffusion Limited Aggregation [29]). An initial particle (seed) is located at some point, a second particle is introduced at a large distance from the seed, and is assumed to perform a random walk. Whenever it contacts the seed before escaping to infinity, it is assumed to irreversibly stick to it. A third particle is then introduced, and so on. When a new particle approaches the cluster, it is more likely to hit a growing part than a part close to the seed, which tends to generate highly irregular patterns. In the context of cell growth, the cluster corresponds to a zone occupied by living cells that grow at a rate depending of the local concentration of some nutrient. The latter substance diffuses outside the cluster, so that, like in the DLA process, nutrients are likely to be eaten before reaching the "fjords". Thus, the tips receive more nutrients than underdeveloped zone, which may tend again to develop branches. Yet, the actual growth of a branching cluster necessitates to account for cell motion. Linear diffusion for the cells does not lead to dendritic growth, but more sophisticated models (Kitsunezaki [13] or Kawaski [12], see details below) reproduce branching patterns at different scales.
We aim here at investigating the possibility that the motion of cells, and thereby the emergence of branching patterns, may result from a very simple mechanism, namely congestion. We shall focus on the specify situation of Bacillus subtilis colonies, but the principles are quite general: we shall consider a population of entities that grow at a rate which depends on the local concentration of a nutrient.
The nutrient propagates by passive diffusion (linear Fick's Law), whereas the main entity does not. The latter is simply subject to a maximal density constraint. When the maximal density is attained, entities push against each other, which induces a motion that ensures that the density constraint is not violated. We shall simply assume here that the correction velocity field that advects the entities is the minimal one (in the sense that it minimizes the L 2 norm) among those that preserve the density constraint. Before detailing this approach, we first give an overview of the different models that have been proposed in the literature to describe the growth of such cells.
Bacterial colonies of Bacillus subtilis develop very different patterns, depending on the conditions they grow into. Many experiments have been conducted to determine the morphological changes arising when the culture parameters are varied (see for instance the work of Matsushita and Fujikawa [9,16,10] and Ben-Jacob et al. [4]). In these experiments, bacteria are inoculated in the center of agar plates of thickness small enough to ensure a two-dimensional evolution of the system. Only two parameters are varied: the concentration of nutrients, denoted by C n , and the hardness of the medium, controlled by the agar concentration C a . Other parameters such as temperature are kept constant.
Five distinct morphological types have been identified (see Figure 1). In region A, for low concentrations of nutrients and hard media, bacterial colonies are growing through DLA processes ( [29]). It appears that bacteria cannot move at these concentrations of agar, and that their evolution is only due to cell-division, thanks to the diffusion of the nutrients towards the colony (see Matsushita et al. [17]). As the nutrient concentration increases, colony branches become thicker and finally form compact patterns (region B). The growth is also faster as C n increases. For intermediate hardness of agar, and high levels of nutrients (region C), the evolution consists of alternative expanding and consolidating phases, forming concentric ring patterns. Decreasing again the softness of the medium (region D), bacteria are able to move on the agar surface. Their evolution is the conjunction of cell-division and diffusion, and the growth consists of a disk expanding rapidly. Finally, for intermediate values of C a and low levels of nutrients, bacterial colonies develop dense branches, the envelope of which is smooth and circular (region E).
Many models have been proposed to reproduce these morphological changes. Some are microscopic and describe the behavior of each bacteria (see for instance the communicating walkers models of Ben-Jacob et al. [5] or the cellular-automaton model described in [3]), but most of them are macroscopic and study the twodimensional density of the bacterial colony via a reaction-diffusion equation. A very simple model to account for the time evolution of the densities of bacteria and nutrients b and n would be Unfortunately, system (1) cannot reproduce the emergence of branching patterns. In [13], Kitsunezaki proposes a model of non-linear diffusion, where two different types of cells are considered: active ones (which move actively and perform celldivisions) of density b, and inactive ones, their density being denoted by a. The reaction-diffusion equation becomes Branching patterns can be obtained for suitable parameters. Kawaski et al. propose another non-linear diffusion model, where the diffusion coefficient is proportional to nb, pointing out that bacteria are active mostly in the regions where nutrients are abundant (see [12]). This nutrient-dependent diffusion ensures the creation of branches without the need of adding a death term, as it is for instance the case in the models proposed in [18,19,20] by Mimura et al. More precisely, their model include a death term in the bacterial population of the form The numerical simulations reproduce different morphological patterns. Moreover, the authors prove in [8,23] the existence and uniqueness of the solution of the reaction-diffusion system, and characterize the asymptotic behavior of the total bacterial population. More complicated models include specific features, such as the presence of some lubricating fluid produced by the cells (see [11]), the effects of chemotaxis described for example in [11,15], or different states of mobility for bacteria, as in [14]. We focus here on the case of hard media (high concentrations C a ), where bacteria cannot move on the agar surface. It has been pointed out in [17] that cell-divisions are in this case the only cause of growth. However, as explained before, the simple reaction-diffusion system cannot produce branching patterns, which must be a consequence of other features which are still to be added to this system. In the spirit of the models studied in the context of crowd movements (see [21,22]) and cell-migration (see [7]), we propose to take into account the congestion of the bacterial population throught the constraint b ≤ 1.
This constraint can be seen as a formal limit of the non-linear diffusion term when m tends to infinity. As a matter of fact, this diffusion term tends to penalize high bacterial densities, and can be seen as a form of soft congestion, opposite to imposing a maximal density constraint, which is usually called hard congestion (see [22] for more details on the different ways of dealing with congestion). Despite its very simple formulation, this model contains many difficulties, both theoretical and numerical. In particular, classical methods are inefficient to prove the wellposedness of this type of problems. In [21], we turned to the theory of optimal transport and gradient flows (see [1,2,27,28]) to prove an existence result. We detail in Section 2 the rigorous formulation of the model. We prove an existence result in section 3, adapting the results obtained in [21] to the case where the total population is growing. Numerical simulations are presented in section 4.

Description of a model driven by congestion
2.1. The congestion model. We propose in this work to take into account the congestion of the bacterial population through the constraint b ≤ 1 which has to be satisfied almost everywhere in the bounded domain Ω. Since the experimental conditions ensure a two-dimensional evolution for the bacteria, this maximal density constraint makes sense. Here the limit density has been fixed to 1 for simplicity, but it can be changed to any other value without loss of generality. The reaction-diffusion system describing our model writes (4) ∂ t b + ∇ · (bu) = nb in Ω, where u is the minimal correction velocity (in a mean square sense) ensuring that the density b does not exceed 1. This velocity has to counterbalance the growth of the bacterial population, given by the term nb. More precisely, we choose the velocity which minimizes the L 2 norm among the admissible velocities, which are roughly speaking the vector fields v such that ∇ · v ≥ n. More precisely, we define the set C b of admissible velocities by Here H 1 b denotes a set of functions acting supported in the saturated area [b = 1] : Since the system evolves in a closed medium (the agar plate), we impose no-flux conditions for the nutrients on the boundary (7) ∇n · ν = 0 on ∂Ω, where ν denotes the outward normal of Ω. Note that this model expresses the basic principles of bacterial growth. The nutrients diffuse through the medium, and are consumed by the bacteria. The spatial expansion of the bacterial colony is due to the correction velocity only. This correction velocity is non-zero whenever there is a sufficient supply of nutrient to induce a growth that activates the maximal density constraint.
Definition 2.1 (Weak solutions). We shall say that (b, n) is a weak solution to (4)-(7) for initial data (b 0 , n 0 ) whenever for all The set C b defined by (5) is closed and convex, therefore there is a unique u with minimal L 2 norm. We shall use in the proof the saddle-point characterization of this minimizer: and such that (complementary slackness condition) Proof. The second equation of (8) ensures that u ∈ C b . Now, from the first equation of (8), u minimizes the quadratic functional Thanks to the slackness condition, it holds that For any v ∈ C b , the second term of the right-hand side above is non-positive, thus ||u|| ≤ ||v||.
i.e. it holds Despite its simple formulation, this model arises several difficulties. In particular, the correction velocity has no a priori regularity (else than L 2 ), and depends on b through a highly non-linear equation. This prevents us from using classical methods both in the theoretical and numerical aspects. We encountered the same difficulties in [21] in the context of crowd motion. The evolution equation studied there was very similar to the present one, since the crowd had also to respect a maximal density constraint. In this work, the problem had been reformulated as a gradient flow in the space of Wasserstein. The Wasserstein distance is indeed well-fitted to consider displacements between measures of same total mass. It is defined (see for more details [27,28]) as the minimal transport cost from a measure µ to a measure ν. More precisely, if µ and ν are absolutely continuous measures, it writes where t # µ denotes the push-forward measure of µ by t. In this context, it was possible to obtain an existence result for the gradient flow problem, and prove that this implies a solution to the original system. We cannot apply straightforwardly these methods to our system. First, the total mass of the bacterial density grows as time runs by, and this prevents us to consider the Wasserstein distance between them. Moreover, since System (4) is a coupled reaction-diffusion system, we have to deal with an evolving density of nutrients. We present in the next section a splitting discrete scheme which makes it possible to overcome these two major difficulties.

2.2.
Description of an approximation scheme. In order to prove the existence of a solution, we consider an approximating scheme of system (4). At each time step, we first let the bacterial density grow thanks to the provision of nutrients, disregarding the congestion constraint. Then we take the nearest bacterial density among the densities which obey the constraint b ≤ 1. More precisely, for a time step τ > 0, starting from initial conditions (more precisely: on each interval [(k − 1)τ, kτ ] the density b k−1 τ is considered to be given, as well as the initial value of n τ at time t = (k − 1)τ , then n τ is computed as the solution of the above equation on [(k − 1)τ, kτ ], which makes it possible to update the value of b k τ and n τ (kτ, ·) for the next interval). Here, P Km denotes the projection for the Wasserstein distance (defined in (9)) on the set of densities of total mass m, which satisfies the congestion constraint: and m k is the total mass ofb k τ (and b k τ as well). Here M + (Ω) stands for the set of positive finite measures on Ω. The constraint b ≤ 1 means in particular that they should be absolutely continuous, and they are therefore confused with their density.
Indeed, it can be established (for details see [22], Proposition 2) that the distance between any given density of mass m k and the set K m attains its minimum at a unique point.
Note thatb k τ is the solution at time kτ of and that the equations on n τ and b k τ are no more coupled, since the equation on n τ only involves b k−1 τ .
Our goal is to prove that these discrete densities converge to a solution of System (4) as τ tends to 0. With this aim in view, we define two families of interpolated curves. The first one will be used to prove that the congestion constraint is satisfied, and is defined by b τ (t, .) = b k τ , t ∈](k − 1)τ, kτ ], whereas n τ is simply given as the solution of (11). We also define discrete velocities for the projection step as where t k τ is the optimal transport from b k τ tob k τ . These velocities are interpolated through v τ (t, .) = v k τ t ∈](k − 1)τ, kτ ], and we also define the quantities The second family of interpolated curves will ensure that the limit densities are solutions of the reaction-diffusion system. We denote (see Fig. 2) k τ (t) (notice that both sides of this equation are functions -densities of non-negative measures of mass m k , actually -of the variable x, but the dependence will be often omitted in the sequel for similar objects), where T k t andb k τ (t) are defined by (15) and r k τ is the optimal transport fromb k τ to b k τ (note that r k τ = (t k τ ) −1 ). The nutrient density is simply given by and we denote byṽ τ the following interpolated velocitỹ We will prove in section 3 that these curves are somehow solutions of the reactiondiffusion equation we consider at a discrete level.

Existence result
We prove in this section the existence of a solution for the model of bacterial growth including congestion. More precisely, we prove the following theorem.
Theorem 3.1. Let Ω be a convex bounded set of R d , b 0 and n 0 given initial admissible densities, and (b τ , n τ ) constructed thanks to the approximating scheme defined above in Section 2.2. Then there exists a family of densities (b(t), n(t)) t>0 , and a family of velocities for a.e. t. Moreover, (b, n, u) satisfies the equation: in Ω, and u satisfies where C b(t) is defined in (5).
To prove this result, the strategy will be the following (we briefly sketch it here since the following subsections will be quite technical): we will use the approximation scheme developed in the previous section and • we prove thatb τ ,ñ τ andṽ τ satisfy the continuity equation; • we prove that v τ is the minimal norm vector field that we consider in the model thanks to a saddle point argument with a pressure function p τ ; • we prove a priori bounds on b τ , E τ ,b τ ,Ẽ τ and stronger bounds (suitable H 1 bounds) on n τ and p τ ; • we use the a priori bounds to get the existence of converging subsequences and we prove that the limits of b τ andb τ , E τ andẼ τ , n τ b τ andñ τbτ are the same; • we show, by using the H 1 bound on n τ (which implies strong L 2 convergence), that the limit of n τ b τ is the product of the two limits; • we show that the minimality of v τ passes to the limit, by taking the limit in the saddle point formulation. The only difficulty is the nonlinear condition p τ (1 − b τ ) = 0, which passes also to the limit due to the stronger bounds on p τ , even if this requires some effort since we only have a bound in L 2 ([0, T ]; H 1 (Ω)).
3.1. Technical lemmas. We prove in this section several lemmas that will be needed in the following. We first focus on the densities b τ and n τ .
Proof. Property (i) is obvious thanks to the definition of b k τ at (10). To prove (ii), the maximum principle ensures that n τ is positive (since b k τ ∈ [0, 1]). We also prove that n τ ≤ N 0 by applying the maximum principle to N 0 − n τ .
The following lemma is the discrete equivalent of the saddle point formulation (8).
Proof. Let us recall that the velocity v k τ is given by where t k τ is the optimal transport from b k τ tob k τ . Moreover, b k τ is the projection ofb k τ on the set K m k of densities of same total mass, which satisfy the congestion constraint. This means that b k τ solves min One can write (following lemma 3.2 p.13 of [21]) the optimality conditions characterizing the optimal b. If one consider another admissible measure h ∈ K m k and defines b(ε) : Yet, this derivative can be computed in terms of the Kantorovitch potential ψ in the optimal transport from b tob k τ (see also [25]), thus getting We recall that ψ is a Lipschitz function satisfying t k τ = id − ∇ψ. This implies for all other densities h ∈ K m k , which means that b solves a linearized problem min ψ dh, h ∈ K m k . Yet, this solution is easy to compute, since one only needs to put as much density as he can on a level set {ψ < l}. More precisely, we must have It is easy to check that the function p : This gives the required result by taking p k τ := p/τ . We define again an interpolation of the quantities p k τ thanks to p τ (t) = p k τ t ∈](k − 1)τ, kτ ].
We switch now to some a priori estimates.
since t k τ is the optimal transport from b k τ tob k τ . Moreover, as n τ is bounded by a 1] a.e., which implies that there exists a universal constant C such that (ii) We have p τ = 0 when b τ is stricly below the maximal density, so thanks to the discrete velocity decomposition (lemma 3.3), and (i).
(iii) Thanks to Cauchy-Schwarz inequality, , and therefore, thanks to inequality (18), We also prove an priori estimate for the nutrient discrete density.
Proof. We will assume that every density b k τ is regular enough just in order to perform the computations; indeed, all the estimates that we obtain could be extended by approximating less regular functions, and only depend on L ∞ bound on b k τ (which are all smaller than 1) and on the H 1 norm of the initial datum. Multiplying the equation [, times n τ , we obtain that the L 2 -norm in space of n τ is decreasing as time runs by, and hence is uniformly bounded by ||n 0 || L 2 . If we multiply the same equation times ∂ t n τ , we obtain moreover This gives bounds both on the L 2 norm of ∂ t n τ and of ∇n τ : and finally implies We now focus on the densitiesb τ andñ τ . We first prove that they are solutions of the reaction-diffusion equation.
Lemma 3.6. The densitiesb τ andñ τ are solutions of The first term writes The second one satisfies We finally obtain which corresponds to the weak formulation -given in Definition 2.1 -of the above equation.
Again, we can prove some a priori estimates.

Convergence of the interpolated curves.
We can now prove that the curves defined at section 2.2 converge to a solution of our problem. Notice that, unless differently stated, all the convergences we use in the sequel will be weak-* convergence in the space of finite measures in the duality with continuous functions over [0, T ] × Ω. We will also refer to this convergence as "weak convergence" and write ⇀ without using any * and without precising any more the space of test functions. When extra bounds are available, this convergence will turn into stronger ones: in particular bothb τ and b τ will actually converge weakly-* in L ∞ (due to the fact that these densities are uniformly bounded by 1) and n τ weakly converges in H 1 (and hence strongly in L 2 ). Let us start our analysis of the limits of the interpolated curves. Thanks to the a priori estimations given in lemmas 3.4 and 3.7, we know that (b τ , E τ ) and (b τ ,Ẽ τ ) converge (up to subsequences) as τ tends to 0. We prove in the following that their limits are the same.
Proposition 1 (Limit of b τ andb τ ). The curves b τ andb τ converge to the same limit b as τ tends to 0.

Proof. For every
The first term writes On the one hand, we have, thanks to Lemma A.1 and to the estimation on W 2 (b k τ ,b k τ ) proved at Inequality (18), On the other hand, since n τ ∈ [0, N 0 ] a.e., The second term also tends to 0, since We can also prove that E τ andẼ τ converge to the same limit.
Proposition 2 (Limit of E τ andẼ τ ). The curves E τ andẼ τ converge to the same limit E as τ tends to 0. Proof.
The first term writes 0, thanks to Lemma 3.7 (i). The second term also tends to 0 since We finally prove that n τ b τ andñ τbτ also converge to the same limit.
Proposition 3 (Limit of n τ b τ andñ τbτ ). The curves n τ b τ andñ τbτ converge to the same limit as τ tends to 0. Proof.
Since n τ is uniformly bounded with respect to τ , we can prove that these terms tend to 0 in the same way as in proposition 1.
3.3. Properties of the limit. We finally prove that the limits (b, n) of the curves are solutions of our problem. First, the a priori estimates given at lemmas 3.5, 3.4, and 3.7 prove that there exist the limits (in the weak sense of measures) ofẼ τ , b τ andñ τbτ , and that this last limit is given by nb, where n is the limit of n τ in H 1 ([0, T ] × Ω)) and b the limit of b τ , the main difficulty is that the product of the limits should be equal to the limit of the product, but n τ converges weakly in H 1 and hence strongly in L 2 . We can then pass to the limit in the linear equations satisfied by the discrete curves, and obtain that (b, n) are weak solutions of We have now to prove that the limit E is of the form bu, where u is the minimal velocities among the admissible ones, i.e. it satisfies the saddle point formulation (8).
To do so, we will use the limit of E τ rather thanẼ τ . Indeed, let us notice that we have Lemma 3.4 gives a bound on p τ , hence we can suppose that we have p τ ⇀ p (up to subsequences and in the weak sense). We need to show that p ∈ H 1 b , which would allow us to write E = −∇p = bu with u = −∇p. Let us suppose for a while that we have done this, and conclude.
Indeed, the limit equation would give, for a.e. t 0 , every h > 0, and every q ∈ Moreover, since b(t 0 , x) = 1 wherever q does not vanish, and since b( Letting h tend to 0, we get and the reverse inequality can be obtained with h < 0. Thanks to the saddle point formulation (8), this allows to conclude that u is the desired vector field, the minimal one in C b . Let us point out that the property we obtained is stronger than expected : we have proved that the complementarity relation holds for every test function in H 1 b(t0,.) . The proof is complete once we establish the following lemma. b . Proof. We obviously have p ≥ 0, and we must prove that p = 0 on [b < 1]. Following the method adopted in [21], let us consider the mean functions The advantage of using p a,b τ is that it is a bounded sequence in H 1 , converging (weakly in H 1 , strongly in L 2 ) to p a,b . Since As in [21], for a.e. a, if we take in the first term first the limit as τ → 0 and then as b → a, we get convergence to Ω p(a, x)(1 − b(a, x)) dx.
We have therefore, uniformly in τ , Thanks to lemma A.1, the second term satisfies b a kt ka Ω This implies that The limit is therefore given by

Numerical simulations
The numerical study of Problem (4), which we recall is given by arises several difficulties. First of all, the correction velocity u is given as the solution of a minimization problem on the space of admissible velocities C b defined at equation (5), which depends on the evolving density b. Its estimation seems to be out of hand, and even if it could be calculated, it would surely lack regularity. A different approach is therefore needed. Following the idea exposed in [22], which we already used in the theoretical work, we propose to use a splitting scheme to estimate the bacterial density. More precisely, for a time step τ > 0, we define the densities b k thanks to (20) bk where P K m k is the projection for the Wasserstein distance on the space of densities of same total mass asb k which respect the congestion constraint b ≤ 1. Notice that the exponential formula in (10), which was suitable for the theoretical proof of convergence, is now replaced with a simple first order expansion. The density n k is given as an approximated solution at time kτ of equation The spatial discretization also requires special attention. This kind of problems, which involve unstable interfaces, are indeed very sensitive to anisotropy. Instead of using cartesian meshes, we therefore turned to randomly generated Voronoi meshes. An example of such a mesh is presented in Figure 3 on the left.
We chose to discretize equation (21) thanks to a semi-implicit finite volume method. On each cell of the mesh, we write which implies, approximating b k and n k with b k i and n k i on the cell I, ).
x i Here, x i denotes the center of cell I, σ ij the length of the edge between cells I and J, and a i the area of cell I (see Figure 3). Concerning the bacterial density, the prediction step of equation (20) is easily computed through the following schemê However, the correction step is not straightforward. We propose to turn to a stochastic method, which has proved to be efficient in the crowd motion context (see [22,24] for more details). We consider the cells where the density is beyond the congestion constraint. For each of these cells, we take the additional mass, and start a random walk which travels through the neighbor cells. As soon as this random walk encounters non-saturated cells, it discharges as much mass as it can, and this until all the additional mass has been redistributed. This method asymptotically approaches the Wasserstein projection, in a sense which is precised in [24], and gives good numerical results on several examples. Note that we perform a random walk on a non regular mesh. The transition probability from one cell to another must therefore depend on the geometry of the mesh. The natural coefficients appear to be the ones that correspond to the discrete Laplacian on the mesh, which write Note that the transition probabilities are consistent to the finite volume discretization of the Laplace operator proposed in [26].
Numerical tests have been conducted with different parameters values. We present in Figure 4 three results for a fixed initial bacterial density b 0 = 0.7 localized in a small area at the center of the domain, a fixed diffusion coefficient D n = 10 −6 , and three different initial nutrient densities, n 0 = 2, n 0 = 0.5 and n 0 = 0.3. On the first case, the evolution is fast and concentric. For lower nutrient densities, however, branches appear and the evolution is slower. This example shows that branching patterns can be obtained without the need of adding non-linear diffusion or death terms.
It is known that the fractal dimension of patterns obtained thanks to diffusionlimited aggregation is of order 1.71 (see [29]). In the present context, the created patterns exhibit some fractal character down to a cut off scale (see the end of this section). In particular we can estimate some sort of fractal dimension of our  Figure 4. We obtain a slope of 1.96 in the first case where n 0 = 2, which corresponds to a concentric evolution, and a slope of 1.7 in the case where n 0 = 0.3, which corresponds to a branching pattern. This fits correctly with the estimations of DLA patterns.
Let us finally underline that the branching patterns we obtain are not fractal at every scale. It seems that there exists a limit scale for the creation of instabilities. More precisely, let us consider a typical situation where the border of the bacterial colony has irregularities, and where nutrients are less concentrated in the holes than on the peaks (see Figure 6). In this case, the irregularities will become unstable if the characteristic time of the nutrients' diffusion is slower than the time of growth of the peaks. The diffusion characteristic time can be estimated by L 2 /D n , where L is the length of the peaks, whereas the growth time is given by L/U . This finally gives the following condition for apparition of instabilities where U is the propagation speed of the bacteria, which is of order 1 if the nutrients concentration is of order 1. Figure 7 gives the evolution of initial irregularities for two different values of D n . The initial nutrient distribution has been taken uniform with n 0 = 1 on the area free of bacteria. This numerical test confirms the apparition of instability only for small diffusion coefficients. We have also plotted on Figure 8 the evolution of the amplitude of the instabilities for different diffusion coefficients. It appears that the limit scale is indeed attained for diffusion coefficients of order the length of the irregularities (here 0.5). Finally, Figure 9 gives the bacterial densities obtained for two different diffusion coefficients D n = 10 −6 and D n = 10 −5 . It appears that the branches are thicker in the second case.

Appendix A. A technical lemma
We recall here a result, first proved in [21] (Lemma 3.4). We give its proof for the sake of completeness, since we used it several times.
Lemma A.1. Assume that µ and ν are absolutely continuous measures, whose densities are bounded by a same constant C. Then, for all function f ∈ H 1 (Ω), we have the following inequality: Proof. Let µ t be the constant speed geodesic between µ and ν, and w t the velocity field such that (µ, w t ) satisfies the continuity equation ∂ t µ t + ∇ · (w t µ t ) = 0, and ||w t || L 2 (µt) = W 2 (µ, ν). For all t, µ t is absolutely continuous, and its density is bounded by the same constant C a.e.. Therefore: Notice that this estimate makes a connection between the distance W 2 and the dual norm of H 1 , exactly as the distance W 1 is connected to the dual of W 1,∞ . On the other hand, it is only valid provided some bound on the densities is supposed, and it is easy to check that (A.1) could not hold without some bounds. In particular, should one be allowed to use atomic measures µ, ν, the existence of an estimate like (A.1) would imply the injection H 1 ֒→ C 0 , which is false in higher dimension.

Appendix B. Estimates in the Wasserstein projection
A useful estimate that we needed throughout the paper was the following: considerρ be a probability density on Ω ⊂ R d satisfyingρ(x) ≤ 1 + τ ≤ 2 a.e. and K = {ρ ∈ P(Ω), : ρ(x) ≤ 1 a.e.}; then there exists a universal constant such that This fact is quite easy to prove when Ω = R d . Indeed, consider the following transport The transported density t #ρ then belongs to K, since it has the same total mass asρ, and its density is given by Therefore, we have which gives min{W 2 (ρ, ρ) : ρ ∈ K} ≤ Cτ, for a constant C only depending on the second moment ofρ.
If on the one hand this last constant would become universal in the case of a bounded domain Ω, on the other hand the same proof could not work, since the dilatation t could fall out of Ω. This is why we develop a different proof in the following theorem.
Proof. We will select a particular density ρ through a minimization problem, prove that ρ ∈ K, and give an estimate on W 2 (ρ, ρ). We will suppose thatρ has some extra regularity properties, namely that it is Hölder continuous and bounded from below by a positive constant. We will remove this assumption later.
Let us solve the following minimization problem min F (ρ) + ε f (ρ(x))dx + 1 2(1 − c) −1 τ W 2 2 (ρ,ρ) , where ε is a small parameter that we will send to 0 later on, c < 1 is such that c|Ω| = 1, f is a given convex function of the form f (t) = 1 t + t p for p > d (so that, in particular, f (ρ) < ∞ only if ρ > 0 a.e.) and It is easy to see that a solution exists, by semicontinuity, and we can assume that ρ ∈ L p , ρ > 0 a.e. and φ ∈ H 1 .
We want now to prove ρ ≤ 1. To do that we first analyze its regularity. Notice that ρ ∈ L p , then φ ∈ W 2,p ⊂ C 1,α . Since ψ is Lipschitz continuous, we deduce that f ′ (ρ) also is Lipschitz, and hence it is bounded. Since f ′ (t) = −1/t 2 + pt p−1 explodes near zero, this shows that ρ must be bounded away from 0. Moreover, since (f ′ ) −1 is Lipschitz (just compute f ′′ and see that it is bounded from below by a positive constant), then ρ is also Lipschitz. Hence we can apply Caffarelli's regularity theory since both ρ andρ are bounded from above and below and Hölder continuous, and we get that ψ is C 2,α . Also φ is C 2,α since ρ is C 0,α (for the regularity on φ, we must apply standard regularity theory with Neumann boundary conditions, which could require some extra regularity on ∂Ω).
We can now look for the maximum point x 0 of ρ, which is also a maximum point for εf ′ (ρ), and hence a minimum point for φ + First we want to send ε → 0. What we have proven so far is that, for every ε > 0 and for fixedρ (Hölder continuous and bounded away from 0), there is a measure ρ = ρ ε , satisfying the inequality (25). Since ε f (ρ(x))dx → 0 (ρ is fixed and it has been chosen so that f (ρ(x))dx < +∞) and ε f (ρ ε (x))dx ≥ 0, if ρ ε ⇀ ρ (and it is always possible to assume convergence, up to subsequences), one gets where we used the semicontinuity of F and of W 2 . Now we want to get rid of the regularity assumption onρ: if we take an arbitraryρ ≤ 1 + τ and a sequence of densitiesρ k ⇀ρ that are Hölder continuous and bounded away from 0, weakly converging in L 2 and satisfying allρ k ≤ 1 + τ (the convergence could be taken in L ∞ , but we do not need it) to our original measure, we can infer F (ρ k ) → F (ρ) (actually, weak convergence in L ∞ , forρ k means H 2 weak convergence for the corresponding functionsφ k , and hence H 1 strong convergence, which gives |∇φ k | 2 → |∇φ| 2 ). Moreover, we can also assume that the corresponding densities ρ k converge to some ρ and, by semicontinuity, we still have where we used both the inequality that we proved in Lemma A.1 (see also Lemma 3.4 in [21] for details), which is valid since both measures are bounded by 2, and the fact that the H 1 norm of φ is uniformly bounded as soon as the L ∞ norm of ρ is bounded. Finally we get W 2 2 (ρ,ρ) ≤ Cτ W 2 (ρ,ρ), which implies the thesis.