Incompressible limit of a continuum model of tissue growth with segregation for two cell populations

This paper proposes a model for the growth two interacting populations of cells that do not mix. The dynamics is driven by pressure and cohesion forces on the one hand and proliferation on the other hand. Following earlier works on the single population case, we show that the model approximates a free boundary Hele Shaw type model that we characterise using both analytical and numerical arguments.


Introduction
During development, tissues and organs grow while generating diverse cell types.Thus, different cell populations co-exist during growth.For example, in a developing limb, prospective muscles, bone and epidermis become distinct during development and as a result grow at different rates.As they grow, these cells types contribute mass gain for the whole structure.However, since the different cell populations grow at different rates, stress arises and must be relieved to ensure that they contribute harmoniously to the final structure.How differential growth rates within a structure are accommodated is therefore an important question in developmental biology.To approach this question from a theoretical point of view, we consider a model whereby an idealised tissue is composed of one contiguous cell population located within a wider area occupied by another cell population.
Our approach relies on a new continuum model for two populations of cells which includes the following biological features.First, we impose a constraint that cells do not overlap which gives rise to maximal packing tissue density.This is ensured in the model by an appropriate pressure-density relationship which becomes singular at the packing density.Second, we model cell-cell contact inhibition by implementing cell motion in the direction opposite to the cell density gradient.These features are believed to play important roles in the development of mono-layered epithelial tissues such as pseudo-stratified epithelia.Third, the model incorporates features that are specific to two non mixing cell populations.The model favors segregation by penalizing the mixing of cells of different populations.Such segregation is observed in various tissues, such as developing tissues when there are two populations of cells which are genetically distinct, or cancerous tissues composed of proliferative and healthy cells.Our model is set in an optimal transport framework.We describe the cell populations by means of continuum densities and assume that they satisfy coupled gradient flow equations which tend to decrease an appropriate free energy.The free energy decay encompasses important properties of the biological system and the use of a gradient flow ensures that this decay is built into the model.The choice of the metric used to compute the free energy gradient is also a critical modelling choice.In this paper, we choose the Wasserstein metric as it is the natural choice to measure the distances between densities and many different types of partial differential equation (PDE) have been derived as gradient flows in this metric (examples include the porous medium equation [36], parabolic of convection-diffusion equations [33], the Fokker Planck equation [30], etc).
After describing the model, we investigate its incompressible limit.In this limit, the the possible values of the cell densities are reduced to either zero or their respective maximal value corresponding to the packing density.We show that this limit model is a free boundary Hele-Shaw model (HSM), which allows one to focus on the geometric evolution of the boundaries of the domain occupied by the two species.
Mathematical models have been widely used to study tissue development or tumour growth.Among these, we distinguish two ways of representing the cells.On the one hand, discrete models consider each cell as an individual entity whose position and other attributes evolve in time [16,23].This provides a high level of accuracy but also results in large computational costs.On the other hand, continuum models consider local averages such as the cell number density, as function of space and time, which evolve according to a suitable PDE [1,7].This description is appropriate when the number of cells is large as it dramatically reduces the computational cost.However it only gives access to the large scale features of the system as the small scale features are averaged out.As the goal of this paper is to study the evolution of the whole tissue, we have chosen to use a continuum model.Continuum models roughly fall into two categories.The first category comprises models which describe the dynamics of the cell density through convection and diffusion [9,14,40].In the second category, we find models which describe the motion of the geometric boundary between the tissue or the tumour and its environment through geometric evolution equations [15,20,27,35].The latter share similarities with Stefan's free boundary problem in solidification [29].These two types of models are related to one another through asymptotic limits.In particular, some tumour growth models of the first type have been related to Hele-Shaw free boundary models in [38,37,39].
In this paper, we propose a two cell population model as an extension of earlier models for single cell populations introduced in [10,38] in the context of tumour growth.In these works, the authors consider a single category of cells whose density is denoted by n(x, t) and depends on time t ≥ 0 and position x ∈ R d .The diffusion of the density is triggered by a mechanical pressure p = p(n) which is a given by a non linear function of the density n.Cell proliferation is modelled by a growth function G = G(p) depending of the pressure.The displacement of the cells occurs with a velocity v = v(x, t) related to the pressure gradient through Darcy's law.The model is written as follows, In [38,37,39] the pressure is expressed as Inserting (3) into (1), (2) leads to the porous medium equation which has been widely studied [43].This model can be expressed as the gradient flow for the Wasserstein metric of the following energy, where P is a primitive of p, i.e. ∂P ∂n = p.Indeed the functional derivative δE δn of E with respect to the density n acting on a density increment δn is given by where < ., .> is the L 2 duality bracket, so that we have δE δn = p.So we can rewrite (1) as the following equation The incompressible limit of this model corresponds to γ → +∞.It has been shown in [28] that this incompressible limit is a Hele-Shaw free boundary model which, classically, used to describe the pattern of tumor growth.In two dimensions, the classical Hele-Shaw problem models an incompressible viscous fluid squeezed between two parallel flat plates.
As more fluid is injected, the region occupied by the fluid expands.It has been shown that the incompressible limits of many PDEs converge towards Hele-Shaw type models [2,3,18,21,24].This incompressible limit of the corresponding Hele-Shaw model have been shown to be particularly relevant to tumor growth modelling [38,37].
In this paper, we present a new Gradient Flow Model (GFM) for two cells populations which is built upon the gradient flow framework presented above.We introduce a free energy E(n 1 , n 2 ) that depends on the cell densities of each cell population n 1 and n 2 .The free energy encompasses a term like (4) which depends on the total cell density n = n 1 +n 2 and models both cell contact inhibition and packing.In addition we introduce active repulsion between unlike cells in order to enforce the segregation property expressed as r = n 1 n 2 = 0 i.e. the two cell densities cannot be simultaneously non-zero.As this term induces repulsion forces and triggers instabilities, we also include regularising terms depending on the gradients ∇n 1 , ∇n 2 .We pursue two goals.The first one is the formal derivation of the incompressible limit model which takes the form of a two species Hele-Shaw Model (HSM).The second goal is to develop a numerical method for the GFM which enables us to illustrate the validity of the limit HSM.
The pressure law (3) previously used in the literature does not prevent cells from overlapping.Indeed, with this expression, the cell density can take a value greater than n = 1, where the value n = 1 is supposed to be the maximal allowed cell density, corresponding to complete packing.In this paper we rather use the expression (instead of (3)) With this expression, the pressure has a singularity at n = 1 which prevents the density to take values above n = 1.Similar pressure laws have been used in [28].The limit γ → ∞ is now replaced by → 0. Systems with multiple populations are studied in many different areas.In chemistry, reaction-diffusion systems are used to model reacting chemical substances [17].In population dynamics, these model are generalised into cross-diffusion systems in which the movement of one species can be induced by the gradient of the population of another species.In biology, Keller-Segel models [31] are used to model bacterial chemotaxis.Another classical example of cross-diffusion in biology is the Lotka-Volterra model [34], which describes the dynamics of a predator-prey system.These have been extended to nonlinear diffusion Lotka-Volterra systems to model cell populations [4,8].In the context of tumor growth, systems with different types of cells have been studied (such as healthy/tumor cells, proliferative/quiescent cells [6,5,42]).Among these models, some conserve the segregation property [4,8,22], which is one of the characteristic we focus on.On the other hand, some models impose the species non mixing property.This is the case of the Cahn-Hilliard equation, which describes the process of phase separation [19].In this model, each phase tends to regroup in one domain and can erase the other phase.In our model, the total number of cells of each species remains constant (in the absence of growth terms and up to possible boundary effects).Though, there are some similarities between our model and the Cahn-Hilliard equation.
The paper is divided into the following five sections.In Section 2, the main results are exposed: the two population model is described; the formal incompressible limit theorem is stated; numerical simulations are shown in support and a discussion is provided.Then, Section 3 contains the derivation of the model from the free energy.The formal proof of the convergence to the Hele Shaw free boundary problem is in Section 4. Finally, Section 5 is devoted to the description of the numerical scheme.A short conclusion is given in Section 6.

