Does indirectness of signal production reduce the explosion-supporting potential in chemotaxis-haptotaxis systems? Global classical solvability in a class of models for cancer invasion (and more)

TL;DR

This paper introduces a class of chemotaxis-haptotaxis models with indirect signal production, proving global classical solvability in 2D and showing that indirectness suppresses blow-up phenomena.

Contribution

It presents the first rigorous proof of global classical solutions for a broad class of indirect signal production chemotaxis models in two dimensions.

Findings

Global classical solutions exist in 2D under mild conditions.

Indirect signal production reduces blow-up potential in chemotaxis systems.

The model applies to cancer invasion and other biological processes.

Abstract

We propose and study a class of parabolic-ODE models involving chemotaxis and haptotaxis of a species following signals indirectly produced by another, non-motile one. The setting is motivated by cancer invasion mediated by interactions with the tumor microenvironment, but has much wider applicability, being able to comprise descriptions of biologically quite different problems. As a main mathematical feature consituting a core difference to both classical Keller-Segel chemotaxis systems and Chaplain-Lolas type chemotaxis-haptotaxis systems, the considered model accounts for certain types of indirect signal production mechanisms. The main results assert unique global classical solvability under suitably mild assumptions on the system parameter functions in associated spatially two-dimensional initial-boundary value problems. In particular, this rigorously confirms that at least in two-dimensional settings, the considered indirectness in signal production induces a significant blow-up suppressing tendency also in taxis systems substantially more general than some particular examples for which corresponding effects have recently been observed.

Figures

Equations

