EARLY AND LATE STAGE PROFILES FOR A CHEMOTAXIS MODEL WITH DENSITY-DEPENDENT JUMP PROBABILITY

In this paper, we derive a chemotaxis model with degenerate diffusion and density-dependent chemotactic sensitivity, and we provide a more realistic description of cell migration process for its early and late stages. Different from the existing studies focusing on the case of non-degenerate diffusion, this model with degenerate diffusion causes us some essential difficulty on the boundedness estimates and the propagation behavior of its compact support. In the presence of logistic damping, for the early stage before tumour cells spread to the whole domain, we first estimate the expanding speed of tumour region as O(tβ) for 0 < β < 1 2 . Then, for the late stage of cell migration, we further prove that the asymptotic profile of the original system is just its corresponding steady state. The global convergence of the original weak solution to the steady state with exponential rate O(e−ct) for some c > 0 is also obtained. 2010 Mathematics Subject Classification. Primary: 35K51, 35K65; Secondary: 92C17.


(Communicated by Yangjin Kim)
Abstract. In this paper, we derive a chemotaxis model with degenerate diffusion and density-dependent chemotactic sensitivity, and we provide a more realistic description of cell migration process for its early and late stages. Different from the existing studies focusing on the case of non-degenerate diffusion, this model with degenerate diffusion causes us some essential difficulty on the boundedness estimates and the propagation behavior of its compact support. In the presence of logistic damping, for the early stage before tumour cells spread to the whole domain, we first estimate the expanding speed of tumour region as O(t β ) for 0 < β < 1 2 . Then, for the late stage of cell migration, we further prove that the asymptotic profile of the original system is just its corresponding steady state. The global convergence of the original weak solution to the steady state with exponential rate O(e −ct ) for some c > 0 is also obtained.
1. Introduction. The motion of cells moving towards the higher concentration of a chemical signal is called chemotaxis. For example, bacteria moves toward the highest concentration of food molecules to find food. A well-known chemotaxis model was initially proposed by Keller and Segel [15] in 1971, subsequently, a number of variations of the Keller-Segel system were proposed and have been extensively studied during the past four decades, for example, see the survey papers [1,12] and the references therein. Especially, chemotaxis models also appear in medical mathematics. Many factors effect the migration mechanisms of tumour cells. For example, the extracellular matrix (ECM), to which the tumour cell to be attached, inhibits the cell polarizes and elongates to migrate. ECM-degrading enzymes (MDE) cleave ECM fibers into smaller chemotactic fragments to facilitate cell-migration [6]. In [4], Chaplain and Anderson introduced a model for tumour invasion mechanism, which describes tumour invasion phenomenon in accounting for the role of chemotactic ECM fragments named ECM*, produced by a biological reaction between ECM and MDE. In these models, the tumour cell random motility is assumed to be a constant, which leads to linear isotropic diffusion. However, in realistic situation, it is emphasized that migration of the tumour cells through the ECM fibers should rather be regarded like movement in a porous medium with degenerate diffusion from a physical point of view [38]. Compared with the classical tumour invasion model with linear diffusion, the mathematical analysis of the nonlinear diffusion system has to cope with considerable additional challenges and is much less understood. Several chemotaxis models with nonlinear diffusion have been recently proposed and analyzed, e.g. [18,38,45,46], where the nonlinear diffusions in these studies were still assumed to be non-degenerate. For tumour angiogenesis model and relevant mathematical analysis with or without degenerate diffusion, we refer to [14,22,23,48,49,50,52,54] and the references therein.
Tumour cells can modify their migration mechanisms in response to different conditions [6]. There are two potentially important factors: (i) the effect of cell-density on the probability of cell movement; (ii) the effect of signal-mediated cell-density sensing mechanisms on movement [28]. For interacting cell population, Painter and Sherratt [29] further presented four different sensing strategies: strictly local, neighbour based, local average and gradient. Cell movement involves the processing of multiple signals, each of them may act on the cells in different ways. For neighbor-based and gradient-based rules, Painter and Hillen [28] proposed volume filling approach, that is, the movement of cells is inhibited by the neighboring site where the cells are densely packed. Inspired by the idea of Painter et al. [28,29] and recently in Xu et al. [57], we extend Chaplain and Anderson's model [4] to a new one with density-dependent jump probability of tumour cells as follows, which is concerned with the competition between the following several biological mechanisms: degenerate diffusion, density-dependent chemotaxis, and general logistic growth. That is, = ∇ · (q(u)∇u) − ∇ · (q(u)u∇v) + µu δ (1 − ru), x ∈ Ω, t > 0, ∂v ∂t = ∆v + wz, x ∈ Ω, t > 0, x ∈ Ω, t > 0, The detailed derivation of the model (1) will be carried out in the Appendix. Here, Ω is a bounded domain in R n with smooth boundary. The four variables u, w, z and v represent the cancer cell density, ECM concentration, the MDE concentration and the ECM* concentration, respectively. q(u) denotes the jump probability of a cell depending on the population pressure at its present location, which is increasing with respect to u with q(0) = 0, q(1) = 1, namely, the jump probability is 1 when the cell density exceeds maximum and it is zero when the cell density is zero, and f (u) = µu δ (1 − ru) is the logistic growth term, where µ > 0 and r > 0 are the proliferation rate and reciprocal of carrying capacity, δ ≥ 1 is a constant. The unbiased cell movement modelled by linear diffusion motility mechanism has been used extensively to study a variety of cell biology problems. However, when cells are close enough for regular contacts, they will inevitably interact [29]. Linear diffusion of each cell type is inappropriate for the close-packed cell populations involved in early tumour growth. The degenerate nonlinear diffusion can represent "population pressure" in cell invasion models [29], which arises from the ecology dispersal literature [9,10,24,58]. A high cell density results in increased probability of a cell being "pushed" from a site. In this case, large dispersal takes place in highly populated regions, but low mobility occurs in the regions of low cell density. The cell invasions described by nonlinear degenerate systems with the density-dependent nonlinear diffusivity function q(u) = Du m−1 with m > 1 in the diffusion term ∇ · (q(u)∇u) have been paid more attention in recent years [13,39,53,57]. These tumour invasions models with porous media diffusion is degenerate at u = 0, that is, when the population density is zero, the diffusion coefficient is zero. In fact, biological evidence suggests that no cell migration (in particular no diffusivity) occurs in noncellular regions [20,59].
Some studies found that degenerate nonlinear diffusion model related to the porous media equation (PME) provides a better match to experimental cell density profiles [34]. Sengers and coworkers [30] developed a set of in vitro cell invasion experiments and image analysis to quantify the migration and proliferation of two different skeletal cell types, including human osteosarcoma MG63 cells and human bone marrow stromal cells (HBMSCs). Comparison of experimental and simulated cell distribution are shown in Fig. 1 in [30], where the cell density considerably increased and simultaneously spread outwards from the centre of the cell circle, producing a new cell migration front every day. Their results show that the MG63 migration with sharp front is best described by a degenerate diffusion model with the diffusivity q(u) = u m−1 with m = 2 [ Fig. 1(a)], while the HBMSC migration with smooth front corresponds to the solution of a linear diffusion equation [ Fig. 1(b)]. Similarly, Sherratt and Murray's work provides a physical connection between epidermal wound healing experimental data and the solutions of either the linear diffusion equation or the porous media equation to represent cell density profiles [33]. They showed that the solution of degenerate diffusion model with the diffusivity function q(u) = u 3 compare well with the experimental data in [42]. Mathematically, the PME raises the possibility of sharp-front waves, whereas the smooth-front waves arise in linear diffusion equations. The difference between these front types is that the sharp-front waves have distinct boundaries, and the population density decreases to zero at a finite point in space, rather than tending to zero asymptotically [31,43,55]. . This diagram was redrawn from the one in Ref. [30].
An interesting work related to the chemotaxis model mentioned above is [7], in which they considered the following chemotaxis system with linear diffusion It is proved the existence of global solutions and the asymptotic behaviors of global solutions as time goes to infinity by using the properties of the Neumann heat semigroup e t∆ in Ω. Recently, Li et al. [17] study the the quasilinear chemotaxis system (2) with the effect of the nonlinear diffusion q(u) ≥ Cu m−1 with C > 0 and the nonlinear chemotactic sensitivity function S(x, u, v) with some structural conditions for the above coupled tumor invasion system. They obtained the boundedness and large time behavior for this system. Apart from the diffusive motility, another important mechanism in cell invasion is cell proliferation. In [35], an assay using gut organ culture validates that proliferation at the invading front is the critical mechanism driving apparently directed invasion. Cells at the invasive front are proliferative and migrate into previously unoccupied tissue. It also has important implications for carcinoma invasion. Tumour invasion systems with proliferative cells have been studied extensively [8,21,26]. Logistic growth is one of important models of proliferation to a carrying capacity limit [24]. Von Bertalanffy derived a general logistic growth law for avascular tumour growth [44], and suggested that where γ, λ, δ, µ > 0 and µ > λ. In our tumour invasion model, the cell proliferation also plays an important role in biological modelling and theoretical study for the evolution of tumour boundary. Based on the structure of degenerate diffusion equation with the Von Bertalanffy's growth law, we compare its solution with the weak upper and lower Barenblatt-type self-similar solutions and we obtain the upper and lower bounds of the expanding rate of its support. These results provide mathematical predictions of the evolution of the tumour invasion boundary. In the proof of the lower bound of the expanding rate of its support, we utilize the combination of the degenerate diffusion and the proliferation from the logistic growth to balance the possible aggregation effect due to the chemotaxis, since this chemotaxis may cause backward diffusion and negative effects on the expanding of the support. Without this logistic growth, we find that the degenerate diffusion alone is insufficient to govern the possible aggregation effect. We note that the upper bound of the expanding rate of the support (i.e. the finite speed propagation property) is also valid for the system without logistic growth, whereas the lower bound of the support or the expanding property is insufficient in this case.
Compared to the linear cases, the chemotactic system with degenerate diffusion and chemotactic sensitivity is more complex and challenging. Since the first equation of (1) is degenerate at any point where u(x, t) = 0, there is no classical solution in general. The spatial derivatives of u may not exist in classical sense, and may even do not belong to the class of locally integrable generalized functions, that is, there might hold u ∈ W 2,1 loc . In this paper, we provide a more realistic description of cell migration process for early and late stages. It is worth to mention that our stability results of the model (1) give a certain estimate for the speed of the expanding speed of tumour region. We construct suitable subsolutions and supersolutions to show the position of the free boundary for the tumour region. Then, we prove that there exist t 0 and two families of monotone increasing open sets {A 1 (t)} t>0 , {A 2 (t)} t∈(0,t0) such that ∂A 1 (t) and ∂A 2 (t) have finite derivatives with respect to t, namely, {A 1 (t)} t>0 and {A 2 (t)} t∈(0,t0) both expand at finite speeds. This indicates the finite speed propagation property of our chemotaxis model. As shown late in Remarks 1 and 2, in the porous media diffusion case, we estimate that, at the early stage the expanding speed of tumour region is somehow like the algebraic rate of (1 + t) β for some β ∈ (0, 1 2 ). This resembles the case of the Stefan problem with a linear diffusion term.
As we all know, for linear diffusion equations with initial data u 0 , the solution u(x, t) > 0 for t > 0 and any x ∈ R N , thus a linear diffusion process predicts a non-zero u for arbitrarily large displacements at arbitrarily small time, namely, the underlying propagation speed is infinite [43]. This means that the initial tumour cells moving into regions of unoccupied tissue immediately in this biological system. However, the spatial support of the solution to the degenerate diffusion equation remains bounded for all time t > 0 [5]. There are distinct boundaries, called interfaces, beyond which the population density is zero. Our stability results of Figure 2. The growth curve of HEPA-1 spheroids. The solid line represents the position of the outer tumour boundary. Dimensional diameters are shown in µm. This diagram was redrawn from the one in Ref. [27]. days, decreased, and approached a horizontal asymptote. It can be difficult to decide what type of model is best suited to a particular biological problem. Different approaches in mathematics can reproduce the same experimental results [3]. Our theoretical results also provide a good fit to the experimental results in [16]. The shape of the growth curve of the cell spheroids is similar to the graph of power function R(t) = (1 + t) γ with 0 < γ < 1. Note that, our estimation of the expanding speed of tumour region with the algebraic rate of (1 + t) β for some β ∈ (0, 1 2 ) compares well with this experimental data. It indicates that the tumour cell model described by the degenerate nonlinear diffusion motility mechanism can describe the progress of the very early stage of tumor growth mathematically. In combination with experiments, this type of tumour model may prove useful in predicting the evolution of tumour cell migration, investigating subsequent stages of tumour progression and testing therapeutic strategies.
In contrast with the well known linear cases, the degenerate diffusion is endowed with the interesting feature of slow diffusion, that is, the compact support of solutions propagates at a finite speed. The slow diffusion feature has some advantages and accuracy for describing specified biological processes in the point of view of the physical reality, and it also leads to more challenges in the mathematical studies. For example, in order to investigate the asymptotic behavior of solutions, one must appropriately describe the propagation behavior of its support, which is more likely to be a compact subset of the prescribed domain for some time interval if the initial data are given so. We mention that the Neumann heat semigroup theory and functional transform methods have been proved to be effective in studying the global boundedness and large time behavior for the linear diffusion equations, but they are all inapplicable in the degenerate diffusion case due to the nonlinearity. We establish the global existence of bounded weak solutions to this model by energy estimate technique and methods based on Moser-type iteration. Then we prove that, at the late stage of the tumour migration, the original weak solution time-asymptotically converges to its steady state, even if the initial perturbation is large, namely, the global stability of the steady state. The adopted approach is the technical compactness analysis with the help of the comparison principle deduced by the approximate Hohmgren's approach and two kinds of lower solutions showing the expanding support and the exponentially convergence. The one is a self similar weak lower solution of Barenblatt type and the other kind is an ODE solution.
This paper is organized as follows. In Section 2, we state our main results. We leave the global existence of weak solutions to the corresponding chemotaxis system and their regularity into Section 3 as preliminaries. Section 4 is devoted to the study of compact support property of the tumour cells at early stage and the large time behavior at late stage, showing the exponentially convergence of solutions.
2. Main results. In this section, we first state our main results on the study of expanding compact support of the tumour cells at early stage and the asymptotic behavior at late stage. We leave the detailed derivation on the new chemotaxis model (1) with density-dependent jump probability in the Appendix. We estimate the upper bound and lower bound for expanding speed of tumour cell region at early stage (before the tumour cells spread to the whole body) and show the exponentially convergence of solutions for large time.
We consider the following system (3) with degenerate diffusion where δ ≥ 1, µ > 0, u 0 , v 0 , w 0 , z 0 are nonnegative functions, ν is the unit outer normal vector, and q(u) ≥ 0 with q(0) = 0. Here and after, the IBVP (3) will be our main target equations. Since degenerate diffusion equations may not have classical solutions in general, we need to formulate the following definition of generalized solutions for the initial boundary value problem (3).
As preliminaries, we leave the global existence and regularity results into Section 3. Our main results concerned with the description of cell invasion processes are as follows. First, we show that the evolution of tumour invasion in the very early stage.
Theorem 2.2 (Early stage profile -upper bound). Let (u, v, w, z) be a globally bounded weak solution of (3) with the initial data for some r 0 > 0 and x 0 ∈ Ω. Then there exists a time t 1 > 0 and a family of monotone increasing open sets {A(t)} t∈(0,t1) such that supp u(·, t) ⊂ A(t) ⊂ Ω, t ∈ (0, t 1 ), and ∂A(t) has a finite derivative with respect to t. More precisely, we can choose with some appropriate η, τ > 0.
Remark 1. As a typical finite propagating model, the Barenblatt solution of the porous medium equation is with k = 1/(m − 1 + 2/n) < n/2 for m > 1, and its support is expanding at the rate (1 + t) k/n . Here we have proved the tumour cells are located within a ball expanding at the rate (1 + t) 1/2 . We note that the upper bound of the expanding rate of the support is also valid for the system without logistic growth.
Next, we show the propagating property of the tumour cells at the early stage.

