A Review on the Qualitative Behavior of Solutions in Some Chemotaxis–Haptotaxis Models of Cancer Invasion

Abstract

Chemotaxis is an oriented movement of cells and organisms in response to chemical signals, and plays an important role in the life of many cells and microorganisms, such as the transport of embryonic cells to developing tissues and immune cells to infection sites. Since the pioneering works of Keller and Segel, there has been a great deal of literature on the qualitative analysis of chemotaxis systems. As an important extension of the Keller–Segel system, a variety of chemotaxis–haptotaxis models have been proposed in order to gain a comprehensive understanding of the invasion–metastasis cascade. From a mathematical point of view, the rigorous analysis thereof is a nontrivial issue due to the fact that partial differential equations (PDEs) for the quantities on the macroscale are strongly coupled with ordinary differential equations (ODEs) modeling the subcellular events. It is the goal of this paper to describe recent results of some chemotaxis–haptotaxis models, inter alia macro cancer invasion models proposed by Chaplain et al., and multiscale cancer invasion models by Stinner et al., and also to introduce some open problems.

1. Introduction

Chemotaxis is an oriented movement of individual cells in response to some signaling chemical (chemoattractant), and is regarded as a universal migration mechanism in a wide range of biological processes such as the migration of embryonic cells to developing tissues, and immune cells to infection sites. Accordingly in the past several decades, chemotaxis models have received a great deal of attention in the academic literature due to their potential to generate aggregation patterns in several relevant situations [1,2,3,4]. In this regard, the most intensively studied chemotaxis model is the celebrated Keller–Segel system of the form