Introduction of the gradient flow model
In this paper, we consider two densities of cells denoted by n 1 (t, x) and n 2 (t, x) and the corresponding pressures p 1 (t, x) and p 2 (t, x) that depend on time t ≥ 0 and position x ∈ R d .We derive our model from the single cell model exposed in the introduction and define the free energy E(n 1 , n 2 ) expressed as follows: where α is a diffusion parameter and , m are parameters of the pressure laws.The first term corresponds to the pressure building up from the volume exclusion constraint and the second term is a repulsion pressure between the two different categories of cells.The last term represents cohesive energy penalising strong gradients of either cell densities.We assume that the two categories of cells have identical geometrical characteristics, so that the volume exclusion pressure resulting from either category of cells is similar and the total volume exclusion pressure is just a function of the total cell density.The expression of the pressure p is given by ( 5) and a parameter supposed to be small.The repulsion pressure is a novel aspect of the model and is defined by with q m = Q m .This expression imposes segregation when m is going to infinity.This can be seen thanks to the equality (1 + r)q m (r) = (1 Passing to the limit m → ∞, we obtain r ∞ q ∞ = 0. Since q ∞ > 1, this implies that r ∞ = 0, which expresses the segregation property (the cell densities cannot be simultaneously non-zero).
For the sake of simplicity, we omit the parameters , m, α in the notations of unknown function n 1 , n 2 , p 1 and p 2 .From the free energy (6), in Section 3, we derive the following system of equations, where G 1 and G 2 are growth functions depending of the pressures p 1 and p 2 , respectively.The particular case q m = 0 and α = 0 has been studied in [4,8,22] where it has been shown that there exists a segregation property between the species, provided that the initial conditions are segregated [6,5,12,32].The present paper treats a different case, as initially the two species may be mixed and the model drives them to a segregated state after some time, except for a thin interface depending on α and m.This is consistant with the biological observation that some mixing between the cell species occurs across the interface.The aim of this paper is to investigate the incompressible limit of this model ( 8)- (11), which consists in letting and α going to 0 and m going to infinity in the system.

Formal limit
In this section, we obtain convergence results of the model when , α → 0 and m → ∞.
First we list some assumptions.As for the growth function, we assume: These assumptions stem from biological considerations.As the pressure increases in the tissue, cell division occurs less frequently, until eventually the pressure reaches a critical value, which either stops the growth or starts to trigger cell death.As for the initial conditions, we assume that there exists 0 > 0 such that, for all ∈ (0, 0 ), Thanks to these assumptions we can establish a priori estimates on n 1 and n 2 .With a Stampaccchia method, we show the positivity of the densities.The L 1 bounds on n 1 and n 2 follow from this result.The initial L ∞ bounds ensure that the densities stay below the critical value 1.However, in order to pass to the limit, we need more a priori estimates.
In the following theorem, we only provide a formal limit as proof rigorous would require obtaining a priori estimates on the two densities, the pressure and on their derivatives in time and in space.The emergence of terms of the type ∇n∇r in the space derivatives of the densities makes it difficult to obtain required estimates since we are not able to control these terms.The same kind of problem is encountered for the time derivative.The main result is the stated following theorem.
satisfy assumptions (12) and (13).Suppose the limits of the densities n 1 , n 2 , of the pressure p and of q m as → 0, m → +∞ and α → 0 exist and are denoted by n ∞ 1 , n ∞ 2 , p ∞ and q ∞ .If the convergence is strong enough then these limits satisfy In addition, we have and the complementary relation Eq. ( 17) suggests to decompose the domain in two parts.We consider the domain Ω(t) = {x | p ∞ (x, t) > 0} and the complementary domain, where the pressure is equal to 0. Notice that Ω(t) ⊂ {x | n ∞ (x, t) = 1}.Moreover, the two domains coincide almost everywhere.Indeed if n ∞ = 1 and p ∞ = 0, the density will continue to grow and the total density will become greater than the maximum packing value.Thanks to the segregation property (18) we can decompose Ω(t) in two subdomain, Ω and Ω 2 (t) are disjoint and that their reunion is Ω(t).It is interesting to remark that with q ∞ = 1.
The equation for the pressure inside Ω(t) is deduced from the complementary relation (19): Since the limit pressure is continuous, the pressure is equal to zero on the boundary ∂Ω(t).
However the normal derivative of the pressure is non-zero and controls the movement of the domain boundary.Knowing that in the system ( 8) -( 11) the densities n 1 and n 2 diffuse with velocities ∇p 1 and ∇p 2 respectively and given that at the limit, according to ( 20) and ( 21), these two pressures are equal to p ∞ , we can deduce that the domain boundary and the interface are moving with a normal velocity, where υ is the outward normal vector.

Numerical validation 2.3.1 Analytical solution of the Hele-Shaw model
The HSM is characterised by the complementary relation (19) and the velocity of the boundary given by ( 22).In the one-dimensional case (1D), the solution of this problem can be computed.We will consider G(p) = g(p * − p), with p * the maximum pressure and g the growth rate, since it is one of the most common growth term in literature.In the case of one population, the complementary relation (19) in 1D can be rewritten as The solutions of this problem are of the form with constants A, B depending on the boundary conditions.We compute the exact solution in two different cases which are going to be used for the numerical validation of the model.

Example 1
We first consider the case where the two species have the same growth term with one species surrounded by the other one.The density n ini 2 is defined as the indicator function of [x ini Γ − ; x ini Γ + ] and the density n ini 1 is defined as the indicator function of where x ini Γ − , x ini Γ + represents the interfaces between the two populations and x ini Σ − , x ini Σ + represent the exterior boundaries of the total density.More specifically, where 1 S is the indicator function of the set S. Since the two populations have the same growth function, the complementary relation can be treated as that of a single population of cells (23) with the boundary conditions p(x Σ − ) = 0 and p(x Σ + ) = 0 at any time.A simple computation shows that where cosh and sinh stand for the hyperbolic cosine and sine.The velocities of the exterior boundaries and of the interfaces can be easily computed: ) , and ) .
Then x Σ ± and x Γ ± evolve respectively according to d dt x so the interface is moving more slowly than the exterior boundary.This means that the density n 1 is not only transported but also spreads.The density n 2 spreads and simultaneously pushes n 1 .
Example 2 We now consider two species having only one contact point, with different growth terms The initial densities are defined as indicator functions of [x ini Σ − ; x ini Γ ] and [x ini Γ ; x ini Σ + ] where x ini Σ − , x ini Σ + define the boundary of the total density and x ini Γ defines the interface between the two densities.More specifically, The pressure follows the complementary relation (19), which can be rewritten as with the additional conditions p(x Σ − ) = 0, p(x Σ + ) = 0, and p and p are continuous at x Γ .After some computations we find, and From this expression we can compute the speed of the exterior boundary and of the interface by computing the derivative of the pressure and evaluating it at x ini Σ − , x ini Σ + and x ini Γ .

