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

This paper investigates the incompressible limit of a system modelling the growth of two cells population. The model describes the dynamics of cell densities, driven by pressure exclusion and cell proliferation. It has been shown that solutions to this system of partial differential equations have the segregation property, meaning that two population initially segregated remain segregated. This work is devoted to the incompressible limit of such system towards a free boundary Hele Shaw type model for two cell populations.


Introduction
Diversity is key in biology. It appears at all kind of level from the human scale to the microscopic scale, with million of cells types; each scales impacting on the others. During development, the coexistence of different cell types following different rules impact on the growth of tissue and then on the global structures. In a more specific case, this can be observed in cancerous tissue with the invasion of tumour cells in an healthy tissue creating a abnormal growth. Furthermore, cancerous cells are not playing all the same roles. They can be proliferative or quiescent depending of their positions, ages, . . . To study the influence of these diverse cells on each others from a theoretical view, we introduce mathematical model for multiple populations. In this paper we are interesting in the global dynamics and interactions of the two populations, meaning that we focus specifically on continuous models.
In the already existing literature on macroscopic model, we distinguish two categories. The most common ones involved partial differential equations (PDE) in which cells are represented by densities. These models have been widely used to model growth of tissue [10,29], in particular for tumor growth [1,7,9,15]. Another way to model tissue growth is by considering free boundary models [16,17,20]. In these models the tissue is described by a domain and its growth and movement are driven by the motion of the boundary. The link between these two types of model has been been made via an incompressible limit in [21,22,24,26,27,28]. This link is interesting as both models have their advantages. On the one hand PDE relying models, also called mechanical models, are widely studied with many numerical and analytical tools. On the other hand free boundary models are closer to the biologic vision of the tissue and allow to study motion and dynamics of the tissue. This paper aims to extend the link between the mechanical and the free boundary models, in the case of multiple populations system.
In the specific case of multiple populations, several mathematical models have been already introduced. In particular in population dynamics, the famous Lotka-Volterra system [23] models the dynamics of a predator-prey system. This model has been extended to nonlinear diffusion Lotka-Volterra systems [3,4,5,8]. For the tumor growth modelling (see e.g. [13]), some models focus on mechanical property of tissues such as contact inhibition [6,2,19] and mutation [18]. They have been extended to multiple populations [18,30]. Solutions to these models may have some interesting spatial pattern known as segregation [5,11,25,30].
The two cell populations system under investigation in this paper is an extension on a simplest cell population model proposed in [10,27]. Let n(x, t) be the density of a single category of cell depending on the position x ∈ R d and the time t > 0, and let p(x, t) be the mechanical pressure of the system. The pressure is generated by the cell density and is defined via a pressure law p = P (n). This pressure exerted on cells induces a motion with a velocity field v = v(x, t) related to the pressure through the Darcy's law. The proliferation is modelled by a growth term G(p) which is pressure dependent. With this assumption, the mathematical model reads ∂ t n + ∇ · (nv) = nG(p), on R d × R + , v = −∇p, p = P (n).
In [22,24,26,27,28], the pressure law is given by P (n) = γ γ−1 n γ−1 which allow to recover the porous medium equation. However, in many tissues, cells may not overlap, implying that the maximal packing density should be bounded by 1. To take into account this non-overlapping constraint, the pressure law P (n) = ǫ n 1−n has been taken in [21]. This latter choice of pressure law has also been taken in the present paper. For this one population model, it has been showed in [21], that at the incompressible limit, ǫ → 0 (or γ → +∞ depending of the pressure expression), the model converges towards a Hele-Shaw type free boundary problem.
The previous model has the particularity to derive from the free energy as a gradient flow for the Wasserstein metric. Using this property we derive a model for two species of cells. Let us denote n 1 (x, t) and n 2 (x, t) the two cell densities depending of the position x ∈ R d and the time t > 0. We assume that the pressure depends on the total density n = n 1 + n 2 . As the pressure depends on a parameter ǫ, we introduce this dependancy in the notation. We define the free energy for the two cell populations by, Restricting to the one dimensional case, the system of equation deriving from this free energy is then defined by, with G 1 , G 2 the growth functions, and p ǫ the pressure. The system is complemented with the initial data n 1ǫ (t = 0) = n 1 ini ǫ ; n 2ǫ (t = 0) = n 2 ini ǫ . The existence of solution for system (1)-(4) has been proven in [6,2]. Moreover, under a set of assumptions which will be defined latter in this paper, solutions to this system satisfies the segregation property, i.e. if the two densities are initially segregated, they remain segregated : We define the segregation property by The aim of this paper is to study the incompressible limit ǫ → 0 for this system. We firstly remark that by adding (1) and (2), we recover the classical equation of the one species model on the total density, Multiplying by P ′ (n ǫ ) we find an equation for the pressure, Formally, passing at the limit ǫ → 0, we expect the relation, In addition, passing formally to the limit ǫ → 0 into (3), it appears clearly that (1 − n 0 )p 0 = 0. We consider the domain Ω(t) = {x, p 0 (x, t) > 0}. It is clear that on Ω(t), n 0 is equal to 1. Moreover, from the segregation property, we have n 1ǫ n 2ǫ = 0 when the two densities are initially segregated. Passing to the limit ǫ → 0 into this relation implies n 10 n 20 = 0. Then we may split Ω(t) in two disjoint sets Ω 1 (t) = {(x, t), n 10 (x, t) = 1} and Ω 2 (t) = {(x, t), n 20 (x, t) = 1}. Formally, it is not difficult to deduce from (6) that when ǫ → 0, we expect to have the relation Then we obtain a free boundary problem of Hele-Shaw type: On Ω 1 (t), we have n 10 = 1 and −∂ xx p 0 = G 1 (p 0 ), on Ω 2 (t), we have n 20 = 1 and −∂ xx p 0 = G 2 (p 0 ), whereas p 0 = 0 on R \Ω 1 (t) Ω 2 (t) .
The outline of the paper is the following. In Section 2 we expose the main results of this paper, which are the convergence of the continuous model (1)-(4) when ǫ → 0 to a Hele-Shaw free boundary model, and uniqueness for this limiting model. Section 3 is devoted to the proof of these main results. The proof on the convergence relies on some a priori estimate and compactness techniques. We use Hilbert duality method to establish uniqueness of solution to the limiting system. Finally in Section 4, we present some numerical simulations of the system (1)-(4) when ǫ is going to 0 and simulations of a specific application on tumor spheroid growth.

