## Abstract

Various models of tumour growth are available in the literature. The first type describe the evolution of the cell number density when considered as a continuous visco-elastic material with growth. The second type describe the tumour as a set, and rules for the free boundary are given related to the classical Hele-Shaw model of fluid dynamics. Following previous papers where the material is described by a purely elastic material, or when active cell motion is included, we make the link between the two types of description considering the ‘stiff pressure law’ limit. Even though viscosity is a regularizing effect, new mathematical difficulties arise in the visco-elastic case because estimates on the pressure field are weaker and do not immediately imply compactness. For instance, travelling wave solutions and numerical simulations show that the pressure is discontinuous in space, which is not the case for an elastic material.

## 1. The visco-elastic cell model

From the mechanical viewpoint, at present living tissues are considered as fluids [1]. For tumour development, cell division is another major aspect of the description. Here, we consider models of tumour growth describing the dynamics under mechanical pressure and cell growth. In the simplest description, mechanical pressure is incorporated in the fluid velocity through Darcy's law as in [2–9], and the tissue is considered to be a porous medium. Then, only friction with the cell surroundings (the extracellular matrix for instance) is considered. In this paper, we propose to take into account viscosity, which is a way to also represent friction between cells themselves, considering the tissue as a Newtonian fluid. Then Brinkman's law is used in place of Darcy's law and has been derived rigorously for inhomogeneous materials [10]. Such visco-elastic models for tumour growth, based on Stokes' law or Brinkman's law, have also been used in the context of tumour growth in [11–13]. Besides the fluid velocity field, another concern is tissue compressibility. In most of the models mentioned previously, incompressibility is assumed, leading to a free boundary problem. However, Stokes' law or Brinkman's law are also used by considering the tissue as ‘compressible’ [14,15]. Note that the theory of mixtures allows for a general formalism containing both Darcy's law and Brinkman's law [16–18].

More precisely, we denote the number density of tumour cells by *n*(*x*,*t*) and the pressure by *p*(*x*,*t*), which is given by a pressure law *p*=*Π*(*n*), where *Π* is an increasing function. We choose the following ‘compressible’ pressure law, parametrized by *k*:
1.1
The ‘incompressible’ limit is obtained by letting the parameter *k* go to . In fact, we formally deduce from (1.1) that the limiting graph is given by for *n*∈[0,1) and for *n*=1. We assume a Brinkman flow, which means that the macroscopic velocity field is given by ∇*W* for a potential *W* closely related to the pressure. With these assumptions, the model, indexed by *k*, for tumour growth is
1.2
and
1.3
The second term of equation (1.2) models the influence of the mechanical pressure. Following [13,19], we assume that growth is directly related to the pressure through a function *G*(⋅) satisfying
1.4
The pressure *P*_{M} is sometimes called the *homeostatic pressure*. We complete equations (1.2) and (1.3) with a family of initial data satisfying (for some constant, *C* independent of *k*)
1.5
The viscosity coefficient, *ν*>0, is supposed to be constant; when viscosity is neglected, that is, equation (1.3) has *ν*=0, we recover Darcy's law. To use the Laplacian in (1.3), rather than the Stokes viscosity terms, is to simplify the mathematical presentation. The case with the Stokes viscosity terms is left as an open question here. Existence of solutions to system (1.1) and (1.2) enters into the class of transport equations coupled to elliptic equations such as the aggregation equation studied, for example, in [20]. However, in this work, we are in the repulsive framework, avoiding blow-up, and therefore consider bounded solutions such as in [21].

The aim of this work is to explain the derivation of ‘incompressible’ models from the ‘compressible’ equations for a mechanical model of tumour growth where the tumour is considered as a visco-elastic medium. This effect arises for the ‘stiff pressure law’ limit of this ‘compressible’ model towards a free boundary model which generalizes the classical Hele-Shaw equation, a limit which has been studied recently in the inviscid case in [8] and extended to the case where cells can undergo active motion in [9]. However, a major mathematical difference occurs in the visco-elastic case; the type of equation changes drastically from parabolic to transport and most of the functional analytic arguments are completely different. A major consequence is that the limiting pressure is discontinuous on the free boundary.

