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

1.   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 cell types contribute to mass gain for the whole structure. However, since the different cell populations grow at different rates, stresses arise 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. 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. Some aspects of the derivation of our model use an optimal transport framework. Indeed, we describe the cell populations by means of continuum densities satisfying gradient flow equations which decrease a free energy encompassing some important biological properties of the system. Optimal transport has proved a fertile concept for many types of Partial Differential Equation (PDE) such as the porous medium equation ^[1], convection-diffusion equations ^[2], the Fokker Planck equation ^[3], etc. However, other aspects of the system such as the inclusion of reaction terms which model tissue growth depart from a strict optimal transport framework. So, optimal transport will not be used in this paper beyond the derivation of the model.

After describing the model, we first present numerical simulations in order to analyse the roles of the various parameters of the model. Then, we investigate its incompressible limit. In this limit, the cell densities can only take two values: either zero or their respective maximal value corresponding to the packing density. We show that the limit model is a free boundary Hele-Shaw model (HSM), which allows us to focus on the geometric evolution of the boundaries of the domains 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 cells. On the one hand, discrete models consider each cell as an individual entity whose position and other attributes evolve in time ^[4,5]. 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 functions of space and time, which evolve according to suitable PDEs ^[6,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 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 ^[8,9,10]. Models of the second category describe the motion of the geometric boundary between the tissue or the tumour and its environment through geometric evolution equations ^{[11,12,13,14]}. The latter share similarities with Hele-Shaw models in fluid mechanics and Stefan's free boundary problem in solidification ^[15]. 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 ^[16,17,18].

In this paper, we propose a two cell population model as an extension of earlier models for single cell populations introduced in ^[19,16] 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\geq 0$ and position $x\in \mathbb{R}^d$. The diffusion of the density is triggered by a mechanical pressure $p = p(n)$ which is a given non linear function of the density $n$. Cell proliferation is modelled by a growth function $G = G(p)$ dependent on the pressure. The displacement of the cells occurs with a velocity $u = u(x, t)$ related to the pressure gradient through Darcy's law. The model is written as follows,

In ^[16,17,18] the pressure is expressed as

where $\gamma > 0$ is a model parameter. Inserting (1.3) into (1.1), (1.2) leads to the porous medium equation which has been widely studied ^[20]. 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. $\frac{\partial P}{\partial n} = p$.

The incompressible limit of this model corresponds to $\gamma \rightarrow +\infty$. It has been shown in ^[21] that this incompressible limit is a Hele-Shaw free boundary model which, classically, is used to describe the pattern of tumor growth ^[12,22]. 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 ^{[23,24,25,22,26]}. This incompressible limit and the corresponding Hele-Shaw model have been shown to be particularly relevant to tumor growth modelling ^[16,17].

In this paper, we present a new Segregation Pressure Model (SPM) for two cells populations which is built upon the gradient flow framework presented above. We introduce a free energy $\mathcal{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 similar to (1.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 cells of different types in order to enforce the segregation property, the latter being expressed as $r = n_1n_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 $\nabla n_1$, $\nabla n_2$ inside the expression of the free energy. 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 SPM which enables us to illustrate the validity of the limit HSM.

The pressure law (1.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 (1.3))

where $\epsilon > 0$ is a modelling parameter which plays a similar role as the parameter $1/\gamma$ in (1.3). With this expression, the pressure has a singularity at $n = 1$ which prevents the density to take values greater than $n = 1$. Similar pressure laws have been used in ^[21]. The limit $\gamma \rightarrow \infty$ is now replaced by $\epsilon \rightarrow 0$.

Systems with multiple populations are studied in many different areas. In chemistry, reaction-diffusion systems are used to model reacting chemical substances ^[27]. 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 ^[28] are used to model bacterial chemotaxis. Another classical example of cross-diffusion in biology is the Lotka-Volterra model ^[29], which describes the dynamics of a predator-prey system. These have been extended to nonlinear diffusion Lotka-Volterra systems to model cell populations ^[30,31]. In the context of tumor growth, systems with different types of cells have been studied (such as healthy/tumor cells, proliferative/quiescent cells ^[32,33,34]). Among these models ^[30,31,35], some preserve the segregation of initially segregated populations. On the other hand, some models generate segregation between initially non-segregated species. The ability of imposing segregation is a distinctive feature of the present work. The Cahn-Hilliard equation gives an example of a model that promotes segregation. It describes the process of phase separation ^[36] where, in the absence of growth terms and up to possible boundary effects, each phase tends to occupy a single connected domain from which it excludes the other phase. Similarities between our model and the Cahn-Hilliard equation will be described below.

The paper is divided into the following five sections. In Section 2 the two populations SPM model is described and numerical simulations which study the influence of the modelling parameters are shown. In Section 3, the main results are exposed: the formal incompressible limit theorem is stated; a formal proof of the convergence to the HSM free boundary problem is given; numerical simulations are shown in support and a discussion is provided. Then, Section 4 contains the description of the numerical scheme. A short conclusion is given in Section 5. Finally, Appendix 6 is devoted to the the derivation of the model from the free energy.

2.   The Segregation Pressure model

2.1.   Introduction of the continuum 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\geq 0$ and position $x\in \mathbb{R}^d$. We derive our model from the single cell model exposed in the introduction and define the free energy $\mathcal{E}(n_1, n_2)$ expressed as follows:

where $\alpha > 0$ is a diffusion parameter and $\epsilon > 0$, $m > 0$ are parameters which set the strengths of the congestion and segregation effects in the pressure laws. The functions $P_{\epsilon}$ and $Q_m$ are primitives of the functions $p_{\epsilon}$ and $q_m$ defined respectively by (1.5) and

From now, the total density will be denoted by $n = n_1+n_2$ and the product of the densities will be denoted by $r = n_1 n_2$. 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. The effect of these different terms are detailed below.

For the sake of simplicity, we omit the parameters $\epsilon, m, \alpha$ in the notations of the unknown functions $n_1$, $n_2$, $p_1$ and $p_2$. From the free energy (2.1), in Section 6, 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$. The pressures $p_1$ and $p_2$ are obtain thanks to the formula

We first comment on the first term of (2.1). We assume that the two categories of cells have identical volume, 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 parameter $\epsilon$ is supposed to be small. So, this first term penalizes configurations where the total density $n$ is close to the packing density (here assumed to be equal to $1$ according to (1.5)).

We now comment on the second term of (2.1). The segregation pressure is a novel aspect of the model. To minimise the second term of the free energy, it is necessary to make $n_1n_2$ as small as possible. In the particular case $\alpha = 0$, for a given density $n = n_1+n_2$, the minimiser of the free energy will be $(n_1^*, n_2^*)$ such that $n_1^*n_2^* = 0$ and $n = n_1^*+n_2^*$. This justifies why the repulsion term tends to increase the segregation of the populations as time evolves. Moreover this term imposes segregation when $m$ is going to infinity. This can be seen thanks to the equality $(1+r) (\frac{m-1}{m}q_m(r) +1) = (\frac{m-1}{m}q_m(r) +1)^{\frac{m}{m-1}}$. Passing to the limit $m \rightarrow \infty$, we obtain

meaning that $r^{\infty}(q^{\infty}+1) = 0$. Since $q^{\infty}\geq 0$, this implies that $r^{\infty} = 0$, which expresses the segregation property (the cell densities cannot be simultaneously non-zero).

The third term of the free energy is a cohesive energy. Indeed this term penalises the gradients of either densities, meaning that it tends to regroup each species in a single cluster. This term also allows the two populations to mix over a small width at their interface. The diffusion coefficient $\alpha$ is related to the width of the mixing region. To understand it, we make a rough order-of-magnitude calculation. We assume that the densities vary linearly in the mixing region supposed to occupy the region $[0, \ell]$. The density $n_1$ varies from a value close to $1$ to the value $0$ and $n_2$ varies in the opposite way, so: $n_1 \approx (1-\nu)(1-\frac{x}{l})$ and $n_2 \approx (1-\nu)\frac{x}{l}$ with $\nu \ll1$. We suppose that at the minimum of the free energy (2.1), the two last terms have the same order of magnitude in the interface zone. The first term of the free energy has been already taken into account by assuming that $n = n_1+n_2 = 1-\nu < 1$. We have,

So by integrating on $[0, \ell]$ we get by a change of variables $z = (1-y)y$:

When $m \gg 1$, the leading order term is $\frac{2}{m-1} \int_0^{\frac{1}{4}} (1+z)^m \frac{dz}{\sqrt{1-4z}}$ that we can estimate by neglecting $\frac{1}{\sqrt{1-4z}}$. We then get

where $\lambda$ is a given constant. The fourth order term can be estimated by

We deduce that

At this stage it is interesting to see the analogy between our model and the Cahn-Hilliard equation. We recall that the Cahn-Hilliard equation is a gradient flow (for the $L^2$ distance) of the free energy

where $d = d(x)$ is the unknown function and $\gamma > 0$ is a modelling coefficient. In order to see the similarity with our model, we rewrite (2.1) in terms of the total density $n$ and the difference $d = n_1-n_2$ and let $m = 2$,

with $\tilde{ \mathcal{E}}(n) = \int (1+P_{\epsilon})dx + \frac{\alpha}{ 4} \int |\nabla n|^2 dx$. Now, fixing $n$ and considering $d$ as the only variable, we can ignore the term $\tilde{ \mathcal{E}}(n)$. Then, the similarity of the first two terms of $\mathcal{E}(n, d)$ in (2.9) with (2.8) becomes clear. The only difference is that the two stable states $u = \pm 1$ of the Cahn-Hilliard energy are replaced by the states $d = \pm n$, which depend on $n$. Hence, we can view our model as a coupling between the Cahn-Hilliard equation for $d$, with a non linear parabolic equation for $n$. If we make $\alpha = 0$ in (2.1) and write the resulting equations in the unknowns $(n, d)$, we obtain an unstable diffusion system, the repulsion pressure giving rise to negative diffusion. The role of the terms in factor of $\alpha$ in (2.1) is to counteract this instability by introducing a stable fourth order diffusion.

Note that the particular case $q_m = 0$ and $\alpha = 0$ has been studied in ^[30,31,35]. There, it has been shown that the system exhibits species segregation, provided that the initial conditions are segregated ^{[32,33,37,38]}. 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 $\alpha$ and $m$. This is consistent 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 (2.3)-(2.6), which consists of letting $\epsilon$ going to 0 and $m$ going to infinity in the system.

2.2.   Numerical simulations

In this section we present numerical simulation for both segregated and mixed initial populations. This section aims to illustrate the dynamics of the system and the role of the parameters $\epsilon, m, \alpha$ in the volume exclusion pressure, the segregation pressure and the cohesive term. In order to facilitate the visualisation of the results, the simulation are performed in one dimension. The scheme used is described in Section 4.

2.2.1.   Numerical simulation with initially segregated and mixed populations

We illustrate the evolution in time of the SPM on two initial configurations. For the first example, we consider initial densities which are segregated and distant, given by

In the second example, the initial densities are mixed and given by

With these two examples, we investigate cases where the initial densities are either segregated or mixed. For both examples, the growth functions are given by

and the numerical parameters are $\epsilon = 0.01$, $m = 10$, $\alpha = 0.01$. The numerical simulation of the SPM model will involve the use of an approximate model, the Relaxation Segregation Pressure Model (RSPM), detailed in Section 4. This relaxation model is introduced to facilitate the numerical treatment of the fourth order term in the SPM. The RSPM model involves a relaxation parameter $\nu$ which is given the value $\nu = 0.001$. The results of the simulations of Examples 1 and 2 are respectively plotted in Figure 1 and Figure 2. On these Figures, the red line represents the species $n_1$ and the blue line, the species $n_2$ on panels (a)-(d) while the orange and green lines respectively represent the pressures $p$ and $q$ on panels (e)-(h). For both figures, the initial densities and pressures are plotted on panel (a) and (e).

In Figure 1 (a)-(b), we observe the dynamics of the two populations when they are not yet in contact with each other. The dynamics is then similar to that of the one species model (1.1)-(1.2) (since the density of one species is equal to zero on the support of the other species). When the densities are smaller than 1, the pressure is small (cf Figure 1(e)) and the reaction terms control the dynamics, resulting in the growth of the densities. When the densities $n_i$ reach the critical values $n^*_i = p^{-1}(p^*_i) = \frac{p^*_i}{\epsilon+p^*_i}$, with $p^*_1 = 20$ and $p^*_2 = 10$ ($p^*_i$ is the value for which the growth functions vanish), the pressure $p$ become significant and creates a moving front encompassing the domain where $n_i\approx n^*_i$ (cf Figure 1(c)-(d)). In panel Figure 1(c) the two populations meet, which creates an interface between them. In panel Figure 1(h), we observe a pressure ditch at the interface which is due to a small variation in the total density. This phenomenon will be investigated later when we will study the influence of the parameters $m$ and $\alpha$. Omitting the pressure ditch, we can consider that the gradient of the pressure is negative at the interface, which creates a movement of the interface toward the right. Indeed the red species, which has reached its critical density, is expanding and pushes the blue species to have more space to grow (cf Figure 1(d) and (h)).

In Figure 2, we observe the dynamics of two cell populations when they are initially mixed. When the overlapping between the two species is strong the dynamics is driven by the segregation pressure as $q$ is nonnegative (cf Figure 2(ⅰ)). This creates domains where the density of one species is close to its maximal value while the other species density is close to 0. Because the parameter $m$ which is the exponent of the segregation pressure is finite, some mixing is allowed between the two cell species. At the interfaces between the domains occupied by either species, we observe ditches in the pressure $p$, which are correlated with peaks in the segregation pressure $q$. In Figure 2(ⅰ-b)-(ⅰ-d), as the densities are close to their maximal value, the interfaces are moving in the directions given by the signs of the pressure gradients $\partial_x p$ (with $p$ plotted in Figure 2(ⅰ-f)-(ⅰ-h)). In this case, the red species is growing faster than the blue one, which explains why the inner species is pushing the outer one to have more space to grow. In Figure 2 (ⅱ), which corresponds to the same simulation as in Figure 2 (ⅰ) but for larger times, we observe on the left side that the red species pushes the blue species outside the domain, until only the red species remains (see Figure 2 (ⅱ-d)). In addition, in Figure 2 (ⅱ) we observe that a small pouch of the blue species surrounded by the red one gradually disappears. Because $20 = p_1^* > p_2^* = 10$, the pressure at the location of the blue species pouch becomes larger than $p_2^* = 10$ (as observed in Figure 2 (ⅱ-e) and Figure 2 (ⅱ-f)) and this triggers the decay of the blue population as its growth term becomes negative. At large times, only the faster growing species remains.

2.2.2.   Influence of the parameters

In Figure 1 and Figure 2, we observe some mixing between the two populations: first, at the interface between their respective domains the species mix in a small interval and create a ditch in the pressure $p$; second, in the central domain occupied by the red species in Figure 2, the density of the blue species is not exactly zero. To understand these phenomena we study the influence of the parameters $\epsilon, m, \alpha$ and plot the simulation results respectively in Figure 3 (ⅰ), (ⅱ), (ⅲ) for a fixed time $t = 0.15$. The initial densities and growth functions are the same as in Example 2 of the previous subsection and are defined by Eqs. (2.11) - (2.12). We vary the parameters about pivot values given by $\epsilon = 0.01$, $m = 10$, $\alpha = 0.01$, $\nu = 0.001$.

In Figure 3 (ⅰ-a), (ⅰ-b), (ⅰ-c) the densities are plotted for the cases $\epsilon = 1$, $\epsilon = 0.1$ and $\epsilon = 0.01$ respectively. We first remark that the parameter $\epsilon$ does not influence the amount of mixing of the two populations or the number of domains each species is divided in. Instead $\epsilon$ influences the densities at the boundaries of the domains occupied by each species. As $\epsilon$ decreases, the gradients of the densities become sharper at the exterior boundaries and at the interfaces between the two populations. In addition, the upper bound of the densities increases and becomes closer to 1.

In Figure 3 (ⅱ), the influence of the exponent of the segregation pressure $m$ is studied. The densities are plotted for cases $m = 5$, $m = 10$ and $m = 50$ in Figure 3 (ⅱ-a), (ⅱ-b), (ⅱ-c) respectively. In Figure 3(ⅱ-c), we can observe that there is less mixing than in Figures 3(ⅱ-a) and (ⅱ-b). In Figure 3(ⅱ-c) the central domain of the red species is smaller than in Figure 3(ⅱ-b) because the blue species has grown in the middle as a result of a larger segregation pressure. Then in Figure 3(ⅱ-c), we do not observe the small dimple in the blue cell density found in Figure 3(ⅱ-a) and (ⅱ-b). To summarise, the parameter $m$ controls the strength of the segregation between the two populations. However, even when $m$ is large, we still observe some mixing at the interfaces.

In Figure 3 (ⅲ), we study the influence of the parameter $\alpha$ in factor of the fourth order diffusion term. The densities are plotted for the cases $\alpha = 0.1, \nu = 0.001$; $\alpha = 0.01, \nu = 0.001$ and $\alpha = 0.001, \nu = 0.0001$ in panels Figure 3 (ⅲ-a), (ⅲ-b), (ⅲ-c) respectively (we have observed that the parameter $\nu$ needs to be smaller than $\alpha$). We can observe that as $\alpha$ becomes smaller, the width of the interface region between the two populations gets smaller. We can also observe in Figure 3(ⅲ-c) that on the right hand side, thin stripes with alternating populations appear. However, even when such stripes appear, there is still some mixing between the two populations, by contrast with the simulations where $m$ was higher in Figure 3(ⅱ-c). The same kind of features, namely the formation of many thin strips alternating $n_1 = 1$, $n_2 = 0$ and $n_1 = 0$, $n_2 = 1$, were observed in simulations performed with $\alpha = 0$. This was observed both with the same numerical scheme as that used in this paper, and with schemes using the gradient flow structure of the system. Without the fourth order term, a species surrounded by an other one will prefer to split, with some proportion of the population jumping across the other species instead of pushing it. For these reasons, the fourth order term is essential to produce realistic dynamics. What this analysis shows is that the choice of the parameter $\alpha$ is critical to prevent this jumping behaviour and maintain the connectivity of the domains occupied by the species.

3.   Incompressible limit $\epsilon \rightarrow 0, m \rightarrow +\infty$ and $\alpha \rightarrow 0$

3.1.   Formal limit

In this Section, we obtain formal convergence results of the model when $\epsilon \rightarrow 0$, $m \rightarrow \infty$ and $\alpha \rightarrow 0$. 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. In what follows, the growth functions will be chosen as $G(p) = g(p^*-p)$. As for the initial conditions, we assume that there exists $\epsilon_0 > 0$ such that, for all $\epsilon\in (0, \epsilon_0)$,

Thanks to these assumptions, it is possible to establish a priori estimates on ${n_1}$ and ${n_2}$ such as positivity and $L^1$ bounds. However, in order to pass to the limits $\epsilon \to 0$, $m \to \infty$ and $\alpha \rightarrow 0$, we need more a priori estimates. In the following theorem, we only provide a formal limit. Indeed, the appearance of terms of the type $\nabla n \cdot \nabla r$ which we are not able to control prevent us from obtaining such theorem. The fourth order term is also difficult to control, as it does not enable us to use maximal principle arguments.

The main result is stated in the following theorem.

Theorem 1. (Formal limit) Let $T > 0$, $Q_T = (0, T)\times \mathbb{R}^d$. Let $G_1$, $G_2$ and $n^{\rm ini}$, $n_1^{\rm ini}$, $n_2^{\rm sini}$ satisfy assumptions (3.1) and (3.2). Suppose the limits of the densities ${n_1}(x, t)$, ${n_2}(x, t)$, of the pressure $p_{\epsilon}(x, t) = p_\epsilon ({n_1}(x, t), {n_2}(x, t))$ and of $q_m(x, t) = q_m ({n_1}(x, t), {n_2}(x, t))$ as $\epsilon \rightarrow 0$, $m \rightarrow +\infty$ and $\alpha \rightarrow 0$ exist and are denoted by $n_1^{\infty}, n_2^{\infty}$, $p^{\infty}$ and $q^{\infty}$. If the convergence is strong enough then these limits satisfy

As a consequence of (3.3) and (3.4) we have,

In addition, we have

with $n^{\infty} = n_1^{\infty} +n_2^{\infty}$,

and the complementary relation

Remark 1. The assumption $n_1^{\rm ini}n_2^{\rm ini} = 0$ is not needed in Theorem 1. Since in the limit, the segregation is imposed by (3.7), the system reorganises instantaneously and the initial density of the limit model must verify ${n_1^{\infty}}^{\rm ini}{n_2^{\infty}}^{\rm ini} = 0$. The densities ${n_1^{\infty}}^{\rm ini}$ and ${n_2^{\infty}}^{\rm ini}$ could be determined through an initial-layer analysis.

Firstly, it is important to notice that, in the limit, $p^{\infty}$ and $q^{\infty}$ become independent of the densities $n_1^{\infty}$ and $n_2^{\infty}$. Indeed, looking at the expression of the pressures (1.5)-(2.2) in the limit, we obtain:

The pressure $p^{\infty}$ can be viewed as a volume exclusion pressure, while the pressures $n_1^{\infty} q^{\infty}$ and $n_2^{\infty} q^{\infty}$ are segregation pressures. Eq. (3.6) suggests to decompose the domain in two parts. We consider the domain $\Omega(t) = \{ x \; |\; {p}^{\infty} (x, t) > 0 \}$ and the complementary domain, where the pressure is equal to 0. Notice that, because of Eq. (3.6), $\Omega(t) \subset \{ x \; |\; {n}^{\infty}(x, t) = 1 \}$. Moreover, the two domains coincide almost everywhere. Indeed, if ${n}^{\infty} = 1$ and ${p}^{\infty} = 0$, because of the reaction term in (3.5) the density would grow following

and the total density would become greater than the maximum density. If initially the pressure ${p}^{\infty} = 0$ and the total density ${n}^{\infty}$ is strictly smaller than 1, the total density will follow Eq. (3.10). It will grow until the value 1 is reached, and then ${p}^{\infty}$ will become nonnegative. Thanks to the segregation property (3.7) we can decompose $\Omega(t)$ in two subdomain, $\Omega_1(t) = \Omega(t) \cap \{ x \; |\; {n_1}^{\infty}(x, t) = 1 \}$ and $\Omega_2(t) = \Omega(t) \cap \{ x \; |\; {n_2}^{\infty}(x, t) = 1 \}$. Then thanks to (3.7), $\Omega_1(t)$ and $\Omega_2(t)$ are disjoint and their union is $\Omega(t)$.

Remark 2. It is interesting to remark that

The pressure $p_1^{\infty}$ is driving the movement of the density $n_1^{\infty}$ which is equal to $0$ on $\Omega_2(t)$. So $q^{\infty}$ does not appear in the limit equation of $n_1^{\infty}$. The same reasoning can be applied to $n_2^{\infty}$. It shows that the pressure $q^{\infty}$ does not influence the limit model. Then at the limit the dynamics of the densities $n_1^{\infty}$ and $n_2^{\infty}$ are driven by $\nabla p^{\infty}$ even though $p_1^{\infty} \neq p_2^{\infty}$.

The equation for the pressure $p^{\infty}$ inside $\Omega(t)$ is deduced from the complementary relation (3.8) and reads:

Thanks to (3.6), the pressure $p^{\infty}$ is equal to zero on the boundary $\partial \Omega(t)$. However the normal derivative of the pressure is non-zero and controls the movement of the domain boundary. The exterior boundary is denoted by $\Sigma$ and the interface between the two populations is called $\Gamma$. Knowing that in the system (2.3)-(2.6) the densities $n_1$ and $n_2$ diffuse with velocities $\nabla p$ (see Remark 2) which in the limit, is equal to $\nabla p^{\infty}$ respectively on $\Omega_1$ and $\Omega_2$, we deduce that the boundaries $\Sigma$ and $\Gamma$ both move in the normal direction with velocities

where $\upsilon$ is any of the two opposite unit normal vectors to $\Gamma$. The velocity $V_{\Sigma}$ is normal to $\Sigma$. Indeed, $p^\infty$ is constant equal to zero along $\Sigma$ from (3.6), so, its gradient must be in the normal direction to $\Sigma$. On the other hand, $\Gamma$ does not need to be a level surface of $p^\infty$, so $\nabla p^{\infty}$ must be projected along the normal direction to $\Gamma$.

To help understand the expression of the velocity (3.12) of the limit model, we now consider a particular solution of the HSM. We consider the case where the initial densities are defined by $n_1^{\rm ini} = {\mathbb{1}}_{B_{R_1^{\rm ini}}}$ and $n_2^{\rm ini}(t, x) = {\mathbb{1}}_{B_{R_2^{\rm ini}} \backslash B_{R_1^{\rm ini}}}$ with $B_R$ the ball of center 0 and radius $R$ and ${\mathbb{1}}_S$ is the indicator function of the set $S$. We suppose that $R_1^{\rm ini} < R_2^{\rm ini}$, then $n^{\rm ini} = {\mathbb{1}}_{B_{R_2^{\rm ini}}}$. We are looking for solutions of the HSM of the form $n_1^{\infty}(x, t) = {\mathbb{1}}_{B_{R_1(t)}(x)}$ and $n_2^{\infty}(x, t) = {\mathbb{1}}_{{B_{R_2(t)}}\backslash{B_{R_1(t)}}(x)}$ with $R_1(0) = R_1^{\rm ini}$ and $R_2(0) = R_2^{\rm ini}$. Thus we have $n(x, t) = {\mathbb{1}}_{B_{R_2(t)}}$ and so, we get:

Using Eq. (3.5) and knowing that $n_1^{\infty}$, $n_2^{\infty}$ are constant equal to 0 or 1, we obtain that $n^{\infty}$ is solution in the sense of the distribution of

Then for all $\varphi \in C_c^{\infty}(\mathbb{R}^d)$

Since $p^{\infty}$ verifies the complementary relation (3.8), applying the Green formula twice we get

where $d\sigma$ is the surface element and $\upsilon$ is the normal directed outward the domain $B(R_2(t))$. Using radial symmetry and knowing that the pressure $p^{\infty}$ is equal to $0$ at the interface, we obtain

Since $\partial B(R_2(t))$ is a level surface of the pressure $p^{\infty}$, its gradient is normal to $\partial B(R_2(t))$. Then identifying Eqs. (3.13) and (3.14), we have

where, by abuse of notation, we have denoted $p^{\infty}(r, t)$ the constant value of $p^{\infty}(., t)$ on the surface $\{ |x| = r\}$. Similarly, applying the same result on the species $n_1^{\infty}$, we obtain

The limit HSM is graphically represented in Figure 4.

3.2.   Formal proof of Theorem 1

This section is dedicated to the formal proof of Theorem 1. First, we prove Eqs. (3.6) and (3.7) 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 $\epsilon \rightarrow 0, m \rightarrow +\infty$. Third, we prove Eq. (3.8), 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.

Proof. Eqs. (3.6) and (3.7) follow from the pressure laws (1.5) and (2.2). Indeed, since $(1-n)p_{\epsilon} = \epsilon n$, it follows in the limit $\epsilon \rightarrow 0$ that

Working out the repulsion pressure law (2.2) leads to the formula

Taking the limit $m \rightarrow +\infty$, we have

which implies that $r^{\infty} (q^{\infty}+1) = 0$. Since $q_m \geq 0$ for all $m$, we deduce $q^{\infty}+1 > 0$ and

Eqs. (3.3) and (3.4) follow directly from taking the limit $\epsilon \rightarrow 0, m \rightarrow +\infty$ and $\alpha \rightarrow 0$ in (2.3) and (2.4). To recover (3.5), it is first interesting to remark that the solution of (2.3)-(2.6) satisfies

where

Indeed, by adding (2.3) and (2.4), we obtain

Then given formula (1.5),

and given formula (2.2),

and inserting (3.18), (3.19) into (3.17) gives (3.15), (3.16). Let us denote

and $h^{\infty}$ its limit when $\epsilon \rightarrow 0, m\rightarrow\infty$. Since $(r+1)^{m} = (1+r) (\frac{m-1}{m} q_m(r)+1)$ and $r^{\infty} = 0$, passing to the limit as $\epsilon \rightarrow 0, m\rightarrow\infty$, we obtain

So at the limit, we have:

To recover the complementary relation (3.8), we need to compute an evolution equation for the pressure. To do so, we multiply the density equation (3.15) by $p_{\epsilon}'(n)$ and obtain,

Multiplying Eq. (3.21) by $\epsilon$ and taking into account that $p_{\epsilon}'(n) = \frac{1}{\epsilon} (p_{\epsilon}+\epsilon)^2$ yields

We now take the limit $\epsilon \rightarrow 0$, $m\rightarrow \infty$ and $\alpha \rightarrow 0$, and using the expression (3.20) for the limit of $H_{\epsilon, m}$, we get,

It leads to the complementary relation (3.8).This concludes the formal proof of the theorem.

Now, we prove Formula (3.12). We first focus on the speed of the exterior boundary $\Sigma$. Thanks to Eq. (3.5), for all $\varphi \in C_c^{\infty}(\mathbb{R}^d)$ and using Green's formula

Since $p^{\infty}$ is equal to $0$ outside $\Omega(t)$ and because $n_1^{\infty}$ and $n_2^{\infty}$ are both constant equal to $1$ on $\Omega_1(t)$ and $\Omega_2(t)$ respectively and equal to $0$ outside, we have

Therefore applying the Green formula twice, yields

where $\upsilon$ is the unit normal vector. For $x \in \Gamma$, we denote

and

On $\partial \Omega(t)$, $\partial \Omega_1(t)$, $\partial \Omega_2(t)$ the unit normals $\upsilon$ are directed outwards to the respective domains. On $\Gamma(t)$, the normal $\upsilon$ is directed from $\Omega_2(t)$ towards $\Omega_1(t)$. From Eq. (3.8), the first integral in (3.23) is equal to 0. In addition, the pressure is equal to zero on the boundary $\partial \Omega(t)$ so that

Moreover since $n^{\infty} = 1$ in the domain $\Omega(t)$, we have

where $V_{\partial\Omega(t)}$ is the speed of the boundary $\partial\Omega(t)$ and is directed along the normal $\upsilon$ in the outwards direction. Since Eqs. (3.24) and (3.25) are verified for all $\varphi \in C_c^{\infty}(\mathbb{R}^d)$, we deduce that

and at the interface $\Gamma$,

Then, both $p^{\infty}$ and its normal derivative are continuous at the interface $\Gamma$. Since $p^{\infty}$ is equal to zero along $\partial\Omega$, its gradient is normal to $\partial\Omega$ and the velocity on the outer domain given by (3.26) can be rewritten

To find the velocity at the interface the same method needs to be applied on either $n_1^{\infty}$ or $n_2^{\infty}$ and it leads to

where $V_{\Gamma(t)}$ is the speed of the boundary $\Gamma(t)$. This gives (3.12).

3.3.   Numerical validation

3.3.1.   Analytical solution of the Hele-Shaw model

The HSM is characterised by the complementary relation (3.8) and the velocity of the boundary given by (3.12). In the one-dimensional case (1D), the solution of this problem can be computed explicitely. 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 terms in the literature ^[16,17,18]. In the case of a single population, the complementary relation (3.8) 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_2^{\mbox{ ini}}$ is defined as the indicator function of $[x^{\mbox{ ini}}_{\Gamma^-}; x^{\mbox{ ini}}_{\Gamma^+} ]$ and the density $n_1^{\mbox{ ini}}$ is defined as the indicator function of $[x^{\mbox{ ini}}_{\Sigma^-}; x^{\mbox{ ini}}_{\Sigma^+} ] \setminus{ [x^{\mbox{ ini}}_{\Gamma^-}; x^{\mbox{ ini}}_{\Gamma^+} ]}$ where $x^{\mbox{ ini}}_{\Gamma^-}, x^{\mbox{ ini}}_{\Gamma^+}$ represents the interfaces between the two populations and $x^{\mbox{ ini}}_{\Sigma^-}, x^{\mbox{ ini}}_{\Sigma^+}$ represent the exterior boundaries of the support of the total density. More specifically,

Since the two populations have the same growth function, the complementary relation can be treated as that of a single population of cells (3.29) with the boundary conditions $p(x_{\Sigma^-}) = 0$ and $p(x_{\Sigma^+}) = 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_{\Sigma^{\pm}}$ and $x_{\Gamma^\pm}$ evolve respectively according to $\frac{d }{dt} x_{\Sigma^{\pm}} = V_{\Sigma^\pm}$ and $\frac{d }{dt} x_{\Gamma^{\pm}} = V_{\Gamma^\pm}$. Since initially $x^{\mbox{ ini}}_{\Sigma^-} < x^{\mbox{ ini}}_{\Gamma^-} < x^{\mbox{ ini}}_{\Gamma^+} < x^{\mbox{ ini}}_{\Sigma^+}$, the velocities $V_{\Sigma^-}$, $V_{\Gamma^-}$ are nonpositive and the velocities $V_{\Sigma^+}$, $V_{\Gamma^+}$ are nonnegative. It means that the densities $n_1$ and $n_2$ spread. However $|V_{\Gamma^-}| \leq V_{\Sigma^-}$ and $|V_{\Gamma^+}| \leq V_{\Sigma^-}$ 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^{\mbox{ ini}}_{\Sigma^-}; x^{\mbox{ ini}}_{\Gamma} ]$ and $[x^{\mbox{ ini}}_{\Gamma}; x^{\mbox{ ini}}_{\Sigma^+} ]$ where $x^{\mbox{ ini}}_{\Sigma^-}, x^{\mbox{ ini}}_{\Sigma^+}$ define the boundaries of the support of the total density and $x^{\mbox{ ini}}_{\Gamma}$ defines the interface between the supports of the two densities. More specifically,

The pressure follows the complementary relation (3.8), which can be rewritten as

with the additional conditions $p(x_{\Sigma^-}) = 0$, $p(x_{\Sigma^+}) = 0$, and $p$ and $p'$ are continuous at $x_{\Gamma}$. After some computations we find,

with

and

From this expressions we can compute the speeds of the exterior boundary and of the interface by computing the derivative of the pressure and evaluating it at $x^{\mbox{ ini}}_{\Sigma^-}, x^{\mbox{ ini}}_{\Sigma^+}$ and $x^{\mbox{ ini}}_{\Gamma}$.

3.3.2.   Numerical validation

In this section we use the numerical scheme presented in Section 4 to illustrate that the solution of the SPM converges to the corresponding solution of the HSM described in Section 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 Section 3.2. For Example 1, we plot the simulation for one set of parameters: $\epsilon = 0.1$, $m = 10$, $\alpha = 0.01$ and $\nu = 0.001$ (where $\nu$ is introduced for the RSPM in Section 4). For Example 2, we plot the simulation for a set of parameters: $(\epsilon, m, \alpha, \nu) = (1, 2, 0.1, 0.01)$, $(\epsilon, m, \alpha, \nu) = (0.1, 10, 0.01, 0.001)$, $(\epsilon, m, \alpha, \nu) = (0.02, 50, 0.001, 0.0001)$.

The first case defined as Example 1 in Section 3.2 is illustrated in Figure 5. We use the growth function (3.30) with the parameters $g = 1$ and $p^* = 10$. The boundary and interface are taken as $x^{\mbox{ ini}}_{\Sigma^\pm} = \pm 1.4$ and $x^{\mbox{ ini}}_{\Gamma^\pm} = \pm 0.6$. Then the initial density of the HSM are defined by Eq. (3.31). As initial conditions of the SPM, we take

We take $n_1 = n_2 = 0.98$ initially as it is close to the singular value 1. We plot the numerical solutions of the SPM (solid line) and of the HSM (dashed line) on the same figure. The species $n_1$ and $n_2$ are plotted in red and blue, respectively. In this example the solution of the SPM is expected to be close to the particular solution of the HSM given in Example 1.

In Figures 5, panel (a) is for the initial configuration and panels (b), (c), (d) are for times $t = 0.001$, $t = 0.01$, $t = 0.15$ 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 SPM exhibits sharp ditches at the interfaces between the two populations. These ditches appear because of the fourth order term in the model (2.3)-(2.6). Indeed, because of the fourth order term, $n_1$ and $n_2$ are continuous. Since $n_1$ is zero in $\Omega_1$ and $n_2$ is zero in $\Omega_1$, both are zero at the interface and since they are continuous, they are both very small in a small vicinity of the interface. Therefore the repulsion pressure $p$ is also small at the interface. Since the width of the the interface is of the order of alpha and is very small, this implies that the region of small pressure is very narrow and generates a sharp ditch near the interface. Away from the interface, the pressure does not feel the influence of the fourth order term because all densities are smooth and therefore, it approximates very well the pressure given by the HSM model. Away from the interfaces, the pressures given by the two models are similar except at small times. In Figure 5(b), the pressure of the SPM is larger than that of the HSM. However, we do not observe any influence of this difference on the densities. On Figures 5(c) and 5(d), the two pressures are almost the same for the two models. Overall, the dynamics of the SPM and the HSM are quite similar. In Figure 5, we observe that the blue (inner) species pushes the red (outer) species in order to be allowed to grow.

The second case is illustrated in Figure 6. The growth functions are defined by (3.32) with the parameters $g_1 = g_2 = 1$, $p_1^* = 20$ and $p_2^* = 10$. As initial conditions of the SPM we take

where $x^{\mbox{ ini}}_{\Sigma^\pm} = \pm 1.4$ and $x^{\mbox{ ini}}_{\Gamma} = 0$. The initial density of the HSM are defined with formula (3.33). In this example, we plot the solution of the HSM (dashed line) in Figure 6 (ⅰ) and the solution of the SPM (solid line) in Figure 6 (ⅱ)-(ⅳ) (respectively for the parameters $(\epsilon, m, \alpha, \nu) = (1, 2, 0.1, 0.01)$, $(\epsilon, m, \alpha, \nu) = (0.1, 10, 0.01, 0.001)$, $(\epsilon, m, \alpha, \nu) = (0.02, 50, 0.001, 0.0001)$). For each line, panels (a), (b), (c) are for times $t = 0.001$, $t = 0.01$, $t = 0.15$. Panel (d) is also a plot at time $t = 0.15$, but only of the densities. The species $n_1$ and $n_2$ are plotted in red and blue, respectively and the pressure $p$ is plotted in black.

In Figure 6 (ⅰ)-(ⅳ), the red species (leftmost) which has the biggest growth rate pushes the blue species (rightmost) toward the right. Therefore by growing more rapidly, the red species exerts a bigger pressure on the blue species than the blue species exerts on the red species. This pressure imbalance which is visible on Figure 6 (a)-(c), triggers the motion of the interface towards regions of lower pressure.

In Figures 6, while going from line (ⅱ) to (ⅳ) (which corresponds to having $\epsilon \rightarrow 0$ and $m \rightarrow \infty$), we observe three important phenomena. As seen previously in Figures 5, for small times the pressure $p$ of the SPM and the HSM are different. However this difference is much bigger in Figures 6 (ⅱ-a) when $(\epsilon, m, \alpha, \nu) = (1, 2, 0.1, 0.01)$ than for the two other cases. So as the parameters $\epsilon$, $m$ converge, the time for the two models to provide similar solutions reduces. This corresponds to an initial layer converging faster to the initial state of the HSM. The other main difference is the pressure ditch at the interface as shown in Figures 6 (a), (b), (c). As $\epsilon$ and $m$ converge, the ditch becomes deeper. Indeed, as $\epsilon$ and $m$ converge, the variations of $n_1$ and $n_2$ near the interface become sharper. Since the pressure near the interface is a function of the densities, its variations become sharper as well, which is consistent with the observation that the pressure ditch becomes more pronounced as these parameters converge. Except near the interface, the pressure $p$ of the SPM and the HSM almost coincide in Figures 6 (b), (c). The last phenomenon we observe is the convergence of the densities of the SPM to the Heaviside functions which characterize the solutions of hte HSM in Figure 6 (d). Indeed as $\epsilon$ gets smaller (and $m$ gets bigger), the fronts of the densities become steeper and the upper bounds of the densities get closer to 1.

3.3.3.   Discussion

The results of the simulations presented in Figures 5 and 6 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 SPM and of the HSM. This is because the initial conditions of the SPM are not well-prepared for the limit model, which generates an initial layer whose function is to relax the initial conditions to well-prepared ones. Indeed, initially in the SPM, the densities are taken as indicator functions, so at the start the pressure $p_{\epsilon}$ of the SPM is also an indicator function whereas that of the HSM is continuous. It takes an initial transient for the model to smooth it out, and then the two pressures coincide. In addition, the initial value of the SPM in these simulations is fixed to $98\%$ of the packing density 1. At the beginning of the simulation, the dynamics of the SPM is driven by the growth terms which increase the densities to their upper bound values $n_M = p_{\epsilon}^{-1}(p^*) = \frac{p^*}{\epsilon +p^*} < 1$ with $p^* = \max(p_1^*, p_2^*)$ (see (3.2)). Until the densities reach their upper bound value, the SPM and the HSM are not close to each other. After this short period, we observe that the pressures of the two models fit well. It also seems that this time becomes smaller when the parameters are closer to convergence, which is consistent with the fact that initial transients are faster when the perturbation parameters are closer to their limit values.

The major difference between the two models is at the interfaces between the two populations. As a result of the fourth order terms, the densities of the SPM are not perfectly segregated and overlap only in a narrow spatial interval. As long $\epsilon$ and $m$ are finite, the densities $n_1$ and $n_2$, which are continuous, are small in the vicinity of the interface. This creates sharp ditches in the pressures (which in this regime are function of the densities). In the limit $\epsilon \rightarrow 0$, $m \rightarrow +\infty$ and $\alpha \rightarrow 0$, the densities become discontinuous and the pressure ditches become narrower and narrower until they disappear completely. For finite $\epsilon$ and $m$ close to the limit, despite the pressure ditches, the interface and boundary speeds are very close in the two models. These are quite remarkable results, considering that the parameters used in the tests were not yet in the asymptotic ranges where the two models should be identical.

However, it is important to remark that the comparisons between the SPM and the HSM can only be made in the case of initially segregated populations. This is necessary to initialize the free boundary model. In the case of mixed population, simulations of the SPM have been made and can be found in Section 2. They confirm the convergence of the SPM towards a free boundary model since we observe that the system self-organizes into separated domains containing the different populations. However we cannot study the convergence towards the HSM since we do not know beforehand which initial conditions for the HSM will correspond to the converged SPM.

4.   Numerical method

This section is devoted to the derivation of the numerical scheme used to perform the simulations of the SPM presented in Sections 2 and 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 derive a CFL stability condition.

4.1.   Relaxation model

In order to simplify the computation of system (2.3)-(2.6), we lower its order through a suitable relaxation approximation. The relaxation system, depending on the relaxation parameter $\nu$, is written as follows

with $p_1$ and $p_2$ given by the constitutive laws (2.5)-(2.6) and where $c_1$ and $c_2$ are two new variables. The fourth order terms in (2.3)-(2.6) are now replaced by second order terms supplemented with Poisson equations which define $c_1$ and $c_2$. This relaxation method has been introduced for the numerical approximation of the Navier-Stokes-Korteweg equations ^[39,40]. This system is referred to as Relaxation Segregation Pressure Model (RSPM). It formally converges toward the SPM when $\nu \rightarrow 0$. Here we give supporting elements for this statement. We suppose that $(n_1^{\nu})$ and $(n_2^{\nu})$ converge when $\nu$ goes to 0. Since the two species play the same role, we only consider $n_1$. Inserting (4.3) in (4.1) we have,

Then to prove the convergence we need to show that $c_1^{\nu} \rightarrow n_1$ when $\nu \rightarrow 0$. We can rewrite (4.3) as

So that

Using Fourier transform, we obtain:

And then,

With the same computation we can obtain that

This formally shows the convergence of the RSPM toward the SPM.

4.2.   The scheme

We aim to numerically solve the RSPM given by (4.1)-(4.4) 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 $\nu$. For the sake of simplicity, we only consider the 1D case on a finite interval $\Omega = (a, b)$.

We divide the computational domain into finite-volume cells $C_j = [x_{j-1/2}, x_{j+1/2} ]$ of uniform size $\Delta x$ with $x_j = j\Delta x, j\in \{1, ..., M_x\}$, and $x_{j} = \frac{x_{j-1/2} +x_{j+1/2} }{2}$ so that

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 (4.1)-(4.4) over $C_j$ and is given by

Here, $c_{\beta, {j}}(t) \approx c_\beta(t, x_j), \ \beta \in \{1, 2\}$ and $F_{\beta, j+1/2}$ are numerical fluxes approximating $-n_{\beta} u_{\beta}: = -n_{\beta} \partial_x(p_{\beta}-\frac{1}{\nu}(c_{\beta}-n_{\beta}))$ and defined by the following upwind scheme:

where

and $(u_{{\beta}, {j+1/2}})^+ = \max(u_{{\beta}, {j+1/2}}, 0)$ and $(u_{{\beta}, {j+1/2}})^- = \min(u_{{\beta}, {j+1/2}}, 0)$ are its positive and negative part, respectively. The numerical pressure $p_{\beta, {j+1}}$ for $\beta = 1, 2$ are computed applying (2.5)-(2.6) to the $\bar{n}_{\beta, j}$.

It is easy to see that scheme (4.5)-(4.7) 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. ^[41].

4.3.   CFL condition and properties of the scheme

In this section, we present some properties of the scheme. We prove that for all times $t\ge 0$, it preserves the non-negativity of the computed densities ${\bar{n}}_{1}$ and ${\bar{n}}_{2}$ and also ensures that the total density $\bar{n}$ stays below $1$. The former property guarantees physically meaningful values of the computed densities while the second one is needed to make sense of the pressure law which blows up when ${\bar{n}} \rightarrow1$.

We consider the semi-discrete finite volume scheme (4.5)-(4.7) and evolve it in time by the forward Euler method. We denote by ${\bar{n}}^{k}_{1, j}$, ${\bar{n}}^{k}_{2, j}$, $p^{k}_{1, j}$ and $p^{k}_{2, j}$ the computed densities and pressures obtained at time $t^k = t^{k-1}+\Delta t^k$, i.e., ${\bar{n}}^{k}_{1, j}: = {\bar{n}}_{1, j}(t^k)$, ${\bar{n}}^{k}_{2, j}: = {\bar{n}}_{2, 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.

Proposition 1. (Positivity of the density) Consider the semi-discrete finite volume scheme (4.5)-(4.7) that is evolved in time by the forward Euler method. Provided that the initial densities ${\bar{n}}^{0}_{1, j} \geq 0$ and ${\bar{n}}^{0}_{2, j}\geq 0$ for all $j \in \{1, ..., M_x\}$, and that the growth terms $G_{\beta}$ are nonnegative for $\beta = \{1, 2\}$, a sufficient CFL condition for the cell averages ${\bar{n}}^{k}_{1_j}$ and ${\bar{n}}^{k}_{2_j}$, with $k \in [0, N_{it}]$, to be positive is

where $({u^{k}_{\beta, {j+1/2}}})^+ : = (u_{\beta, {j+1/2}})^+ (t^k)$ and $({u^{k}_{\beta, {j+1/2}}})^- : = (u_{\beta, {j+1/2}})^- (t^k)$.

Proof. Assume that at a given time $t^k$, ${\bar{n}}^{k}_{1, j} \geq 0$ and ${\bar{n}}^{k}_{2, j}\geq 0$ for all $j \in \{1, ..., M_x\}$. Then, the new densities are given by the general formula

where $F^{k}_{\beta, j+1/2} : = F_{\beta, j+1/2}(t^k)$. Taking into account formula (4.6) for fluxes $F^{k}_{\beta, j+1/2}$ and the fact that the growth terms $G_{\beta}$ are nonnegative for $\beta = 1, 2$, we obtain

By definition, $(u^{k}_{\beta, {j+1/2}})^-\leq 0$ and $(u^{k}_{\beta, {j+1/2}})^+\geq 0$ for all $j \in \{1, ..., M_x\}$. Then, provided that the CFL condition (4.8) is satisfied,

Since $({\bar{n}}^{k}_{1, j})_{j \in \{1, ..., M_x\}}, ({\bar{n}}^{k}_{2, j})_{j \in \{1, ..., M_x\}}$ are non-negative, we conclude that ${\bar{n}}^{k+1}_{1, j} \geq 0$ and ${\bar{n}}^{k+1}_{2, j}\geq 0$ for all $j \in \{1, ..., M_x\}$.

Proposition 2. (Maximum total density) Consider the semi-discrete finite volume scheme (4.5)-(4.7) that is evolved in time by the forward Euler method. Provided that the initial total density ${\bar{n}}^{0}_j = {\bar{n}}^{0}_{1, j} +{\bar{n}}^{0}_{2, j} < 1$ for all $j \in \{1, ..., M_x\}$, a sufficient CFL condition for the average total densities ${\bar{n}}^{k}_j = {\bar{n}}^{k}_{1, j}+{\bar{n}}^{k}_{2, j}$, to be below 1 for all $k \in [0, N_{it}]$ is

Proof. Assume that at a given time $t^k$, ${\bar{n}}^{k}_{1, j} \leq 1$ and ${\bar{n}}^{k}_{2, j}\leq 1$ for all $j \in \{1, ..., M_x\}$. The densities ${\bar{n}}^{k+1}_{\beta, j}, \ \beta = \{1, 2\}$, at the following time step are given by

Adding the last inequality for the two densities we get that

Then a sufficient condition for the total density to stay below 1 is (4.10).

Remark 3. Eqs. (4.8) and (4.10) give conditions enabling us to choose the time step in the numerical simulations. However, in some computations we need to reduce the time step in order to avoid oscillations in the pressure when the density is close to the singular density $n = 1$. These oscillations originate from small oscillation in the densities. However despite these varations the total density still verifies the condition $0\leq n < 1$. To enforce stability in the simulations, $\Delta x$ is replaced by $\Delta x^2$ in (4.8)-(4.10) which give rise to more stringent stability constraints. This kind of condition is in fact expected for a standard parabolic CFL.

Remark 4. Similar theorems can also be proven if the second-order upwind spatial scheme from ^[41] 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., ^[42,43].

5.   Conclusion

In this paper, we have presented a continuum 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.

Through a combination of analytical and numerical arguments, this paper demonstrates that the continuum 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 continuum 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 continuum model to the Hele-Shaw one seems within reach in the case $\alpha = 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.

Acknowledgements

The work of AC was supported in part by NSF Grant DMS-1521051, DMS-1818684 and RNMS Grant DMS-1107444 (KI-Net). SH and JPV acknowledge support from the Francis Crick Institute which receives its core funding from Cancer Research UK (FC001204), the UK Medical Research Council (FC001204), and the Wellcome Trust (FC001204). PD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1, EP/N014529/1 and EP/P013651/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. PD and SH would like to thank Andrea Bertozzi for stimulating discussions.

Data Statement

No new data were collected in the course of this research.

Conflict of interest

All authors declare no conflicts of interest in this paper.

AMS subject classifications

35K55; 35R35; 65M08; 92C15; 92C10

Appendix

A.   Derivation of the continuum model

Eqs. (2.3) and (2.4) can be derived from the gradient flow associated with the free energy (2.1). Indeed the functional derivatives $\frac{\delta \mathcal{E}}{\delta n_1}$ and $\frac{\delta \mathcal{E}}{\delta n_2}$ of $\mathcal{E}$ with respect to $n_1$ and $n_2$, acting on a density increment $\tilde{n}_1(x)$ and $\tilde{n}_2(x)$ are given by

Indeed, the computation for the first derivative is given by

The expression of $\frac{\delta \mathcal{E}}{\delta n_2}$ follows from a similar computation. Therefore the gradient flow associated to the free energy (2.1) according to the Wasserstein metric can be written as follows ^[1],

The metric used here is not the traditional distance on $L^2(\mathbb{R}^d)$, but the Wasserstein distance. This distance is defined on the set of probability distributions on $\mathbb{R}^d$. It is used in a wide variety of physically meaningful equations like the porous medium equation or degenerate parabolic equations [2,3,1]. We recover Eqs. (2.3) and (2.4) of the SPM by adding growth terms $G_1$ and $G_2$ to (A.3) and (A.4).

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. (A.3) and (A.4), 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 on three terms, representing the levels of volume exclusion, segregation and cohesion. Therefore the model seeks to minimize the pressures, which means moving away from situations where the volume exclusion or the segregation constraints are violated.