Main results
In this paper we aim to prove the convergence the incompressible limit ǫ → 0 of the two populations model with non overlapping constraint (1)-(4). We first introduce a list of assumptions on the growth terms and the initial conditions. For the growth, we consider the following set of assumptions: The set of assumptions on the growth rate is standard and similar to the one in [21]. The parameters P 1 M and P 2 M are called homeostatic pressures which represent the maximal pressure that the tissue can handle before starting dying. For the initial datas, we assume that there exists ǫ 0 > 0 such that, for all ǫ ∈ (0, ǫ 0 ), These initial conditions imply that n 1 ini ǫ and n 2 ini ǫ are uniformly bounded in W 1,1 (R). Notice also that the existence of ζ 0 being the interface between the two species imply that the two population are initially segregated.
We recognise the Hele-Shaw model. Noticing also that taking t = 0 into the relation (17), we recover the expected relation p ini 0 n ini 0 = p ini 0 . The proof of this convergence result is given in Section 3. It is straightforward to observe that adding (1) and (2) provides an equation on the total density similar to the one found in the one species case [21,27]. Then we use a similar strategy for the proof relying on a compatness method. However the presence of the two populations generate some technical difficulties. To overcome them, we use the segregation property. Notice that this paper is written in the specific case where the two species are separated by one interface, but could be generalised to many interfaces. Using the segregation of the species we are able to obtain a priori estimate on the densities, the pressure and their spatial derivatives. Compactness in time is deduced thanks to the Aubin-Lions theorem. The proof of convergence follows from these new estimates. However, the lack of estimates on the time derivative makes obtaining the complementary relation difficult, then we are not able to recover the usual relation but the one stated in (17) which may be seen as an integral in time of the usual one as explained in the above remark.
To complete the result on the asymptotic limit of the model, an uniqueness result for the Hele-Shaw free boundary model for two populations is provided in Proposition 1 in §3.4. The proof of this uniqueness result for the limiting problem is based on Hilbert's duality method.

