Incompressible limit for a two-species tumour model with coupling through Brinkman's law in one dimension

We present a two-species model with applications in tumour modelling. The main novelty is the coupling of both species through the so-called Brinkman law which is typically used in the context of visco-elastic media, where the velocity field is linked to the total population pressure via an elliptic equation. The same model for only one species has been studied by Perthame and Vauchelet in the past. The first part of this paper is dedicated to establishing existence of solutions to the problem, while the second part deals with the incompressible limit as the stiffness of the pressure law tends to infinity. Here we present a novel approach in one spatial dimension that differs from the kinetic reformulation used in the aforementioned study and, instead, relies on uniform BV-estimates.


INTRODUCTION
In recent years there has been an increasing interest in multi-phase models applied to tumour growth. Traditionally, tumour growth was modelled using a single equation describing the evolution of the abnormal cell density. This paper is dedicated to studying the two-species where n (i) represents the normal (resp. abnormal) cells, for i = 1, 2, and k ∈ N is a given constant modelling the stiffness of the total population pressure, p k , which is generated by both species, i.e., In addition, ν > 0 is a fixed positive constant that is understood as a measure of viscosity. The elliptic equation linking the macroscopic velocity, W k , with the pressure p k is typically referred to as Brinkman's law, for instance cf. [1]. The growth of the two densities is assumed to be modulated by two functions G (i) , for i = 1, 2, that are assumed to be decreasing in their variable, p k , similar to [7,20]. Throughout, we shall use the shorthand notation n k := n (1) k + n (2) k , in order to denote the total population. Upon adding up the two equations for the individual species, we obtain an equation for the total population density, n k , i.e., ∂n k ∂t − ∇ · (n k ∇W k ) = n k r k G (1) (p k ) + (1 − r k )G (2) where r k is the population fraction r k := n (1) k /n k . Related models have been extensively studied in the past. We refer to [17,19], and references therein, for a treatise of the incompressible limit for a single-species visco-elastic tumour model. As above, the velocity field is given by an elliptic equation involving the pressure that, in their case, is just given by a power of the sole species. Introducing the coupling of the two equations for the individual species drastically changes the behaviour and the same tool employed in [19] cannot be applied, at least not in a straightforward manner, and a different strategy has to be found. Even in the case ν = 0 corresponding to the inviscid case, the system nature of the problem gives rise to a whole range of difficulties, cf. [6,8,13]. At first glance, the pressure gains in regularity, however, it gains just enough regularity to obtain compactness of its gradient, requiring a minute derivation of suitable estimates. Let us stress that the same type of difficulties are also encountered when the pressure is not given as a power law, cf. [11,12,14]. A key tool in obtaining existence results and stable (with respect to the parameter k) estimates is to devise and manipulate the equation satisfied by the (joint) population pressure, cf. [6, 8, 11-14, 16, 18, 19]. In this work we shall follow this path. An easy application of the chain rule in conjunction with Eq. (1) leads to where the population fraction r k satisfies The change to these new variables was first introduced in [2][3][4] in the context of a two-species system where the two species avoid overcrowding. In a way, their works paved the way for more modern approaches to tumour models linked through Darcy's law, cf. [5,6,8,13].
The rest of this paper is organised as follows. In the subsequent section we set up precisely the problem and state our assumptions. In Section 3 we establish existence of solutions to the main system under consideration, Eq. (1), and discuss their regularity necessary for our purposes. Section 4 is dedicated to establishing a range of a priori estimates necessary in the analysis of the incompressible limit. Section 5 is devoted to establishing the strong compactness of the pressure, which is key in passing to the stiff limit. Finally, with all information at hand, we pass to the incompressible limit in the pressure equation and derive the so-called complementarity relation in Section 6. We round off the analytical results in Section 7 by presenting some numerical simulations for different parameter choices.