Numerical validation
In this section we use the numerical scheme presented in Section 4 to illustrate that the solution of the GFM converges to the corresponding solution of the HSM described in the previous section 2.3.1.In order to facilitate the comparison with the analytical solutions of the HSM shown in the previous section, the simulations are performed in one dimension.To compare the two models we initialize them with the same initial configuration: the densities are taken as segregated indicator functions.We consider the two different cases with different growth functions and different initial conditions, corresponding to Examples 1 and 2 of the previous section.
The parameters are = 0.1, m = 10, α = 0.01.Another parameter which is introduced for the Relaxation Gradient Flow Model (RGFM) in Section 4 is given the value ν = 0.001.We plot on the same figure the numerical solutions of the GFM (solid line) and the HSM (dashed line).The species n 1 and n 2 are plotted in red and blue, respectively.
The first case defined as Example 1 in Section 2.3 is illustrated in Fig. 1.We use the growth function (24) with the parameters g = 1 and p * = 10.The boundary and interface are taken as x ini Σ ± = ±1.4 and x ini Γ ± = ±0.6.Then the initial density of the HSM are defined by the formula (25).As initial conditions of the GMF, we take We take n 1 = n 2 = 0.98 initially as it is close to the singular value 1.In this example the solution of the GFM is expected to be close to the particular solution of the HSM given in Example 1.
The second case is illustrated in Fig. 2. The growth functions are defined by ( 26) with the parameters g 1 = g 2 = 1, p * 1 = 20 and p * 2 = 10.As initial conditions of the GMF we take where x ini Σ ± = ±1.4 and x ini Γ = 0.The initial density of the HSM are defined with formula (27).In this example the solution of the GFM is expected to be close to the the particular solution of the HSM given in Example 2.
We can draw the same conclusions about the approximation of the HSM by the GFM in the two sets of simulations displayed in Figs. 1 and 2 (where panel (a) is for the initial configuration and panels (b), (c), (d) are for times t = 0.001, t = 0.01, t = 0.15 or t = 0.25 respectively).First of all, we notice that the approximation is excellent except for the pressure at small times, and at the interface between the two cell population.While the pressure of the HSM is smooth wherever the total density is positive, the pressure of the GFM exhibits sharp ditches at the interfaces between the two populations.Otherwise the pressures given by the two models are similar except at small times.In Fig. 1(b), the pressure of the GFM is larger than that of the HSM, whereas on Fig. 2(b) it is the opposite.However, we do not observe any impact of this difference on the densities.On Figs. 2(c) and 2(d), the two pressures are almost the same for the two models.Overall, the dynamics of the GFM and the HSM are quite similar.In Fig. 1, we observe that the blue (inner) species pushes the red (outer) species in order to be allowed to grow.In Fig. 2, the red species (leftmost) which has the biggest growth rate pushes the blue  28), density of the HSM from Eq. ( 25), growth functions from Eq. (24).species (rightmost) toward the right.Therefore by growing more rapidly, the red species exert a bigger pressure on the blue species than the blue species exerts on the red species.This pressure imbalance which is visible on Fig. 2, triggers the motion of the interface towards regions of lower pressure.