We first explain formally what can be expected. The limit strongly uses the equation satisfied by the pressure. Multiplying equation (1.2) by *Π*_{k}′(*n*) and using the chain rule, we deduce
From our choice for the law of state (1.1), we deduce that *nΠ*_{k}′(*n*)=*kn*^{k−1}=(*k*−1)*Π*_{k}(*n*). Injecting this expression into the above equation, we deduce that
1.6
where we have defined the function *H*, as well as some of its properties, as
1.7
Indeed, *G* is non-increasing and thus (*I*−*νG*) is invertible on [0,*P*_{M}] onto [−*νG*(0),*P*_{M}]. Furthermore, note that (*I*−*νG*)′>1.

Back to the limit : at least when *n*_{k} and *p*_{k} converge strongly, from (1.1), we deduce that , then we find formally the relation
1.8
Letting and assuming we can pass to the limit in all terms of (1.6), we formally deduce

Therefore, at the limit we can distinguish between two different regions. The first region is defined by the set
1.9
on which we have the system
1.10
and
1.11
Thus, the latter system reduces to, for *x*∈*Ω*(*t*),

On the second region, , the limiting system becomes

To establish this limit rigorously, we need some additional assumptions on the initial data. Namely, we need that the family is ‘well prepared’. By this, we mean that, for some open set *Ω*^{0},
1.12
With the notation in (1.6), these assumptions imply that and as , which is needed to identify the kinetic defect measure in our proof (see (2.26)). Regarding the first assumption, which prevents an initial layer, it is needed because regularizing effects are not known for transport equations, unlike for porous media equations, when Darcy's law is considered, which yield lower estimates on ∂_{t}*n*_{k}; it is a major difficulty here to prove time compactness. The latter assumption can however be slightly relaxed to for all *A*>0 in . With our present proof, we need to avoid the existence of a domain where remains strictly between 0 and 1, a case which we leave open at this stage and when our conclusion in equation (1.13) is certainly wrong.

Our goal is to prove the following.

### Theorem 1.1

*Under assumptions* (*1.4*), (*1.5*) *and* (*1.12*), *consider a solution of the system* (*1.2*) *and* (*1.1*). *After extraction of subsequences, both the density n*_{k} *and the pressure p*_{k} *converge strongly in* *for all T*>0, *as* *towards, respectively,* *and* *belonging to* *; up to a subsequence, W*_{k} *converges strongly in* *for all q*≥1, *towards* *. Moreover, these functions satisfy*
1.13
1.14
1.15
1.16

The first relation in (1.15) is equivalent to the statement (1.16) and replaces the usual ‘complementary relation’ in Hele-Shaw flow, [8,9,22].