Proof of the main results
This section is devoted to the proof of Theorem 1, whereas in Section 3.4 the uniqueness of the solution to the Hele Shaw system is established. We first establish some a priori estimates.

Nonnegativity principle
The following Lemma establishes the nonnegativity of the densities.
Proof. To show the nonnegativity we use the Stampaccchia method. We multiply (1) by 1 n 1ǫ <0 and denote |n| − = max(0, −n) for the negative part, we get With the above notation, it reads We integrate in space, using assumption (7), we deduce Then we integrate in time, With the initial condition n 1 ini ǫ > 0 we deduce n 1ǫ > 0. With the same method we can show that if n 2 ini ǫ > 0 we have n 2ǫ > 0.
Remark 3. We notice that the positivity gives a formal proof of the segregation of any solution of (1)-(4). Indeed, defining r ǫ = n 1ǫ n 2ǫ and multiplying (1) by n 2ǫ , (2) by n 1ǫ and adding, we obtain the following equation for r ǫ , Given that r ini = 0, we get that r = 0 at all time.

A priori estimates
To show the compactness result we need to find a priori estimate on the densities, pressure and their derivatives. We first compute the equation on the total density. As show earlier n 1ǫ and n 2ǫ are respectively weak solution of (1) and (2). By summing the two equations we deduce that n ǫ is a weak solution of (5). Notice that this equation can be rewritten as, with We establish the following a priori estimates Lemma 2. Let us assume that (7) and (8) hold. Let (n 1ǫ , n 2ǫ , p ǫ ) be a solution to (1)-(4). Then, for all T > 0, we have the uniform bounds in ǫ ∈ (0, ǫ 0 ), Moreover, we have that (n 1ǫ ) ǫ and (n 2ǫ ) ǫ are uniformly bounded in L ∞ ([0, T ], W 1,1 (R)) and Proof. Comparison principle.
The usual comparison principle is not true for this system of equation. However we are able to show some comparison between the total density and n M defined by where we use the monotonicity of G 1 and G 2 from assumption (7). Notice that, since the function H is nondecreasing, the sign of n ǫ − m ǫ is the same as the sign of H(n ǫ ) − H(m ǫ ). Moreover, so for y = H(n ǫ ) − H(n M ) and f (y) = y + the positive part, the so-called Kato inequality reads ∂ xx f (y) ≥ f ′ (y)∂ xx y. Thus multiplying the latter equation by 1 nǫ−n M >0 and given (9) we obtain Since the function P is increasing and G 1 and G 2 are decreasing (see (7)), we deduce that the last term is nonpositive. Then, integrating on R, we deduce Then, integrating in time, we deduce With the maximum principle, we conclude that n ǫ ≤ n M . We deduce easily with the non-negativity principle (1) that 0 ≤ p ǫ ≤ P M , 0 ≤ n 1ǫ ≤ n M and 0 ≤ n 2ǫ ≤ n M . L 1 bounds of n ǫ , n 1ǫ , n 2ǫ and p ǫ . Integrating (18) on R and using the nonnegativity of the densities from Lemma 1, we deduce Integrating in time, we deduce Since n 1ǫ ≥ 0 and n 2ǫ ≥ 0, we deduce the uniform bounds on n 1ǫ L 1 (R) and on n 2ǫ L 1 (R) . From the relation (3), we deduce p ǫ = n ǫ (ǫ + p ǫ ). Moreover, the bound p ǫ ≤ P M := max(P 1 M , P 2 M ) implies L 1 estimates on the x derivatives.
Recalling (9), we can refomulate (5) by The space derivative of this growth function is given by, We derive (19) with respect to x, We multiply by sign(∂ x n ǫ ) = sign(∂ x p ǫ ) and use the Kato inequality, We integrate in space on R. Using the fact that max [0, Using Gronwall's lemma and the uniform bound on n ǫ and G 1 and G 2 (see (7)), we deduce that This conclude the proof for the estimate on ∂ x n ǫ . Then, We split the latter integral in two: either n ǫ ≤ 1/2 and then Then, we integrate in time and we deduce using (20) Hence we have an uniform bound on ∂ x p ǫ in L 1 (Q T ). To recover the estimate on ∂ x n 1ǫ and ∂ x n 2ǫ we derive in time (9), and This concludes the proof.
Proof. From the equation on p ǫ , we deduce Let us definep Therefore, for all t ∈ [0, θ], meaning thatp is a supersolution of the pressure equation (21). Let us defineñ = N (p) =p ǫ+p . Multiplying (22) Moreover, from (1) and (2) we have In other words,ñ is a supersolution of the density equation and by the comparison principle we have n ini ǫ ≤ñ ini ⇒ n ǫ ≤ñ. Thus, given a C large enough so that p ini ǫ (x) ≤p ini (x) we have p ǫ ≤p for all t ∈ [0, θ] and p ǫ (t) is compactly supported in [−R θ (t), R θ (t)]. Moreover p ǫ is uniformly bounded in L ∞ (R), so the process can be iterated to reach the time T in a finite number of step.

3.1.4
L 2 estimate for ∂ x p Lemma 4 (L 2 estimate for ∂ x p). Let us assume that (7) and (8) hold. Let (n 1ǫ , n 2ǫ , p ǫ ) be a solution to (1)- (4). Then, for all T > 0 we have a uniform bound on ∂ x p ǫ in L 2 (Q T ).
Proof. For a given function ψ we have, multiplying (4) by ψ(n ǫ ), Integrating on R, we have where Ψ is an antiderivative of ψ. We choose ψ(n) = ǫ(ln(n) − ln(1 − n) + 1 1−n ) so that n ǫ ψ ′ (n ǫ ) = p ′ (n ǫ ). Inserting the expression of ψ, we get After integrating in time and using the expression of the pressure (3), we have Then, to prove that ∂ x p ǫ ∈ L 2 (Q T ), we are left to find a uniform bound on R ǫn ǫ | ln( pǫ ǫ )|dx. Using the expression of p ǫ in (3), we have Since n ǫ is bounded in L 1 , the second term of the right hand side is uniformly bounded with respect to ǫ. Moreover given that 0 ≤ p ǫ ≤ P M and x → x| ln x| is uniformly bounded on [0, P M ], we get As the Lemma 3 provides a uniform control on the support of p ǫ , we conclude the proof.
3.2 Proof of theorem 1

Convergence
In the last paragraph we have found a priori estimates for the densities and their space derivatives. However we still need to obtain estimates on the time derivative. To overcome this lack of estimates, we are going to use the Aubin Lions theorem [31]. According to Lemma 4, n 1ǫ ∂ x p ǫ and n 2ǫ ∂ x p ǫ are in L 2 (Q T ). Moreover thanks to Lemma 2, we have that n 1ǫ G 1 (p ǫ ) and n 2ǫ G 2 (p ǫ ) are uniformly bounded in L ∞ ([0, T ]; L 1 ∩ L ∞ (R)), so ∂ t n 1ǫ and ∂ t n 2ǫ are uniformly bounded in L 2 ([0; T ], W −1,2 loc (R)). We also have n 1ǫ and n 2ǫ bounded in L 1 ([0; T ], W 1,1 (R)). Since we are working in one dimension The Aubin Lions theorem implies that {u ∈ L 1 ([0; T ], W 1,1 loc (R));u ∈ L 2 ([0; T ], W −1,2 loc (R))} is compactly embedded in L 1 loc ([0; T ], L 1 (R)). So we can extract strongly converging subsequences n 1ǫ and n 2ǫ in L 1 loc (Q T ). In order to extend this local convergence to a global convergence in L 1 (Q T ) we need to prove that we can control the mass in an initial strip and in the tail. Indeed, let ǫ, ǫ ′ > 0, R > 0, Since we have strong convergence of n ǫ in L 1 loc (Q T ), Then we have to control the two other terms in the right hand side. The control of the initial strip comes from the L 1 estimate of n, where the notation C stand for a generic nonnegative constant. Moreover, using equation (18), Then, integrating on [0, T ], we get By assumption (8), since the initial data is uniformly compactly supported, we deduce that the right hand side tends to 0 as R goes to +∞ and ǫ goes to 0. Thanks to the segregation property the control of the mass in an initial strip and in the tail is verified for n 1ǫ and n 2ǫ . Then (n 1ǫ ) ǫ and (n 2ǫ ) ǫ are Cauchy sequences in L 1 (Q T ). It implies their convergence in L 1 (Q T ). The convergence of the pressure follows from the same kind of computation. The only difference is for the control of the tail which is directly given by the estimate

Limit model
From the above results, up to extraction of subsequences, (n 1ǫ ) ǫ , (n 2ǫ ) ǫ , and (p ǫ ) ǫ converge strongly in L 1 (Q T ) towards some limits denoted n 10 , n 20 , and p 0 , respectively. Moreover, (∂ x p ǫ ) ǫ converges weakly in L 2 (Q T ) towards ∂ x p 0 thanks to Lemma 4. Passing to the limit in the uniform estimate of Lemma (2) gives (11) and n 10 , n 20 , n 0 , p 0 belongs to BV (Q T ).
Thus, the term in the Laplacian converges strongly to p 0 . Then, thanks to the strong convergence of n ǫ and p ǫ , we deduce that in the sense of distributions ∂ t n 0 − ∂ xx p 0 = n 10 G 1 (p 0 ) + n 20 G 2 (p 0 ).
Moreover, due to the uniform estimate on ∂ x p ǫ in L 2 (Q T ) from Lemma 4, we can show that, by passing into the limit in a product of a weak-strong convergence, in the sense of distributions, (n 10 , n 20 , p 0 ) satisfies (13)- (14).
Passing into the limit in the relation (1 − n ǫ )p ǫ = ǫn ǫ implies We can also pass to the limit for the segregation and deduce n 10 n 20 = 0. To conclude the proof of Theorem 1, we are left to establish the relation (17).

Complementary relation
In this section we want to pass to the limit in the equation for the pressure (6). However, this task can not be performed easily since we only have uniform estimates on the gradient of n and p, whereas we need strong convergence of the gradient to pass to the limit in (6). Then we propose to work on the time antiderivative. Let us denote q ǫ = p ǫ − ǫ ln(p ǫ + ǫ). Then, we have proved above that q ǫ → p 0 strongly as ǫ → 0, and Let us introduce Q ǫ a time antiderivative of q ǫ , Q ǫ (t, x) := t 0 q ǫ (s, x) ds. From the strong convergence of q ǫ , we deduce that Q ǫ → P 0 := t 0 p 0 (s, x) ds as ǫ → 0. By a simple time integration of (23), we have From Lemma 2, we deduce that ∂ xx Q ǫ is uniformly bounded in L 1 ∩ L ∞ ([0, T ] × R). Moreover, using the relation q ǫ = p ǫ − ǫ ln(p ǫ + ǫ), we get From the uniform bound on ∂ x p ǫ in L 2 (Q T ) in Lemma 4, we deduce that the sequence (∂ t ∂ x Q ǫ ) ǫ is uniformly bounded in L 2 ([0, T ] × R). Thus we have obtained that the sequence (∂ x Q ǫ ) ǫ is uniformly bounded in H 1 ([0, T ] × R). We deduce from the compact embedding of H 1 loc (Q T ) into L 2 loc (Q T ) that we can extract a subsequence, still denoted (∂ x Q ǫ ) ǫ , converging strongly in L 2 loc (Q T ) and weakly in H 1 (Q T ) towards a limit denoted ∂Q. Since Q ǫ → P 0 as ǫ → 0, we deduce ∂Q = ∂ x P 0 .
Thus, we can pass to the limit ǫ → 0 into the equation (24). We obtain Multiplying by p 0 and using the relation p 0 n 0 = p 0 , we deduce the complementary relation (17). This concludes the proof of Theorem 1.

Uniqueness of solutions
In this section, we focus on the uniqueness of solutions to the limiting problem (12)- (16). We first observe that from (12) and (16), we have Since we have the segregation property given by (16), we deduce that the support of n 10 and of n 20 are disjoints. Then, by taking test functions with support included in the support of n 10 or of n 20 in the weak formulation of (25), we deduce that We are going to prove that system (26)- (27) complemented with the segregation property (16) and the relation (15)  Proof. We follow the idea developped in [27] and adapt the Hilbert's duality method. Consider two solutions (n 10 , n 20 , p 0 ) and ( n 10 , n 20 , p 0 ) of the system (26)-(27)-(15)- (16). Making the difference and denoting q i = n i0 p 0 and q i = n i0 p 0 , for i = 1, 2, we have We first observe that on the set {n 10 > 0}∩ {p 0 > 0}, we have q 1 = p 0 from (15). Hence we have n 10 G 1 (p 0 ) = n 10 G 1 (q 1 ). The same observation holds for the other terms in the right hand side of these latter equations. Let Ω be a bounded domain containing the supports of all solutions, we denote Ω T = [0, T ] × Ω. For any suitable test functions ψ 1 and ψ 2 , we have, for i = 1, 2, This can be rewritten as, for i = 1, 2, where and we define A i = 0 as soon as n i0 = n i0 and B i = 0 as soon as q i = q i , whatever is the value of their denominators. It is shown in Lemma 5 below that, for i = 1, 2, we have 0 The idea of the Hilbert's duality method consists in solving the dual problem, which is defined here by, for any smooth function Φ i , i = 1, 2, in Ω T , If such a system admits a smooth solution, then, by choosing ψ i as a test function in (29), we get From the expression of A i , we deduce for any smooth function Φ i , i = 1, 2. It is obvious to deduce the uniqueness for the density. Uniqueness for the pressure will follow from (28). However, the dual problem (30) is not uniformly parabolic and its coefficients are not smooth. Then, in order to make this step rigorous, a regularization procedure is required. It can be done exactly as in [27, p 109-110]. For the sake of completeness of this paper, this regularizing procedure is recalled in Appendix A.
Proof. We observe that, for i = 1, 2, n i0 > n i0 implies q i ≥ q i . Indeed, either n i0 = 0 and then q i = 0 ≤ q i , or 0 < n i0 < 1 and then from the segregation property (16) we have n 0 = n i0 and from the relation (1 − n 0 ) p 0 = 0 we deduce that p 0 = 0, thus q i = 0 ≤ q i . Similarly, for i = 1, 2, n i0 > n i0 implies q i ≥ q i . By setting A i = 0 whenever n i0 = n i0 , we conclude that 0 ≤ A i ≤ 1.
Finally, the bound on C i is a direct consequence of the fact that G i is nonincreasing and Lipschitz (see (7)) and that 0 ≤ n i0 ≤ 1.

Numerical simulations 4.1 Numerical scheme
The numerical simulations are performed using a finite volume method similar as the one proposed in [12,14]. The scheme used for the conservative part is a classical explicit upwind scheme. To facilitate the reading of this paper, we recall here the scheme used. We divide the computational domain into finite-volume cells C j = [x j−1/2 , x j+1/2 ] of uniform size ∆x with and define the cell average of functions n 1 (t, x) and n 2 (t, x) on the cell C j bȳ The scheme is obtained by integrating system (1)-(2) over C j and is given bȳ where F k β,j+1/2 are numerical fluxes approximating −n k β u k β := −n k β ∂ x (p k β ) and defined by: with the discretized pressure We use the usual notation (u) + = max(u, 0) and (u) − = min(u, 0) for the positive part and, respectively, the negative part of u.
We recall that we have defined the parameters P 1 M and P 2 M as the values of the pressure for which the growth functions vanishe (see (7)). In this case their numerical values are given by P 1 M = 2 and P 2 M = 1. Then, we define Since the growth functions are different, clearly N 2 M ǫ < N 1 M ǫ . In Fig 1, the red and blue species are initially separated. As the time evolves the two species grow and reach their respective maximal packing values N 1 M ǫ and N 2 M ǫ . Once these values are reached, the species start to expand symmetrically in both directions until the two species collide. The collision creates an interface between the two species. From the figure, we can easily get that the total density is the envelop of the two densities n 1 and n 2 . We can observe this collision at t = 0.6 on Fig. 1 (d). The total densities at the interface then increase from 0 to the maximum packing value of the blue species (which is smaller than the one of the red species) while staying continuous. Once this value is reached (t = 1, 2 on both panel (e) and (f)), we observe two phenomena. First a bump is created on the left side of the interface, in the domain of n 2 . This bump help the total densities to stay continuous, as it joins the two maximal densities. It also means that, at the interface, the pressure is going to be higher than the limit pressure P 2 M . Then the derivative of the pressure at the interface is positive, which induces a motion of the interface representing the fact that the the red species n 1 pushes the blue species n 2 . This motion of the interface is the second phenomenon which is observed.

Influence of the parameter ǫ
In order to illustrate our main result on the limit ǫ → 0, we show, in this section, some numerical simulations of the model (1) We observe in Fig. 2 that the time dynamics of the numerical solutions is similar for each case and follows the dynamics presented above for the case ǫ = 1. The main difference observed is the maximal packing value N 1 M ǫ and N 2 M ǫ . Indeed since the maximal packing values are given by (34), when ǫ → 0, the maximal packing value converges to 1. This is consistent with the numerical results shown in Fig. 2. In addition we observe that as ǫ decreases the stiffness of the densities increases. In overall we observe that as ǫ → 0 densities converge to Heaviside functions with a support that slightly varies with respect to ǫ. Indeed, we observe, for instance, that the position of the right boundary of the species n 2 is smaller than 10 for ǫ = 1, 0.1 while it is bigger than 10 for ǫ = 0.01.

Particular solutions: tumor spheroid
One interested application of this study is tissue development. Since we consider a system with two populations of cells, we can for example consider the case of tumour with proliferative cells, whose density is denoted n 2 , and quiescent cells, whose density is denoted n 1 .
Solution of the limiting Hele-Shaw problem. We assume that initially the tumor is a spheroid centered in 0 and is composed by a spherical core representing the quiescent cells surrounded by a ring representing the proliferative cells. Then, we are looking for particular solution of the limiting Hele-Shaw problem (1)-(2) under the form: The two radius R 1 (t) and R 2 (t), with R 1 (t) < R 2 (t), are computed according to the geometric motion rules Such functions n 1 and n 2 are solutions to the limiting Hele-Shaw problem (1)- (2). Indeed by derivating the densities, in the distributional sense, we get, By applying the same computation on n 2 we get, Analytical solution. As this paper is reduced to the case of dimension 1, we can compute the exact solution of the limiting Hele-Shaw problem (1)-(2) with this initial configuration for some simple expression of the growth terms G 1 and G 2 . For instance, let us suppose that This choice means that as the pressure increases, the tumor will grow more slowly, until the pressure reach a critical value (P 1 M or P 2 M depending of the species) where the growth rate takes negative values, modelling the apoptosis of cells. The solution of the pressure equation is given by, on Ω 1 (t), (P 1 M − P 2 M ) √ g 1 sinh( √ g 2 (x−R 2 (t))) sinh( √ g 1 R 1 (t)) λ +P 2 M (1 − e √ g 2 (R 2 (t)−x) + µ λ e √ g 2 (R 2 (t)−R 1 (t)) sinh( √ g 2 (x − R 2 (t)))) on Ω 2 (t), 0 otherwise. with λ = √ g 1 sinh( √ g 2 (R 1 − R 2 )) sinh( √ g 1 R 1 ) − √ g 2 cosh( √ g 2 (R 1 − R 2 )) cosh( √ g 1 R 1 ), µ = √ g 1 sinh( √ g 1 R 1 ) + √ g 2 cosh( √ g 1 R 1 ). At last, in a third example we display an example where the species n 1 grows and pushes the surrounding species n 2 .
In Fig 3, we display the time dynamics of the densities of these three examples at different time step: (i) t = 0, (ii) t = 0.1, (iii) t = 0.3, (iv) t = 0.6, (v) t = 1. It illustrates the three different behaviours mentionned above by (36), (37) and (38). On Fig 3 (a) the blue species n 1 is initially surrounded by the red species n 2 . In Fig 3 (b) the red species grows and the blue species disappears since the pressure in the domain is bigger that P 1 M . Fig 3 (c) illustrates the second case where the blue species stabilizes. In Fig 3 (d), the blue species pushes the red species and propagates.

A Uniqueness of solutions: Regularized dual problem
In this appendix we prove rigorously Proposition 1 using a regularization procedure for the dual problem 30. We follow closely the ideas in [27, p 109-110] which are recall here for the sake of completness of this paper. Since the coefficients A i , B i are not strictly positive and not smooth, then we need to regularize the problem 30. For i = 1, 2, let A k i , B k i , C k i and G k i be sequences of smooth functions such that, for some constant α i , β i , δ 1,i , δ 2,i , M 1,i , M 2,i , K 1,i , K 2,i . For any smooth function Φ i , i = 1, 2, we consider the following regularised dual system,