Discussion
The results of the simulations exposed in Figs. 1 and 2 show that the dynamics of the two models are almost identical after some time.However, at the beginning of the simulation we observe some differences between the pressure of the GFM and the HSM.We notice that initially in the GFM, the densities are taken as indicator functions, so at the start the pressure p of the HSM is also an indicator function whereas that of the HSM is continuous.It takes some time for the scheme to smooth it out, and then the two pressures coincide.In addition, the initial value of the GFM in these simulations is fixed at 98% of the packing density 1.The growth term produces an increase of the densities to their upper bound value.The densities equate then to n M = p −1 (p * ) = p * +p *   29), density of the HSM from Eq. ( 27), growth functions from Eq. ( 26).
(because the pressure stays below p * = max(p * 1 , p * 2 ), see also (13)).During this period the GFM and the HSM are not synchronized.After a short transient, we observe that the pressures of the two models have the same shape.
The major difference between the two models is at the interfaces between the two populations.As result of the fourth order term, the densities of the GFM are not perfectly segregated and overlap only in a narrow spatial interval.When the parameter α tends to 0, the system amplifies the effect of segregation.However in the numerical simulation, the parameter α = 0.01 allows some mixing of the two populations.This mixing is observed at the interface and creates some small disturbances on the total density (i.e. the total density is not constantly equal to n M ), which are then amplified in the pressure.Despite this, the speeds of the interfaces and of the boundaries seem to be close in the two models.These are quite remarkable results, considering that the parameters are not yet in the asymptotic ranges where the two models should be identical.
However, it is important to remark that the simulations have only been made in the case of initially segregated populations.This was necessary to initialize the free boundary model.In the case of mixed population, simulations of the GFM have been made and can be found in Section 6.They confirm the convergence of the GFM towards a free boundary model since we observe that the system self-organizes into separated domains containing the different populations.However we can't study the convergence towards the HSM since we do not know beforehand which initial condition for the HSM will correspond to the the converged GFM.