PRELIMINARIES AND STATEMENT OF THE MAIN RESULTS
We study the system posed on the whole domain R. It is coupled through the Brinkman law The system is equipped with non-negative initial data for any integer k ≥ 2. Moreover, we assume that there exists a constant, C > 0, such that for i = 1, 2, and every k ≥ 2. As before, the pressure is given in form of a power of the joint population, i.e., Recall that the pressure satisfies with the population fraction, r k := n (1) k /n k , given by Throughout the paper we assume the following regularity and properties of the growth functions G (i) , i = 1, 2, p denotes the derivative of the function G (i) . The pressure p M is often called the homeostatic pressure.
Recall that a solution W k to Brinkman's equation −ν∂ 2 x W k + W k = p k can be written as W k = K p k , where K is the fundamental solution to the equation −ν∂ 2 x K + K = δ 0 , i.e., Then K ≥ 0, K(x) dx = 1 and K, ∂ x K ∈ L q (R) for 1 ≤ q ≤ ∞. By the elliptic regularity theory we have W k (t, ·) ∈ W 2,q (R), for any t ∈ [0, T ], 1 ≤ q ≤ ∞.
Below we formulate the main results of this work.

Theorem 2.1 (Existence of Solutions).
For any initial data satisfying (2), system (1) admits a solution n . We highlight the fact that solutions are essentially bounded since these bounds are not a consequence of the BV-bounds. Rather, they are obtained independently. This may prove useful for an extension to higher dimensions in future works. Theorem 2.2 (Incompressible Limit and Complementarity Relation). We may pass to the limit k → ∞ in the pressure equation, Eq. (2). This yields the so-called complementarity relation in the distributional sense, where n (i) Moreover, the following holds true The subsequent sections are concerned with the proof of the two main theorems.