Because the function *H*(⋅) does not vanish, we conclude from the first relation in (1.15) that is discontinuous. This is a major difference from elastic materials (Darcy's law), where is continuous in space, and this is illustrated by the travelling wave solutions we build in §3. The pressure jump is however related to the potential , which is different from models that include surface tension, where the jump is related to the free boundary curvature (see [23,24] and references therein).

We first prove theorem 1.1 in several steps. In a first step, we derive *a priori* estimates. Because they do not give compactness for the pressure, we analyse possible oscillations using a kinetic formulation [25]. From the properties of solutions of the corresponding kinetic equation, we conclude that strong compactness occurs. All these steps are performed in §2. The one-dimensional travelling wave profiles are presented in §3 with numerical illustrations. The final section is devoted to the conclusion and some perspectives.

## 2. Proof of the Hele-Shaw limit

We divide the proof of our main result theorem 1.1 into several steps. We begin with several bounds which are useful for the sequel. Then, in order to prove strong convergence of the pressure *p*_{k}, we analyse possible oscillations using the kinetic formulation of (1.6) in the spirit of Perthame & Dalibard [21].

### (a) Estimates

### Lemma 2.1 (*A priori* estimates)

*Under the previous assumptions, for all T*>0, *the following uniform bounds with respect to k hold*:
*For some non-negative constant C*(*T*), *independent of k, we have*
2.1

We can draw several conclusions about this lemma. First, after extracting subsequences, it is immediately obvious that the following convergences hold as :
and these limits belong to for all *t*>0. Also, we have
Passing to the limit in (1.3), we get
2.2

The second consequence concerns the backward flow with velocity ∇*W*_{k} defined as
2.3
as well as the forward flow
2.4
even though ∇*W*_{k} is not uniformly Lipschitz continuous but slightly less, and according to DiPerna–Lions theory [26] these flows are well defined a.e. and, after extraction of subsequences as in lemma 2.1, it converges a.e. to the limiting flows defined by (2.24) for the backward flow and by (2.11) for the forward flow.

The third conclusion uses a combination of the above forward flow with equation (1.6). We have 2.5

### Proof.

*First step. A priori bounds in .* Clearly, *n*_{k} is non-negative provided *n*_{k}(*t*=0)≥0. Integrating equation (1.2), we deduce a bound for *n*_{k} in , uniformly with respect to *k*.

By definition of *p*_{k} in (1.1), we clearly have that *Π*_{k}′(*n*_{k})≥0, when *k*>1. We can apply the maximum principle of [27], lemma 2.1 to obtain the uniform bound
Therefore, still using relation (1.1), we have *n*_{k}=(((*k*−1)/*k*)*p*_{k})^{1/(k−1)} and *n*_{k} is uniformly bounded in . Then, writing , we deduce a uniform bound of (*p*_{k})_{k} in .

*Second step. Representation of W_{k}.* Using elliptic regularity on (1.3), we conclude that, for all

*t*∈[0,

*T*],

*W*

_{k}(

*t*,⋅) is bounded in , . Moreover, denoting by

*K*the fundamental solution of −

*ν*Δ

*K*+

*K*=

*δ*

_{0}, we have 2.6 We recall that and that

*K*≥0, , which we use below.

Taking the convolution of (1.6), we deduce 2.7

*Third step. Bounds on* *Q*_{k}. Then, by definition of *Q*_{k} and using (1.6), we compute
Therefore, from a standard computation, we deduce
We may integrate in *x* and *t*. Because *p*_{k} and *W*_{k} are uniformly bounded in , and |*G*′|≥*α* from (1.4), we find
The three first terms on the right-hand side are all controlled uniformly and, to conclude the bound (2.1), we have to estimate the last two terms. Using (1.3), the first term is
and this term is controlled, for *k* large enough, by the *αk* term on the left-hand side. The second term is
Using the uniform bounds on *p*_{k}, we have that *p*_{k}∂_{i}*W*_{k}=*p*_{k}∂_{i}*K*★*p*_{k} is uniformly bounded, with respect to *k*, in , , and, thus, is also uniformly bounded in , . Finally, *p*_{k}Δ*W*_{k} is also uniformly bounded in , . This immediately concludes the proof of estimates (2.1).

*Fourth step. Estimate on* ∂_{t}*W*_{k}. Finally, using the above estimates and equation (2.7), we deduce that ∂_{t}*W*_{k} is uniformly bounded with respect to *k* in , .

For the estimate for ∂_{t}∇*W*_{k}, we can again use the above calculation and write, for *i*=1,…,*d*,
As ∂_{ij}*K* is a bounded operator in *L*^{1}, we conclude the last bound in lemma 2.1. ▪

### (b) Which oscillations for the pressure?

We deduce from lemma 2.1 that, up to a subsequence, the sequence (*W*_{k})_{k} converges strongly in . However, we only get weak-★ convergence for the pressure (*p*_{k})_{k} and the density (*n*_{k})_{k}. Here, we give an argument showing that the only obstruction to strong compactness is oscillations of *p*_{k} between the values *p*_{k}≈0 and .

### Lemma 2.2

*Let T*>0 *and let H be defined in* (*1.7*) *with the assumptions* (*1.4*). *Consider real numbers β*_{1}>0, *β*_{2}>0 *small enough, and let p*_{k} *be as in lemma 2.1, then we have*
and
*where meas denotes the Lebesgue measure*.

### Proof.

Let 0<*β*_{1}<*β*_{2}<*p*_{m}, *p*_{m} being defined in (1.7); therefore, we have for all
2.8
From assumption (1.4), the function *I*−*νG* is increasing and, by definition (1.7), (the non-negativity is because is a solution of (2.2)). Therefore, on the set , we have, for some *ω*(*β*_{2})>0,
and
Thus, we can estimate

Hence, using estimate (2.1), and the strong convergence of *W*_{k}, we deduce that
2.9
Thus, estimates (2.8) and (2.9) prove the first statement of lemma 2.2.

The second statement can be proved in the same way. ▪

### Remark 2.3

We note, for future use, that, in the same spirit, we also have that 2.10

### (c) Strong convergence of the pressure

However, we need strong convergence to recover the asymptotic limit, in particular the equation satisfied by . A difficulty here is that we do not have estimates on the derivatives on *p*_{k}, unlike in [8,9]. Therefore, we develop another strategy, based on estimate (2.1), to obtain the following strong convergence result.

### Lemma 2.4 (Strong convergence of *p*_{k})

*Up to a subsequence, p*_{k} *converges strongly locally in* *towards* . *Moreover*, *a.e*.

*Furthermore, we have that*
*is the image of Ω*^{0} *by the limiting flow Y* _{(x)}(*t*), *defined by*
2.11

*Finally, we have, for all T*>0,
2.12

### Proof.

The strategy is to pass to the limit in equation (1.6) for *p*_{k} and to combine this information with the possible oscillations of *p*_{k} as described by lemma 2.2. For that, we need a representation of the weak limit of *p*_{k}, which we can obtain thanks to a kinetic representation.

*First step. Representation of nonlinear weak limits.* Our first result is that there is a measurable function 0≤*f*(*x*,*t*)≤1 such that, for all smooth functions , we have, up to a subsequence,
2.13
and
2.14
Interpreted in terms of Young measures, this means that *p*_{k} oscillates between the values 0 and with the weights 1−*f*(*x*,*t*) and *f*(*x*,*t*). Note that, for *S*(*p*)=*p*, we find
2.15

To prove these results, we define
and we write
2.16
We can extract a subsequence, still denoted (*p*_{k})_{k}, such that **1**_{{0<ξ<pk}} converges in towards a function *χ*(*x*,*ξ*,*t*), which satisfies 0≤*χ*(*x*,*ξ*,*t*)≤1. Then, *S*(*p*_{k}) converges weakly to .

We define by *f* the weak-★ limit of **1**_{{pk(x,t)≥pm/2}},
where we recall that *p*_{m} is defined in (1.7). As , we may use lemma 2.2 to conclude (2.13) and (2.14).

*Second step. Equation satisfied by* *χ*_{k}. We use equation (1.6),
For any function , multiplying it by *S*′(*p*_{k}) leads to
Denoting by *δ* the Dirac mass, we can rewrite the latter equation as
2.17
and
2.18
Eliminating the test function *S*′(⋅), this is equivalent to writing
2.19

However, this formula is not enough to pass to the limit and we need the divergence form, Therefore, using (2.16) and the fact that , we have 2.20

Because , and integrating by parts, we have Therefore, (2.20) is equivalent to our final formulation 2.21

One can simplify this relation and write
Finally, (2.20) is equivalent to
In particular, integrating in *ξ* we recover the expected formula

*Third step. Equation satisfied by* *f*. We may pass to the limit in (2.21). For all *T*>0, the sequence *μ*_{k} is uniformly bounded in owing to estimate (2.1). Thus, we can extract a subsequence converging, in the weak sense of measures, towards a measure denoted *μ* in . Because *Q*_{k}(*x*,*ξ*,*t*)=*W*_{k}−*ξ*+*νG*(*ξ*) is positive for *ξ*≤*p*_{m}, we have

Therefore, passing the limit into (2.21), in the sense of distributions, This last equation can also be written with (2.14) and thus 2.22

Using assumption (1.12), this equation is complemented with the initial condition and

It is useful to keep in mind the equivalent form of this equation, and thus, using (2.15), 2.23

We can also integrate (2.22) and recover

*Fourth step. The set* {*g*(*x*,*t*)=1 and *ξ*<*p*_{m}}. It is useful to consider the function
with the characteristics defined by
2.24
This function *g* is the solution to the transport equation

Using (2.23) and 0≤*f*≤1, we find
2.25
From the comparison principle, we conclude that *f*(*x*,*t*)≥*g*(*x*,*t*) and that
2.26

*Fifth step. Strong convergence of* *p*_{k}. Another wording for step 4 is that
with *Y* _{(x)}(*t*) the limiting flow of defined in (2.4). Indeed, from (2.5) and the strong convergence of the flow, we infer that
Then, we have . We recall that, by definition, .

We show that it implies the strong convergence locally in of *p*_{k} towards . Let *U* be an open bounded subset of ; therefore, we have
2.27
with
For the first term *I*_{k}, we have that
Using estimate (2.1), we deduce that the first term on the right-hand side goes to 0 as . From the local strong convergence of *W*_{k} towards , the second term on the right-hand side converges to 0 too. We conclude that . Moreover, it has been proved in lemma 2.2, see equation (2.10), that . For the last term, we have, using the fact that is bounded in , that, for some non-negative constant *C*,
We have shown in the fourth step above that **1**_{{pk≥pm/2}} converges weakly towards . Then passing to the limit in the latter inequality, we deduce that . We conclude from (2.27) that, for any open bounded subset *U*,
By uniqueness of the weak limit, we deduce that a.e.

*Sixth step. Derivation of* (2.12). From definition (2.18), this limit is now a consequence of
But *μ*_{k} vanishes for because, from (2.25), we infer that *μ*=0 both when *f*=1 and when *f*=0. Therefore, we find (2.12). ▪

### (d) Proof of theorem 1.1

The proof of theorem 1.1 can now be easily deduced from lemma 2.4. First, up to a subsequence, we have that *p*_{k} converges a.e. towards . On the one hand, recalling that the sequence (*p*_{k}) is uniformly bounded in , we use the Lebesgue dominated convergence theorem to show that, for any bounded open *U*,
On the other hand, we have from estimate (2.1) that
We deduce that a.e., that is, (1.16).

We may apply the strong convergence for transport equations, as in [26,28], to conclude that, as the term *G*(*p*_{k}) converges strongly, *n*_{k}, which solves the transport equation (1.2), itself converges strongly. Note in particular that, from assumption (1.12), we have . Passing to the limit in equation (1.2), we recover the limit equation for , (1.13).

Finally, passing to the limit in the relation we deduce that . The relation (1.15) is then a direct consequence of lemma 2.4.

## 3. One-dimensional travelling waves

In order to exemplify theorem 1.1 and to give a simple case, with a solution that can be built analytically, we look for a one-dimensional travelling wave solution to the Hele-Shaw limit.

Because travelling waves are defined up to a translation, we may set, in the moving frame, . Then, the system rewrites as
3.1
and
3.2
Moreover, the jump condition at the interface *x*=0 implies −*σ*[*n*]−[*nW*′]=0, which leads to the travelling velocity *σ*=−*W*′(0). We denote *W*_{0}:=*W*(0). For *x*>0, we have
3.3
from which we deduce that . Then we can rewrite the first equation in (3.1) as
Taking the limit *x*→0 leads to *n*(0)=0. Moreover, as *n*′≤0, we deduce that *n*=0 on .

For *x*<0, we solve the second-order ODE for *W* (3.2) with the boundary condition given by the continuity at the interface: *W*(0)=*W*_{0} and . As an example, we choose for the growth term the function
3.4
Then, replacing *H*(*W*) by its expression in equation (3.2), we deduce that *W* satisfies
The only bounded solution on such that *W*(0)=*W*_{0} is given by
3.5
Moreover, from (3.3), we have that . With (3.5) we also have that . From the continuity of *W*′ at the point *x*=0, we deduce the value for *W*_{0}
Then, we conclude that, for *x*<0,
The pressure is then given by
and the travelling velocity by
We note that the pressure is non-negative and has a jump at the interface *x*=0. The height of the jump is given by . We observe moreover that *σ* is a decreasing function of *ν*. Letting *ν*→0, we recover the result for the Hele-Shaw model for purely elastic tumours [27,29].

*Numerical simulations.* Finally, we present numerical simulations of the system (1.2) and (1.1) in one dimension. We use a discretization owing to a Cartesian grid of a bounded domain [−*L*,*L*] of the real line. Equation (1.2) is discretized by a finite volume upwind scheme. Equation (1.3) is discretized owing to a finite difference scheme. As we focus on the case where *k* is large, we use *k*=100 in the numerical computation. For the initial data, we choose *n*^{0}=**1**_{[−0.2,0.2]}. The growth function *G* is chosen as in (3.4) with *P*_{M}=1.

In figure 1, we display the shape of the density *n*, the pressure *p* and *W* obtained by the numerical simulation. Figure 1*a* displays the result with a viscosity coefficient *ν*=1. For comparison, we plot, in figure 1*b*, the shape in the case without viscosity (*ν*=0). Comparing figure 1*a* and 1*b*, we observe that, in the case *ν*=1, we have a jump in pressure at the interface of the solid tumour, whereas in the case *ν*=0 the pressure is continuous at the interface.

We display in figure 2 the first steps in the formation of the propagating front with the initial data *n*^{0}=**1**_{[−0.2,0.2]}. For this simulation, we take *ν*=1 and *k*=100. The dynamics is represented at four successive times of the density *n*, pressure *p* and *W*. After a transitory regime during which the pressure increases until it reaches its maximal value *P*_{M}=1, the shape of the travelling waves is obtained and the front of the tumour invades the healthy tissue.

## 4. Conclusion and perspectives

A geometric model, also called incompressible, has been derived from a cell density model (also called compressible) when the pressure law is stiff. The limiting problem is a free boundary problem for the set *Ω*(*t*) of non-zero pressure. Because viscosity is considered here, the limiting system for the pressure consists in an algebraic relation between the pressure and the limiting velocity potential (1.16), coupled with an elliptic equation for the potential set in the whole space (1.14).

This is a major difference from the case where viscosity is neglected, the so-called Hele-Shaw system [4,7]; in the latter, the pressure is given by an elliptic equation for the pressure in the moving domain *Ω*(*t*). A paradox is that the effect of retaining the viscosity generates a jump in the pressure at the interface of the region defining the tumour, unlike in [8,9], where the Hele-Shaw problem is complemented with Dirichlet boundary conditions and therefore the pressure is continuous. This point is also observed in the numerical simulations in §3. The velocity of the propagating front of the tumour is given by the equation satisfied by the density (1.13). Because the pressure is discontinuous, it has weaker regularity than in the inviscid case treated in [8,9], and we need to develop a new strategy of proof to derive the incompressible limit. Our approach is based on a kinetic formulation of the equation satisfied by the pressure.

This work also opens several additional questions. First, the case with general initial data is not treated here because we assume that *n*^{0} vanishes outside *Ω*^{0}. Then, it would be interesting to consider the case with active motion as in [9]. In such a case, equation (1.2) is replaced by a parabolic equation. Then, the structure of the problem is different but the limiting system should be the same, except the equation for the density, which implies a faster propagation of the region *Ω*(*t*). Moreover, it is formally clear from (1.10) and (1.11) that, letting *ν*→0, we recover the Hele-Shaw system. However, a rigorous proof of this fact requires compactness of the sequence which was not directly available with the method developed here. Finally, uniqueness for the limit problem is still open.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the French ‘ANR blanche’ project Kibord: ANR-13-BS01-0004.

## Acknowledgements

While this work was in press, we learned about recent work by K. Trivisa, D. Donatelli and F. Weber on visco-elastic tissues and numerical methods (http://arXiv.org/abs/1408.4794, http://arXiv.org/abs/1408.4606, http://arXiv.org/abs/1504.05982).

## Footnotes

One contribution of 15 to a theme issue ‘Free boundary problems and related topics’.

- Accepted January 8, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.