Derivation of the gradient flow model
Eqs. ( 8) and ( 9) can be derived from the gradient flow associated with the free energy (6).Indeed the functional derivatives δE δn 1 and δE δn 2 of E with respect to n 1 and n 2 , acting on a density increment δn 1 (x) and δn 2 (x) are given by Indeed, the computation for the first derivative is given by The expression of δE δn 2 follows from a similar computation.Therefore the gradient flow associated to the free energy (6) according to the Wasserstein metric can be written as follows [36], The metric used here is not the traditional distance on L 2 (R d ), but the Wasserstein distance.This distance is defined on the set of probability distributions on R d .It is used in a wide variety of physically meaningful equations like the porous medium equation or degenerate parabolic equations [33,30,36].We recover Eqs. ( 8) and ( 9) of the GFM by adding growth terms G 1 and G 2 to ( 32) and (33).
It is worth noting that the free energy decreases along the trajectories of the equation in the absence of growth terms.Indeed, using the Green formula and the Eqs.( 32) and (33), we have Therefore, when the growth rates are set to 0, the model evolves in a way that minimize the free energy.This free energy depends of the pressures, representing the levels of volume exclusion and the segregation.Therefore the model seeks to minimize the pressures, which means moving away from situations where the volume exclusion or the segregation constraints are violated.