$$\left\{\begin{array}{c}{n}_t=▵n-\nabla ·(n\nabla c),\hfill \\ c_t=\Delta c-c+n,\hfill \end{array}\right.$$

(1)

which was introduced in [5,6] to model the aggregation phenomenon undergone by the slime mold Dictyostelium discoideum. Indeed, the most striking feature of (1) is the occurrence of the critical mass blow-up phenomena in two-dimensional domains (see [7,8,9,10,11] for the analogue addressing its parabolic–elliptic version).

In contrast to chemotaxis, haptotaxis is understood as the migration of individual cells in response to the gradient of an immovable signal. Many biochemical mechanisms, inter alia chemotaxis and haptotaxis, play an important role in a plethora of biochemical processes and thereby influence cancer invasion and metastasis [12,13]. Indeed, in the process of the cancer cell invasion of surrounding healthy tissue, apart from random motion, cancer cells migrate toward increasing concentrations of a diffusible enzyme as well as toward higher densities of non-diffusible tissue by detecting matrix molecules such as vitronectin adhered therein. The latter biased migration of cancer cells is usually called chemotaxis, whereas the former is referred to as haptotaxis [1]. In order to gain a comprehensive understanding of the invasion–metastasis cascade, a number of mathematical models have been introduced for various aspects of cancer invasion and metastasis [12,13,14,15,16,17,18,19,20]. For example, a reaction–diffusion system was introduced in [17] to describe the interaction between the density of normal cells, cancer cells, and the concentration of $H^{+}$ ions produced by the latter. Particularly, it is recognized that cancer cells can up-regulate certain mechanisms that allow for extrusion of excessive protons and accordingly acidify the peritumoral region. The high acidity triggers the apoptosis of normal cells and then allows tumor cells to proliferate and invade into the surrounding tissue [21]. Furthermore, taking into account the microscopic dynamics of intracellular protons and their exchange with extracellular counterparts, a population-based micro–macro model for acid-mediated tumor invasion was proposed by Meral et al. in [18]. These continuum micro–macro models explicitly involving subcellular events are rather novel, especially in the context of cancer cell migration [22,23].

From a mathematical point of view, one substantial obstacle to any qualitative analysis of models of cancer invasion and metastasis consists of the coupling of partial differential equations (PDEs) for the quantities on the macroscale with ordinary differential equations (ODEs) modeling the subcellular events. In fact, the considerable difficulty in the context of the rigorous analysis stems from the lack of smoothing action on the spatial regularity of ODE. To the best of our knowledge, the mathematical well-posedness of various models of cancer invasion has been receiving increased interest in the literature [1,22,23,24,25,26,27,28,29,30,31,32,33]. Without the pretension of exhaustiveness, this paper provides a short review of the global bounded results on some cancer invasion models and sketches necessary proofs thereof.

The rest of paper is organized as follows: Section 2 provides the global existence and large time behavior to the macroscopic cancer invasion models, particularly in the case that the extracellular matrix (ECM) remodeling is taken into account. Section 3 shows the global existence of weak solutions to multiscale cancer invasion models. Finally, a brief summary and open questions are presented in Section 4.

2. Macroscopic Cancer Invasion Models

As an important extension of the Keller–Segel system (1), the following system was proposed by Chaplain and Lolas [12,13] to simulate the cancer invasion of surrounding normal tissue

$$\left\{\begin{array}{c}{u}_t=\Delta u-\chi \nabla ·(u\nabla v)-\xi \nabla ·(u\nabla w)+\mu u(r-u-w),x\in \mathsf{\Omega },t>0,\hfill \\ {\displaystyle \sigma v_t=\Delta v-v+u,x\in \mathsf{\Omega },t>0,}\hfill \\ {\displaystyle w_t=-vw+\eta w(1-w-u),x\in \mathsf{\Omega },t>0,}\hfill \\ {\displaystyle \frac{\partial u}{\partial \nu }-\chi u\frac{\partial v}{\partial \nu }-\xi u\frac{\partial w}{\partial \nu }=\frac{\partial v}{\partial \nu }=0,x\in \partial \mathsf{\Omega },t>0,}\hfill \\ {\displaystyle u(x,0)=u_0\left(x\right),\sigma v(x,0)=\sigma v_0\left(x\right),w(x,0)=w_0\left(x\right),x\in \mathsf{\Omega }}\hfill \end{array}\right.$$

(2)

where $\mathsf{\Omega }\subset {\mathbb{R}}^n$ is a bounded domain, $\partial /\partial \nu$ denotes the outward normal derivative on $\partial \mathsf{\Omega }$, and the model variables are u, the density of cancer cells; v, the concentration of the matrix-degrading enzyme (MDE); and w, the concentration of the extracellular matrix (ECM), respectively. $\chi$ and $\xi$ measure the chemotactic and haptotactic sensitivities, respectively. The term $\mu u(r-u-w)$ in the first equation of (2) implies that in the absence of the ECM, cancer cells proliferate according to the standard logistic law, and $\eta >0$ embodies the ability of the ECM to remodel back to a healthy level. Here the parameter $\sigma =0$ or 1, especially the underlying mechanism of $\sigma =0$ is that the diffusion of the enzyme is much faster in comparison to that of cancer cells [12], which may also follow an approach of the quasi-steady-state approximation frequently used to study minimal chemotaxis systems [10].

When $\chi =0$, the system (2) reduces to the haptotaxis-only system and receives some attention in the literature. For example, the existence and uniqueness of local classical solutions to system (2) with $\sigma =1,\chi =\mu =\eta =0$ was proved in [34]. Further, the existence and large time behavior of global weak solutions was investigated in [26,31,35] in the case of $\eta =0$, and the existence and uniqueness of global classical solutions was considered in [36] when $\eta >0$, respectively.

In this section we focus on the global existence and large time behavior of classical solutions to the full parabolic system (2). To state these results more precisely, the initial data $(u_0,v_0,w_0)$ are assumed that for some $\vartheta \in (0,1)$,

$$\left\{\begin{array}{c}{\displaystyle u_0\in C^1\left(\overline{\mathsf{\Omega }}\right)\mathrm{with}{u}_0\ge 0\mathrm{in}\mathsf{\Omega },u_0\not\equiv 0,}\hfill \\ {\displaystyle \sigma v_0\in W^{1,\infty }\left(\mathsf{\Omega }\right)\mathrm{with}\sigma v_0\ge 0\mathrm{in}\mathsf{\Omega },}\hfill \\ {\displaystyle w_0\in C^{2+\vartheta }\left(\overline{\mathsf{\Omega }}\right)\mathrm{with}{w}_0\ge 0\mathrm{in}\overline{\mathsf{\Omega }}\mathrm{and}\frac{\partial w_0}{\partial \nu }=0on\partial \mathsf{\Omega }.}\hfill \end{array}\right.$$

(3)

For the full parabolic system (2) with $\eta =0$, we have the following theorem [32,37,38,39].

 Theorem 1. 

(i) 

(Tao) Let $n=2$, $\xi >0,\chi >0$, $r>0$, and $\mu >0$. Then, for each $(u_0,v_0,w_0)$ fulfilling (3), system (2) admits a unique global classical solution which is bounded in $\mathsf{\Omega }\times (0,\infty )$.

(ii) 

(Cao) Let $n=3$. Then, for each $(u_0,v_0,w_0)$ fulfilling (3), there exists ${\mu }_{*}:={\mu }_{*}(\chi ,\xi ,\parallel w_0{\parallel }_{W^{2,\infty }\left(\mathsf{\Omega }\right)})$ such that the conclusion of (i) holds whenever $\mu >{\mu }_{*}$.

(iii) 

(Ke and Wang) Let $3\le n\le 8$. Then, for each $(u_0,v_0,w_0)$ fulfilling (3), there exists ${\theta }_{*}:={\theta }_{*}\left(n\right)$ such that the conclusion of (i) remains valid whenever $\mu >\chi {\theta }_{*}$.

(iv) 

(Ke and Zheng) Let $n\ge 3$. Then, for each $(u_0,v_0,w_0)$ fulfilling (3), there exists ${\mu }_{*}:={\mu }_{*}(\chi ,\xi ,n,\parallel w_0{\parallel }_{W^{2,\infty }\left(\mathsf{\Omega }\right)})$ such that whenever $\mu >{\mu }_{*}$, the conclusion of (i) remains valid.

Proof. 

Thanks to the one-sided estimate for $-\Delta w$ of the form:

$$\begin{array}{cc}\hfill -\Delta w(x,t)\le & {\displaystyle -\Delta w_0\left(x\right)e^{-{\int }_0^{t}v(x,s)ds}+2e^{-{\int }_0^{t}v(x,s)ds}\nabla w_0\left(x\right)·{\int }_0^t\nabla v(x,s)ds}\hfill \\ \hfill & {\displaystyle +\frac{1}{e}{w}_0\left(x\right)+w_0\left(x\right)v(x,t)e^{-{\int }_0^{t}v(x,s)ds}}\hfill \end{array}$$

(4)

for all $x\in \mathsf{\Omega }$ and $t\in (0,T)$, which can be gained by directly solving the third equation in (2) with $\eta =0$, in the spatial two-dimensional setting one shall track the time evolution of the quantity ${\int }_{\mathsf{\Omega }}{u}^2+{\int }_{\mathsf{\Omega }}{|\nabla v|}^4$. Along with a bootstrap argument, the latter can provide a bound for ${\parallel u(·,t)\parallel }_{L^p\left(\mathsf{\Omega }\right)}$ with any $p>2$, which thus leads the boundedness of u in norm of $L^{\infty }\left(\mathsf{\Omega }\right)$ by means of a Moser-type iteration procedure and thereby completes the proof of (i).

As for the higher-dimensional case, in addition to taking advantage of (4), one utilizes the maximal Sobolev regularity properties of the heat equation to derive a bound of ${\parallel u(·,t)\parallel }_{L^p\left(\mathsf{\Omega }\right)}$ for $p>2$. Indeed, multiplying the first equation in (2) by $u^{p-1}(p>2)$ and applying the Young inequality, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{p}\frac{d}{dt}{\int }_{\mathsf{\Omega }}{u}^p+(p-1){\int }_{\mathsf{\Omega }}{u}^{p-2}{|\nabla u|}^2+\mu {\int }_{\mathsf{\Omega }}{u}^{p+1}}\hfill \\ \hfill \le & {\displaystyle (p-1)\chi {\int }_{\mathsf{\Omega }}{u}^{p-1}\nabla u·\nabla v+(p-1)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}\nabla u·\nabla w+\mu r{\int }_{\mathsf{\Omega }}{u}^p}\hfill \\ \hfill =& {\displaystyle -\frac{p-1}{p}\chi {\int }_{\mathsf{\Omega }}{u}^p\Delta v+(p-1)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}\nabla u·\nabla w+\mu r{\int }_{\mathsf{\Omega }}{u}^p}\hfill \\ \hfill \le & {\displaystyle {\epsilon }_1{\int }_{\mathsf{\Omega }}{u}^{p+1}+{\epsilon }_1^{-p}{\chi }^{p+1}{\int }_{\mathsf{\Omega }}{|\Delta v|}^{p+1}+(p-1)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}\nabla u·\nabla w+\mu r{\int }_{\mathsf{\Omega }}{u}^p.}\hfill \end{array}$$