EXISTENCE OF SOLUTIONS AND REGULARITY
This section is dedicated to proving the existence of solutions to the (p, r)-system. The proof is based on an application of Banach's fixed point theorem. Let k ≥ 2 be fixed throughout this section. Further, assume for now, that the initial data u (i) 0 are Lipschitz continuous. For given functions p, r ∈ L ∞ (0, T ; L ∞ (R)) we construct solutions u (i) to the linearised system, i = 1, 2, ∂u ∂t where f 1 (p, r) = K p − p + νrG (1) For the fixed p from above, we may construct the backward flow We readily observe that , r(τ, x)) dτ, i = 1, 2, solve the linearised system (3). Now, considering another element (p,r) in L ∞ (0, T ; L ∞ (R)), we observe that Thus, upon passing to the supremum, we obtain the following stability estimate for two solutions In particular, for T 1 > 0 small enough the estimate gives rise to a contraction in the Banach space L ∞ (0, T 1 ; L ∞ (R)), which is sufficient to infer the existence of a unique fixed point, by an application of Banach's fixed point theorem. Since the supremum norm of the solution does not blow up, a finite number of iterations of the above argument leads to existence of solutions for all times T > 0.
For the subsequent analysis, let us call this fixed point (u . It remains to prove the expected BV-regularity of solutions. This is an easy consequence of the "transport nature" of the system, i.e., Multiplying by sign(u (i) * ) and adding the two equations, for i = 1, 2, we obtain, after integrating where the constant C > 0 depends only on the Lipschitz constants of the functions f (i) and the L ∞ -bounds on the fixed point. In particular, from Gronwall's inequality we deduce a control on the BV-seminorm and, more importantly, the existence of solutions even in cases where u (i) 0 is not Lipschitz continuous but only of bounded variation. Using the fact that the existence result transfers to the original system for n (i) Remark 3.1 (Extension to Higher Dimensions). Let us remark here that the same strategy can be easily extended to higher dimensions since the transport nature is the same in any dimension. In fact, the only "problematic" point in our strategy is the contraction argument which depends on ∂ x K L 1 . However, this norm is finite in any dimension, and therefore our existence result holds in any dimension.

A PRIORI ESTIMATES
In this section we derive some bounds for the main quantities of interests, uniformly in k. These will be vital when passing to the limit with k → ∞.

Lemma 4.1 (A priori estimates I).
The following hold uniformly in k for any T > 0.
Proof. Clearly when n k (t = 0) ≥ 0, then n k stays non-negative at all times. Integrating Eq. (1) in space and time we deduce that n k ∈ L ∞ (0, T ; L 1 (R)) uniformly in k. By the maximum principle we have the bound 0 ≤ p k ≤ p M . Then using n k p 1 k−1 we deduce n k ∈ L ∞ (0, T ; L ∞ (R)) uniformly. Writing p k ≤ n k n k k−2 ∞ we see that p k ∈ L ∞ (0, T ; L 1 (R)). Finally, we use that n (1) k = r k n k and 0 ≤ r k ≤ 1 to deduce the last bounds.
Using the above Lemma and the boundedness of W k , we have the following result.
Proof. Here and henceforth we shall employ the notation The supremum is taken only up to p M , because in principle the functions G (i) can decrease arbitrarily. The uniform bound obtained in the previous proof shows however, that only the range 0 ≤ p k ≤ p M is relevant.
Using the equation for the population fraction and boundedness of the growth functions having used the a priori bounds on the pressure, p k .
The following lemma establishes an L 1 -bound on the right-hand side of the pressure equation.

Lemma 4.3 (A priori estimates II).
The following estimate holds for any T > 0 for a constant C(T ) > 0, independent of k. Furthermore, the following bounds hold uniformly in k Proof. Let us introduce the following notation and follow the strategy of [19]. Using and consequently, Integrating in space-time and using the assumption that |G The first three terms on the right-hand side are controlled uniformly, as is the very last term. For the remaining two terms we write which, for k large enough, is controlled by the left-hand side of the last inequality, and Using Lemma 4.1, we see that the right-hand side is uniformly bounded in L ∞ (0, T ; L q (R)), 1 ≤ q ≤ ∞. It follows that as desired. Now, using Eq. (4) and the above computations, it is clear that ∂ t W k is uniformly bounded in and use the definition of K, cf. Eq. (2), to conclude the proof.

Remark 4.4.
All the results of this section remain valid in any spatial dimension d ≥ 1, see for example [19] for the a priori estimates, and the L 1 -bound on the quantity kp k Q k .

STRONG COMPACTNESS OF THE PRESSURE
This section is solely dedicated to the derivation suitable estimates in order to obtain strong compactness of the pressure, p k . A key step in this pursuit is the following BV-estimate on the individual species as well as the total population. Proof. For i = 1, 2, we consider Upon differentiating in space, we obtain for i = 1, 2. Upon adding up both equations we get ∂ ∂t Multiplying the equation for the individual species by σ (i) := sign(n (i) k ) and the equation for the total population by σ := sign(n k ), we get, upon adding the three equations (5), (5) and integrating in space First we notice that the exact derivatives vanish and the estimate simplifies to Next, we notice that the all the terms involving the pressure gradient cancel due to opposite signs, whence Thus we are left with Using the fact that the integrand of the last line of Eq. (5) may be simplified to Using the fact that |σ (i) |, |σ| ≤ 1 and exploiting the bounds we may bound the terms of Eq. (5), and the last line of Eq. (5) becomes The last integral in Eq. (5) vanishes due to the fact that n k = n (1) k +n (2) k . Thus, substituting Eq. (5) into Eq. (5), an application of Gronwall's lemma yields the BV-estimate in space.

Corollary 5.2. From the proof of the preceding lemma we deduce
where C > 0 is independent of k.
Proof. Let us revisit the equation for ∂ t n k , i.e., Now we use the bounds G (i) p ≤ −α < 0, for i = 1, 2, and integrate in time to see that Thus we infer that n k ∂ x p k is uniformly bounded in L 1 (0, T ; L 1 (R)).

Lemma 5.3 (Strong Compactness of the Pressure). There exists a function
such that there holds up to a subsequence, as k → ∞ in any L p loc (0, T ; L q (R)), for 2 ≤ p, q < ∞. In addition, the convergence also holds in the pointwise almost everywhere sense.
Proof. Let us write the quantity n k |∂ x p k | as a spatial derivative of a non-decreasing function of the pressure. We compute as follows i.e., ∂ x φ k (p k ) ∈ L 1 (0, T ; L 1 (R)), uniformly in k. Moreover, we have the same L 1 -bound for the time derivative of φ k (p k ). Indeed where we have used Lemma 4.3.
We conclude that the sequence (φ k (p k )) k converges strongly in L 2 ((0, T ) × R). On the other hand, from the uniform bounds on p k we infer that p k p ∞ , weakly in L 2 ((0, T ) × R), up to the subsequence. We can therefore apply Lemma 8.1 to conclude that φ k (p k ) → p ∞ , strongly in L 2 loc ((0, T )×R). We claim that this in fact implies strong convergence of the sequence of pressures (p k ) k itself. Indeed, using the triangle inequality yields with the right-hand side of the last line converging to zero. We conclude that p k → p ∞ , strongly in L 2 loc ((0, T ) × R). In combination with the L ∞ -bounds, we deduce that this convergence holds strongly in L p loc (0, T ; L q (R)), for any 2 ≤ p, q < ∞, using the dominated convergence theorem. Moreover, the convergence is also true pointwise almost everywhere.

INCOMPRESSIBLE LIMIT AND COMPLEMENTARITY RELATION
We have garnered all information necessary to pass to the incompressible limit in the pressure equation (2) and prove Theorem 2.2.
Proof of Theorem 2.2. Having established strong convergence of the sequence (p k ) k , and weak convergence of (n k ) k due to the a priori estimates, we can pass to the limit in the relation (14) n to deduce the relation (1−n ∞ )p ∞ = 0, almost everywhere. For a test function ϕ ∈ C 1 c ((0, T )×R), let us recall the weak formulation of the equation for the pressure Due to the uniform bounds on the right-hand side, cf. Lemma 4.3, we may divide by k − 1 to obtain Note that, writing n (1) k = n k r k and expressing n k in terms of p k , in a fashion similar to Eq. (6), we may readily pass to the limit in all of these terms due to the strong convergence of the pressure and the a priori bounds of Lemma 4.1. We thus obtain for i = 1, 2. Indeed these equations follow by passing to the limit in the weak formulation of (1) where ϕ ∈ C 1 c ((0, T ) × R). Remark 6.1. In fact, using the strategy of the previous section, i.e., the BV-bounds in space, in conjunction with a control on the time derivative obtained from bounding the right-hand side of the equation for the individual species, one can deduce also strong convergence of the sequence (n k ) k . As a consequence, the limit functions n ∞ , n

NUMERICAL INVESTIGATIONS
In this section, we revisit the results from the preceding sections and showcase certain properties of the system. The numerical simulations are performed using the positivity-preserving upwind finite volume scheme proposed for a system of two interacting species in [9,10] where the reaction terms are computed on each finite volume cell as simple ODEs. The implementation hinges on the fact that the elliptic Brinkman law (1) can be solved using the integral representation (2). Figure 1 displays the role of the viscosity parameter, ν. The same initial data n (1) k,0 (x) = m(x − 4.5)(6.5 − x), and n (2) k,0 (x) = m(x − 8.5)(10.5 − x), are used in both cases and m > 0 is chosen to normalise the initial mass to 1. In both cases we used k = 100, as we are interested in the incompressible regime. In addition, we chose G (i) (p) = 1 − p, for i = 1, 2, corresponding to a homeostatic pressure of p M = 1, cf. Eq. (2). In  both cases we observe the propagation of segregation in agreement with Lemma 4.2. Moreover, we observe a drastic drop in the pressure in Figure 1(b). This was already observed in the one species case, cf. [19], where the fact was exploited that the limiting pressure has an integral representation formula. In stark contrast, Figure 1(d) shows an almost smooth transition of the pressure indicating a much higher regularity. This is in perfect alignment with the findings of [6], as the case ν = 0 yields, at least formally, the system studied in the latter. As a matter of fact, the pressure gradient was shown to be square-integrable in the Darcy case, i.e., ν = 0. We conclude, by remarking that the front propagation is much faster in the regime of small ν, another fact that was already observed in the single-species case.
In Figure 2 we present the effect of different growth terms of the tumour cells and healthy tissue.
To be more precise, we choose the same initial condition as above but use as growth terms for the two species. We see that the first species, n k , proliferates much faster compared to the second one. More interestingly, we see that the pressure not only has a jump at the boundaries of the support of the total population, but also at the internal layer.  We run the simulation for the same initial data for two different growth functions, G (i) (p). In both cases, we chose k = 100 and the individual species are represented by solid lines in red and blue, the pressure is superimposed as a black dotted line, as before. The pressure drops not only at the boundary of its support. We also observe jumps in internal layers. Figure 3 shows the evolution of system (1) for initial data representing a regime where healthy tissue has already been intruded by cancerous cells, i.e., n k,0 (x) = m(x − 6.5)(8.5 − x), and n (2) k,0 (x) = m(x − 6)(9 − x), where, again, m > 0 normalises the mass. In addition, we choose the same unequal growth functions, G (i) , as before, thus promoting the tumour growth compared to the normal tissue.

CONCLUSIONS
The goal of the paper was twofold. We extended an established tumour growth model to an interaction system of two cell populations, i.e., normal and abnormal cells. The interaction is given through the Brinkman flow, an elliptic equation that yields the velocity fields for each cell population. In the first part of this paper we proved the existence of solutions to the interaction system, cf. Theorem 2.1. Building upon this result, we passed to the "incompressible" limit in the pressure equation, Eq. (2), and obtained the limiting equation, also referred to as complementarity relation, cf. Theorem 2.2. This way we were able to derive a geometric model from the cell-density model we presented. Note that both the existence result and the incompressible limit rely on strategies different from the ones adapted for related models (either in the parabolic two-species case when Brinkman's law is replaced by Darcy's law (ν = 0) or the one-species model with Brinkman flow). The  results are complemented with a numerical investigation showcasing the segregation result, the discontinuities in the pressure and the two individual population densities which is why we do not expect better regularity than bounded variation. In summary, this paper extends known results in the literature to two species. As far as the existence of solutions is concerned, no additional difficulties are expected in the multi-dimensional case. However, when it comes to the stiff limit not only our method fails but also the kinetic reformulation that was employed in the one-species case, cf. [19], would need a serious makeover that is, at this stage, far from clear -even in one dimension. New singularities appear at internal layers when the two species meet and it appears different tools are required, such as the extension of the kinetic reformulation to systems, which, to our knowledge, does not exist.
The exploration of such a technique is left for future works. In addition, the rigorous inviscid limit, ν → 0, remains an open question that is left for future work.

APPENDIX
For the readers' convenience we shall recall here the compactness method invoked in [15] in the context of the fast reaction limit in a cross-diffusion system with growth and death processes. Roughly speaking, it allows to identify the limit of the composition of a uniformly compact nonlinear function and a weakly convergent sequence.

Lemma 8.1 ("Lemma A").
Let Ω ⊂ R be a compact domain and set Q T := (0, T ) × Ω. Furthermore, let {u n } ⊂ L ∞ (Q T ) and {f n } ⊂ C(R) be sequences with the properties (i) u n u, weakly in L 2 (Q T ), (ii) f n is nondecreasing, (iii) f n → f , uniformly on compact subsets of R, and (iv) f n (u n ) → χ, strongly in L 2 (Q T ). Then For the sake of exposition, we recall here that the assumptions of the above lemma are indeed met in our case.
Remark 8.2 (The assumptions are met). The first assumption is the easiest to check as it follows directly from the uniform L ∞ -bounds on the pressure. Similarly, it is readily verified that each element of the sequence of functions, in our context given by φ k (x) = x k/k−1 , is indeed nondecreasing. Moreover, the uniform convergence towards the identity is straightforward. Thus the only requirement that needs a more minute argument is (iv) which we present in the first part of the proof of Lemma 5.3.