Formal proof of Theorem 1
This section is dedicated to the formal proof of Theorem 1. First, we prove Eqs. ( 17) and ( 18) that lead to the definition of the evolving domains of the HSM.Secondly, we compute the equations for the densities n 1 , n 2 and n in the limit → 0, m → +∞, α → 0. Third, we prove Eq. ( 19), also named complementary relation, which gives the evolution of the pressure inside the domain for the HSM.Finally, we compute the speeds of the boundaries of the domains.Eqs. ( 17) and ( 18) follow from the pressure laws ( 5) and (7).Indeed, (1 − n)p = n, it follows in the limit → 0 that Working out the repulsion pressure law (7) leads to the formula Taking the limit m → +∞, we recover the segregation property, i.e. ( Since q m > m m−1 for all m, we deduce q ∞ > 1 and Eqs. ( 14) and ( 15) follow directly from taking the limit → 0, m → +∞, α → 0 in ( 8) and (9).To recover (16), it is first interesting to remark that the solution of ( 8)-( 11) satisfies where Indeed, by adding ( 8) and ( 9), we obtain Then given formula (5), and given formula (7), and inserting (38), ( 39) into (37) gives (36).Since (r + 1) m = m−1 m (1 + r)q m (r) and r ∞ = 0, passing to the limit as → 0, m → ∞, we obtain So, taking also the limit α → 0, we have: To recover the complementary relation (19), we need to compute an evolution equation for the pressure.To do so, first we pass to the limit α → 0 in (35) but keep the notations of n 1 , n 2 , p 1 , p 2 , p , q m for the limit.We multiply the resulting equation by p (n) and obtain, Multiplying the Eq. ( 41) by and taking into account that p (n) = 1 (p + ) 2 yields ).
We now take the limit → 0, and denoting the pressure p at the limit by p ∞ but keep denoting the limit of the densities and the repulsion pressure by n 1 , n 2 , p 1 , p 2 ,q m , we get Then passing to the limit m → ∞, and using the expression (40) for the limit of H, we recover the complementary relation (19).This concludes the formal proof of the theorem.
In order to recover the speed of the boundary of the Hele-Shaw model, we first focus on the speed of the exterior boundary.Thanks to Eq. ( 16), for all Therefore applying the Green formula twice, yields where Γ(t) = Ω 1 (t) ∩ Ω 2 (t) is the interface between the two domains Ω 1 (t) and Ω 2 (t).On ∂Ω(t), ∂Ω 1 (t), ∂Ω 2 (t) the unit normals υ are directed outwards to the respective domains.On Γ(t), the normal υ is directed from Ω 2 (t) towards Ω 1 (t).For x ∈ Γ, we denote From Eq. ( 19), the first integral is equal to 0. In addition, the pressure is equal to zero on the boundary ∂Ω(t) so that Moreover since n ∞ = 1 in the domain Ω(t), we have where V ∂Ω(t) is the normal speed of the boundary ∂Ω(t) in the outwards direction.Since Eqs. ( 42) and ( 43) are verified for all ϕ ∈ C ∞ c (R d ), we deduce that and at the interface Γ, Then, both p ∞ and its normal derivative are continuous at the interface Γ.To find the velocity at the interface the same method needs to be applied on either n ∞ 1 or n ∞ 2 and it leads to where V Γ(t) is the speed of the boundary Γ(t) measured in the direction of υ.This gives (22).
If we consider the case where the initial densities of the HSM are defined by with B R the ball of center 0 and radius R and Using that the pressure and its normal derivative are continuous through the interface Eq. ( 42), and using radial symmetry, we obtain where, by abuse of notation, we have denoted p ∞ (r, t) the constant value of p ∞ (., t) on the surface {|x| = r}.Similarly, in applying the same result on the species n ∞ 1 , we obtain