(5)

On the other hand, in view of (4), one can find $c_1>0,c_2>0$ only depending upon $w_0$ such that

$$\begin{array}{cc}\hfill & {\displaystyle \left(p-1\right)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}\nabla u·\nabla w}\hfill \\ \hfill =& {\displaystyle -\frac{p-1}{p}\xi {\int }_{\mathsf{\Omega }}{u}^p\Delta w}\hfill \\ \hfill \le & {\displaystyle -\frac{p-1}{p}\xi {\int }_{\mathsf{\Omega }}{u}^p\left(\Delta w_0\left(x\right)e^{-{\int }_0^{t}v(x,s)ds}-2e^{-{\int }_0^{t}v(x,s)ds}\nabla w_0\left(x\right)·{\int }_0^t\nabla v(x,s)ds-\right.}\hfill \\ \hfill & \left.{\displaystyle \frac{1}{e}{w}_0\left(x\right)-w_0\left(x\right)v(x,t)e^{-{\int }_0^{t}v(x,s)ds}}\right)\hfill \\ \hfill \le & {\displaystyle c_1\xi {\int }_{\mathsf{\Omega }}{u}^p+c_1\xi {\int }_{\mathsf{\Omega }}{u}^{p}v-\frac{2(p-1)}{p}\xi {\int }_{\mathsf{\Omega }}{u}^p\nabla e^{-{\int }_0^{t}v(x,s)ds}·\nabla w_0\left(x\right)}\hfill \\ \hfill =& {\displaystyle c_1\xi {\int }_{\mathsf{\Omega }}{u}^p+c_1\xi {\int }_{\mathsf{\Omega }}{u}^{p}v+\frac{2(p-1)}{p}\xi {\int }_{\mathsf{\Omega }}(u^p{e}^{-{\int }_0^{t}v(x,s)ds}▵w_0\left(x\right)+\nabla u^p·\nabla w_0\left(x\right)e^{-{\int }_0^{t}v(x,s)ds})}\hfill \\ \hfill \le & {\displaystyle c_2\xi {\int }_{\mathsf{\Omega }}{u}^p+c_1\xi {\int }_{\mathsf{\Omega }}{u}^{p}v+c_2(p-1)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}|\nabla u|.}\hfill \end{array}$$

(6)

Inserting (6) into (5) yields

$$\begin{array}{cc}\hfill & {\displaystyle \frac{1}{p}\frac{d}{dt}{\int }_{\mathsf{\Omega }}{u}^p+(p-1){\int }_{\mathsf{\Omega }}{u}^{p-2}{|\nabla u|}^2+\mu {\int }_{\mathsf{\Omega }}{u}^{p+1}}\hfill \\ \hfill \le & {\displaystyle {\displaystyle {\epsilon }_1{\int }_{\mathsf{\Omega }}{u}^{p+1}+{\epsilon }_1^{-p}{\chi }^{p+1}{\int }_{\mathsf{\Omega }}{|\Delta v|}^{p+1}+(1+c_2\xi ){\int }_{\mathsf{\Omega }}{u}^p+c_1\xi {\int }_{\mathsf{\Omega }}{u}^{p}v+c_2(p-1)\xi {\int }_{\mathsf{\Omega }}{u}^{p-1}|\nabla u|}}\hfill \\ \hfill \le & {\displaystyle ({\epsilon }_1+{\epsilon }_2){\int }_{\mathsf{\Omega }}{u}^{p+1}+{\epsilon }_1^{-p}{\chi }^{p+1}{\int }_{\mathsf{\Omega }}{|\Delta v|}^{p+1}+c_3({\epsilon }_2,p,\xi )+{\epsilon }_2^{-p}{\left(c_1\xi \right)}^{p+1}{\int }_{\mathsf{\Omega }}{v}^{p+1}+\frac{p-1}{2}{\int }_{\mathsf{\Omega }}{u}^{p-2}{|\nabla u|}^2}\hfill \end{array}$$

which together with the Young inequality leads to

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d}{dt}{\int }_{\mathsf{\Omega }}{u}^p+(p+1){\int }_{\mathsf{\Omega }}{u}^p+\frac{p(p-1)}{2}{\int }_{\mathsf{\Omega }}{u}^{p-2}{|\nabla u|}^2+p(\mu -2{\epsilon }_2-{\epsilon }_1){\int }_{\mathsf{\Omega }}{u}^{p+1}}\hfill \\ \hfill \le & {\displaystyle {\epsilon }_1^{-p}p{\chi }^{p+1}{\int }_{\mathsf{\Omega }}{|\Delta v|}^{p+1}+c_4({\epsilon }_2,p,\xi )+{\epsilon }_2^{-p}p{\left(c_1\xi \right)}^{p+1}{\parallel v(·,t)\parallel }_{L^{p+1}\left(\mathsf{\Omega }\right)}^{p+1}.}\hfill \end{array}$$