Theorem 2.3 (Early stage profile -lower bound)
. Let (u, v, w, z) be a globally bounded weak solution of (3). Assume that 1 ≤ δ < m, Ω is convex and u 0 ≡ 0. Then there exists a time t 2 > 0 such that the support of u expands to the whole Ω when t ≥ t 2 . Precisely speaking, there exist a family of monotone increasing open sets and A(t) = Ω for t ≥ t 2 , ∂A(t) has a finite derivative with respect to t.
Remark 2. For this chemotaxis system, we proved that the tumour cells will expand to the whole body when the time t increases. Compared with the porous medium equation, whose Barenblatt solution B(x, t) in (4) is expanding at the rate (1 + t) 2k/n , the tumour cells of (3) migrate to at least a ball expanding at the rate (1 + t) β . Here in the proof we have selected β > 0 sufficiently small, which means the support is expanding with a much slower rate.
Under the hypotheses of Theorem 2.2 and Theorem 2.3, we see that there exist t 0 and two family of monotone increasing open sets ∂A 1 (t) and ∂A 2 (t) have finite derivatives with respect to t, which means that {A 1 (t)} t>0 and {A 2 (t)} t∈(0,t0) both expand at finite speeds. This indicates immediately the finite speed propagation property of this chemotaxis model. After the tumour cells spread to the whole domain, we can investigate the large time behavior. We show that the solution converges to its steady state exponentially.
Theorem 2.4 (Late stage profile). Let (u, v, w, z) be a globally bounded weak solution of (3). Assume that the hypothesis in Theorem 2.3 is valid. Then there exist C and c > 0 such that The main difficulty lies in proving the expanding property of the support of the first component u. We first prove the comparison principle by the approximate Hohmgren's approach, and then construct two kinds of lower solutions. The one is a self similar weak lower solution with much slower expanding support and slightly faster decaying maximum compared with the Barenblatt solution to the porous medium equation, the other kind is an ODE solution. After showing the expanding property, we formulate several upper and lower solutions that converge to steady state exponentially by utilizing the exponential decay of other components.
3. Preliminaries: Global existence, boundedness and regularity. As preliminaries, we prove the existence, boundedness and regularity of a global weak solution in this section. The main preliminary results are as follows. Theorem 3.2 (Regularity). Let (u, v, w, z) be a globally bounded weak solution of (3). Then there exist α ∈ (0, 1) and C(p) > 0 such that , for any p > 1 and t ≥ 1.
We first use the artificial viscosity method to get smooth approximate solutions. Despite the absence of comparison principle, we can prove a special case compared with a lower solution, which is helpful for establishing the regularity estimates. By making use of the special structure of dispersion, we carry on the estimates on u m in W 1,2 (Q T ), instead of u. These energy estimates ensure the global existence of weak solution.
Consider the following corresponding regularized problem and ∂w0ε ∂ν = 0 on ∂Ω. The local existence of the regularized problem (5) is trivial and we denote the unique solution by (u ε , v ε , w ε , z ε ). Let (0, T max ) be its maximal existence interval.
As usual, there is no comparison principle for the system, because the system is strongly coupled. However, we have the following lemma.
Since u ε ≥ 0, the first equation of (5) is equivalent to Now we present some energy estimates independent of time t and the parameter ε.
Lemma 3.4. The first solution component u ε satisfies Proof. We denote u ε by u in this proof for the sake of simplicity. Since u is nonnegative and ∂u ∂ν = ∂v ∂ν = 0 on ∂Ω, integration of the first equation of (5) over Ω yields for all t ∈ (0, T max ). We note that and
Here we recall some lemmas about the L p -L q type estimates for the components of the solution, and we refer the readers to [7] for details.
Then for any T ∈ (0, T max ], there exists a constant C v (q, r) such that Lemma 3.7. There holds and Proof. Since both w ε and z ε are nonnegative, it is clear from the third equation of We add the third to the second equation of (5) and integrate over Ω to obtain Thus, for all t ∈ (0, T max ).
For any r ≥ 1, there exists a constant C(r) independent of t and ε such that Proof. According to Lemma 3.4, u ε L 1 (Ω) is uniformly bounded. Since n ≤ 3, we can apply Lemma 3.5 and 3.6 to complete this proof.
The following Gagliardo-Nirenberg inequality (see [46,51]) will be used in deriving the L p estimates of u ε .
We present the following L p estimate of u ε .
Lemma 3.10. Let 1 ≤ n ≤ 3. For any given p ≥ 1, there exists a constant C(p) > 0 independent of t and ε such that Proof. We denote u ε , v ε by u, v in this proof for the sake of simplicity. By a straightforward computation, testing the first equation in (5) by u r for r > 0 and integrating by parts we find that We note that and where C 1 and C 2 are constants independent of t. Then by Young's inequality, we see that
We now employ the following Moser-type iteration to get the L ∞ (Ω) estimate of u.
Lemma 3.12. Let 1 ≤ n ≤ 3. There exists a constant C > 0 independent of T max and ε such that Proof. We denote u ε , v ε by u, v in this proof for the sake of simplicity. We test the first equation in (5) by u r for r > 0 and integrating by parts we find that Similar to the proof of Lemma 3.10, using Young's inequality we can estimate and Ω u m ∇v · ∇u r dx ≤ r where according to Lemma 3.11 ∇v L ∞ (Ω) is uniformly bounded. Now we apply the Gagliardo-Nirenberg inequality Lemma 3.9 to obtain Ω u m+r dx = u where a = n/(n + 2) ∈ (0, 1) and C 0 is the constant in the Gagliardo-Nirenberg inequality which is independent of r. Therefore, we have where Inserting the above estimates (12), (13) into (11) yields where C 2 (r) = (r + 1) 4 δ+r µ|Ω| + 4 µ r δ+r |Ω| + 2 µ r+1 δ |Ω|.
Now we use the following Moser-type iteration. Let r = r j , with r j = 2 j + m − 2 for j ∈ N + , that is, r 1 = m and We can invoke Lemma 3.10 to find C 0 such that sup t∈(0,Tmax) From (14) and an ODE comparison, we have A simple analysis shows that C 1 (r)(r + 1) ≤ a 1 r b1 and C 2 (r) ≤ a 2 b r 2 for some positive constants a 1 , a 2 and b 1 , b 2 that all are greater than 1 and independent of r. Therefore, we can rewrite the above inequality (15) Let Since boundedness of u in L ∞ (Ω) is evident in the case when M j ≤ max{ Ω (u 0 + 1) rj +1 dx, 1} for infinitely many j ≥ 1, we may assume that M j ≥ max{ Ω (u 0 + 1) rj +1 dx, 1} and thus, according to (16), there holds We note that if M 2 j−1 ≤ a 2 b rj 2 for infinitely many j ≥ 1, then for j sufficiently large, which shows the boundedness of u in L ∞ (Ω). Otherwise, Therefore, we can rewrite (17) into for all j ≥ j 0 with a constant D independent of j, whence upon enlarge D if necessary we can achieve that (18) actually holds for all j ≥ 1. By introduction, this yields for all j ≥ 1, and hence that for all j ≥ 1. This implies that u indeed belongs to L ∞ (Ω × (0, T max )). Now we turn to the regularity estimates.
And the third solution component w ε fulfills for all t ∈ (0, T max ).
Proof. According to Lemma 3.12, Lemma 3.5, Lemma 3.6, we see that u ε L ∞ (Ω , z ε L ∞ (Ω , ∇v ε L ∞ (Ω are uniformly bounded in (0, T max ). The standard L p − L q type estimates also shows the boundedness of v ε L ∞ (Ω and ∇z ε L ∞ (Ω . We denote v ε , w ε , z ε by v, w, z in this proof for the sake of simplicity. Since both w and z are nonnegative according to the third and fourth equation in (5) and the initial data, we have Therefore, This completes the proof.
Lemma 3.14. There exists a constant C > 0 independent of ε and T , such that Proof. We denote v ε , w ε , z ε by v, w, z in this proof for the sake of simplicity. Multiplying the second equation in (5) by −∆v and integrating over Ω yields since ∇v, z and ∇z are uniformly bounded in L ∞ (Ω) according to Lemma 3.13.
Integrating over (0, T ), we complete this proof.
Lemma 3.15. There exists a constant C > 0 independent of ε and T , such that Proof. We denote u ε , v ε by u, v in this proof for the sake of simplicity. We test the first equation in (5) by (u + ε) m and get According to Lemma 3.11 and Lemma 3.12, ∇v and u are uniformly bounded. Thus, where C 1 is a constant independent of t and ε. Integrating (19) on (0, T ) yields We note that This completes the proof.
Lemma 3.16. There exists a constant C > 0 independent of ε and T , such that for all T ∈ (0, T max ).
Proof. We denote u ε , v ε by u, v in this proof for the sake of simplicity. We multiply the first equation in (5) by [(u + ε) m ] t and then we have We note that u L ∞ (Ω) is uniformly bounded and then where C 1 , C 2 , C 3 are constants independent of t and ε. We also have and

TIANYUAN XU, SHANMING JI, CHUNHUA JIN, MING MEI AND JINGXUE YIN
There holds where C 4 and C 5 are constants independent of t and ε, since the uniform boundedness of ∇v L ∞ (Ω) . Inserting the above inequalities into (21), and noticing the inequality (20) in the proof of Lemma 3.15, we find a constant C independent of t and ε such that The proof is completed.
Proof of Theorem 3.1. According to the estimates, for any ε, the approximation solution (u ε , v ε , w ε , z ε ) exists globally. The regularity estimates of v ε , w ε , z ε are trivial. For any T ∈ (0, ∞), we see that u m ε ∈ L ∞ (Q T ), ∇u m ε ∈ L 2 (Q T ), and ∂u m ε /∂t ∈ L 2 (Q T ), Thus, there exists a function u ∈ W 1,2 (Q T ), such that u m ε weakly in W 1,2 (Q T ) and strongly in L 2 (Q T ) converges to u. We denote u = u 1/m since u ≥ 0. Thus, u m ε converges almost everywhere to u m , and u ε converges almost everywhere to u. We can verify the integral identities in the definition of weak solutions. By taking a sequence of T ∈ (0, ∞) and the diagonal subsequence procedure, we can find the existence of a global weak solution.

Now we show the regularity of the globally bounded weak solution.
Lemma 3.17. Let (u, v, w, z) be a globally bounded weak solution of (3) such that u(·, t) L γ (Ω) is uniformly bounded with γ = max{1, n/3}. Then there exists a constant C such that Proof. Since u(·, t) L n 3 (Ω) is uniformly bounded, for any r ≥ 1 we can apply Lemma 3.5 and 3.6 to find a constant C(r) independent of t such that ∇v(·, t) L r (Ω) ≤ C(r) for all t > 0. The estimates in the proof of Theorem 3.1 in section 3 can be carried on to complete this proof here.
In Lemma 3.13, we have proved ∇w(·, t) L ∞ (Ω) ≤ C(1 + t) (same as ∇w ε ) for some constant C > 0. However that is an estimate depending on time t. Employing the method in the proof of Lemma 4.4 in next Section and iteration technique, we can prove the following uniform estimate.
Lemma 3.18. Let (u, v, w, z) be a globally bounded weak solution of (3). Then for any p ≥ 1 there holds for some constant C(p) independent of time t.
Proof. This proof proceeds along the idea of the arguments of Lemma 4.3 in [47] and Lemma 4.1 in [40]. Since We see that

≤C.
Using the same method, we have ∇z(x, s)ds according to the proof of the previous estimate on ∇w L 2 (Ω) and the boundedness of z L ∞ (Ω) . Repeating this process for ∇w L k (Ω) with k = 6, 8, . . . , we complete this proof by iteration.
Proof of Theorem 3.2. Lemma 3.18 shows the uniform bound of ∇w L n+2 (Ω) . According to the third equation of (3), we see that w t L ∞ (Ω) = wz L ∞ (Ω) ≤ C. Therefore, w ∈ W 1,n+2 (Ω × (t, t + 1)) and its norm is uniformly bounded for any t > 0. Sobolev embedding theorem implies the existence of α ∈ (0, 1) and C > 0 such that Since u is uniformly bounded, the strong solution theory of parabolic equation applied to the fourth equation in (3) shows z t L p (Ω×(t,t+1)) + ∆z L p (Ω×(t,t+1)) ≤ C(p), t > 0, for some constant C(p) > 0. Taking p > 1 + n/2, we see that for some α ∈ (0, 1) Thus, This can also be deduced by ∇(wz) L p (Ω) + (wz) t L p (Ω×(t,t+1)) ≤ C, t > 0, with p > n + 1. Using bootstrap arguments involving the standard parabolic regularity theory, we can verify that v C 2+α,1+α/2 (Ω×[t,t+1]) + z W 2,1 p (Ω×(t,t+1)) ≤ C(p). The proof is completed. 4. Propagating properties and large time behavier. This section is devoted to the study of the propagating properties of the tumour cells and the large time behavior of the weak solution (u, v, w, z) to the problem (3). In contrast with the heat equation, it is known that the porous medium equation has the property of finite speed of propagation. Therefore, the first component u may not have positive minimum for some time t > 0. We use the comparison principle together with two kinds of weak lower solutions, one is decaying but its support is expanding with finite speed of propagation, the other one is an increasing function of time t, to overcome the difficulty of degenerate dispersion.
We first present the following comparison principle of the first component.
Proof. The following inequality holds for all 0 ≤ ϕ ∈ L 2 ((0, T ); W 1,2 (Ω)) ∩ W 1,2 ((0, T ); L 2 (Ω)) with ϕ(x, T ) = 0. Let Since ∇v, u 1 , u 2 are bounded, there exists a constant C > 0 such that |b| ≤ Ca and |c| ≤ C. Henceforth, a generic positive constant (possibly changing from line to line) is denoted by C. However, c is not bounded by Ca and we have no estimate on ∇c since we only assume that δ ≥ 1. Then for all 0 ≤ ϕ ∈ L 2 ((0, T ); W 1,2 (Ω)) ∩ W 1,2 ((0, T ); L 2 (Ω)) with ϕ(x, T ) = 0 on Ω and ∂ϕ ∂ν = 0 on ∂Ω × (0, T ), there holds We employ the standard duality proof method or the approximate Hohmgren's approach to complete this proof (see Theorem 6.5 in [43], Chapter 1.3 and 3.2 in [56]). For any smooth function ψ(x, t) ≥ 0, we solve the inverse-time problem where κ > 0, θ > 0, a ε is a smooth approximation of a, a ε ≥ a, and This definition of c θ allows us to find a constant C(θ) such that We may also need to replace b(x, t) and c θ (x, t) by their smooth approximation functions b ε (x, t) and c θ,ε (x, t) respectively in (22). For the sake of simplicity we omit this procedure. Here we note that (22) is a standard parabolic problem as the initial data is imposed at the end time t = T . Therefore, it has a smooth solution ϕ ≥ 0. Maximum principle shows the boundedness of ϕ. Then we get the estimate Now we need the a priori estimate on a ε |∆ϕ| 2 . We can assume that T is appropriately small, otherwise we can prove step by step on each time interval. We multiply the equation (22) by η(t)∆ϕ where 1/2 ≤ η(t) ≤ 1 is a smooth function with η (t) ≥ M > 0 for t ∈ (0, T ). Since T is small, we can choose M appropriately large. Integrating over Q T yields Using ϕ(x, T ) = 0, we have Therefore, It follows that which converges to zero if we let ε → 0. For any fixed γ > 0, denote Then there exists a constant C(γ) such that a(x, t) ≥ C(γ) on F γ and , which converges to zero if we first let κ → 0 and then let γ → 0. We also have which converges to zero if we let θ → 0. Now we conclude that for any given ψ ≥ 0 and then u 1 ≥ u 2 almost everywhere on Q T .
Here we recall some lemmas about the asymptotic behavior of solutions to evolutionary equations. In particular, ∇v(·, t) L ∞ (Ω) → 0, as t → ∞.
where u(x, t) ≥ 0 is given. Then there exist constants C 1 and C 2 > 0 only depend on diamΩ and sup τ <t u L 1 (Ω) respectively, such that with u ≥ 0 on Ω × R + and ∂w0 ∂ν = 0 on ∂Ω, then for all t > 0. Now we construct a self similar weak lower solution with expanding support.
We denote C 1 = ∇v L ∞ (Ω×R + ) and C 2 = ∆v L ∞ (Ω×R + ) for convenience, since they are bounded according to Theorem 3.2. A sufficient condition of inequality (25) is As we have chosen d = 1/(m − 1) and κ = (1 − β)/(m − 1), we rewrite (26) into Let ε, β and η be chosen such that Thus, for (24) and (28), it suffices to choose η = r 2 , and then β = 4ε m−1 m/(m − 1). Now, we find a weak lower solution with expanding support and comparison principle Lemma 4.1 implies There exists a t 0 such that and thus Next, we construct another constant lower solution with 0 < ε 0 ≤ ε(1 + t 0 ) −κ (η/2) d ≤ 1/2 to be determined. Clearly, ∂u ∂ν = 0 on ∂Ω. We only need to check the following differential inequality since δ < m and ∆v is uniformly bounded according to Theorem 3.2. Applying the comparison principle Lemma 4.1 again, we find This completes the proof. with k = 1/(m − 1 + 2/n). The Barenblatt solution B(x, t) is decaying at the rate (1 + t) −1/(m−1+2/n) in L ∞ (R n ) and the support is expanding at the rate (1 + t) 2k/n . While the self similar weak lower solution g(x, t) is decaying at the rate (1 + t) −(1−β)/(m−1) and its support is expanding at the rate (1 + t) β . Here in the proof we have selected β > 0 sufficiently small, which means the support of g is expanding with a much slower rate and the maximum of g is decaying at a slightly faster rate.
Proof of Theorem 2.3. This has been proved in Lemma 4.5.

Proof. Applying Lemma 4.3, we see that
since u(x, t) ≥ ε 0 for x ∈ Ω and t > t 0 according to Lemma 4.5. Therefore, This is also valid for t ∈ (0, t 0 ) upon enlarging C 1 if necessary and hereafter we only need to prove this lemma for t > t 0 . We also have with 0 < c 1 < c 1 . We may write C 1 and c 1 as C 1 and c 1 for simplicity. Therefore, It follows form the second equation in (3) where λ 1 > 0 is the first nonzero eigenvalue of −∆ with homogeneous Neumann boundary condition. The L ∞ estimate of ∇v can be deduced in a similar way. In the proof of Lemma 3.7, we have obtained which is the same as the estimate of v ε + w ε . It follows from (29) that w(x, t) is decaying to zero exponentially. This implies that The proof is completed.
Lemma 4.7. For constants C, c > 0 and m > 1, the local solution g of the following ODE Integrating over (0, t) shows .
A simple analysis completes this proof.
Lemma 4.8. Let (u, v, w, z) be a globally bounded weak solution of (3) with the first component initial data u 0 ≥ 0, u 0 ≡ 0 and 1 ≤ δ < m. Then there exist constants C 3 > 0 and c 3 > 0 independent of t such that for all t > 0.
Proof. Lemma 4.5 implies that u(x, t) ≥ ε 0 for x ∈ Ω and t > t 0 . It suffices to prove this lemma for t ≥ t 1 with some fixed t 1 ≥ t 0 to be determined. We use upper and lower solution method to achieve this. Let u 1 (t) and u 2 (t) be one pair of the solutions of the following ODE Lemma 4.1 shows that

TIANYUAN XU, SHANMING JI, CHUNHUA JIN, MING MEI AND JINGXUE YIN
We only need to find one pair of (u 1 , u 2 ) such that u 1 and u 2 both converge to 1 exponentially. A sufficient condition of (30) is since ∆v(·, t) L ∞ (Ω) ≤ C 2 e −c2t according to Lemma 4.6. We note that we can choose t 1 sufficiently large such that Lemma 4.7 implies that u 1 (t) is uniformly bounded by some constant C. And a simple ODE comparison shows that u 1 (t) > 1 for all t > t 1 . Therefore, We see that u 1 (t) is an upper solution of u(x, t) and an upper solution of u 1 (t) is u 1 (t) such that which can be solved as On the other hand, the lower solution of u(x, t) satisfies We note that we can choose t 1 sufficiently large that ). An ODE comparison shows that ε 0 ≤ u 2 (t) < 1 for all t > t 1 and We see that u 2 (t) is a lower solution of u(x, t) and a lower solution of u 2 (t) is u 2 (t) such that This can also be solved as Thus, we conclude and u 2 (t), u 1 (t) converge to 1 exponentially.
Lemma 4.9. Let (u, v, w, z) be a globally bounded weak solution of (3) with the first component initial data u 0 ≥ 0, u 0 ≡ 0 and 1 ≤ δ < m. Then there exist constants C 4 > 0 and c 4 > 0 independent of t such that for all t > 0.
Proof. From the fourth equation in (3), we have We note that which can be deduced by solving the ODE z = −z + 1 with z(0) = 0. Therefore, The proof is completed.
Proof of Theorem 2.4. This is proved by collecting Lemma 4.5, Lemma 4.6, Lemma 4.8 and Lemma 4.9.
Finally, we construct a self similar upper solution with expanding support to prove Theorem 2.2. We note that for constructing a weak upper solution for the heat equation, one should replace the cut-off composite function (·) + by (·) − . But here for the degenerate porous medium type equation and the self similar function of the form g = [(1 − |x| 2 ) + ] d with md > 1, we can check that ∇g m is continuous and ∆g m ∈ L q (Ω) for some q > 1. This shows that the differential inequality for an upper solution only need to be valid almost everywhere, without the possible Radon measures on the boundary of its support. Lemma 4.10. Let (u, v, w, z) be a globally bounded weak solution of (3). We further assume that supp u 0 ⊂ B r0 (x 0 ) ⊂ Ω, for some r 0 > 0 and x 0 ∈ Ω. Define a function . Then by appropriately selecting β, σ ε, η and τ , the support of g(x, t) is contained in Ω for t ∈ (0, t 0 ) with some t 0 > 0 and the function g(x, t) can be an upper solution of the first equation in (3) on Ω × (0, t 0 ), that is, in the sense that the following inequality holds for all test functions 0 ≤ ϕ ∈ L 2 ((0, t 0 ); W 1,2 (Ω)) ∩ W 1,2 ((0, t 0 ); L 2 (Ω)) with ϕ(x, t 0 ) = 0 on Ω and g(x, 0) ≥ u 0 (x) ≥ 0 on Ω. Therefore, u(x, t) ≤ g(x, t) and there exist a family of monotone increasing open sets {A(t)} t∈(0,t0) such that supp u(·, t) ⊂ A(t) ⊂ Ω, t ∈ (0, t 0 ), and ∂A(t) has a finite derivative with respect to t.
The comparison principle Lemma 4.1 implies that u(x, t) ≤ g(x, t) for all x ∈ Ω and t ∈ (0, t 0 ). Thus, supp u(·, t) ⊂ A(t) = {x ∈ Ω; |x − x 0 | 2 < η(τ + t) β }, t ∈ (0, t 0 ), and ∂A(t) = {x ∈ Ω; |x − x 0 | = η 1 2 (τ + t) β 2 }, t ∈ (0, t 0 ), which has finite derivative with respect to t. with k = 1/(m − 1 + 2/n). The Barenblatt solution B(x, t) is decaying at the rate (1 + t) −1/(m−1+2/n) in L ∞ (R n ) and the support is expanding at the rate (1 + t) 2k/n . As we have shown the support of the lower solution in Lemma 4.5 is expanding with a much slower rate and decaying at a slightly faster rate. Here, the upper solution is increasing at the rate (τ + t) σ and its support is expanding at the rate (τ + t) β . The increasing of g(x, t) makes it possible to be an upper solution, which can be seen from the proof.
Remark 5. From the proof of Lemma 4.10, we can choose β > 0 to be as small as we want. But we note that supp u 0 ⊂⊂ supp g(·, 0) and if we choose a smaller β > 0, then the parameters τ and t 0 are also smaller. This shows if we let the upper solution expands slower, then it may only be an upper solution for a smaller time interval. Thus, the slower expanding upper solution g(x, t) on a smaller time interval does not contradict to the possible feature that the solution u(x, t) expands at a fixed rate since supp u 0 ⊂⊂ supp g(·, 0) at the initial time.
Proof of Theorem 2.2. This has been proved in Lemma 4.10.
Appendix. In this section, we extend the derivation of the classical taxis models in [36]. The derivation of the model begins with a master equation for a continuoustime and discrete-space random walk where T ± i (·) denote the transitional-probabilities per unit time of a one-step jump to i ± 1 and u i denotes the cell density at i.
Painter and Hillen [11,28] proposed volume filling approach. In this model, the transitional probability then takes the form where q(u) denotes the probability of a cell finding space at its neighboring location, constant α is the intrinsic dispersion coefficient, constant β the coefficient signal detection, v i the signal concentration, and τ the mechanism of tactic responses in cell populations, such as chemotaxis, haptotaxis or phototaxis. Substituting (40) to the master equation (39), in the PDE limits they derives ∂u ∂t = ∇ · (d 1 (q(u) − q (u))∇u − χ(v)q(u)u∇u) where d 1 = kα, χ(v) = 2kβ dτ (v) dv , k is a scaling constant. Note that q(u) is a nonincreasing function in this model, which says that the probability of a cell finding space at its neighboring site decreases in the cell density at that site.
Since a different combination of the above strategies may be necessary to reflect cell movement, we combine the local and gradient-based strategies and assume the transitional probability of the form where q(u) represents the jump probability of a cell due to the population pressure at present site. At the microscopic level, a high cell density results in increased probability of a cell being "pushed" from departure site [19,25,29], for example due to the pressure exerted by neighboring cells. We shall assume that only a finite number of cells, U max , can be accommodated at any site. We study the relative densityũ = u/U max , (and drop the symbol˜for simplicity). Moreover, the jump probability is 1 when the cell density exceeds U max and it is zero when the cell density is zero. Thus we stipulate the following conditions on q: q(0) = 0, q(1) = 1 and q(u) ≥ 0, for all 0 ≤ u ≤ 1.
A natural choice for q(u) is which states that the probability of a jump leaving one site increases with the cell density at that site [24,37].
Substituting (41) into the Master Equation (39) gives: We set x = kh, interpret x as a continuous variable and extend the definition of u i accordingly. The transitional probabilities of jumping to a neighboring location depend on the spatial scale h. Thus we assume that T ± h = k h 2 T ± for some scaling constant k. Expanding the right-hand side with respect to h, we obtain for the cell density u(x, t): By taking the limit of h → 0, we arrive at the following model where D u = kα, χ(v) = 2k dτ (v) dv . The function χ(v) is commonly referred as the tactic sensitivity function. The simplest form is χ(v) = χ 0 with χ 0 being a constant.
Apart from that, we consider a modification of the Verhulst logistic growth term to model organ size evolution introduced by Blumberg [2] and Turner [41], which is called hyper-logistic function, accordingly f (u) = ru δ (1 − µu).