Numerical method
This section is devoted to the derivation of a numerical scheme used to perform the simulations of the GFM presented in Section 1.3.Since the equations for the densities are of fourth order, we first introduce a relaxation model, in which the fourth order terms are replaced by a coupled system of second order equations.Then, we present the numerical scheme together with some of its properties and CFL stability condition.The last part contains some simulations without initial segregation.

The scheme
We aim to numerically solve the RGFM given by ( 47)-( 50) with Neumann boundary conditions by a finite-volume method.The choice of the boundary condition is arbitrary because we are interested in the dynamics of the system whilst the densities have not reached the boundary.To facilitate the comprehension, we omit the indices ν.For the sake of simplicity, we only consider the 1D case on a finite interval Ω = (a, b).
We divide the computational domain into finite-volume cells C j = [x j−1/2 , x j+1/2 ] of uniform size ∆x with x j = j∆x, j ∈ {1, ..., M x }, and x j = and define the cell average of functions n 1 (t, x) and n 2 (t, x) on the cell C j by A semi-discrete finite-volume scheme is obtained by integrating system (47)-(50) over C j and is given by Here, c β j (t) ≈ c β (t, x j ), β ∈ {1, 2} and F β,j+1/2 are numerical fluxes approximating ) and defined by: where and (u β j+1/2 ) + = max(u β j+1/2 , 0) and (u β j+1/2 ) − = min(u β j+1/2 , 0) are its positive and negative part, respectively.It is easy to see that scheme (51)-( 53) is first order in space and if one uses an explicit Euler method for time evolution, then the scheme is also first order in time.Higher order approximations can also be obtained, see, e.g., [11].

CFL condition and properties of the scheme
In this section, we present some fundamental properties the scheme is endowed with.In particular, we prove that for all times t ≥ 0, scheme (51)-( 53) preserves the nonnegativity of the computed densities n1 and n2 and also ensures that the total density n stays below 1.The former property guarantees physically meaningful values of the computed densities while the second one is enforced to make sense of the pressure law for which 1 is a singular point.
We consider the semi-discrete finite volume scheme (51)-( 53) and assume that is evolved in time by the forward Euler method.We denote by nk 1 j , nk 2 j , p k 1 j and p k 2 j the computed densities and pressures obtained at time t k = k∆t, i.e., nk 1 j := n1 j (t k ), nk 2 j := n2 j (t k ), p k 1 j := p 1 j (t k ) and p k 2 j := p 2 j (t k ), and prove the following two propositions.
Remark 1. Eqs. ( 54) and (56) give the conditions to choose the time step in the numerical simulations.It should be however observed that in some computations we needed to reduce the time step in order to avoid oscillations the pressure develops when the density is close to the singular point n = 1.
Remark 2. Similar theorem can also be proven if the second-order upwind spatial scheme from [11] is used and the forward Euler method is replaced by a higher-order SSP ODE solver since a time step in such solvers can be written as a convex combination of several forward Euler steps, see, e.g., [25,26].