(7)

According to the maximal Sobolev regularity properties of the Neumann heat equation, one can see that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^t{\int }_{\mathsf{\Omega }}{e}^{(p+1)s}{\left(\right|\Delta v(x,s)|}^{p+1}+{\left|v(x,s)\right|}^{p+1})dxds}\\ \hfill {\displaystyle \le c_3{\int }_0^t{\int }_{\mathsf{\Omega }}{e}^{(p+1)s}{\left|u(x,s)\right|}^{p+1}dxds+c_3{\parallel v_0\parallel }_{W^{2,p+1}\left(\mathsf{\Omega }\right)}.}\end{array}$$

(8)

Therefore, combining (7) with (8), one can establish a bound for ${\int }_{\mathsf{\Omega }}{u}^p$ for some $p>\frac{n}{2}$, which turns out to be a starting point for a bootstrap procedure to achieve the global boundedness of solutions. □

Beyond the global boundedness stated above, the large time behavior of solutions to (2) have been achieved under some conditions on the model parameters thereof, which implies that although haptotaxis mechanisms may have some important influence on the properties of (2), they will become extinct asymptotically [32,39,40,41]. In particular, under an explicit condition which is independent of all further model ingredients such as $\xi$, the spatial domain or the initial data, the asymptotic behavior of all solution components in (2) is derived in [41].

 Theorem 2. 

(Tao and Winkler) Let $\mathsf{\Omega }\subset {\mathbb{R}}^n(n\le 3)$ be a bounded convex domain, $\xi >0$, $\chi >0$, $\mu >\frac{{\chi }^2}{8}$, and assume that $(u,v,w)$ is a bounded global classical solution of (2) with initial data satisfying (3) as well as

$$\sqrt{w_0}\in W^{1,\infty }\left(\mathsf{\Omega }\right),w_0\left(x\right)\le rforallx\in \mathsf{\Omega }.$$

Then, there exist $\lambda >0$ and $c>0$ such that

$${\parallel u(·,t)-r\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}{+\parallel v(·,t)-r\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}+{\parallel w(·,t)\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}\le c{e}^{-\lambda t}$$

(9)

for all $t>0$.

Proof. 

The proof can be divided into two steps. Firstly, under the condition $w_0\left(x\right)\le r$, there exist $c_1>0$ and $\kappa >0$ such that

$${\parallel w(·,t)\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}\le c_1{e}^{-\kappa t}$$

for all $t>0$. To achieve this, by a straightforward testing procedure, one finds some $p\in (0,1)$ and $c_2>0$ such that

$${\int }_{\mathsf{\Omega }}{u}^p(x,t)dx>c_2\mathrm{for}\mathrm{all}t>1,$$

which yields a lower bound for the total mass ${\int }_{\mathsf{\Omega }}u(x,t)dx$.

At this position, the proof of (9) shows that the hypothesis $\mu >\frac{{\chi }^2}{8}$ warrants that

$$\mathcal{F}\left(t\right):={\int }_{\mathsf{\Omega }}(u-r-rln\frac{u}{r})+\alpha {\int }_{\mathsf{\Omega }}{(v-r)}^2+\beta {\int }_{\mathsf{\Omega }}{|\nabla w|}^2+\gamma {\int }_{\mathsf{\Omega }}{w}^2+\delta {\int }_{\mathsf{\Omega }}w$$

for all $t>0$, and acts as a Lyapunov functional for (2) under appropriate choices of the positive constants $\alpha ,\beta ,\gamma$, and $\delta$. By means of an analysis of the corresponding energy inequality, one can first establish the mere convergence of $(u,v)$ to $(r,r)$ with respect to the spatial $L^{\infty }\left(\mathsf{\Omega }\right)$ norm. Accordingly, making essential use of interpolation argument based on the latter and respective higher regularity properties of the solution, we thereby prove that the convergence actually takes place at an exponential rate. □

Theorem 2 indicates that, although the behavior of solutions to (2) can be affected by two taxis mechanisms on intermediate time scales, the destabilizing effects thereof are substantially overbalanced by the zero-order dissipative action of logistic damping when $\mu$ is suitably large as related to $\chi >0$.

Toward the understanding of possible effects that tissue remodeling may have on the qualitative behavior of solutions to (2), we turn to consider (2) with $\eta >0$.

Compared to the case of $\eta =0$, the additional mathematical challenges stem from the coupling between w and the crucial quantity u in the third equation of system (2) when $\eta >0$. Recently, the global solvability for the two-dimensional system (2) with $\sigma =0,\eta >0$ was addressed in [30].

 Theorem 3. 

(Tao and Winkler) Let $\mathsf{\Omega }\subset {\mathbb{R}}^2$, $\xi >0$, $\chi >0$, $\mu >0$, $\eta >0$ and $\sigma =0$. Then, for each $(u_0,w_0)$ fulfilling