$$\left\{\begin{array}[]{l}u_{t}=\Delta u-\chi\nabla\cdot(u\nabla h)-\xi\nabla\cdot(u\nabla v)+\mu u(1-u-v-w),\\[2.84526pt] h_{t}=\Delta h-h+\alpha w,\\[2.84526pt] v_{t}=-hv+\eta v(1-u-v)+\beta\frac{w}{1+w},\\[2.84526pt] w_{t}=u,\end{array}\right.$$

$$\left\{\begin{array}[]{l}u_{t}=D_{u}\Delta u+\chi\nabla\cdot(u\nabla h)-\xi\nabla\cdot(u\nabla v)-k_{1}u+k_{2}\frac{hw}{1+hw},\\[2.84526pt] h_{t}=D_{h}\Delta h+k_{3}w-k_{4}h,\\[2.84526pt] v_{t}=-k_{5}(h+u+w)v+k_{6}v(1-v)_{+}^{2},\\[2.84526pt] w_{t}=k_{7}u-k_{8}hw+k_{9}w(1-u-v-w)_{+}^{2},\end{array}\right.$$

$$\left\{\begin{array}[]{ll}u_{t}=D_{u}\Delta u-\chi\nabla\cdot(u\nabla h)-\xi\nabla\cdot(u\nabla v)+f(u,v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] h_{t}=D_{h}\Delta h+g(u,v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] v_{t}=-\alpha uv+v\phi(u,v,w,h)+\Phi(w),&x\in\Omega,\ t>0,\\[2.84526pt] w_{t}=\beta u+w\psi(u,v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] D_{u}\frac{\partial u}{\partial\nu}-\xi u\frac{\partial v}{\partial\nu}=\frac{\partial h}{\partial\nu}=0,&x\in\partial\Omega,\ t>0,\\[2.84526pt] u(x,0)=u_{0}(x),\quad h(x,0)=h_{0}(x),\quad v(x,0)=v_{0}(x),\quad w(x,0)=w_{0}(x),&x\in\Omega,\end{array}\right.$$

$$\left\{\begin{array}[]{l}u_{0}\in C^{2+\vartheta}(\overline{\Omega})\quad\mbox{is nonnegative with }u_{0}\not\equiv 0,\\[2.84526pt] h_{0}\in C^{2+\vartheta}(\overline{\Omega})\quad\mbox{is nonnegative},\\[2.84526pt] v_{0}\in C^{2+\vartheta}(\overline{\Omega})\quad\mbox{is positive in $\overline{\Omega}$ with $\frac{\partial v_{0}}{\partial\nu}=0$ on $\partial\Omega$, and that}\\[2.84526pt] w_{0}\in C^{2+\vartheta}(\overline{\Omega})\quad\mbox{is positive in $\overline{\Omega}$.}\end{array}\right.$$

$$f,g,\varphi \text{mbox}and\psi \text{mbox}belongto{C}^1([0,\infty {)}^4)\text{mbox}and\Phi \in C^1([0,\infty )),$$

$$\left\{\begin{array}[]{l}f_{0}(u)\leq f(u,v,w,h)\leq C_{f}(v)\cdot(u+w+1)\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\\[2.84526pt] \mbox{with some $f_{0}\in C^{1}([0,\infty))$ such that $f_{0}(0)\geq 0$, and some nondecreasing $C_{f}:[0,\infty)\to(0,\infty)$}\end{array}\right.$$

$$\left\{\begin{array}[]{l}|g(u,v,w,h)|\leq C_{g}(v)\cdot(w+h+1)\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\\[2.84526pt] \mbox{with some nondecreasing $C_{g}:[0,\infty)\to(0,\infty)$ \quad and}\\[2.84526pt] g(0,v,w,0)\geq 0\quad\mbox{for all }(v,w)\in[0,\infty)^{2},\end{array}\right.$$

$$\left\{\begin{array}[]{l}\phi(u,v,w,h)\leq-c_{\phi}\cdot w+C_{\phi}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\phi_{u}(u,v,w,h)|\leq\frac{C_{\phi}}{\sqrt{uv+1}}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\phi_{v}(u,v,w,h)|\leq\frac{C_{\phi}}{v+1}+C_{\phi}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\phi_{w}(u,v,w,h)|\leq\frac{C_{\phi}}{\sqrt{v+1}}+C_{\phi}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\quad\mbox{and}\\[2.84526pt] |\phi_{h}(u,v,w,h)|\leq\frac{C_{\phi}}{\sqrt{v+1}}+C_{\phi}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\\[2.84526pt] \mbox{with some positive constants $c_{\phi}$ and $C_{\phi}$,}\end{array}\right.$$

$$\left\{\begin{array}[]{l}0\leq\Phi(w)\leq C_{\Phi}\quad\mbox{for all }w\geq 0\quad\mbox{and}\\[2.84526pt] w\Phi^{\prime 2}(w)\leq C_{\Phi}\Phi(w)\quad\mbox{for all }w\geq 0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\psi(u,v,w,h)\leq C_{\psi}(v)\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\psi_{u}(u,v,w,h)|\leq\frac{C_{\psi}(v)}{\sqrt{uw+1}}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\psi_{v}(u,v,w,h)|\leq C_{\psi}(v)\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4},\\[2.84526pt] |\psi_{w}(u,v,w,h)|\leq\frac{C_{\psi}(v)}{w+1}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\quad\mbox{and}\\[2.84526pt] |\psi_{h}(u,v,w,h)|\leq\frac{C_{\psi}(v)}{(w+1)^{\gamma}}\quad\mbox{for all }(u,v,w,h)\in[0,\infty)^{4}\\[2.84526pt] \mbox{with some nondecreasing $C_{\psi}:[0,\infty)\to(0,\infty)$ and some $\gamma\in(0,\frac{1}{2})$.}\end{array}\right.$$

$$\left\{\begin{array}[]{l}u_{t}=\Delta u-\chi\nabla\cdot(u\nabla h)-\xi\nabla\cdot(u\nabla v)+\mu u(1-u-v),\\[2.84526pt] h_{t}=\Delta h-h+u,\\[2.84526pt] v_{t}=-hv+\eta v(1-u-v),\end{array}\right.$$

$$z:=u{e}^{-\lambda v}\text{mbox}with\lambda :=\frac{\xi }{D_u}>0,$$

$$\left\{\begin{array}[]{ll}z_{t}=D_{u}e^{-\lambda v}\nabla\cdot(e^{\lambda v}\nabla z)-\chi e^{-\lambda v}\nabla\cdot(ze^{\lambda v}\nabla h)+e^{-\lambda v}f(ze^{\lambda v},v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] h_{t}=D_{h}\Delta h+g(ze^{\lambda v},v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] v_{t}=-\alpha ve^{\lambda v}z+v\phi(ze^{\lambda v},v,w,h)+\Phi(w),&x\in\Omega,\ t>0,\\[2.84526pt] w_{t}=\beta e^{\lambda v}z+w\psi(ze^{\lambda v},v,w,h),&x\in\Omega,\ t>0,\\[2.84526pt] \frac{\partial z}{\partial\nu}=\frac{\partial h}{\partial\nu}=0,&x\in\partial\Omega,\ t>0,\\[2.84526pt] z(x,0)=u_{0}(x)e^{\lambda v_{0}(x)},\ h(x,0)=h_{0}(x),\ v(x,0)=v_{0}(x),\ w(x,0)=w_{0}(x),&x\in\Omega,\end{array}\right.$$

$$\mbox{either $T_{max}=\infty$, \quad or}\quad\limsup_{t\nearrow T_{max}}\Big{\{}\|z(\cdot,t)\|_{L^{\infty}(\Omega)}+\|v(\cdot,t)\|_{W^{1,5}(\Omega)}\Big{\}}=\infty.$$

$$v(x,t)\leq\Big{\{}\|v_{0}\|_{L^{\infty}(\Omega)}+\frac{C_{\phi}}{C_{\Phi}}\Big{\}}\cdot e^{C_{\phi}t}\qquad\mbox{for all $x\in\Omega$ and }t\in(0,T_{max}).$$

$$v_t\le C_{\varphi }v+C_{\Phi }\text{mbox}in\Omega \times (0,T_{max}),$$

$$\Vert v(\cdot ,t){\Vert }_{L^{\infty }(\Omega )}$$

$${\int }_{\Omega }u(\cdot ,t)\le C(T)\text{mbox}and{\int }_{\Omega }w(\cdot ,t)\le C(T)\text{mbox}forallt\in (0,\hat{T})$$

$${\int }_{\Omega }h(\cdot ,t)\le C(T)\text{mbox}forallt\in (0,\hat{T}),$$

$$\displaystyle\frac{d}{dt}\bigg{\{}\int_{\Omega}u+\int_{\Omega}w\bigg{\}}$$

$$\displaystyle\frac{d}{dt}\int_{\Omega}h=\int_{\Omega}g(u,v,w,h)\leq C_{g}(c_{1}(T))\cdot\bigg{\{}\int_{\Omega}w+\int_{\Omega}h+|\Omega|\bigg{\}}$$

$$\frac{d}{dt}{\int }_{\Omega }u\ln u+D_u{\int }_{\Omega }\frac{|\nabla u{|}^2}{u}$$

$$\frac{d}{dt}{\int }_{\Omega }u\ln u+D_u{\int }_{\Omega }\frac{|\nabla u{|}^2}{u}=\chi {\int }_{\Omega }\nabla u\cdot \nabla h+\xi {\int }_{\Omega }\nabla u\cdot \nabla v+{\int }_{\Omega }f(u,v,w,h)\ln u+{\int }_{\Omega }f(u,v,w,h),$$

$$\chi {\int }_{\Omega }\nabla u\cdot \nabla h=-\chi {\int }_{\Omega }u\Delta h\le \frac{\eta }{2}{\int }_{\Omega }{u}^2+\frac{{\chi }^2}{2\eta }{\int }_{\Omega }|\Delta h{|}^2\text{mbox}forallt\in (0,T_{max}).$$

$$\int_{\Omega}f(u,v,w,h)\leq C_{f}(\|v\|_{L^{\infty}(\Omega)})\cdot\bigg{\{}\int_{\Omega}u+\int_{\Omega}w+|\Omega|\bigg{\}}$$

$${\int }_{\Omega }f(u,v,w,h)\ln u$$

$$\frac{d}{dt}{\int }_{\Omega }\frac{|\nabla v{|}^2}{v}+\frac{c_{\varphi }}{2}{\int }_{\Omega }\frac{w}{v}|\nabla v{|}^2$$

$$\frac{d}{dt}{\int }_{\Omega }\frac{|\nabla v{|}^2}{v}$$

$${\int }_{\Omega }\frac{\varphi (u,v,w,h)}{v}|\nabla v{|}^2\le -c_{\varphi }\cdot {\int }_{\Omega }\frac{w}{v}|\nabla v{|}^2+C_{\varphi }\cdot {\int }_{\Omega }\frac{|\nabla v{|}^2}{v}\text{mbox}forallt\in (0,T_{max}),$$

$$2{\int }_{\Omega }{\varphi }_u(u,v,w,h)\nabla u\cdot \nabla v$$

Full text

The main results assert unique global classical solvability under suitably mild assumptions on the system parameter functions in associated spatially two-dimensional initial-boundary value problems. In particular, this rigorously confirms that at least in two-dimensional settings, the considered indirectness in signal production induces a significant blow-up suppressing tendency also in taxis systems substantially more general than some particular examples for which corresponding effects have recently been observed.

Key words: chemotaxis; haptotaxis; indirect signal production; global existence

MSC: 35Q92, 35B44, 35K55, 92C17, 35A01

1 Introduction

We study here a general class of parabolic-parabolic-ODE-ODE systems (see (1.3) below) containing the following model of cancer invasion with chemotaxis and haptotaxis:

[TABLE]

supplemented with adequate initial conditions and no-flux boundary conditions. The model variables are: $u$: density of tumor cells, $h$: concentration of matrix metalloproteinases (MMP), $v$: density of tissue fibers (extracellular matrix, ECM), $w$: density of cancer associated fibroblasts (CAFs). The chemotactic bias of the cells is in the direction of the MMP gradient, while haptotaxis means as usual following the gradient of tissue density. In most of the previous chemotaxis and chemotaxis-haptotaxis models the chemoattractant is directly produced by the population performing diffusion and taxis, thereby involving ([38]) or not ([31], [43], [45]) the ECM in this production. In our present model the chemoattractant is generated in an indirect way, by the CAFs, which are activated by the tumor cells. In [5] was introduced a complex model for the evolution of a population of tumor cells interacting with two chemoattractants and also performing haptotaxis. One of the chemoattractants therein is directly produced by the cells, while the other’s production is mediated by another substance, the latter being in turn produced by the cells under the influence of the first chemoattractant. All model variables (except tissue) are diffusing in a linear way; in particular, the producers of all chemoattractants are diffusive. This is not the case in our model (1.1); moreover, both signals are completely ($h$) or partially ($v$) obtained from the non-diffusing, indirect producer $w$.

Model (1.1) is motivated by the problem of investigating CAF-mediated cancer invasion into the surrounding tissue. CAFs are major components of the neoplasm microenvironment. They secrete a variety of extracellular matrix components and are involved in the formation of the desmoplastic stroma characterizing many advanced carcinomas ([15]). For a long time the general belief was that tumor development, invasion, and metastasis occur as a result of cancer progression. Recent studies revealed, however, that CAFs contribute instead of tumor cells to these processes, via expression of various growth factors, cytokines, chemokines, and degradation of ECM ([7], [15], [34], [36]), but also by restructuring the latter to facilitate migration ([22]). Fibronectin (Fn) assembled by CAFs mediates cell association and directional migration. Compared with normal fibroblasts, CAFs produce an Fn-rich ECM with anisotropic fiber orientation, along which the tumor cells preferentially migrate ([10]). The origin of CAFs is not completely elucidated; we refer e.g. to [6, 20] for a couple of reviews. There is evidence that they can arise among others from carcinoma cells through epithelial-mesenchymal transition (EMT) ([23]), thus allowing the cancer cells to adopt a mesenchymal phenotype associated to enhanced migratory capacity and invasiveness ([9], [30]). It has also been shown (cf. e.g. [29]) that cancer cells can reprogram resident tissue fibroblasts to become CAFs through the actions of miRNAs. MMPs are primarily derived from CAFs in various types of tumor ([23], [26], [42]). In particular, it has been shown e.g., that MMP-9, an endoproteinase involved in ECM degradation and implicated as a prerequisite of metastasis, has very limited or no expression in various cancer cell lines. Instead, MMP-9 is well-known to be secreted from cancer stromal fibroblasts and endothelial cells ([39]). For further information about CAFs and their MMP production we refer, e.g., to the review [26]. The main features of CAF-mediated tumor invasion mentioned above are captured in our model (1.1) below.

Closely related from a mathematical viewpoint is the chemotaxis-haptotaxis model

[TABLE]

which describes the evolution of two cancer cell subpopulations either proliferating ($w$) or migrating ($u$), with the corresponding transitions between the two phenotypes, and in interaction with tissue ($v$) and acidity ($h$). The latter is mainly produced by the highly glycolytic proliferating tumor cells – hence indirect signal production. The proliferation/migration (also known as go-or-grow) dichotomy asserting that moving cells defer their proliferation seems to be a relevant feature for some types of cancer ([13], [21],[49]) and can lead to interesting mathematical problems and qualitative behavior of the corresponding heterogeneous tumor ([37], [51]). The migrating cells in (1.2) perform haptotaxis and pH-taxis: they move towards increasing ECM gradient and away from acidic (thus hypoxic) areas. The conversion from proliferating to migrating cells depends on the concentration $h$ of protons in the peritumoral region and infers limitations, as only a rather small part of the tumor becomes motile -usually cells situated at the tumor margins. Moreover, the protons are buffered by the environment (e.g., uptake by vasculature), contribute to tissue degradation, and restrict tumor proliferation. Supplementary to hypoxia the tissue can be degraded by chemical or biological agents directly or indirectly produced or activated by the tumor cells, see e.g. [27] and [28]. The ECM remodeling also involves limited, logistic-like growth. The latter is also used in (1.2) to describe tumor cell proliferation. We note here that taking squares of the positive parts of the proliferation terms is motivated by the mainly technical need of satisfying the regularity assumptions in Theorem 1.1 and could perhaps be relaxed; in fact, the difference between functions of this type just involving the positive parts and those taking their squares is rather small. System (1.2) fits in the theoretical framework proposed and analyzed below.

Thus the general model class to be introduced and studied here includes descriptions of biologically quite different problems, as exemplified above. It therefore provides a comprehensive mathematical structure for several issues related to cancer cell migration under the influence of chemo- and haptotactic effects, including some variants of the often addressed Chaplain-Lolas model from [5]. Other problems characterizing tactic cell migration, e.g. in wound healing and/or angiogenesis or interacting microbial populations in biofilm formation and persistence, can potentially be cast into this mathematical framework. Models describing the dynamics of a species performing chemotaxis towards the gradient of a signal produced by another, non-tactic species (see e.g. [19] and [47]) form a subclass of this setting.

Approaching a mathematical core feature: indirectness of chemoattractant production.  From a purely mathematical perspective, a common feature distinguishing both (1.2) and (1.1) from classical chemotaxis or chemotaxis-haptotaxis systems of Keller-Segel or Chaplain-Lolas type consists in the circumstance that the respective mechanisms of chemoattractant production are indirect in the sense that the corresponding signal is not produced directly by individuals of the cell population, but rather through a third agent. Possible implications on the system dynamics, however, have apparently been detected only in some quite particular examples of chemotaxis-only models: Indeed, the only findings in this direction we are aware of concentrate on associated derivates of the classical Keller-Segel system, for which substantial blow-up preventing effects of such indirect signal production mechanisms have recently been revealed. More precisely, in sharp contrast to what is known for standard Keller-Segel systems ([17], [2]), regardless of the size of the initial data some corresponding spatially two-dimensional initial-boundary value problems always possess globally defined classical solutions ([47], [25]).

Complementary to its biological motivation, the main mathematical purpose of the present work is to rigorously confirm that such relaxing effects of indirect chemoattractant production are actually not restricted to cases of chemotaxis-only systems, but rather seem to form a much more general and robust feature of chemotactic interaction also in significantly more contexts involving haptotactic migration mechanisms with all their potentially regularity-limiting properties due to lack of haptoattractant diffusion.

This will subsequently be examined in the framework of the problem

[TABLE]

in a bounded domain $\Omega\subset\mathbb{R}^{2}$ with smooth boundary, where the parameters $\chi,\xi,D_{u}$, $D_{h}$, $\alpha$ and $\beta$ are assumed to be positive, and where for simplicity we suppose throughout this paper that with some $\vartheta\in(0,1)$,

[TABLE]

As for the parameter functions in (1.3), in order to create a setup sufficiently general so as to include both (1.2) and (1.1), we shall require that

[TABLE]

and are such that

[TABLE]

and

[TABLE]

that

[TABLE]

and that

[TABLE]

and finally

[TABLE]

As can readily be verified, indeed both models (1.2) and (1.1) then become special cases of the PDE system in (1.3) whenever the parameters $D_{u},D_{h},\chi,\xi,k_{4}$ and $k_{6}$ therein are positive, whereas $k_{i}$ for $i\in\{1,2,3,5,7,8\}$ and $\eta$ and, in particular, $\mu$ is merely required to be nonnegative.

Our main results in this context then read as follows.

Theorem 1.1

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary, assume that $\chi,\xi,\alpha,\beta,D_{u}$ and $D_{h}$ are positive, and suppose that $f$, $g$, $\phi$, $\Phi$ and $\psi$ satisfy (1.5), (1.6), (1.7), (1.8), (1.9) and (1.10). Then for any choice of $(u_{0},h_{0},v_{0},w_{0})$ fulfilling (1.4) with some $\vartheta\in(0,1)$, the problem (1.3) possesses a uniquely determined globally defined classical solution $(u,h,v,w)\in(C^{2,1}(\overline{\Omega}\times[0,\infty))^{4}$ for which $u,h,v$ and $w$ are nonnegative.

In particular, Theorem 1.1 asserts that indeed no finite-time blow-up occurs in (1.3) under the above assumptions, and that in this regard the solution behavior in (1.3) quite drastically differs from that in the corresponding variants of (1.2) and (1.1) in which the second equation is replaced with e.g. $h_{t}=D_{h}\Delta h+u$, and in which already in the semi-trivial case when $v\equiv 0$, known results on the actually resulting two-component Keller-Segel system for $(u,h)$ assert finite-time blow-up for some solutions ([16]).

In line with this, all known results even on global existence, but also on qualitative properties, in the original Chaplain-Lolas model ([5])

[TABLE]

seem to strongly rely on the assumption that $\mu$ be positive, thus guaranteeing the presence of a logistic-type quadratic growth restriction on the cell density ([4], [32], [31], [40], [41], [46], [48]). Our results show that in the context of the variant (1.1) of (1.11) involving indirect chemoattractant production, no such additional dampening is necessary: Indeed, even when tissue remodeling is included by supposing that $\eta>0$ in (1.1), Theorem 1.1 asserts global classical solvability in (1.1) for actually any nonnegative value of $\mu$.

The paper is structured as follows: After stating a result on local existence and extensibility as well as some preliminary estimates in Section 2, in Section 3 we shall construct a quasi-energy functional for (1.3) and draw some immediate conclusions concerning regularity of solutions. The accordingly obtained estimates are used as a starting point for a Moser-type iterative argument yielding $L^{\infty}$ bounds for the key solution component $u$ in Section 4, and hence paving an essential part of the way toward our proof of Theorem 1.1 in Section 5. Finally, in Section 6 we illustrate the theoretical findings by numerical simulations of (1.1) and the corresponding model with direct signal production and provide some comments about the obtained results.

2 Local existence and basic estimates

Following several precedents in the literature ([11], [12], [31]), in order to establish a preliminary result on local existence, but also to prepare our subsequent estimation procedure, we note that on substituting

[TABLE]

the problem (1.3) is equivalently transformed to

[TABLE]

In this formulation, with respect to the construction of local-in-time solutions the problem (1.3) indeed becomes accessible to appropriate fixed point frameworks; by straightforward and minor adaptations of the corresponding arguments detailed e.g. in [31], it is thereby possible to establish the following basic statement on unique solvability and extensibility.

Lemma 2.1

Let $\chi,\xi,\alpha,\beta,D_{u}$ and $D_{h}$ be positive, let $f$, $g$, $\phi$, $\Phi$ and $\psi$ comply with (1.5), (1.6), (1.7), (1.8), (1.9) and (1.10), and suppose that $(u_{0},h_{0},v_{0},w_{0})$ satisfies (1.4) with some $\vartheta\in(0,1)$. Then there exist $T_{max}\in(0,\infty]$ and a uniquely determined classical solution $(z,h,v,w)\in(C^{2,1}(\overline{\Omega}\times[0,\infty)))^{4}$ of (2.2) such that

[TABLE]

Moreover, we have $z>0,h\geq 0,v>0$ and $w>0$ in $\overline{\Omega}\times(0,T_{max})$.

Without any further explicit mentioning, throughout the sequel we shall suppose that the assumptions of Lemma 2.1 are satisfied, and that $(z,v,w,h)$ and $T_{max}\in(0,\infty]$ are as provided by the latter. Moreover, we shall tactitly switch between these variables and the quadruple $(u,v,w,h)$ solving (1.3) classically in $\Omega\times(0,T_{max})$, as thereby defined through (2.1).

A first boundedness property of this solution is immediate.

Lemma 2.2

The solution of (1.3) satisfies

[TABLE]

Proof.  As from (1.3), (1.8) and (1.9) we know that

[TABLE]

by means of a simple comparison argument we conclude that

[TABLE]

Herein estimating $\int_{0}^{t}e^{C_{\phi}(t-s)}ds\leq\frac{1}{C_{\phi}}e^{C_{\phi}t}$ for $t\geq 0$, from this we readily obtain (2.4). $\Box$ Next, thanks to (1.6) and (1.7) the first and fourth solution components in (1.3) can at least controlled with respect to their norm in $L^{1}(\Omega)$ in a fairly simple manner.

Lemma 2.3

Let $T>0$. Then there exists $C(T)>0$ such that

[TABLE]

as well as

[TABLE]

where $\widehat{T}:=\min\{T,T_{max}\}$.

Proof.  Using that v\leq c_{1}(T):=\Big{\{}\|v_{0}\|_{L^{\infty}(\Omega)}+\frac{C_{\Phi}}{C_{\phi}}\Big{\}}\cdot e^{C_{\phi}T} in $\Omega\times(0,\widehat{T})$ due to Lemma 2.2, after spatial integration in the first and the fourth equation in (1.3) and adding the respective results we may rely on (1.6) and (1.10) in estimating

[TABLE]

for $t\in(0,\widehat{T})$. A time integration of this linear ODI for $\int_{\Omega}u+\int_{\Omega}w$ directly yields (2.5). Similarly, (1.7) implies that

[TABLE]

for all $t\in(0,\widehat{T})$, so that (2.6) becomes a consequence of (2.5). $\Box$

3 A quasi-energy inequality

The purpose of this section consists in the construction of an Lyapunov-like functional which through a corresponding energy-dissipation inequality will provide some fundamental regularity information that will form the starting point of a series of a priori estimates which in the presently considered spatially two-dimensional setting will finally allow for the conclusion that $(u,v,w,h)$ is actually global in time.

As our first step in this direction, let us perform a standard testing procedure by which the crucial haptotactic contribution to the first equation in (1.3) is reduced to an $L^{2}$ inner product of gradients:

Lemma 3.1

Let $\eta>0$. Then

[TABLE]

for all $t\in(0,T_{max})$.

Proof.  In the identity

[TABLE]

valid for all $t\in(0,T_{max})$ due to (1.3), we use Young’s inequality to estimate

[TABLE]

Moreover, by means of (1.6) we see that for all $t\in(0,T_{max})$,

[TABLE]

and

[TABLE]

because $\ln y\leq y$ for all $y\geq 1$. In conjunction with (3.2)-(3.4), this establishes (3.1). $\Box$ Now in order to achieve a precise cancelation of the integral in (3.1) stemming from haptotactic cross-diffusion, inspired by several precedent works on similar types of interaction (see e.g. [8], [45] and [38]), we track the evolution of the Dirichlet integral associated with $\sqrt{v}$.

Lemma 3.2

The solution of (1.3) has the property that

[TABLE]

Proof.  According to the third equation in (1.3), we have

[TABLE]

Here from (1.8) we know that

[TABLE]

and that due to Young’s inequality, for all $t\in(0,T_{max})$ we have

[TABLE]

and

[TABLE]

as well as

[TABLE]

and

[TABLE]

Moreover, combining Young’s inequality with (1.9) we can estimate

[TABLE]

As the second summand on the right of (3.7) is nonpositive, collecting (3.8)-(3.13) we thus infer (3.5) from (3.7). $\Box$ Next, several expressions on the right-hand sides of (3.1) and (3.5) need to be controlled in modulus. Here the second integral on the right of (3.1) can in fact be absorbed by the dissipation rate appearing in the following inequality gained by means of a standard procedure.

Lemma 3.3

For all $t\in(0,T_{max})$,

[TABLE]

Proof.  On testing the second equation in (1.3) by $\Delta h$ and using Young’s inequality in a standard manner, we obtain

[TABLE]

which implies (3.14) due to the fact that by (1.7) and again Young’s inequality,

[TABLE]

in $\Omega\times(0,T_{max})$. $\Box$ The second last summand in (3.5), referring to the component $w$ with evolution governed by an ODE only, apparently cannot be expected to be absorbed by some suitable dissipation rate. The following lemma indicates that at least some exponential control thereof will eventually be possible.

Lemma 3.4

We have

[TABLE]

Proof.  Using the fourth equation in (1.3), we compute

[TABLE]

where by Young’s inequality,

[TABLE]

and where by (1.10),

[TABLE]

Apart from that, combining (1.10) with Young’s inequality we see that

[TABLE]

and

[TABLE]

as well as

[TABLE]

and

[TABLE]

for all $t\in(0,T_{max})$. Inserting (3.17)-(3.22) into (3.16) directly yields (3.15). $\Box$ A final minor ingredient to our quasi-energy inequality is addressed in the following.

Lemma 3.5

If $\eta>0$, then

[TABLE]

Proof.  As

[TABLE]

by (1.3), this follows by observing that

[TABLE]

due to Young’s inequality, and that

[TABLE]

according to (1.10). $\Box$ As a last preparation, let us make use of appropriate parabolic regularization features to estimate terms of the form appearing in the last integral from (3.15).

Lemma 3.6

Let $q\in[1,4)$. Then there exists $\theta=\theta(q)\in(0,1)$ with the property that for all $T>0$ one can find $C(q,T)>0$ such that

[TABLE]

where again $\widehat{T}:=\min\{T,T_{max}\}$.

Proof.  We first note that according to (1.7), Lemma 2.2 and Lemma 2.3 we can fix positive constants $c_{1}(T),c_{2}(T)$ and $c_{3}(T)$ such that

[TABLE]

and that

[TABLE]

and observe that since $q<4$ it is possible to choose $r\equiv r(q)\in(1,2)$ suitably close to $2$ satisfying $r>q-2$. In the Duhamel representation

[TABLE]

we may then use standard $L^{p}$-$L^{q}$ estimates for the Neumann heat semigroup $(e^{\sigma\Delta})_{\sigma\geq 0}$ ([50, Lemma 1.3]) to find $c_{4}(T)>0$ such that for all $t\in(0,\widehat{T})$,

[TABLE]

As (3.25) and (3.26) warrant that

[TABLE]

this readily entails that

[TABLE]

with $c_{6}(T)$ being finite due to our restriction $r<2$. Now since a combination of the Gagliardo-Nirenberg inequality with elliptic regularity theory yields $c_{7}>0$ such that

[TABLE]

from (3.27) we immediately obtain (3.24) with $\theta\equiv\theta(q):=\frac{q-r}{2}$ fulfilling $\theta\in(0,1)$ according to the inequality $r>q-2$. $\Box$ We can now proceed to our detection of an energy-like structure in (1.3), as expressed in the following lemma.

Lemma 3.7

Let $T>0$. Then there exist $a=a(T)>0,b>0$ and $C=C(T)>0$ such that for

[TABLE]

and

[TABLE]

we have

[TABLE]

where $\widehat{T}:=\min\{T,T_{max}\}$.

Proof.  Thanks to Lemma 2.3 and Lemma 2.2, we can fix positive constants $c_{1}(T),c_{2}(T)$ and $c_{3}(T)$ such that

[TABLE]

as well as

[TABLE]

Then in accordance with the Gagliardo-Nirenberg inequality, choosing $c_{4}>0$ such that

[TABLE]

we define

[TABLE]

and take $a\equiv a(T)>0$ large enough fulfilling

[TABLE]

Finally picking $b>0$ small such that both

[TABLE]

and

[TABLE]

hold, we let ${\mathcal{F}}$ and ${\mathcal{D}}$ be as determined through (3.28) and (3.29) and claim that then (3.30) is valid with some suitably large $C(T)>0$.

To this end, we first take an appropriate linear combination of the inequalities provided by Lemma 3.1, Lemma 3.3, Lemma 3.2, Lemma 3.4 and Lemma 3.5, which when applied to our particular value of $\eta$ namely show that

[TABLE]

[TABLE]

with obvious definitions of $c_{i}(T)$ for $i\in\{5,...,12\}$, where we have made use of a favorable cancelation in some contributions stemming from the haptotactic interaction.

Now employing (3.33) we see that due to (3.31) and (3.34) we have

[TABLE]

whereas Lemma 3.6 says that as $\frac{4}{1+2\gamma}<4$, there exist $\theta\in(0,1)$ and $c_{13}(T)>0$ such that

[TABLE]

whence by Young’s inequality,

[TABLE]

with some $c_{14}(T)>0$. Since moreover the Poincaré inequality provides $c_{15}(T)>0$ satisfying

[TABLE]

thanks to the restrictions in (3.35), (3.36) and (3.37) we conclude from (3.38), (3.39) and (3.40) that there exists $c_{16}(T)>c_{6}(T)$ such that

[TABLE]

for all $t\in(0,\widehat{T})$. Using that $y\ln y\geq-\frac{1}{e}$ for all $y>0$ and thus $\int_{\Omega}u\ln u\geq-\frac{|\Omega|}{e}$ for all $t\in(0,T_{max})$, from this we finally obtain that for all $t\in(0,\widehat{T})$,

[TABLE]

as intended. $\Box$ Upon integration, the latter implies several a priori estimates, significantly going beyond those from Lemma 2.2 and Lemma 2.3, among which we explicitly state only those three inequalities that will be referred to later on.

Lemma 3.8

For all $T>0$ there exists $C(T)>0$ such that again writing $\widehat{T}:=\min\{T,T_{max}\}$ we have

[TABLE]

and

[TABLE]

as well as

[TABLE]

Proof.  Upon integrating (3.30) in time, we can find $c_{1}(T)>0$ such that with ${\mathcal{F}}$ and ${\mathcal{D}}$ taken from (3.28) and (3.29) we have

[TABLE]

Once more using that $y\ln y\geq-\frac{1}{e}$ for $y>0$, from this we readily obtain the claimed inequalities as particular consequences. $\Box$ As an immediate consequence, let us add the following observation about regularity of $h$.

Corollary 3.9

Let $p\geq 1$ and $T>0$. Then there exists $C(p,T)>0$ such that with $\widehat{T}:=\min\{T,T_{max}\}$,

[TABLE]

Proof.  As $W^{1,2}(\Omega)\hookrightarrow L^{p}(\Omega)$, combining (2.6) with (3.42) immediately yields (3.44). $\Box$

4 $L^{\infty}$ bounds for $u$

In this section we intend to make use of the information from Lemma 3.8 in order to finally achieve an a priori bound for the quantity $z$ from (2.1), and hence for $u$, with respect to the norm in $L^{\infty}(\Omega)$. Here a key role will be played by the following implication of the estimate (3.41) on a Gagliardo-Nirenberg-type interpolation, as expressed in the following.

Lemma 4.1

Let $p>1$ and $T>0$. Then for all $\eta>0$ one can find $C(\eta,p,T)>0$ such that

[TABLE]

with $\widehat{T}:=\min\{T,T_{max}\}$.

Proof.  This follows in a standard manner from the boundedness property of $z$ in $L\log L(\Omega)$ as implied by Lemma 3.8 and Lemma 2.2, by means of a refined interpolation inequality of Gagliardo-Nirenberg type, the latter going back to [3] and extended to a version applicable to the present setting in [45, Lemma A.5] (for details, see e.g. [45, p.800]). $\Box$ In order to make appropriate use of this, let us perform another well-established testing procedure to (2.2), a basic outcome of which is the following.

Lemma 4.2

Let $T>0$. Then there exists $C(T)>0$ such that for all $p\geq 2$ and all $t\in(0,\widehat{T})$ with $\widehat{T}:=\min\{T,T_{max}\}$,

[TABLE]

Proof.  By means of (2.2), we obtain

[TABLE]

for all $t\in(0,T_{max})$, where by Young’s inequality,

[TABLE]

for all $t\in(0,\widehat{T})$, with $c_{1}(T):=\|v\|_{L^{\infty}(\Omega\times(0,\widehat{T}))}$ being finite according to Lemma 2.2. Furthermore, (1.6) and Young’s inequality warrant that for all $t\in(0,\widehat{T})$,

[TABLE]

while from (1.8) and (1.9) we know that

[TABLE]

and that

[TABLE]

As clearly $e^{\lambda v}\geq 1$, by nonnegativity of $\alpha$ we therefore infer (4.2) from (4.3) when combined with (4.4), (4.5), (4.6) and (4.7). $\Box$ Here a suitable control of the crucial rightmost summand in (4.2) will rely, besides on Lemma 4.1, also on the following statement which can be viewed as partially generalizing Lemma 3.5.

Lemma 4.3

Let $p\geq 1$ and $T>0$. Then writing $\widehat{T}:=\min\{T,T_{max}\}$ and $\tau:=\frac{1}{2}\widehat{T}$, with some $C(p,T)>0$ we have

[TABLE]

Proof.  Choosing $c_{1}(T)>0$ large enough fulfilling $v\leq c_{1}(T)$ in $\Omega\times(0,\hat{T})$ according to Lemma 2.2, equation (1.3), and due to Young’s inequality and (1.10) we can estimate

[TABLE]

so that with some $c_{2}(p,T)>0$ we have

[TABLE]

An integration thereof shows that

[TABLE]

and hence implies (4.8). $\Box$ We can thereby use Lemma 4.2 along with Lemma 4.1 to derive the following $L^{p}$ estimate for $z$, at this stage yet involving bounds possibly depending on the finite number $p\geq 2$.

Lemma 4.4

Let $p\geq 2$ and $T>0$. Then there exists $C(p,T)>0$ such that

[TABLE]

where again $\widehat{T}:=\min\{T,T_{max}\}$.

Proof.  From Lemma 4.2 we obtain $c_{1}(p)>0$ and $c_{2}(p,T)>0$ such that

[TABLE]

again because $e^{\lambda v}\geq 1$. Here since Lemma 3.8 provides $c_{3}(T)>0$ such that

[TABLE]

by means of Young’s inequality and the Gagliardo-Nirenberg inequality we infer that with some $c_{4}(p,T)>0$ we have

[TABLE]

Therefore, (4.10) shows that there exists $c_{5}(p,T)>0$ such that $y(t):=\int_{\Omega}e^{\lambda v(\cdot,t)}z^{p}(\cdot,t)$, $t\in[0,T_{max})$, satisfies

[TABLE]

for all $t\in(0,\widehat{T})$ and thus, upon integration,

[TABLE]

for all $t\in(\tau,\widehat{T})$, where again we have set $\tau:=\frac{1}{2}\widehat{T}$. As thus $\tau$ is positive, a well-known result on maximal Sobolev regularity in the parabolic subproblem of (1.3) satisfied by $h$ ([14]) becomes applicable so as to yield $c_{6}(p,T)>0$ satisfying

[TABLE]

for all $t\in(\tau,\widehat{T})$ because of (1.7). Since Corollary 3.9 and Lemma 4.3 provide $c_{7}(p,T)>0$ and $c_{8}(p,T)>0$ such that

[TABLE]

and that

[TABLE]

this means that with some $c_{9}(p,T)>0$ we have

[TABLE]

so that from (4.12) we infer the existence of $c_{10}(p,T)>0$ satisfying

[TABLE]

We now employ Lemma 4.1 to see that with some $c_{11}(p,T)>0$,

[TABLE]

whence (4.13) ensures that

[TABLE]

and that thus, clearly, (4.9) holds. $\Box$ By adapting a well-established Moser-type iteration ([1], [44]) to the present context, however, one can readily turn the latter into estimates in $L^{\infty}$.

Lemma 4.5

Given $T>0$, one can find $C(T)>0$ such that with $\widehat{T}:=\min\{T,T_{max}\}$ we have

[TABLE]

Proof.  We first fix any $p_{\star}>2$ and then obtain upon combining Lemma 4.4 with Lemma 4.3, Lemma 2.2, (1.7) and Corollary 3.9 that

[TABLE]

with $\tau:=\frac{1}{2}\widehat{T}$ and some $c_{1}(T)>0$. As a consequence thereof, standard regularization features of the Neumann heat semigroup ([50, Lemma 1.3], [18, Lemma 4.1]) entail boundedness of $\nabla h$ in $\Omega\times(\tau,\widehat{T})$, so that from Lemma 4.2 we infer the existence of $c_{2}(T)>0$ and $c_{3}(T)>0$ such that for $p:=2^{k+1}$ and any nonnegative integer $k$,

[TABLE]

On integrating and recalling Lemma 2.2 and Lemma 4.3, we see that with some $c_{4}(T)>0$, for any such $p$ and arbitrary $t\in(\tau,\widehat{T})$ we have

[TABLE]

The remaining part now follows a well-established reasoning: According to the Gagliardo-Nirenberg inequality and Young’s inequality, we can find $c_{5}(T)>0$ such that introducing the numbers

[TABLE]

all finite due to Lemma 4.4, for $p=p_{k}$ and each $k\in\{1,2,3,...\}$ we have

[TABLE]

whence (4.15) entails that for some $c_{6}(T)>0$,

[TABLE]

By means of a standard recursive argument, both when $p_{k}^{4}M_{k-1}^{2}\leq\|z(\cdot,\tau)\|_{L^{\infty}(\Omega)}^{p_{k}}$ for infinitely many $k\geq 1$, and as well in the opposite case this can readily be seen to imply the existence of $c_{7}(T)>0$ such that

[TABLE]

from which the claim immediately follows. $\Box$

5 Global extensibility. Proof of Theorem 1.1

Having thus ruled out blow-up of the first quantity appearing in the second alternative from (2.3), it hence remains to derive appropriate bounds for the haptotactic gradient. A first observation relates the latter to some spatio-temporal regularity properties of $z$ and $h$.

Lemma 5.1

Let $T>0$ and $q\geq 1$. Then there exists $C(q,T)>0$ fulfilling

[TABLE]

for all $t\in(\tau,\widehat{T})$, with $\widehat{T}:=\min\{T,T_{max}\}$ and $\tau:=\frac{1}{2}\widehat{T}$.

Proof.  Differentiating in (1.3) and recalling (2.1), we see that throughout $\Omega\times(0,T_{max})$,

[TABLE]

and

[TABLE]

where we have suppressed the argument $(u,v,w,h)$ in $\phi$, $\psi$ and the derivatives thereof. Now as a consequence of Lemma 2.2, Lemma 4.5 and (2.1), our requirements (1.8), (1.9) and (1.10) guarantee that herein all the functions $-\alpha v+v\phi_{u}e^{\lambda v}$, $-\alpha\lambda uv-\alpha u+\lambda uv\phi_{u}+v\phi_{v}+\phi$, $v\phi_{w}+\Phi^{\prime}(w)$, $v\phi_{h}$, $\beta e^{\lambda v}+we^{\lambda v}\psi_{u}$, $\beta\lambda u+\lambda uw+w\psi_{v}$, $w\psi_{w}+\psi$ and $w\psi_{h}$ are bounded in $\Omega\times(0,\widehat{T})$, so that (5.2) and (5.3) imply that with some $c_{1}(T)>0$ we have

[TABLE]

and, similarly,

[TABLE]

for all $t\in(\tau,\widehat{T})$. Adding these inequalities and invoking Gronwall’s lemma readily leads to (5.1). $\Box$ Here the crucial ingredient containing $\nabla z$ will ultimately be controlled by using the following result which extends a corresponding statement from [45, Lemma 3.14] to the present more complex system (2.2).

Lemma 5.2

Let $T>0$. Then there exists $C(T)>0$ with the property that

[TABLE]

where again $\widehat{T}:=\min\{T,T_{max}\}$ and $\tau:=\frac{1}{2}\widehat{T}$.

Proof.  Defining

[TABLE]

for $(x,t)\in\Omega\times(0,T_{max})$, we see that the first equation in (2.2) simplifies to the identity

[TABLE]

which we multiply by $-\Delta z$ and integrate to find using Young’s inequality that for all $t_{0}\in[0,T_{max})$ and any $t\in(t_{0},T_{max})$,

[TABLE]

Here we note that given $T>0$, in view of the boundedness property of $z$ in $\Omega\times(0,\widehat{T})$ asserted by Lemma 4.5 we can fix $c_{1}(T)>0$ such that due to Young’s inequality,

[TABLE]

Moreover, as the Gagliardo-Nirenberg inequality together with elliptic estimates and Lemma 4.5 provides $c_{2}>0$ and $c_{3}(T)>0$ such that

[TABLE]

by combining the Cauchy-Schwarz inequality with Young’s inequality we obtain

[TABLE]

which along with (5.6) shows that (5.5) implies the inequality

[TABLE]

with some appropriately large $c_{4}(T)>0$.

To proceed from this, we observe that as a consequence of Lemma 5.1 we can find $c_{5}(T)>0$ such that

[TABLE]

and that another application of the Gagliardo-Nirenberg inequality in conjunction with Lemma 3.8 provides $c_{6}>0$ and $c_{7}(T)>0$ fulfilling

[TABLE]

Therefore, again through Lemma 3.8, (5.8) reduces to

[TABLE]

where by the Fubini theorem and, again, (5.7),

[TABLE]

We now let $\delta>0$ be small enough such that $c_{3}(T)c_{8}(T)\delta\leq\frac{D_{u}}{8}$, and choose $N\in\mathbb{N}$ and $(t_{i})_{i\in\{1,...,N\}}\subset[\tau,\widehat{T}]$ such that $t_{1}=\tau,t_{N}=\widehat{T}$ and $0<t_{i+1}-t_{i}\leq\delta$ for all $i\in\{1,...,N-1\}$. Then inserting (5.10) into (5.9) shows that

[TABLE]

along with $I_{0}:=\int_{\Omega}|\nabla z(\cdot,\tau)|^{2}$ and $J_{0}:=0$ satisfy

[TABLE]

whence if we let

[TABLE]

then

[TABLE]

with $c_{9}(T):=\max\{1+(N-2)c_{3}(T)c_{8}(T)T\,,\,I_{0}+(N-1)c_{8}(T)+1\}$. Therefore, $K_{i}\leq K_{0}\cdot c_{9}^{N-1}(T)$ for all $i\in\{1,...,N-1\}$, which in view of the definitions of $(I_{i})_{i\in\{1,...,N-1\}}$ and $(J_{i})_{i\in\{1,...,N-1\}}$ yields (5.4). $\Box$ Combining this with Lemma 5.1 readily implies the following.

Lemma 5.3

For all $q\geq 1$ and $T>0$ one can fix $C(q,T)>0$ such that with $\widehat{T}:=\min\{T,T_{max}\}$ and $\tau:=\frac{1}{2}\widehat{T}$ we have

[TABLE]

Proof.  From Lemma 5.2 in conjunction with elliptic regularity theory we obtain $c_{1}(T)>0$ such that

[TABLE]

which we combine with the continuity of the embedding $W^{2,2}(\Omega)\hookrightarrow W^{1,q}(\Omega)$ and the Cauchy-Schwarz inequality to conclude that with some $c_{2}(q)>0$ we have

[TABLE]

Since Lemma 3.8 in quite a similar fashion yields the existence of $c_{3}(q,T)>0$ such that

[TABLE]

the claimed estimate is a consequence of Lemma 5.1. $\Box$ Now our main result has actually already been established:

Proof  of Theorem 1.1.  Thanks to Lemma 2.1, and in particular the extensibility criterion (2.3) therein, we only need to combine the outcome of Lemma 4.5 with an application of Lemma 5.3 to $q:=5$. $\Box$

6 Simulations and discussion

To illustrate our theoretical results we also performed some numerical simulations of the version

[TABLE]

of the indirect signal production model (1.1) and – for comparison purposes – we also simulated some solutions of the system

[TABLE]

in which the signal (MDEs) is directly produced by the tumor cells. Since for both models there is no blow-up as long as $\mu>0$ (see [24] for the result concerning (6.13)), we only considered here the case with no tumor cell proliferation, i.e. $\mu=0$.

The simulations were performed using a discontinuous Galerkin FEM method. Thereby, the diffusion was discretized in space by using a symmetric interior penalty Galerkin (SIPG) method (see [35]), while for the drift term we did an upwind discretization. The time was discretized with an IMEX procedure handling the diffusion implicitly and the reaction and taxis terms explicitly. The computational domain was $\Omega=[0,1]^{2}$ and the initial conditions for tumor cells and MDEs were chosen in the form $u_{0}(x)=\exp(-\frac{|x|^{2}}{2\varepsilon_{u}})$, $h_{0}(x)=\exp(-\frac{|x|^{2}}{2\varepsilon_{h}})$, with $\varepsilon_{u}=0.05$, $\varepsilon_{h}=0.1$. For the initial density of CAFs we considered a radially symetric form $w_{0}(r)=\exp(-\frac{(r-r_{0})^{2}}{2\varepsilon_{w}})$ with $\varepsilon_{w}=0.01$ and $r_{0}=0.5$. Finally, the initial tissue density was characterized by v_{0}(x)=\left\{\begin{array}[]{cc}v_{\text{max}},&x\in\Omega_{v}\\ v_{\text{min}},&x\not\in\Omega_{v}\end{array}\right., with $v_{\text{max}}=1$ and $v_{\text{min}}=0.2$, thus $\Omega_{v}\subset\Omega$ representing the stripes shown in the second columns of Figures 1 and 2. Together these initial conditions describe a heterogeneous tissue structure, in which a Gaussian-shaped tumor is embedded, surrrounded by activated CAFs and featuring higher MDE concentration in the areas with many tumor cells. We considered for both models (where applicable) the following parameters: $\chi=0.6$; $\xi=0.5$; $D_{u}=10^{-10}$; $D_{h}=0.1$; $\eta=10.6$; $\alpha=5$; $\beta=1$; $\gamma=1.0$. The results are shown in Figure 1 for system (6.12) and in Figure 2 for system (6.13). The first rows represent the initial conditions, and the subsequent rows illustrate the solution behavior at several successive time steps.

The simulations of (6.12) depicture an increase in tumor cell density, which is, however, limited. Same applies to the MDE concentration and CAFs density, but with smaller rates. The tissue is correspondingly degraded, and the CAFs spread into the region containing the main tumor mass, at the same time building up some tissue (for a sufficiently large $\beta$, as in these simulations) at the sites where they are abundant enough. In contrast, model (6.13) predicts localized tumor cell aggregates of very high density which are almost three orders of magnitude higher than the initial condition and keep growing in time, thus hinting on blow-up of the solution. Likewise, in this latter setting the MDE concentration is directly produced by the tumor cells and keeps growing as well, although to a much smaller rate than the cell density. The tissue degradation is much more localized and – where it happens – stronger. These results are conform with the theoretical findings in this and previous papers predicting blow-up of solutions when the signal was directly produced by the agents performing chemotaxis ([2]), while the solutions stay bounded in the case of indirect signal production.

From a biological viewpoint the tumor cells use CAFs (which are originally ’harmless’ stroma cells only becoming supporters of tumor invasion upon activation) to produce matrix degrading factors. As mentioned above, the production of the latter seems to be decisively controled by CAFs, hence is rather indirect, as the neoplastic cells first need to activate the CAFs, which then enhance degradation of surrounding tissue and cell motility, including chemotaxis towards the gradient of proteolytic agents. Our mathematical result actually tells that such mediation of invasion leads to avoidance of blow-up, unlike previous models where the direct production of chemoattractant let the solution(s) become unbounded, a rather unrealistic biological scenario. The result is in line with many in vivo and in vitro observations that cancer cells ’hijack’ their environment in order to gain migratory, survival, and proliferative advantages.

We considered here that CAFs were non-diffusing, although they are indeed able to spread [42]. Accounting for diffusion of the chemoattractant producer $w$ does not pose, however, any further challenge to the analysis of our model (1.1); in that case our setting belongs to the same mathematical class as the one in [5], to which it is also biologically related: both describe the evolution of a tumor under chemotaxis and haptotaxis, the chemoattractant(s) – of which MDEs are considered in both models – being supposed to diffuse.

Further chemotaxis-haptotaxis models belonging to the class studied here can be considered, of which (1.2) is just one example. As mentioned above, a diffusing signal producer can be easily accomodated to this model class – if the diffusion is linear. Solution-dependent diffusions of either involved species need further investigation; in [33] one such system extending that from [47] to allow for nonlinear diffusion of the chemotactic species has been studied, however without haptotaxis.

Acknowledgement

The authors thank Aydar Uatay (TU Kaiserslautern) for performing the numerical simulations. The second author acknowledges support by Deutsche Forschungsgemeinschaft in the framework of the project Emergence of structures and advantages in cross-diffusion systems (No. 411007140, GZ: WI 3707/5-1).