Numerical simulation of the GFM with initial mixing
Finally, we illustrate the performance of RGFM in a number of numerical examples when the initial populations are not segregated.We take as growth functions, The numerical parameter values are = 0.01, m = 10, α = 0.01, ν = 0.001, and the initial densities are given by n ini 1 (x) = 0.7e −5x 2 and n ini 2 (x) = 0.5e −5(x−0.5) 2 + 0.6e −5(x+1) 2 . (58) These initial conditions depicted in Fig. 3(a) have been chosen to show the evolution with non segregated and non symmetric initial conditions.Fig. 3 (b), (c), (d) represent the solution at time t = 0.01, t = 0.1 and t = 0.15, respectively.The red line represents the species n 1 while the blue line represents the species n 2 .As expected, the two species are growing and diffusing.Since the initial overlap is quite strong, segregation does not appear immediately.In Fig. 3(b), we see that the red species, which is growing faster than the blue one, also grows in the interval occupied by the blue species.When the species reach the maximum density in Fig. 3(c), interfaces between the two populations are created.We can observe that the red species is split into two groups.In Fig. 3(d), we also observe that the inner species is pushing the other one to have more space to grow.However we still observe some mixing.On the one hand, looking at the central red species, we observe that the density of the blue species is not exactly equal to 0. This is due to the low value of the parameter m which is the exponent of the repulsion pressure.As m becomes larger, the simulation fulfills the segregation property faster.This can be observed on Fig. 4, where we plot results obtained from the simulation with the same initial condition than previously given by Eqs.(57) and (58) at time t = 0.15 with parameters = 0.1, α = 0.01, ν = 0.001.Respectively in Fig. 4 (a), (b), (c) the densities are plotted for the cases m = 5, m = 10 and m = 50.We can observe in Fig. 4(c) that there is more segregation than in Figs.4(a) and (b).If we compare the red central domain, in Fig. 4(c) the domain is smaller than in Fig. 4(b) because the blue species has grown in the middle of it in order to fulfill the non-mixing property.Then in Fig. 4(c), we don't observe the small bump of the blue population found in Fig. 4(b).On the other hand, at the interface, the species mix in a small interval.This is due to the finiteness of the parameter α which multiplies the fourth order diffusion term and of ν which is the relaxation parameter.We observe these phenomena in Fig. 5, where we show the results obtained from simulations done with the same initial conditions (57), (58) than previously, at time t = 0.15 with parameters = 0.1, m = 10.Respectively in panels Fig. 5 (a), (b), (c) the densities are plotted for the cases α = 0.1, ν = 0.001; α = 0.01, ν = 0.001 and α = 0.001, ν = 0.0001 (we have observed that the parameter ν needs to be smaller than α).We can observe that as α becomes smaller, the width of the interface between the two populations gets smaller.We can also observe in Fig. 5(c) that on the right side, thin stripes with alternating populations appear.This is the kind of dynamics we observed in simulations performed with α = 0. Indeed without the fourth order term in the GFM, the numerical simulations do not have the same dynamics as the HSM.This has been verified with numerical schemes of a similar structure to the one exposed in this paper, as well as with schemes derived from the gradient flow in the Wasserstein metric structure of the system.It appears that without the fourth order term, a species surrounded by an other one will prefer to split into two small subdomains crossing through the other species in order to grow instead of pushing it, which does not correspond to what we find with the HSM dynamics in the incompressible limit.For these reasons, the fourth order term is essential to producing realistic dynamics.What this analysis shows is that the choice of the parameter α is critical.In particular, in order to recover the HSM, α has to be made smaller as → 0 and m → ∞.The values of the parameters , m, α must be carefully chosen to match a particular application.

Conclusion
In this paper, we have presented a gradient flow model for two populations which avoid mixing.In addition we have introduced a numerical scheme and used it to study the behaviour of the system with different parameters.This model is a generalisation of a single population case which has been studied in the literature previously.
This paper demonstrates by a combination of analytical and numerical arguments, that the gradient flow model converges to a Hele-Shaw free interface/boundary model, when an appropriate set of parameters is sent to zero (or infinity).The analytical proof is only formal and is supported by numerical simulations.In particular, we observe that the speed of the boundaries and interfaces of the gradient flow model for parameters taken in the asymptotic regime is the same as those computed by the Hele-Shaw model.This is verified with values of the parameters fairly far from the asymptotic regime which means that the convergence is quite fast.
Perspectives for this work are both on the analytical and numerical sides.On the analytical side, a rigorous proof of the convergence of the gradient flow model to the Hele-Shaw one seems within reach in the case α = 0 and q = 0. On the numerical side, simulations of two-dimensional cases will be developed and applied to the modelling of tissue growth and growth termination.
1 j + nk 2 j , to be below 1 for all k ∈ [0, N ∆t ] is