$$\left\{\begin{array}{c}{\displaystyle u_0\in C^{2+\vartheta }\left(\overline{\mathsf{\Omega }}\right)with{u}_0>0in\mathsf{\Omega },}\hfill \\ {\displaystyle w_0\in C^{2+\vartheta }\left(\overline{\mathsf{\Omega }}\right)with{w}_0>0in\overline{\mathsf{\Omega }},}\hfill \\ {\displaystyle \frac{\partial u_0}{\partial \nu }-\xi u_0\frac{\partial w_0}{\partial \nu }=0on\partial \mathsf{\Omega }}\hfill \end{array}\right.$$

with some $\vartheta \in (0,1)$, the system (2) admits a unique global classical solution $(u,v,w)$ with $u>0,v>0$ and $0<w\le max\{1,\parallel w_0{\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}\}$.

Proof. 

The subsential issue consists of taking advantage of the dampening effect of $-\eta vw$ in the w- equation of (2) to derive an energy-like inequality, which yields an a priori estimate of ${\int }_{\mathsf{\Omega }}ulnu$ in bounded time intervals and also provides $c\left(T\right)>0$ such that ${\int }_{\mathsf{\Omega }}{|\nabla v(·,t)|}^2\le c\left(T\right)$ for all $t<T$ with the help of a result from regularity theory of elliptic equations. The latter can act as a starting point for an iterative bootstrap argument used to derive higher regularity estimates. □

As for the global boundedness of the full parabolic model (2) with $\eta >0$, thanks to $L^q-L^p$ estimates for the Neumann heat semigroup, the authors of [28] can deal with the chemotaxis-related integral term ${\int }_0^t{\int }_{\mathsf{\Omega }}{e}^{-(p+1)(t-s)}{u}^p{|\nabla v|}^{2}dxds$, and thereby derive the global boundedness of the first component of solution $(u,v,w)$ when $\mu$ is sufficiently large. It is noticed that very recently, the corresponding results of [28,30] have been improved in [24,42], which can be stated as follows.

 Theorem 4. 

(Jin, Pang and Wang) Let $\mathsf{\Omega }\subset {\mathbb{R}}^2$, $\chi >0,\xi >0$, $r>0$, $\sigma \in \{0,1\}$, and $\eta >0$, and suppose that $(u_0,v_0,w_0)$ fulfils (3). Then for any $\mu >0$, the problem (2) admits a unique global classical solution $(u,v,w)$, where ${\parallel u(·,t)\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}$ is uniformly bounded for $t\in (0,\infty )$.

Proof. 

In the case of $\sigma =1$, the crucial idea of the proof is to discover that the integral of the form ${\int }_{\mathsf{\Omega }}{e}^{\xi w}{a}^2$ with $a=e^{-\xi w}u$ enjoys a certain dissipative property. In fact, it is observed that a certain variant of the latter satisfies

$${\displaystyle \frac{d}{dt}{\int }_{\mathsf{\Omega }}{e}^{\xi w}{a}^2+\frac{1}{\epsilon }{\int }_{\mathsf{\Omega }}{e}^{\xi w}{a}^2\le c_1{(\parallel \Delta v\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+{\parallel a\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2){\int }_{\mathsf{\Omega }}{e}^{\xi w}{a}^2+c_2\left(\epsilon \right)}$$

with some constants $c_1>0,c_2\left(\epsilon \right)>0$ for any $\epsilon >0$. The latter will provide a bound for ${\int }_{\mathsf{\Omega }}{u}^2(·,t)$, which forms the cornerstone for the derivation of higher regularity estimates of solutions, inter alia the global bound for ${\parallel u(·,t)\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}$.

As for the case of $\sigma =0$, the initial but crucial step is to derive the estimate for u in $LlogL$, which is a consequences of the following inequality:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d}{dt}{\int }_{\mathsf{\Omega }}{e}^{\xi w}alna+\frac{\mu }{2}{\int }_{\mathsf{\Omega }}{e}^{\xi w}alna}\hfill \\ \hfill \le & {\displaystyle {\epsilon c(\parallel \Delta v\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+{\parallel \nabla v\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2){\int }_{\mathsf{\Omega }}{e}^{\xi w}alna+{c(\parallel \Delta v\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+{\parallel \nabla v\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2)+c\left(\epsilon \right)}\hfill \end{array}$$

with some $c>0,c\left(\epsilon \right)>0$ for all $\epsilon >0$. □

In the spatial three-dimensional setting, the global solvability of the full parabolic system (2) with tissue remodeling is much more delicate. Indeed, very little is known for this higher-dimensional full parabolic model, and as far as we know, the only result is presented in the survey paper [1], where a certain global weak solution of (2) with $\sigma =1$ was constructed. In this direction, under the smallness restriction on the growth rate r, the global boundedness of solutions was recently addressed in [42]. A natural question is whether $r_0>0$ obtained there is optimal for the global existence of classical solutions, and if the weak solution constructed in [1] is eventually smooth.

 Theorem 5. 

(Pang and Wang) Let $\mathsf{\Omega }\subset {\mathbb{R}}^3$ be a bounded convex domain with smooth boundary and suppose that $\chi >0$, $\xi >0$, $\eta >0$, and $\mu >0$. Then, there exists $r_0>0$ with the property such that for any $r<r_0$, the problem (2) admits a unique global classical solution provided that $\parallel u_0{\parallel }_{L^2\left(\mathsf{\Omega }\right)}$ and $\parallel v_0{\parallel }_{W^{1,4}\left(\mathsf{\Omega }\right)}$ are appropriately small.

Proof. 

On the basis of the mass evolution of solutions to (2), one can verify that the quantity ${\int }_{\mathsf{\Omega }}{a}^2(·,t)+{\int }_{\mathsf{\Omega }}{|\nabla v(·,t)|}^4$ with $a=e^{-\xi w}u$ satisfies a differential inequality. Accordingly, thanks to the comparison argument of the respective ordinary differential equation, u is indeed bounded provided that initial data and $r>0$ are appropriately small. In the context of some straightforward $L^p$ testing procedure, the latter forms a cornerstone for the bootstrap argument to yield a bound for u in $L^{\infty }\left(\mathsf{\Omega }\right)$. □

3. Multiscale Cancer Invasion Models

It is well-established that the macroscopic behavior of tumor cells is influenced by the internal state of cells, hence by microscopic processes taking place on the subcellular scale such as receptor binding to chemoattractants or adhesion molecules initiating intracellular signaling pathways. As far as we know, there are several types of multiscale cancer invasion models which connect the macroscopic evolution of cells with the microscopic dynamics of (some of) their subcellular events, and thus particularly lead to the system of PDEs strongly coupled with ODEs [12,18,22,23,43]. Due to the feature of the multiscale models in which different types of equations are coupled in a highly nonlinear way, the problem of well-posedness thereof seems to be a challenge [22,23,33,44]. Here we briefly review the existence of global weak solutions to a go-or-grow multiscale system for tumor invasion with therapy given by

$$\left\{\begin{array}{c}{\displaystyle m_t=\nabla ·(D(m,q,v)\kappa \nabla m)-\nabla ·(\frac{\kappa vm}{1+v}\nabla v)+\lambda \left(y\right)q-\gamma \left(y\right)m-r_m\left(t\right)m,}\hfill \\ {\displaystyle q_t={\mu }_{q}q(1-\frac{m+q}{K_c}-{\eta }_1\frac{v}{K_v})-\lambda \left(y\right)q+\gamma \left(y\right)m-r_q\left(t\right)q,}\hfill \\ {\displaystyle v_t=-{\delta }_{v}mv+{\mu }_{v}v(1-{\eta }_2\frac{m+q}{K_c}-\frac{v}{K_v})-r_v\left(t\right)v,}\hfill \\ y_t=-k_{-1}\left(d_c\right)y+k_1\left(d_c\right)(1-y)v,\hfill \\ {\kappa }_t=-{\delta }_{\kappa }\kappa +H\left(y(·,t-\tau )\right)\hfill \end{array}\right.$$

(10)

supplemented by initial conditions

$$m(x,0)=m_0\left(x\right),q(x,0)=q_0\left(x\right),v(x,0)=v_0\left(x\right),\kappa (x,0)={\kappa }_0\left(x\right),x\in \mathsf{\Omega },$$

(11)

$$y(x,t)=y_0(x,t),x\in \mathsf{\Omega },t\in [-\tau ,0].$$

(12)

and no-flux boundary conditions

$$D(m,q,v)\kappa {\partial }_{\nu }m-\frac{\kappa v}{1+v}m{\partial }_{\nu }v=0,x\in \partial \mathsf{\Omega },t>0.$$

(13)

where $\mathsf{\Omega }\subset {\mathbb{R}}^n(n\le 3)$ is a bounded domain with smooth boundary and $\nu$ denotes the outward unit normal on $\partial \mathsf{\Omega }$. Here m and q denote the densities of migrating cancer cells and proliferating cancer cells, respectively, v is the density of tissue fibers in the ECM, y denotes the concentration of integrins bound to ECM fibers, and $\kappa$ represents the contractility function of cancer cells. $\lambda \left(y\right)$ denotes the rate at which the proliferating cells cease their proliferation in the cycle and advance to migration phases, while $\gamma \left(y\right)$ is the rate at which the moving cells enter a proliferative state. It is natural to assume that these rates are influenced by subcellular dynamics, featured by the amount of cell surface receptors bound to insoluble ligands in the tumor microenvironment. ${\delta }_v>0$ is the decay rate of ECM due to interactions with motile cells, and ${\mu }_q$ and ${\mu }_v$ are growth rates for the tumor cells and the tissue, respectively.

Concerning the initial data, we suppose that

$$m_0\in C\left(\overline{\mathsf{\Omega }}\right),q_0,v_0\in W^{1,2}\left(\mathsf{\Omega }\right)\cap \left(\overline{\mathsf{\Omega }}\right),{\kappa }_0\in W^{1,4}\left(\mathsf{\Omega }\right),y_0\in C([-\tau ,0];W^{1,4}\left(\mathsf{\Omega }\right))$$

satisfy

$$m_0\ge 0,q_0\ge 0,v_0\ge 0,{\kappa }_0\ge 0\mathrm{in}\mathsf{\Omega }$$

as well as

$$0\le y_0\le 1\mathrm{in}\mathsf{\Omega }\times [-\tau ,0].$$

Furthermore, it is assumed that for any $A>0$ and $L>0$ there exist positive constants $C_1$ and $C_2$ such that

$$D\in C^3\left({[0,\infty )}^3\right)\cap W^{2,\infty }([0,\infty )\times [0,A]\times [0,L]),\lambda \in C^1\left([0,1]\right),H\in C^3\left([0,1]\right),$$

$$r_i\in C^1\left([0,\infty ]\right),i\in \{m,q,v\},k_j\in C^1\left([0,\infty ]\right),j\in \{-1,1\},$$

$$0<C_1\le D(m,q,v)\le C_2\mathrm{for}\mathrm{all}(m,q,v)\in ([0,\infty )\times [0,A]\times [0,L]),$$

$$0<{\lambda }_1\le \lambda \left(y\right)\le {\lambda }_2,\gamma \left(y\right)=\gamma ,H\left(y\right)\ge 0\mathrm{for}\mathrm{all}y\in [0,1],$$

$$0\le r_i\left(t\right)\le C_4,0<C_3\le k_j\left(t\right)\le C_4,$$

with positive constants $C_i$, ${\lambda }_i$.

Under the above assumptions, the global existence of weak solutions to the system (10)–(13) in the following sense was established in [23].

 Definition 1. 

(Weak solution) Let $T\in (0,\infty )$. A weak solution to (10)–(13) consists of nonnegative functions $m\in L^1((0,T);L^2\left(\mathsf{\Omega }\right))$ with $\sqrt{1+m}\in L^2((0,T);W^{1,2}\left(\mathsf{\Omega }\right))$ and $\sqrt{m}\nabla v\in L^2(\mathsf{\Omega }\times (0,T)),v\in L^{\infty }(\mathsf{\Omega }\times (0,T))\cap L^2((0,T);W^{1,2}\left(\mathsf{\Omega }\right))$, $q,\kappa \in L^{\infty }(\mathsf{\Omega }\times (0,T)),y\in L^{\infty }(\mathsf{\Omega }\times (-\tau ,T))$ such that for all $\phi \in C_0^{\infty }(\overline{\mathsf{\Omega }}\times (0,T))$

$$\begin{array}{cc}\hfill {\displaystyle -{\int }_0^T{\int }_{\mathsf{\Omega }}m{\partial }_t\phi -{\int }_0^T{\int }_{\mathsf{\Omega }}{m}_0\phi (·,0)=}& {\displaystyle -2{\int }_0^T{\int }_{\mathsf{\Omega }}D(m,q,v)\kappa \sqrt{1+m}\nabla \sqrt{1+m}\nabla \phi }\hfill \\ \hfill {\displaystyle +{\int }_0^T{\int }_{\mathsf{\Omega }}\frac{\kappa v}{1+v}m\nabla v\nabla \phi }& {\displaystyle +{\int }_0^T{\int }_{\mathsf{\Omega }}(\lambda \left(y\right)q-\gamma m-r_m\left(t\right)m)\phi ,}\hfill \end{array}$$

(14)

$${\displaystyle -{\int }_0^T{\int }_{\mathsf{\Omega }}q{\partial }_t\phi -{\int }_0^T{\int }_{\mathsf{\Omega }}{q}_0\phi (·,0)={\int }_0^T{\int }_{\mathsf{\Omega }}\{{\mu }_{q}q(1-\frac{m+q}{K_c}-{\eta }_1\frac{v}{K_v})-\lambda \left(y\right)q+\gamma m-r_q\left(t\right)q\}\phi ,}$$

(15)

$${\displaystyle -{\int }_0^T{\int }_{\mathsf{\Omega }}v{\partial }_t\phi -{\int }_0^T{\int }_{\mathsf{\Omega }}{v}_0\phi (·,0)={\int }_0^T{\int }_{\mathsf{\Omega }}\{-{\delta }_{v}mv+{\mu }_{v}v(1-{\eta }_2\frac{m+q}{K_c}-\frac{v}{K_v})-r_v\left(t\right)v\}\phi ,}$$

(16)

$${\displaystyle -{\int }_0^T{\int }_{\mathsf{\Omega }}y{\partial }_t\phi -{\int }_0^T{\int }_{\mathsf{\Omega }}{y}_0\phi (·,0)={\int }_0^T{\int }_{\mathsf{\Omega }}\{k_1\left(d_c\right)(1-y)v-k_{-1}\left(d_c\right)y\}\phi }$$

(17)

as well as

$${\displaystyle -{\int }_0^T{\int }_{\mathsf{\Omega }}\kappa {\partial }_t\phi -{\int }_0^T{\int }_{\mathsf{\Omega }}{\kappa }_0\phi (·,0)={\int }_0^T{\int }_{\mathsf{\Omega }}\{-{\delta }_{\kappa }\kappa +H\left(y(·,t-\tau )\right)\}\phi }$$

(18)

are fulfilled. A quintuple of functions $(m,q,v,y,\kappa )$is called the global weak solution of (10)–(13) if it is a weak solution in $\mathsf{\Omega }\times (0,T)$ for all $T>0$.

Resorting to the weak solution concepts above, Stinner et al. [23] obtained the following global solvability result. However, the global solvability of (10)–(13) with the non-constant $\gamma \left(y\right)$ is still an open issue.

 Theorem 6.> 

(Stinner, Surulescu, and Uatay) Let $n\le 3$. Then, (10)–(13) possess a global weak solution in the sense of Definition 1 with additional properties that $m\in L^{\infty }((0,\infty ),L^1\left(\mathsf{\Omega }\right)),q,v,\kappa \in L^{\infty }(\mathsf{\Omega }\times (0,\infty ))$, and $y\in L^{\infty }(\mathsf{\Omega }\times (-\tau ,\infty ))$.

Proof. 

The weak solution is constructed as the limit of global smooth solutions to the regularized problems:

$$\left\{\begin{array}{c}{\displaystyle m_{\epsilon t}=\nabla ·(D(m_{\epsilon },q_{\epsilon },v_{\epsilon }){\kappa }_{\epsilon }\nabla m_{\epsilon })-\nabla ·(\frac{{\kappa }_{\epsilon }{v}_{\epsilon }}{1+v_{\epsilon }}{m}_{\epsilon }\nabla v_{\epsilon })}\hfill \\ +\lambda \left(y_{\epsilon }\right)q_{\epsilon }-\gamma m_{\epsilon }-r_m\left(t\right)m_{\epsilon }-\epsilon m_{\epsilon }^{\theta },x\in \mathsf{\Omega },t>0,\hfill \\ {\displaystyle q_{\epsilon t}=\epsilon \Delta q_{\epsilon }+{\mu }_q{q}_{\epsilon }(1-\frac{m_{\epsilon }+q_{\epsilon }}{K_c}-{\eta }_1\frac{v_{\epsilon }}{K_v})-\lambda \left(y_{\epsilon }\right)q_{\epsilon }+\gamma m_{\epsilon }-r_q\left(t\right)q_{\epsilon },x\in \mathsf{\Omega },t>0,}\hfill \\ {\displaystyle v_{\epsilon t}=\epsilon \Delta v_{\epsilon }-{\delta }_v{m}_{\epsilon }{v}_{\epsilon }+{\mu }_v{v}_{\epsilon }(1-{\eta }_2\frac{m_{\epsilon }+q_{\epsilon }}{K_c}-\frac{v_{\epsilon }}{K_v})-r_v\left(t\right)v_{\epsilon },x\in \mathsf{\Omega },t>0.}\hfill \\ y_{\epsilon t}=k_1\left(d_c\right)(1-y_{\epsilon })v_{\epsilon }-k_{-1}\left(d_c\right)y_{\epsilon },x\in \mathsf{\Omega },t>0\hfill \\ {\kappa }_{\epsilon t}=-{\delta }_{\kappa }{\kappa }_{\epsilon }+H\left(y_{\epsilon }(·,t-\tau )\right),x\in \mathsf{\Omega },t>0\hfill \\ {\partial }_{\nu }{m}_{\epsilon }={\partial }_{\nu }{q}_{\epsilon }={\partial }_{\nu }{v}_{\epsilon }=0,x\in \partial \mathsf{\Omega },t>0,\hfill \\ m_{\epsilon }(x,0)=m_{0\epsilon }\left(x\right),q_{\epsilon }(x,0)=q_{0\epsilon }\left(x\right),v_{\epsilon }(x,0)=v_{0\epsilon }\left(x\right),{\kappa }_{\epsilon }(x,0)={\kappa }_0\left(x\right),x\in \mathsf{\Omega },\hfill \\ y_{\epsilon }(x,t)=y_{0\epsilon }(x,t),x\in \mathsf{\Omega },t\in [-\tau ,0],\hfill \end{array}\right.$$

(19)

with $\epsilon \in (0,1)$ and $\theta >max\{n,2\}$ and $m_{0\epsilon },q_{0\epsilon },v_{0\epsilon },{\kappa }_{0\epsilon }$ and $y_{0\epsilon }$ the regularization of the respective initial data. In view of $\theta >max\{n,2\}$ and well-known results on maximal Sobolev regularity to the third equation in (19), these approximate problems are globally solvable. The essential step toward the existence of global weak solutions is to establish $\epsilon$-independent a priori estimates stemming from tracking the evolution of the functional

$${\mathcal{F}}_{\epsilon }\left(t\right):={\int }_{\mathsf{\Omega }}{m}_{\epsilon }ln{m}_{\epsilon }+{\int }_{\mathsf{\Omega }}\frac{{\kappa }_{\epsilon }{|\nabla v_{\epsilon }|}^2}{1+v_{\epsilon }}+{\int }_{\mathsf{\Omega }}|\nabla q_{\epsilon }{|}^2+{\int }_{\mathsf{\Omega }}{|\nabla y_{\epsilon }|}^2$$

which in fact satisfies the inequality

$$\frac{d}{dt}{\mathcal{F}}_{\epsilon }\left(t\right)+{\mathcal{D}}_{\epsilon }\left(t\right)\le c_1\left(T\right)({\mathcal{F}}_{\epsilon }\left(t\right)+\epsilon h_{\epsilon }(t-\tau ))\mathrm{for}\mathrm{all}t\in (0,T)$$

with

$${\mathcal{D}}_{\epsilon }\left(t\right):=\frac{1}{2}{\int }_{\mathsf{\Omega }}D(m_{\epsilon },q_{\epsilon },v_{\epsilon }){\kappa }_{\epsilon }\frac{|\nabla m_{\epsilon }{|}^2}{m_{\epsilon }}+{\int }_{\mathsf{\Omega }}\frac{m_{\epsilon }{\kappa }_{\epsilon }{|\nabla v_{\epsilon }|}^2}{{(1+v_{\epsilon })}^2}+\frac{\epsilon }{2}{\int }_{\mathsf{\Omega }}{m}_{\epsilon }^{\theta }ln(m_{\epsilon }+2),$$

$$h_{\epsilon }\left(t\right):={\int }_{\mathsf{\Omega }}{|\nabla y_{\epsilon }|}^4$$

satisfying ${\displaystyle \epsilon h_{\epsilon }(t-\epsilon )\le c_2\left(T\right){\mathcal{F}}_{\epsilon }(t-\epsilon )+\underset{s\in [-\tau ,0]}{sup}{\int }_{\mathsf{\Omega }}{|\nabla y_{\epsilon }(·,s)|}^4}$ for suitable $c_1\left(T\right)>0,c_2\left(T\right)>0$. From the entropy-type function ${\mathcal{F}}_{\epsilon }$, one can derive the a priori estimates which provide suitable compactness properties of the approximate solution families and thereby allow for extracting subsequences which convergence to a global weak solution in the desired sense. □

4. Conclusions

This paper describes recent results of some chemotaxis–haptotaxis models, inter alia macro cancer invasion models proposed by Chaplain and Lolas in [12,13] and the multiscale cancer invasion models by Stinner et al. in [20].

It is observed that in the case of $\eta =0$, one-sided pointwise bound for $-▵w$ in the flavour of (4) plays an essential role in the analysis of (2) not only at the stage of the global boundedness of solutions, but also in the description of large time behavior, whereas $\eta >0$ apparently makes (4) inaccessible. Thanks to the variable transformation $a=e^{-\xi w}u$, one obtains the global boundedness of the two-dimensional version of (2) with $\eta >0$, and even its three-dimensional version under a smallness condition on initial data and growth rate. In synopsis of the above results, the naturally arising problem consists of determining whether the initial boundary value problem for the higher-dimensional model (2) possesses a certain generalized solution which becomes eventually smooth and approaches a spatially homogeneous steady state, and to which extent the nonlinear diffusion of cancer cells may influence the solution behavior when the the ECM remodeling is taken into account.

For the multiscale cancer invasion models (10)–(13) with the constant $\gamma \left(y\right)$, Stinner et al. [23] established the global solvability in the framework of weak solutions. However, from a biological point of view, the rate $\gamma \left(y\right)$ at which the moving cells start proliferating should depend on the concentration of integrins bound to ECM fibers, rather than the constant considered in [23]. On the other hand, it is recognized that tumor cell motility is triggered (among others) by cancer cell population growth and acidification due to excessive glycolysis, and a kind of repellent taxis of the form $+\nabla ·\left(\chi \right(m,q,v)m\nabla q)$ should be added into the first equation of (10). Therefore, it is quite interesting to explore the qualitative behavior of solutions to the corresponding initial boundary value problem for (10)–(13).