# Asymptotic behavior of a doubly haptotactic cross-diffusion model for
oncolytic virotherapy

Yifu Wang1,∗, Chi Xu1
School of Mathematics and Statistics, Beijing Institute of Technology
Beijing 100081, P.R. China1

###### Abstract

This paper considers a model for oncolytic virotherapy given by the doubly
haptotactic cross-diffusion system

$\left\\{\begin{array}[]{ll}u_{t}=D_{u}\Delta u-\xi_{u}\nabla\cdot(u\nabla
v)+\mu_{u}u(1-u)-\rho uz,\\\ v_{t}=-(\alpha_{u}u+\alpha_{w}w)v,\\\
w_{t}=D_{w}\Delta w-\xi_{w}\nabla\cdot(w\nabla v)-w+\rho uz,\\\
z_{t}=D_{z}\Delta z-\delta_{z}z-\rho uz+\beta w,\end{array}\right.$

with positive parameters $D_{u},D_{w},D_{z},\xi_{u},\xi_{w},\delta_{z},\rho$,
$\alpha_{u},\alpha_{w},\mu_{u},\beta$. When posed under no-flux boundary
conditions in a smoothly bounded domain $\Omega\subset{\mathbb{R}}^{2}$, and
along with initial conditions involving suitably regular data, the global
existence of classical solution to this system was asserted in Tao and Winkler
(2020). Based on the suitable quasi-Lyapunov functional, it is shown that when
the virus replication rate $\beta<1$, the global classical solution
$(u,v,w,z)$ is uniformly bounded and exponentially stabilizes to the constant
equilibrium $(1,0,0,0)$ in the topology $(L^{\infty}(\Omega))^{4}$ as
$t\rightarrow\infty$.

††∗Corresponding author. E-mail<EMAIL_ADDRESS>(Y.Wang),
<EMAIL_ADDRESS>(C. Xu)

Keywords: Haptotaxis; $LlogL$-estimates; Asymptotic behavior.

2010 Mathematics Subject Classification: 35K57, 35B45, 35Q92, 92C17

##  1 Introduction

As compared to the traditional treatment like chemotherapy or radiotherapy for
cancer diseases, the prominent advantage of virotherapy is that the therapy
can reduce the side-effect on the healthy tissue. In clinical treatments, the
so-call oncolytic viruses (OV) which are either genetically engineered or
naturally occurring can selectively attack the cancer cells and eventually
destroy them without harming normal cells because virus can replicate inside
the infected cells and proceed to infect adjacent cancer cells with the aim to
drive the tumor cells to extinction ([8, 9]). Despite some partial success,
implementation of virotherapy is not in sight. In fact, clinical data reveal
that the efficacy of virotherapy will be reduced by many factors, such as
circulating antibodies, various immune cells or even deposits of extracellular
matrix may essentially decrease ([10, 16]). Therefore, to facilitate the
understanding of the mechanisms that hinder virus spread, the authors of [1]
proposed a mathematical model to describe the interaction between both
uninfected and infected cancer cells, as well as extracellular matrix (ECM)
and oncolytic virus particles, which is given by

$\left\\{\begin{array}[]{lll}u_{t}=D_{u}\Delta u-\xi_{u}\nabla\cdot(u\nabla
v)+\mu_{u}u(1-u)-\rho_{u}uz,&&x\in\Omega,~{}t>0,\\\
v_{t}=-(\alpha_{u}u+\alpha_{w}w)v+\mu_{v}v(1-v),&&x\in\Omega,~{}t>0,\\\
w_{t}=D_{w}\Delta w-\xi_{w}\nabla\cdot(w\nabla
v)-\delta_{w}w+\rho_{w}uz,&&x\in\Omega,~{}t>0,\\\ z_{t}=D_{z}\Delta
z-\delta_{z}z-\rho_{z}uz+\beta w,&&x\in\Omega,~{}t>0,\end{array}\right.$ (1.1)

in a smoothly bounded domain $\Omega\subset\mathbb{R}^{n}$, with positive
parameters $D_{u},D_{w},D_{z},\xi_{u},\xi_{w},\alpha_{u},\alpha_{w}$,
$\mu_{u},\delta_{w},\delta_{z},\beta$ and nonnegative constants
$\mu_{v},\rho_{u},\rho_{w},\rho_{z}$, and with the unknown variables $u,w,z$
and $v$ denoting the population densities of uninfected cancer cells, infected
cancer cells, virus particles and ECM, respectively. Here the crucial modeling
hypothesis underlying (1.1), which accounts for haptotactic motion of cancer
cells and thereby marks a substantial difference between (1.1) and related
more classical reaction-diffusion models for virus dynamics ([13, 19]), is
that apart from its random diffusion, both uninfected and infected cancer
cells bias their motion upward ECM gradients simultaneously due to the
attraction by some macromolecules trapped in the ECM. In addition, the
oncolytic virus particles infect the uninfected cancer cells upon contact with
uninfected tumour cells, and new infectious virus particles are released at
rate $\beta>0$ when infected cells burst (a process known as lysis); beyond
this, (1.1) presupposes that the ECM is degraded upon interacting with both
type of cancer cells, and is possibly remodeled by the normal tissue according
to logistic laws.

Due to its relevance in several biological contexts, inter alia the cancer
invasion ([2, 5]), haptotaxis mechanism has received considerable attention in
the analytical literature ([3, 11, 12, 14, 17, 21, 27, 28, 29, 30, 33, 32]).
The most characteristic ingredient of the model (1.1) is the presence of two
simultaneous haptotaxis processes of cancer cells, and thereby distinguishes
it from most haptotaxis ([11, 29, 33]) and chemotaxis-haptotaxis systems ([3,
17, 27]) studied in the literature, especially the ECM is degraded by both
type of cancer cells in (1.1) directly, rather matrix degrading enzymes (MDEs)
secreted by tumor cells (see [12, 17, 18, 27] for example). It is observed
that the former circumstance seems to widely restrict the accessibility to the
approaches well established in the analysis of related reaction–diffusion
systems, and accordingly the considerable challenges arise for the rigorous
analysis of (1.1), particularly when addressing issues related to qualitative
solution behaviour.

To the best of our knowledge, so far the quantitative comprehension available
for (1.1) is yet mainly limited in some simple setting ([4, 15, 20, 23, 24,
22, 25, 26]). For instance, based on the construction of certain
quasi–Lyapunov functional, Tao and Winkler ([22]) established the global
classical solvability of (1.1) in the two-dimensional case. With respect to
the boundedness of solutions to (1.1), authors in [15] considered some
slightly more comprehensive variants of (1.1), which accounts for the
haptotaxis mechanisms of both cancer cells and virions, in the situation when
zero-order term has suitably strong degradation. Apart from that, existing
analytical works indicate that the virus reproduction rate relative to the
lysis rate of infected cancer cells appears to be critical in determining the
large time behavior of the corresponding solutions at least in some simplified
version of (1.1), inter alia upon neglecting haptotactic cross-diffusion of
infected cancer cells and renew of ECM. Indeed, for the
reaction–diffusion–taxis system

$\left\\{\begin{array}[]{lll}u_{t}=D_{u}\Delta u-\xi_{u}\nabla\cdot(u\nabla
v)+\mu_{u}u(1-u)-\rho uz,&&x\in\Omega,~{}t>0,\\\
v_{t}=-(\alpha_{u}u+\alpha_{w}w)v,&&x\in\Omega,~{}t>0,\\\ w_{t}=D_{w}\Delta
w-w+uz,&&x\in\Omega,~{}t>0,\\\ z_{t}=D_{z}\Delta z-z-uz+\beta
w,&&x\in\Omega,~{}t>0,\end{array}\right.$ (1.2)

it is shown in [24] that if $\beta>1$, then for any reasonably regular initial
data satisfying $\overline{u_{0}}>\frac{1}{\beta-1}$ the global classical
solution of (1.2) with $\rho=0,\mu_{u}=0$ must blow up in infinite time, which
is also implemented by the result on boundedness in the case when
$\overline{u_{0}}<\frac{1}{(\beta-1)_{+}}$ and $v_{0}\equiv 0$ for any
$\beta>0$. Beyond the latter, it was proved that when $\rho>0$ and
$\mu_{u}=0$, the first solution component $u$ of (1.2) possesses a positive
lower bounds whenever $0<\beta<1$ and the initial data $u_{0}\not\equiv 0$
([25]). Furthermore, as a extension of above outcome, the asymptotic behavior
of solution was investigated in [26] if $0<\beta<1$. It is remarked that for
system (1.2) with $\mu_{u}>0$ and $0<\beta<1$, the convergence properties of
the corresponding solutions was also discussed in [4]. We would like to
mention that as the complementing results of [26], the recent paper [23]
reveals that for any prescribed level $\gamma\in(0,\frac{1}{(\beta-1)_{+}})$,
the corresponding solution of (1.2) with $\mu_{u}=0,\rho\geq 0,\beta>0$ will
approach the constant equilibrium $(u_{\infty},0,0,0)$ asymptotically with
some $u_{\infty}>0$ whenever the initial deviation from homogeneous
distribution $(\gamma,0,0,0)$ is suitably small.

The purpose of this work is to investigate the dynamical features of the
models involving the simultaneous haptotactic processes of both uninfected and
infected cancer cells when the virus replication rate $\beta<1$. To this end,
we are concerned with the comprehensive haptotactic cross-diffusion systems of
the form

$\left\\{\begin{array}[]{ll}u_{t}=D_{u}\Delta u-\xi_{u}\nabla\cdot(u\nabla
v)+\mu_{u}u(1-u)-\displaystyle\rho uz,&x\in\Omega,t>0,\\\
\displaystyle{w_{t}=D_{w}\Delta w-\xi_{w}\nabla\cdot(w\nabla
v)-w+\displaystyle\rho uz},&x\in\Omega,t>0,\\\
\displaystyle{v_{t}=-(\alpha_{u}u+\alpha_{w}w)v},&x\in\Omega,t>0,\\\
\displaystyle{z_{t}=D_{z}\Delta z-\delta_{z}z-\displaystyle\rho uz+\beta
w},&x\in\Omega,t>0,\\\ \displaystyle{(D_{u}\nabla u-\xi_{u}u\nabla
v)\cdot\nu=(D_{w}\nabla w-\xi_{w}w\nabla v)\cdot\nu=\nabla
z\cdot\nu=0},&x\in\partial\Omega,t>0,\\\
\displaystyle{u(x,0)=u_{0}(x)},w(x,0)=w_{0}(x),v(x,0)=v_{0}(x),z(x,0)=z_{0}(x),&x\in\Omega\end{array}\right.$
(1.3)

in a smoothly bounded domain $\Omega\subset\mathbb{R}^{2}$. To import the
precise framework underlying the basic theory from [22] we shall henceforth
assume that

$\left\\{\begin{array}[]{ll}\displaystyle{u_{0},w_{0},z_{0}~{}\hbox{and}\,v_{0}~{}\hbox{are
nonnegative functions from}~{}C^{2+\vartheta}(\bar{\Omega})~{}\hbox{for
some}~{}\vartheta\in(0,1),}\\\ \displaystyle{\mbox{with}~{}u_{0}\not\equiv
0,~{}w_{0}\not\equiv 0,~{}z_{0}\not\equiv 0,~{}v_{0}\not\equiv
0~{}\hbox{and}~{}\frac{\partial
w_{0}}{\partial\nu}=0~{}~{}\mbox{on}~{}~{}\partial\Omega.}\\\
\end{array}\right.$ (1.4)

Hence the outcome of [22] asserts the global existence of a unique classical
solution $(u,v,w,z)$ to (1.3). Our main results reveal that whenever
$\beta<1$, $(u,v,w,z)$ is uniformly bounded and exponentially converges to the
constant equilibrium $(1,0,0,0)$ in the topology $(L^{\infty}(\Omega))^{4}$ in
a large time limit, which can be stated as follows

######  Theorem 1.1.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary,
$D_{u},D_{w},D_{z},\xi_{u},\xi_{w}$,
$\mu_{u},\rho,\alpha_{u},\alpha_{w},\delta_{z}$ are positive parameters, and
suppose that $0<\beta<1$. Then system (1.3) admits a unique global classic
positive solution satisfying

$\sup\limits_{t>0}\left\\{\|u(\cdot,t)\|_{L^{\infty}(\Omega)}+\|v(\cdot,t)\|_{L^{\infty}(\Omega)}+\|w(\cdot,t)\|_{L^{\infty}(\Omega)}+\|z(\cdot,t)\|_{L^{\infty}(\Omega)}\right\\}<\infty.$
(1.5)

Moreover there exist positive constants $\eta,\varrho,\gamma_{1}$,
$\gamma_{2}$ and $C>0$ such that

$\|u(\cdot,t)-1\|_{L^{\infty}(\Omega)}\leq Ce^{-\eta t},$ (1.6)
$\|w(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\varrho t},$ (1.7)
$\|z(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\gamma_{1}t}$ (1.8)

as well as

$\|v(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\gamma_{2}t}.$ (1.9)

Since the third equation in (1.3) is merely an ordinary differential equation,
no smoothing action on the spatial regularity of $v$ can be expected. To
overcome the analytical difficulties arising from the latter, inter alia in
the derivation of global boundedness of solutions, we accordingly introduce
the variable transformation $a=e^{-\chi_{u}v}u$ and $b=e^{-\chi_{w}v}w$, and
establish a priori estimate for the solution components $a,b$ of the
corresponding equivalent system (2.2) below in the space $LlogL(\Omega)$
rather than the solution components $u,w$ to (1.3). Note that in the evolution
of density $v$ of ECM fibers the quantity $u,w$ appears via a sink term,
whereas it turns to a genuine superlinear production terms of system (2.2) in
the style of $\chi_{u}a(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v$
and $\chi_{w}b(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v$. Taking
advantage of the exponential decay of $w$ in $L^{1}$ norm in the case
$0<\beta<1$, we shall track the time evolution of

$\mathcal{F}(t)=\int_{\Omega}e^{\chi_{u}v}a(\cdot,t)\log
a(\cdot,t)+\int_{\Omega}e^{\chi_{w}v}b(\cdot,t)\log
b(\cdot,t)+\int_{\Omega}z^{2}(\cdot,t)$

with $a=e^{-\chi_{u}v}u$ and $b=e^{-\chi_{w}v}w$, which is somewhat different
from the quasi-Lyapunov functional (4.11) in [22] where a Dirichlet integral
of $\sqrt{v}$ is involved. Here the quadratic degradation term in the first
equation of (1.3) seems to be necessary. Thereafter applying a variant of the
Gagliardo–Nirenberg inequality involving certain $LlogL$-type norms and
performing a Moser-type iteration, the $L^{\infty}$-bounds of solutions is
derived.

In addition, our result indicates that although haptotaxis mechanism may have
some important influence on the properties of the related system on short or
intermediate time scales, the large time behavior of solution to (1.3) can
essentially be described by the corresponding haptotaxis-free system at least
under the biological meaningful restriction $\beta<1$. In order to prove
Theorem 1.1, a first step is to derive a pointwise lower bound for
$a:=e^{-\chi_{u}v}u$ (Lemma 4.2), which, in turn, amounts to establishing an
exponential decay of $z$ with respect to the norm in $L^{\infty}(\Omega)$
(Lemma 4.2). To achieve the latter, we will make use of $L^{1}(\Omega)$-decay
information of $w,z$ explicitly contained in Lemma 2.3. Scondly, as a
consequence of the former, the exponential decay of $v$ with respect to
$L^{\infty}(\Omega)$ norm is achieved (Lemma 4.3), which along a
$a^{-1}$–testing procedure will provide quite weak convergence information of
$u$, inter alia the integrability property of $\nabla\sqrt{a}$ in
$L^{2}((0,\infty);L^{2}(\Omega))$ (Lemma 4.4). The next step will consist of
verifying the integrability of $\nabla v$ in $L^{2}((0,\infty);L^{2}(\Omega))$
rather than that of $a_{t}$ (Lemma 4.5), which will turn out to be sufficient
a condition in the derivation of exponential decay property of
$\|u(\cdot,t)-1\|_{L^{p}(\Omega)}$. Indeed, this integrability property of
$\nabla v$ enables us to derive an exponential decay of $\int_{\Omega}|\nabla
v|^{2}$ (Lemma 4.6), upon which and through a testing procedure, it is shown
that the convergence property of $u$ actually takes place in the type of
(4.22) (Lemma 4.7). Further, upon the above decay properties, we are able to
verify that $\int_{\Omega}|\nabla v|^{4}$ decays exponentially by means of the
suitable quasi-Lyapunov functional (Lemma 4.9). At this position, thanks to
the integrability exponent in $\int_{\Omega}|\nabla v|^{4}$ exceeding the
considered spatial dimension $n=2$, the desired decay property stated in
Theorem 1.1 can be exactly achieved.

This paper will be organized as follows: Section 2 will introduce an
equivalent system of (1.3) and give out some basic priori estimates of
classical solutions thereof, inter alia the weak decay properties of $w,z$.
Section 3 will focus on the construction of an entropy-type functional, which
entails certain $LlogL$-type norms and thereby allows us to establish the
$L^{\infty}$-bounds. Finally, starting from the exponential decay of
quantities $w,z$ with respect to the norm in $L^{1}(\Omega)$, we established
the exponential convergence properties of the solutions in Section 4.

##  2 Preliminaries

Let us firstly recall the result in [22] which warrants the global smooth
solvability of problem (1.3).

######  Lemma 2.1.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary,
$D_{u},D_{w},D_{z},\xi_{u},\xi_{w}$,
$\mu_{u},\rho,\alpha_{u},\alpha_{w},\delta_{z},\beta$ are positive parameters.
Then for any choice of $(u_{0},v_{0},w_{0},z_{0})$ fulfilling (1.4), the
problem (1.3) (1.4) possesses a uniquely determined classical solution
$(u,v,w,z)\in(C^{2,1}(\overline{\Omega}\times[0,\infty)))^{4}$ for which
$u>0,w>0,z>0$ and $v\geq 0$.

Following the variable of change used in related literature ([6, 15, 26, 28]),
which can conveniently reformulate the haptotactic interaction in (1.3), we
define $\chi_{u}:=\frac{\xi_{u}}{D_{u}}$ and $\chi_{w}:=\frac{\xi_{w}}{D_{w}}$
and set

$a:=e^{-\chi_{u}v}u~{}~{}\hbox{and}~{}~{}b:=e^{-\chi_{w}v}w.$ (2.1)

Then we transform (1.3) into an equivalent system as below

$\left\\{\begin{array}[]{ll}a_{t}=D_{u}e^{-\chi_{u}v}\nabla\cdot(e^{\chi_{u}v}\nabla
a)+f(a,b,v,c),&x\in\Omega,t>0,\\\
b_{t}=D_{w}e^{-\chi_{w}v}\nabla\cdot(e^{\chi_{w}v}\nabla
b)+g(a,b,v,c),&x\in\Omega,t>0,\\\
v_{t}=-(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v,&x\in\Omega,t>0,\\\
z_{t}=D_{z}\Delta z-\delta_{z}z-\rho_{z}uz+\beta w,&x\in\Omega,t>0,\\\
\displaystyle\frac{\partial a}{\partial\nu}=\frac{\partial
b}{\partial\nu}=\frac{\partial z}{\partial\nu}=0,&x\in\partial\Omega,t>0,\\\
a(x,0)=u_{0}(x)e^{-\chi_{u}v_{0}(x)},~{}b(x,0)=w_{0}(x)e^{-\chi_{w}v_{0}(x)},&x\in\Omega,\\\
v_{0}(x,0)=v_{0}(x),~{}~{}z(x,0)=z_{0}(x),&x\in\Omega\end{array}\right.$ (2.2)

with

$f(a,b,v,c):=\mu_{u}a(1-ae^{\chi_{u}v})-\displaystyle\rho
az+\chi_{u}a(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v,$

as well as

$g(a,b,v,c):=-b+\displaystyle\rho
aze^{(\chi_{u}-\chi_{w})v}+\chi_{w}b(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v.$

In our subsequent analysis, unless otherwise stated we shall assume that
$(a,v,b,z)$ is the global classical solution to (2.2) addressed in Lemma 2.1.

The damping effects of quadratic degradation in the first equation in (1.3)
will be important for us to verify the global boundedness of the solutions.
Let we first apply straightforward argument to achieve the following basic
$L^{1}$-bounds for $u,w$ and $z$, which is also valid for the solution
components $a,b$ of (2.2).

######  Lemma 2.2.

For all $t>0$, the solution $(u,w,v,z)$ satisfies

$\int_{\Omega}u(\cdot,t)\leq\max\left\\{\int_{\Omega}u_{0},|\Omega|\right\\}:=m_{u},$
(2.3)

and

$v\leq\|v_{0}\|_{L^{\infty}(\Omega)}:=m_{v}$ (2.4)

and

$\int_{\Omega}w(\cdot,t)\leq\max\left\\{\|u_{0}\|_{L^{1}(\Omega)}+\|w_{0}\|_{L^{1}(\Omega)},\frac{|\Omega|\mu_{u}}{\min\\{1,\mu_{u}\\}}\right\\}:=m_{w},$
(2.5)

as well as

$\int_{\Omega}z(\cdot,t)\leq\max\left\\{\int_{\Omega}z_{0},\frac{\beta
m_{w}}{\delta_{z}}\right\\}:=m_{z}.$ (2.6)

###### Proof.

It is easy to see that (2.3) can be derived through an integration of the
first equation in (1.3) along with Cauchy–Schwarz’s inequality, and (2.4) is a
direct consequence of $(\alpha_{u}u+\alpha_{w}w)v\geq 0$ due to the
nonnegativity of $u,w$ and $v$.

In addition, integrating the $w$-equation as well as $u$-equation respectively
and adding the corresponding results, we then have

$\frac{d}{dt}\left(\int_{\Omega}u+\int_{\Omega}w\right)+\int_{\Omega}w+\mu_{u}\int_{\Omega}u\leq|\Omega|\mu_{u},$
(2.7)

which readily leads to (2.5) upon an ODE comparison. At last, thanks to (2.5),
(2.6) clearly results from the integration of $z$-equation in (1.3). ∎

Beyond that, making use of the restriction $\beta\in(0,1)$, one can derive the
decay properties of the solution components $w$ and $z$ with respect to
$L^{1}(\Omega)$, which will be used later on.

######  Lemma 2.3.

Suppose that $0<\beta<1$, then there exists constant $C>0$ such that

$\int_{\Omega}w(\cdot,t)+\int_{\Omega}z(\cdot,t)\leq Ce^{-\delta
t}~{}~{}\hbox{for all $t>0$.}$ (2.8)

with $\delta=\min\\{1-\beta,\delta_{z}\\}$.

###### Proof.

We use the $z$-equation and $w$-equation to compute

$\frac{d}{dt}\left(\int_{\Omega}w+\int_{\Omega}z\right)+(1-\beta)\int_{\Omega}w+\delta_{z}\int_{\Omega}z=0.$
(2.9)

Due to $0<\beta<1$, this readily implies that

$\left(\int_{\Omega}w(\cdot,t)+\int_{\Omega}z(\cdot,t)\right)\leq\left(\int_{\Omega}w_{0}+\int_{\Omega}z_{0}\right)e^{-\min\\{1-\beta,\delta_{z}\\}t}$

and hence (2.8) is valid with $C=\int_{\Omega}w_{0}+\int_{\Omega}z_{0}$. ∎

##  3 Global boundedness

As in [15, 28], the crucial step in establishing a priori $L^{\infty}$ bounds
for $a,b$ and $z$ is to derive estimates for $a$ and $b$ in $LlogL$, which
turn out to be consequences of a quasi-energy structure associated with the
system (2.2) rather than the system (1.3). Indeed, making appropriate use of
the logistic degradation in the first equation of (2.2) and inter alia the
$L^{1}$-decay property of the solution component $w$, one can verify that
functional

$\mathcal{F}(t):=\int_{\Omega}e^{\chi_{u}v}a(\cdot,t)\log
a(\cdot,t)+\int_{\Omega}e^{\chi_{w}v}b(\cdot,t)\log
b(\cdot,t)+\frac{1}{2}\int_{\Omega}z^{2}(\cdot,t),$

which does not involve the Dirichlet integral of $\sqrt{v}$, actually
possesses a certain quasi-dissipative property for all $t>t_{0}$ with constant
$t_{0}>1$ suitably chosen.

######  Lemma 3.1.

For any $\varepsilon>0$, there exists $C(\varepsilon)>0$ such that

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a\log
a+\int_{\Omega}e^{\chi_{u}v}a\log
a+D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a}+\frac{\mu_{u}}{2}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log
a\\\\[5.69054pt]
\leq&\displaystyle\varepsilon\int_{\Omega}b^{2}+C(\varepsilon)\end{array}$
(3.1)

for all $t>0$.

###### Proof.

From the first equation in (1.3), it follows that

$(e^{\chi_{u}v}a)_{t}=D_{u}\nabla\cdot(e^{\chi_{u}v}\nabla
a)+\mu_{u}e^{\chi_{u}v}a(1-e^{\chi_{u}v}a)-\rho e^{\chi_{u}v}az.$ (3.2)

Hence, a testing procedure on the first equation in (2.2) leads to

$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a\log
a&=\displaystyle\int_{\Omega}(e^{\chi_{u}v}a)_{t}\log
a+\int_{\Omega}e^{\chi_{u}v}a_{t}\\\ &\leq\displaystyle-
D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a}+\mu_{u}\int_{\Omega}e^{\chi_{u}v}a\log
a-\mu_{u}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log a\\\
&\displaystyle-\rho\int_{\Omega}e^{\chi_{u}v}za\log
a+\mu_{u}\int_{\Omega}ae^{\chi_{u}v}-\rho\int_{\Omega}e^{\chi_{u}v}za\\\
&\displaystyle+\chi_{u}\alpha_{u}\int_{\Omega}e^{2\chi_{u}v}va^{2}+\chi_{u}\alpha_{w}\int_{\Omega}e^{(\chi_{w}+\chi_{u})v}vab.\end{array}$

Thanks to Lemma 2.2 and the elementary inequality $a\log a\geq-\frac{1}{e}$
valid for all $a>0$, one can find $c_{1}>0$ such that

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a\log
a+D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a}+\mu_{u}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log a\\\\[5.69054pt]
\leq&\displaystyle\alpha_{u}\chi_{u}e^{2\chi_{u}m_{v}}m_{v}\int_{\Omega}a^{2}+\alpha_{w}\chi_{u}e^{(\chi_{u}+\chi_{w})m_{v}}m_{v}\int_{\Omega}ab+\mu_{u}\int_{\Omega}e^{\chi_{u}v}a\log
a+c_{1},\end{array}$

which along with Young’s inequality implies that for any $\varepsilon>0$

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a\log
a+D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a}+\mu_{u}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log a\\\\[5.69054pt]
\leq&(\displaystyle\alpha_{u}\chi_{u}e^{2\chi_{u}m_{v}}m_{v}+\frac{1}{\varepsilon}\alpha_{w}^{2}\chi_{u}^{2}e^{2(\chi_{u}+\chi_{w})m_{v}}m_{v}^{2})\int_{\Omega}a^{2}+\varepsilon\int_{\Omega}b^{2}+\mu_{u}\int_{\Omega}e^{\chi_{u}v}a\log
a+c_{1}.\end{array}$

Further invoking the inequality $a^{2}\leq\varepsilon_{1}a^{2}\log
a+e^{\frac{2}{\varepsilon_{1}}}$ for any $\varepsilon_{1}>0$, we arrive at

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a\log
a+D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a}+\frac{3\mu_{u}}{4}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log a\\\
\leq&\varepsilon\displaystyle\int_{\Omega}b^{2}+\mu_{u}\int_{\Omega}e^{\chi_{u}v}a\log
a+c_{2}(\varepsilon)\end{array}$ (3.3)

with some $c_{2}(\varepsilon)>0$. Accordingly, (3.1) is a consequence of (3.3)
and the fact that $a\log a\leq\varepsilon_{2}a^{2}\log
a-\varepsilon_{2}^{-1}\ln\varepsilon_{2}$ with
$\varepsilon_{2}=\frac{\mu_{u}}{4(\mu_{u}+1)}$. ∎

For the solution component $b$ of (2.2), we also have

######  Lemma 3.2.

Let $0<\beta<1$. Then one can find $C>0$ and $t_{0}>0$ such that for all
$t>t_{0}$

$\displaystyle\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b\log
b+\int_{\Omega}e^{\chi_{w}v}b\log
b+\frac{D_{w}}{2}\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}$ $\displaystyle\leq$
$\displaystyle\int_{\Omega}a^{2}+C(\int_{\Omega}z^{4})^{\frac{1}{2}}.$ (3.4)

###### Proof.

From the second equation in (1.3), it follows that

$(be^{\chi_{w}v})_{t}=D_{w}\nabla\cdot(e^{\chi_{w}v}\nabla b)-w+\rho uz.$

Relying on $0\leq v\leq m_{v}$ in $\overline{\Omega}\times(0,\infty)$, a
straightforward calculation along with the Young inequality yields

$\displaystyle\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b\log
b+\int_{\Omega}e^{\chi_{w}v}b\log
b+D_{w}\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}$ $\displaystyle=$ $\displaystyle\int_{\Omega}\log(eb)(\rho
uz-w)+\chi_{w}\displaystyle\int_{\Omega}w(\alpha_{u}u+\alpha_{w}w)v+\int_{\Omega}w\log
b$ (3.5) $\displaystyle\leq$
$\displaystyle\alpha_{w}\chi_{w}e^{2\chi_{w}m_{v}}m_{v}\int_{\Omega}b^{2}+\rho\int_{\Omega}aze^{\chi_{u}v}\log
b+\rho
e^{\chi_{u}m_{v}}\int_{\Omega}az+\alpha_{u}\chi_{w}e^{(\chi_{u}+\chi_{w})m_{v}}m_{v}\int_{\Omega}ab$
$\displaystyle\leq$ $\displaystyle
c_{1}\int_{\Omega}b^{2}+\int_{\Omega}a^{2}+c_{1}\int_{\Omega}z^{2}+c_{1}\displaystyle\int_{\\{x\in\Omega;b(x,t)\geq
1\\}}z^{2}\log^{2}b$ $\displaystyle\leq$ $\displaystyle
c_{1}\int_{\Omega}b^{2}+\int_{\Omega}a^{2}+c_{1}\int_{\Omega}z^{2}+c_{1}\displaystyle(\int_{\Omega}z^{4})^{\frac{1}{2}}(\int_{\\{x\in\Omega;b(x,t)\geq
1\\}}|\log b|^{4})^{\frac{1}{2}}$ $\displaystyle\leq$ $\displaystyle
c_{1}\int_{\Omega}b^{2}+\int_{\Omega}a^{2}+c_{1}\int_{\Omega}z^{2}+c_{1}\displaystyle(\int_{\Omega}z^{4})^{\frac{1}{2}}(\int_{\Omega}b+c_{2}|\Omega|)^{\frac{1}{2}}$

with some $c_{1}>0$, where we use the fact that there exits $c_{2}>0$ such
that $\log^{4}s\leq s+c_{2}$ for all $s\geq 1$.

Furthermore, in order to appropriately estimate the first summand on the
right-hand of (3), we apply the two-dimensional Gagliardo–Nirenberg
inequalities

$\|\varphi\|^{4}_{L^{4}(\Omega)}\leq
C_{g}\|\nabla\varphi\|_{L^{2}(\Omega)}^{2}\|\varphi\|^{2}_{L^{2}(\Omega)}+C_{g}\|\varphi\|^{4}_{L^{2}(\Omega)}~{}~{}~{}~{}\hbox{for
some}~{}C_{g}>0~{}\hbox{and all}~{}\varphi\in W^{1,2}(\Omega)$

to get

$\begin{array}[]{rl}c_{1}\displaystyle\int_{\Omega}b^{2}=&c_{1}\|\sqrt{b}\|^{4}_{L^{4}(\Omega)}\\\
\leq&\displaystyle
c_{1}C_{g}\displaystyle\int_{\Omega}b\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}+c_{1}C_{g}(\displaystyle\int_{\Omega}b)^{2}\\\ \leq&\displaystyle
c_{1}C_{g}\displaystyle\int_{\Omega}w\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}+c_{1}C_{g}(\displaystyle\int_{\Omega}w)^{2}\\\ \end{array}$ (3.6)

Therefore combining (3.6) with (3), we arrive at

$\displaystyle\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b\log
b+\int_{\Omega}e^{\chi_{w}v}b\log
b+D_{w}\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}$ (3.7) $\displaystyle\leq$ $\displaystyle\displaystyle
c_{1}C_{g}\displaystyle\int_{\Omega}w\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}+c_{1}C_{g}(\displaystyle\int_{\Omega}w)^{2}+\int_{\Omega}a^{2}+c_{1}\int_{\Omega}z^{2}+c_{1}\displaystyle(\int_{\Omega}z^{4})^{\frac{1}{2}}(\int_{\Omega}w+c_{2}|\Omega|)^{\frac{1}{2}}.$

By Lemma 2.3, we can pick $t_{0}>0$ suitably large such that

$c_{1}C_{g}\displaystyle\int_{\Omega}w\leq\frac{D_{w}}{2}$

and thereby for $t\geq t_{0}$,

$\displaystyle\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b\log
b+\int_{\Omega}e^{\chi_{w}v}b\log
b+\frac{D_{w}}{2}\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}$ $\displaystyle\leq$ $\displaystyle
c_{1}C_{g}(\displaystyle\int_{\Omega}w)^{2}+\int_{\Omega}a^{2}+c_{1}\int_{\Omega}z^{2}+c_{1}\displaystyle(\int_{\Omega}z^{4})^{\frac{1}{2}}(\int_{\Omega}w+c_{2}|\Omega|)^{\frac{1}{2}},$

which along with Lemma 2.2 and the Young inequality completes the proof. ∎

Whereas the expressions $\int_{\Omega}a^{2}$ appearing in (3.2) turns out to
be conveniently digestible through the dissipation rate in (3.1), it remains
to estimate $(\int_{\Omega}z^{4})^{\frac{1}{2}}$ by means of an interpolation
argument.

######  Lemma 3.3.

Let $0<\beta<1$ and define

$\mathcal{F}(t):=\int_{\Omega}e^{\chi_{u}v}a(\cdot,t)\log
a(\cdot,t)+\int_{\Omega}e^{\chi_{w}v}b(\cdot,t)\log
b(\cdot,t)+\int_{\Omega}z^{2}(\cdot,t).$

Then there exists constant $C>0$ such that

$\mathcal{F^{\prime}}(t)+\mathcal{F}(t)\leq C$ (3.8)

for all $t\geq t_{0}$ as given in Lemma 3.2.

###### Proof.

Testing the fourth equation in (2.2) by $z$, we get

$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int_{\Omega}z^{2}+\int_{\Omega}z^{2}+2D_{z}\int_{\Omega}|\nabla
z|^{2}&\leq\displaystyle
2\beta\int_{\Omega}e^{\chi_{w}v}bz+\int_{\Omega}z^{2}\\\\[5.69054pt]
&\leq\displaystyle\varepsilon\int_{\Omega}b^{2}+(1+\frac{e^{2\chi_{w}m_{v}}\beta^{2}}{\varepsilon})\int_{\Omega}z^{2}.\end{array}$

Hence Lemma 3.1 and Lemma 3.2 provide positive constants
$c_{i}(\varepsilon)>0(i=1,2)$ such that

$\begin{array}[]{rl}&\mathcal{F^{\prime}}(t)+\mathcal{F}(t)+\displaystyle
D_{w}\int_{\Omega}e^{\chi_{w}v}\frac{|\nabla
b|^{2}}{b}+D_{z}\int_{\Omega}|\nabla
z|^{2}+\frac{\mu_{u}}{2}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log
a+\int_{\Omega}z^{2}\\\
\leq&\displaystyle\int_{\Omega}a^{2}+c_{1}(\varepsilon)(\int_{\Omega}z^{4})^{\frac{1}{2}}+2\varepsilon\int_{\Omega}b^{2}+c_{2}(\varepsilon)\end{array}$
(3.9)

for all $t\geq t_{0}$ with $t_{0}$ given by Lemma 3.2.

Now again since $a^{2}\leq\varepsilon_{1}a^{2}\log
a+e^{\frac{2}{\varepsilon_{1}}}$ for any $\varepsilon_{1}>0$,

$\displaystyle\int_{\Omega}a^{2}\leq\frac{\mu_{u}}{4}\int_{\Omega}e^{2\chi_{u}v}a^{2}\log
a+c_{3}$ (3.10)

with $c_{3}>0$, whereas according to the two-dimensional Gagliardo–Nirenberg
inequalities,

$\begin{array}[]{rl}2\varepsilon\displaystyle\int_{\Omega}b^{2}=&2\varepsilon\|\sqrt{b}\|^{4}_{L^{4}(\Omega)}\\\
\leq&\displaystyle 2\varepsilon
C_{g}\displaystyle\int_{\Omega}b\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}+2\varepsilon C_{g}(\displaystyle\int_{\Omega}b)^{2}\\\
\leq&\displaystyle 2\varepsilon
m_{w}\displaystyle\int_{\Omega}e^{\chi_{w}v}\displaystyle\frac{|\nabla
b|^{2}}{b}+2\varepsilon m_{w}^{2}C_{g}\\\
\leq&\displaystyle\frac{D_{w}}{2}\int_{\Omega}e^{\chi_{w}v}\frac{|\nabla
b|^{2}}{b}+c_{4}\end{array}$ (3.11)

by the choice of $\varepsilon=\frac{D_{w}}{4m_{w}}$.

In summary, (3.11), (3.10) and (3.9) show that for $t\geq t_{0}$

$\mathcal{F^{\prime}}(t)+\mathcal{F}(t)+D_{z}\int_{\Omega}|\nabla
z|^{2}+\int_{\Omega}z^{2}\leq c_{1}(\int_{\Omega}z^{4})^{\frac{1}{2}}+c_{5}$
(3.12)

with some $c_{5}>0$.

Furthermore to estimate $\left(\int_{\Omega}z^{4}\right)^{\frac{1}{2}}$ on the
right-hand of (3.12), we employ the Gagliardo–Nirenberg inequalities once more
to obtain that

$\begin{array}[]{rl}c_{1}\|z(\cdot,t)\|^{2}_{L^{4}(\Omega)}&\leq
c_{1}c_{g}\|\nabla
z(\cdot,t)\|^{\frac{3}{2}}_{L^{2}(\Omega)}\|z(\cdot,t)\|^{\frac{1}{2}}_{L^{1}(\Omega)}+c_{1}c_{g}\|z(\cdot,t)\|^{2}_{L^{1}(\Omega)}\\\
&\leq\displaystyle\frac{D_{z}}{2}\int_{\Omega}|\nabla z|^{2}+c_{6}\end{array}$
(3.13)

with $c_{6}>0$, which together with (3.12) readily establishes (3.8). ∎

As a consequence of (3.8), the $LlogL$-estimate of quantities $a$ and $b$ is
achieved as follows

######  Lemma 3.4.

Let $0<\beta<1$. Then there exists $C>0$ such that for all $t>t_{0}$,

$\int_{\Omega}a(\cdot,t)|\log a(\cdot,t)|\leq C$ (3.14)

as well as

$\int_{\Omega}b(\cdot,t)|\log b(\cdot,t)|\leq C.$ (3.15)

###### Proof.

As in Lemma 3.4 of [15], by the inequality $a\log a>-e^{-1}$ in
$\Omega\times(0,\infty)$, we have

$\begin{array}[]{rl}\displaystyle\int_{\Omega}a(\cdot,t)|\log
a(\cdot,t)|&\leq\displaystyle\int_{\Omega}e^{\chi_{u}v}a(\cdot,t)\log
a(\cdot,t)-2\int_{a<1}e^{\chi_{u}v}a(\cdot,t)\log a(\cdot,t)\\\\[5.69054pt]
&\leq\displaystyle\int_{\Omega}e^{\chi_{u}v}a(\cdot,t)\log
a(\cdot,t)+\frac{2|\Omega|e^{\chi_{u}m_{v}}}{e}.\end{array}$

Likewise, we can also obtain

$\int_{\Omega}b(\cdot,t)|\log
b(\cdot,t)|\leq\int_{\Omega}e^{\chi_{w}v}b(\cdot,t)\log
b(\cdot,t)+\frac{2|\Omega|e^{\chi_{w}m_{v}}}{e}.$

According to (3.8), there exits $c_{1}>0$ such that

$\mathcal{F}(t)\leq c_{1}.$ (3.16)

Therefore by the definition of $\mathcal{F}(t)$, we can see that

$\displaystyle\int_{\Omega}a(\cdot,t)|\log
a(\cdot,t)|+\int_{\Omega}b(\cdot,t)|\log
b(\cdot,t)|\leq\mathcal{F}(t)+\frac{2|\Omega|e^{\chi_{u}m_{v}}}{e}+\frac{2|\Omega|e^{\chi_{w}m_{v}}}{e}\leq
c_{2}$

with $c_{2}=\frac{2|\Omega|}{e}(e^{\chi_{w}m_{v}}+e^{\chi_{u}m_{v}})+c_{1}$
for all $t>t_{0}$, and thus complete the proof. ∎

The a priori estimates for $a,b$ gained in Lemma 3.4 is the cornerstone to
establish a $L^{\infty}(\Omega)$-bound for solution $(u,v,w,z)$. Indeed, one
can proceed to derive obtain $L^{\infty}(\Omega)$-bound by means of some quite
straightforward $L^{p}$ testing procedures.

######  Lemma 3.5.

Let $0<\beta<1$. Then there exists $C>0$ such that for all $t>0$

$\|a(\cdot,t)\|_{L^{\infty}(\Omega)}+\|b(\cdot,t)\|_{L^{\infty}(\Omega)}+\|z(\cdot,t)\|_{L^{\infty}(\Omega)}\leq
C.$ (3.17)

###### Proof.

Testing the first equation in (2.2) by $e^{\chi_{u}v}a^{r-1}$ with $r\geq 2$,
we obtain $c_{1}(r)>0$

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a^{r}+\frac{4(r-1)D_{u}}{r}\int_{\Omega}e^{\chi_{u}v}|\nabla
a^{\frac{r}{2}}|^{2}+r\mu_{u}\int_{\Omega}e^{2\chi_{v}}a^{r+1}\\\
\leq&\displaystyle\mu_{u}r\int_{\Omega}e^{\chi_{u}v}a^{r}+r\chi_{u}\int_{\Omega}e^{\chi_{u}v}a^{r}(\alpha_{u}u+\alpha_{w}w)v\\\
\leq&\displaystyle\mu_{u}r\int_{\Omega}e^{\chi_{u}v}a^{r}+c_{1}(r)\int_{\Omega}a^{r}(a+b),\end{array}$

which together with the Young inequality, leads to

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}a^{r}+2D_{u}\int_{\Omega}|\nabla
a^{\frac{r}{2}}|^{2}+\int_{\Omega}e^{\chi_{u}v}a^{r}\\\ \leq&\displaystyle
c_{2}(r)\int_{\Omega}a^{r+1}+c_{2}(r)\int_{\Omega}b^{r+1}+c_{2}(r)\end{array}$
(3.18)

with constant $c_{2}(r)>0$. Likewise, there exist $c_{i}(r)>0(i=3,4)$ such
that

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b^{r}+2D_{w}\int_{\Omega}|\nabla
b^{\frac{r}{2}}|^{2}+r\int_{\Omega}e^{\chi_{v}}b^{r}\\\ \leq&\displaystyle
c_{3}(r)\int_{\Omega}b^{r-1}az+c_{3}(r)\int_{\Omega}a^{r+1}+c_{3}(r)\int_{\Omega}b^{r+1}\\\\[5.69054pt]
\leq&\displaystyle
c_{4}(r)\int_{\Omega}a^{r+1}+c_{4}(r)\int_{\Omega}b^{r+1}+c_{4}(r)\int_{\Omega}z^{r+1}\end{array}$
(3.19)

as well as

$\displaystyle\frac{d}{dt}\int_{\Omega}z^{r}+2D_{z}\int_{\Omega}|\nabla
z^{\frac{r}{2}}|^{2}+\frac{r\delta_{z}}{2}\int_{\Omega}z^{r}\leq\frac{2r\beta^{r}}{\delta_{z}}\int_{\Omega}w^{r}.$
(3.20)

Collecting (3.18)–(3.20), we then arrive at

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})+c_{5}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})+c_{5}\int_{\Omega}(|\nabla
a^{\frac{r}{2}}|^{2}+|\nabla b^{\frac{r}{2}}|^{2}+|\nabla
z^{\frac{r}{2}}|^{2})\\\ \leq&\displaystyle
c_{6}(r)(\int_{\Omega}a^{r+1}+\int_{\Omega}b^{r+1}+\int_{\Omega}z^{r+1})+c_{6}(r)\end{array}$
(3.21)

with $c_{5}>0$ and $c_{6}(r)>0$.

Now we invoke the logarithm-type Gagliardo–Nirenberg inequality (we refer to
Lemma A.5 in [28] for details) to obtain that for any $\varepsilon>0$,

$\int_{\Omega}a^{r+1}\leq\varepsilon\|\nabla
a^{\frac{r}{2}}\|^{2}_{L^{2}(\Omega)}\cdot\int_{\Omega}a\cdot|\log|a^{\frac{r}{2}}||+C(\varepsilon,r)\left(\|a^{\frac{r}{2}}\|^{\frac{2(r+1)}{r}}_{L^{\frac{2}{r}}(\Omega)}+1\right),$
(3.22)

and

$\int_{\Omega}b^{r+1}\leq\varepsilon\|\nabla
b^{\frac{r}{2}}\|^{2}_{L^{2}(\Omega)}\cdot\int_{\Omega}b\cdot|\log|b^{\frac{r}{2}}||+C(\varepsilon,r)\left(\|b^{\frac{r}{2}}\|^{\frac{2(r+1)}{r}}_{L^{\frac{2}{r}}(\Omega)}+1\right)$
(3.23)

as well as

$\int_{\Omega}z^{r+1}\leq\varepsilon\|\nabla
z^{\frac{r}{2}}\|^{2}_{L^{2}(\Omega)}\cdot\int_{\Omega}z\cdot|\log|z^{\frac{r}{2}}||+C(\varepsilon,r)\left(\|z^{\frac{r}{2}}\|^{\frac{2(r+1)}{r}}_{L^{\frac{2}{r}}(\Omega)}+1\right).$
(3.24)

Therefore combining (3.21)–(3.24) with Lemma 3.4, one can find $c_{7}(r)>0$
and $c_{8}(\varepsilon,r)>0$ such that

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})+c_{5}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})+c_{5}\int_{\Omega}(|\nabla
a^{\frac{r}{2}}|^{2}+|\nabla b^{\frac{r}{2}}|^{2}+|\nabla
z^{\frac{r}{2}}|^{2})\\\\[5.69054pt] \leq&\displaystyle
c_{7}(r)\varepsilon\int_{\Omega}(|\nabla a^{\frac{r}{2}}|^{2}+|\nabla
b^{\frac{r}{2}}|^{2}+|\nabla
z^{\frac{r}{2}}|^{2})+c_{8}(\varepsilon,r).\end{array}$ (3.25)

Accordingly upon the choice of $\varepsilon=\frac{c_{5}}{2c_{7}(r)}$, (3.25)
shows that

$\displaystyle\displaystyle\frac{d}{dt}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})+c_{5}\int_{\Omega}(e^{\chi_{u}v}a^{r}+e^{\chi_{w}v}b^{r}+z^{r})\leq
c_{9}(r)$

with some $c_{9}(r)>0$ for all $t>t_{0}$, and hence entails

$\int_{\Omega}a^{r}(\cdot,t)+\int_{\Omega}b^{r}(\cdot,t)+\int_{\Omega}z^{r}(\cdot,t)\leq
c_{10}(r)$

with some $c_{10}(r)>0$ by a standard ODE comparison argument. At this
position, one can derive a bound for $a,b,z$ with respect to the norm in
$L^{\infty}(\Omega)$ by means of a Moser-type iteration argument in quite a
standard manner. We omit the proof thereof, and would like refer to [15, 26,
28] for details in a closely related setting. ∎

##  4 Asymptotic behavior

On the basis of the exponential decay of quantities $w,z$ with respect to the
norm in $L^{1}(\Omega)$ and global boundedness of solutions, we will address
the large time asymptotics of the solution $(u,v,w,z)$ to (1.3). To this end,
we first turn the $L^{1}$-decay information explicitly contained in Lemma 2.3
to the decay property of $z$ in $L^{\infty}$-norm by an appropriate
application of the parabolic smoothing estimates in the two-dimensional
domain.

######  Lemma 4.1.

Let $(u,v,w,z)$ be the global solutions of (1.3) with $\beta\in(0,1)$. Then
for all $t>0$, we have

$\|z(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\gamma_{1}t}$ (4.1)

with $\gamma_{1}=\frac{\min\\{1-\beta,\delta_{z}\\}}{2}$ and some $C>0$.

###### Proof.

We invoke Lemma 2.3 along with (3.17) to see that there exist $c_{1}>0$ and
$c_{2}>0$ such that

$\|w(\cdot,t)\|_{L^{2}(\Omega)}^{2}\leq
c_{1}\|w(\cdot,t)\|_{L^{\infty}(\Omega)}e^{-\delta t}\leq c_{2}e^{-\delta t}$
(4.2)

with $\delta=\min\\{1-\beta,\delta_{z}\\}>0$. According to known smoothing
properties of the the Neumann heat semigroup
$\left(e^{\sigma\Delta}\right)_{\sigma>0}$ on the domain
$\Omega\subset\mathbb{R}^{2}$ ([31]), there exists $c_{3}>0$ such that for
each $\varphi\in C^{0}(\Omega)$,

$\left\|e^{\sigma D_{z}\Delta}\varphi\right\|_{L^{\infty}(\Omega)}\leq
c_{3}(1+\sigma^{-\frac{1}{2}})\|\varphi\|_{L^{2}(\Omega)}.$ (4.3)

Due to the nonnegativity of $z$ and the comparison principle, we may use (4.2)
and (4.3) to infer that

$\begin{array}[]{rl}z(\cdot,t)&=e^{t(D_{z}\Delta-\delta_{z})}z_{0}+\displaystyle\int^{t}_{0}e^{(t-s)(D_{z}\Delta-\delta_{z})}(\beta
w-\rho uz)(\cdot,s)ds\\\ &\leq
e^{t(D_{z}\Delta-\delta_{z})}z_{0}+\beta\displaystyle\int^{t}_{0}e^{(t-s)(D_{z}\Delta-\delta_{z})}w(\cdot,s)ds\\\
&\leq\displaystyle e^{-\delta_{z}t}\|z_{0}\|_{L^{\infty}(\Omega)}+\beta
c_{3}\int^{t}_{0}(1+(t-s)^{-\frac{1}{2}})e^{-\delta_{z}(t-s)}\|w(\cdot,s)\|_{L^{2}(\Omega)}ds\\\
&\leq\displaystyle e^{-\delta_{z}t}\|z_{0}\|_{L^{\infty}(\Omega)}+\beta
c_{2}c_{3}\int^{t}_{0}(1+(t-s)^{-\frac{1}{2}})e^{-\delta_{z}(t-s)}e^{-\frac{\delta
s}{2}}ds\\\ &\leq\displaystyle
e^{-\delta_{z}t}\|z_{0}\|_{L^{\infty}(\Omega)}+\beta
c_{2}c_{3}c_{4}e^{-min\\{\delta_{z},\frac{\delta}{2})t}\end{array}$ (4.4)

with some $c_{4}>0$, which along with the nonnegativity of $z$ entails that
(4.1) holds with $C=\|z_{0}\|_{L^{\infty}(\Omega)}+\beta c_{2}c_{3}c_{4}.$ ∎

Now thanks to the uniform decay property of $z$, a pointwise lower bound for
$a=ue^{-\chi_{u}v}$ can be achieved by means of an argument based on
comparison with spatially flat functions, which is documented as follows.

######  Lemma 4.2.

Under the assumption of Lemma 4.1, one can find constant $\gamma>0$ such that

$a(x,t)>\gamma~{}~{}\hbox{for all $(x,t)\in\Omega\times(0,\infty)$}.$ (4.5)

###### Proof.

According to Lemma 4.1, one can pick $t_{0}>0$ sufficiently large such that
for all $t>t_{0}$

$\|z(\cdot,t)\|_{L^{\infty}(\Omega)}\leq\frac{\mu_{u}}{2\rho}.$ (4.6)

Hence by means of a straightforward computation based on (2.2), one can see
that

$\displaystyle a_{t}$ $\displaystyle\geq$ $\displaystyle
D_{u}e^{-\chi_{u}v}\nabla\cdot(e^{\chi_{u}v}\nabla
a)+\mu_{u}a(1-ae^{\chi_{u}v})-\displaystyle\rho az$ $\displaystyle\geq$
$\displaystyle D_{u}e^{-\chi_{u}v}\nabla\cdot(e^{\chi_{u}v}\nabla
a)+a\left(\frac{\mu_{u}}{2}-e^{\chi_{u}\|v_{0}\|_{L^{\infty}(\Omega)}}a\right)$

for all $t>t_{0}$.

Now let $\underline{a}(t)$ be the smooth solution to the initial value
problem:

$\left\\{\begin{array}[]{lll}\underline{a}_{t}=\underline{a}\left(\frac{\mu_{u}}{2}-e^{\chi_{u}\|v_{0}\|_{L^{\infty}(\Omega)}}\underline{a}\right),\\\
\underline{a}(t_{0})=\inf\limits_{x\in{\Omega}}\\{u(x,t_{0})e^{-\chi_{u}\|v_{0}\|_{L^{\infty}(\Omega)}}\\},\end{array}\right.$
(4.7)

then through the explicit solution of above Bernoulli-type ODE, we have

$\underline{a}(t)\geq
c_{1}\coloneqq\min\\{\underline{a}(t_{0}),\frac{\mu_{u}}{2}e^{-\chi_{u}\|v_{0}\|_{L^{\infty}(\Omega)}}\\}>0$
(4.8)

for all $t>t_{0}$. It is observed that

$\displaystyle\underline{a}_{t}=D_{u}e^{-\chi_{u}v}\nabla\cdot(e^{\chi_{u}v}\nabla\underline{a})+\underline{a}\left(\frac{\mu_{u}}{2}-e^{\chi_{u}\|v_{0}\|_{L^{\infty}(\Omega)}}\underline{a}\right)$

and $a(x,t_{0})\geq\underline{a}(t_{0})$. Hence from the comparison principle
of the parabolic equation, one can conclude that

$a(x,t)\geq\underline{a}(t)\geq c_{1}\qquad\hbox{for all
$(x,t)\in\Omega\times(t_{0},\infty)$}.$ (4.9)

On the other hand, for any $t\in[0,t_{0}]$, the continuity of $a$ allows us to
find constant
$c_{2}\coloneqq\inf\limits_{(x,t)\in\overline{\Omega}\times[0,t_{0}]}a(x,t)$
such that

$a(x,t)\geq c_{2}\qquad\hbox{for all $(x,t)\in\Omega\times(0,t_{0}]$}.$ (4.10)

and thereby (4.5) results from (4.9) and (4.10) with
$\gamma\coloneqq\min\\{c_{1},c_{2}\\}$. ∎

In view of the $v$-equation in (2.2), the latter information immediately
entails the exponential decay of $v$ with respect to $L^{\infty}(\Omega)$
norm.

######  Lemma 4.3.

Under the assumption of Lemma 4.1, we have for all $t>0$

$\|v(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\gamma_{2}t}$ (4.11)

with $\gamma_{2}=\alpha_{u}\gamma$ and some $C>0$.

###### Proof.

By recalling the outcomes of Lemma 4.2, we have

$v_{t}=-(\alpha_{u}u+\alpha_{w}w)v\leq-\alpha_{u}\gamma v$

and hence

$v(x,t)\leq v_{0}(x)e^{-\alpha_{u}\gamma t}$

from which (4.11) follows. ∎

Furthermore upon the decay property of $v$ with respect to
$L^{\infty}(\Omega)$, one can derive the following basic stabilization feature
of $a(=e^{-\chi_{u}v}u)$.

######  Lemma 4.4.

Let the assumption of Lemma 4.1 hold. Then we have

$\int^{\infty}_{0}\int_{\Omega}\frac{|\nabla a|^{2}}{a^{2}}<\infty$ (4.12)

as well as

$\int^{\infty}_{0}\int_{\Omega}(u-1)^{2}<\infty.$ (4.13)

###### Proof.

In view of $s-1-\log s>0$ for all $s>0$ and $v_{t}<0$, we can conclude that

$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}(a-1-\log
a)=&\displaystyle\int_{\Omega}e^{\chi_{u}v}(a-1-\log
a)v_{t}+\int_{\Omega}e^{\chi_{u}v}\left(\frac{a-1}{a}\right)a_{t}\\\\[5.69054pt]
\leq&\displaystyle-D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a^{2}}+\mu_{u}\int_{\Omega}e^{\chi_{u}v}(a-1)(1-u)\\\
&+\displaystyle\chi_{u}\int_{\Omega}e^{\chi_{u}v}(a-1)(\alpha_{u}u+\alpha_{w}w)v-\displaystyle\rho\int_{\Omega}e^{\chi_{u}v}z(a-1).\end{array}$
(4.14)

Here by Young’s inequality,

$\begin{array}[]{rl}(1-a)(1-u)&=(1-u)^{2}+(u-a)(1-u)\\\
&\geq\displaystyle\frac{1}{2}(1-u)^{2}-a^{2}(e^{\chi_{u}v}-1)^{2}.\end{array}$
(4.15)

Due to the fact that $e^{s}\leq 1+2s$ for all $s\in[0,log2]$, (4.11) allows us
to fix a $t_{1}>1$ suitable large such that for all $t\geq t_{1}$,

$(e^{\chi_{u}v(\cdot,t)}-1)^{2}\leq 4\chi_{u}^{2}v^{2}(\cdot,t),$

which together with (4.15) entails that for $t\geq t_{1}$

$(1-a)(1-u)\geq\displaystyle\frac{1}{2}(1-u)^{2}-4a^{2}\chi_{u}^{2}v^{2}.$

Therefore we infer from (3.17) and (4.14) that for all $t\geq t_{1}$

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{u}v}(a-1-\log
a)+D_{u}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a^{2}}+\mu_{u}\int_{\Omega}(u-1)^{2}\\\\[5.69054pt] \leq&\displaystyle
c_{1}\int_{\Omega}z+c_{1}\int_{\Omega}v\end{array}$

with some $c_{1}>0$. After a time integration this leads to

$\begin{array}[]{rl}&\displaystyle
D_{u}\int^{t}_{t_{1}}\int_{\Omega}e^{\chi_{u}v}\frac{|\nabla
a|^{2}}{a^{2}}+\mu_{u}\int^{t}_{t_{1}}\int_{\Omega}(u-1)^{2}\\\\[5.69054pt]
\leq&c_{1}\displaystyle\int^{t}_{t_{1}}\int_{\Omega}z+c_{1}\int^{t}_{t_{1}}\int_{\Omega}v+\int_{\Omega}e^{\chi_{u}m_{v}}(a(\cdot,t_{1})-1-\log
a(\cdot,t_{1}))\end{array}$

and thereby implies that both (4.12) and (4.13) is valid thanks to (4.1) and
(4.11). ∎

In order to improve yet quite weak decay information of $u$, we turn to
consider the exponential decay properties of $\int_{\Omega}|\nabla
v(\cdot,t)|^{2}$, rather than the integrability of $a_{t}$ in
$L^{2}((0,\infty);L^{2}(\Omega))$. As the first step toward this, we first
show the convergence of integral $\int^{\infty}_{0}\int_{\Omega}|\nabla
v|^{2}$, which is stated below.

######  Lemma 4.5.

Let the assumption of Lemma 4.1 hold. Then we have

$\int^{\infty}_{0}\int_{\Omega}|\nabla v|^{2}<\infty.$ (4.16)

###### Proof.

Multiplying the second equation in (2.2) by $e^{\chi_{w}v}b$ and integrating
by parts, one can conclude that

$\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b^{2}+2D_{w}\int_{\Omega}e^{\chi_{w}v}|\nabla
b|^{2}+2\int_{\Omega}e^{\chi_{w}v}b^{2}$
$\displaystyle=2\rho\int_{\Omega}abze^{\chi_{w}v}+2\chi_{w}\int_{\Omega}e^{\chi_{w}v}b^{2}(\alpha_{u}ae^{\chi_{u}v}+\alpha_{w}be^{\chi_{w}v})v,$

which together with the global-in-time boundedness property of $a$ and $b$,
implies that

$\displaystyle\frac{d}{dt}\int_{\Omega}e^{\chi_{w}v}b^{2}+2D_{w}\int_{\Omega}e^{\chi_{w}v}|\nabla
b|^{2}+2\int_{\Omega}e^{\chi_{w}v}b^{2}\leq
c_{1}\left(\int_{\Omega}v+\int_{\Omega}z\right)$

for some $c_{1}>0$. Hence according to Lemma 4.3 and Lemma 2.3, we can get

$\int^{\infty}_{0}\int_{\Omega}|\nabla b|^{2}<\infty.$ (4.17)

Now since

$\nabla v_{t}=-(\alpha_{u}\nabla u+\alpha_{w}\nabla
w)v-(\alpha_{u}u+\alpha_{w}w)\nabla v,$

a direct computation shows that

$\begin{array}[]{rl}&\displaystyle\frac{1}{2}\frac{d}{dt}\int_{\Omega}|\nabla
v|^{2}+\int_{\Omega}(u+w)|\nabla v|^{2}\\\
=&-\alpha_{w}\chi_{w}\displaystyle\int_{\Omega}ve^{\chi_{w}v}b|\nabla
v|^{2}-\alpha_{w}\displaystyle\int_{\Omega}ve^{\chi_{w}v}\nabla v\cdot\nabla
b-\alpha_{u}\displaystyle\int_{\Omega}v\nabla v\cdot\nabla u\\\
\leq&-\alpha_{w}\displaystyle\int_{\Omega}ve^{\chi_{w}v}\nabla v\cdot\nabla
b-\alpha_{u}\displaystyle\int_{\Omega}ve^{\chi_{u}v}\nabla v\cdot\nabla
a.\end{array}$

Therefore, recalling the pointwise lower bound in (4.5), (3.17) and by the
Young inequality, we can find find a constant $c_{2}>0$ such that

$\frac{d}{dt}\int_{\Omega}|\nabla v|^{2}+\gamma\int_{\Omega}|\nabla v|^{2}\leq
c_{2}\left(\int_{\Omega}|\nabla b|^{2}+\int_{\Omega}\frac{|\nabla
a|^{2}}{a^{2}}\right).$ (4.18)

Hence combining this with (4.17) and (4.12), (4.16) is actually valid. ∎

Beyond the integrability of $\int_{\Omega}|\nabla v|^{2}$ over $(0,\infty)$,
we make use of the explicit expression of $\nabla v$ together with (4.17) and
(4.12) to identify that $\int_{\Omega}|\nabla v|^{2}$ exponentially decays.

######  Lemma 4.6.

Under the assumption of Lemma 4.1, one can find constant $C>0$ such that

$\int_{\Omega}|\nabla v|^{2}\leq C(t+1)e^{-2\gamma t}~{}~{}\hbox{for all
$t>0$,}$ (4.19)

where $\gamma$ is given by Lemma 4.2.

###### Proof.

On the basis of the $v$-equation in (1.3), we have

$\nabla v(\cdot,t)=\nabla
v(\cdot,0)e^{-\int^{t}_{0}(u+w)(\cdot,s)ds}-v(\cdot,0)e^{-\int^{t}_{0}(u+w)(\cdot,s)ds}\int^{t}_{0}(\nabla
u(\cdot,s)+\nabla w(\cdot,s))ds$

which along with (4.5) and the Young inequality entails that

$\int_{\Omega}|\nabla v|^{2}\leq 2e^{-2\gamma t}\|\nabla
v_{0}\|^{2}_{L^{2}(\Omega)}+4te^{-2\gamma
t}\|v_{0}\|^{2}_{L^{\infty}(\Omega)}(\int^{t}_{0}\int_{\Omega}|\nabla
u|^{2}ds+\int^{t}_{0}\int_{\Omega}|\nabla w|^{2}ds).$ (4.20)

Furthermore observing that

$|\nabla w|\leq\chi_{w}e^{\chi_{w}v}|\nabla v|b+e^{\chi_{w}v}|\nabla b|$

as well as

$|\nabla u|\leq\chi_{u}e^{\chi_{u}v}|\nabla v|a+e^{\chi_{u}v}|\nabla a|,$

we conclude from (4.20) that there exists $c_{1}>0$ such that

$\begin{array}[]{rl}\displaystyle\int_{\Omega}|\nabla v(\cdot,t)|^{2}&\leq
c_{1}e^{-2\gamma t}+c_{1}te^{-2\gamma
t}\displaystyle\int^{t}_{0}\int_{\Omega}(|\nabla b|^{2}+|\nabla a|^{2}+|\nabla
v|^{2})ds\\\\[8.53581pt] &\leq c_{1}e^{-2\gamma t}+c_{1}te^{-2\gamma
t}\displaystyle\int^{\infty}_{0}\int_{\Omega}(|\nabla b|^{2}+|\nabla
a|^{2}+|\nabla v|^{2})ds\end{array}$ (4.21)

and thus

$\displaystyle\int_{\Omega}|\nabla v(\cdot,t)|^{2}\leq c_{2}(t+1)e^{-2\gamma
t}$

with some $c_{2}>0$, thanks to (4.17), (4.12) and (4.16). ∎

On the basis of smoothing estimates for the Neumann heat semigroup on
$\Omega$, and decay information provided by Lemma 4.6, we can make sure that
$u-1$ decays exponentially with respect to $L^{p}(\Omega)$-norm.

######  Lemma 4.7.

Suppose the condition in Lemma 4.1 holds, then there exists $\eta_{1}>0$ such
that for every $p\geq 2$,

$\|u(\cdot,t)-1\|_{L^{p}(\Omega)}\leq C(p)e^{-\eta_{1}t}$ (4.22)

with some $C(p)>0$ for all $t>0$.

###### Proof.

Testing the first equation in (1.3) by $u-1$ and integrating by parts, we have

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}(u-1)^{2}+2D_{u}\int_{\Omega}|\nabla
u|^{2}+2\mu_{u}\int_{\Omega}u(u-1)^{2}\\\\[5.69054pt]
=&2\xi_{u}\displaystyle\int_{\Omega}u\nabla v\cdot\nabla
u-2\rho\int_{\Omega}(u-1)uz.\end{array}$

We thereupon make use of Lemma 4.2, Lemma 3.5 along with the Young inequality
to get

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}(u-1)^{2}+D_{u}\int_{\Omega}|\nabla
u|^{2}+2\mu_{u}\gamma\int_{\Omega}(u-1)^{2}\\\
\leq&\displaystyle\frac{\xi^{2}_{u}}{D_{u}}\int_{\Omega}u^{2}|\nabla
v|^{2}+2\rho\int_{\Omega}uz\\\ \leq&c_{1}\displaystyle\int_{\Omega}|\nabla
v|^{2}+c_{1}\int_{\Omega}z\end{array}$ (4.23)

with some $c_{1}>0$.

According to Lemma 4.6 and Lemma 2.3, (4.23) implies that

$\int_{\Omega}(u-1)^{2}\leq c_{2}e^{-\eta_{1}t}$ (4.24)

with $\eta_{1}:=\min\\{2\mu_{u}\gamma,2\gamma,\delta\\}$ and $c_{2}>0$ for all
$t>0$.

Recalling known smoothing estimates for the Neumann heat semigroup on
$\Omega\subset\mathbb{R}^{2}$ ([31]), there exist $c_{3}=c_{3}(p,q)>0$,
$c_{4}=c_{4}(p,q)>0$ fulfilling

$\left\|e^{\sigma D_{u}\Delta}\varphi\right\|_{L^{p}(\Omega)}\leq
c_{3}\sigma^{-(\frac{1}{q}-\frac{1}{p})}\|\varphi\|_{L^{q}(\Omega)}$ (4.25)

for each $\varphi\in C^{0}(\Omega)$, and for all
$\varphi\in\left(L^{q}\left(\Omega\right)\right)^{2}$,

$\left\|e^{\sigma
D_{u}\Delta}\nabla\cdot\varphi\right\|_{L^{p}\left(\Omega\right)}\leq
c_{4}(1+\sigma^{-\frac{1}{2}-(\frac{1}{q}-\frac{1}{p})})e^{-\lambda_{1}\sigma}\left\|\varphi\right\|_{L^{q}\left(\Omega\right)}$
(4.26)

with $\lambda_{1}>0$ the first nonzero eigenvalue of $-\Delta$ in $\Omega$
under the Neumann boundary condition.

Relying on a variation-of-constants representation of $u$ related to the the
first equation in (1.3), we utilize (4.25) and (4.26) to infer that

$\begin{array}[]{rl}&\|(u-1)(\cdot,t)\|_{L^{p}(\Omega)}\\\
\leq&\|e^{t(D_{u}\Delta-\delta)}(u_{0}-1)\|_{L^{p}(\Omega)}+\xi_{u}\displaystyle\int^{t}_{0}\|e^{(t-s)(D_{u}\Delta-\delta)}\nabla\cdot(u\nabla
v)\|_{L^{p}(\Omega)}ds\\\
&+\displaystyle\int^{t}_{0}\|e^{(t-s)(D_{u}\Delta-\delta)}((\mu_{u}u-\delta)(1-u)-\rho
uz)\|_{L^{p}(\Omega)}ds\\\ \leq&\displaystyle e^{-\delta
t}\|u_{0}-1\|_{L^{p}(\Omega)}+c_{5}(p)\int^{t}_{0}(1+(t-s)^{-1+\frac{1}{p}})e^{-(\delta+\lambda_{1})(t-s)}\|\nabla
v(\cdot,s)\|_{L^{2}(\Omega)}ds\\\
&+c_{5}(p)\displaystyle\int^{t}_{0}(1+(t-s)^{-\frac{1}{2}+\frac{1}{p}})e^{-\delta(t-s)}\|(u-1)(\cdot,s)\|_{L^{2}(\Omega)}ds\\\
&+c_{5}(p)\displaystyle\int^{t}_{0}(1+(t-s)^{-1+\frac{1}{p}})e^{-\delta(t-s)}\|z(\cdot,s)\|_{L^{1}(\Omega)}ds\end{array}$
(4.27)

for some $c_{5}(p)>0$. Therefore by (4.24), (4.19), (2.8) and thanks to the
fact that for $\alpha\in(0,1)$ $\gamma_{1}$ and $\delta_{1}$ positive
constants with $\gamma_{1}\neq\delta_{1}$, there exists $c_{6}>0$ such that

$\int_{0}^{t}(1+(t-s)^{-\alpha})e^{-\gamma_{1}s}e^{-\delta_{1}(t-s)}ds\leq
c_{6}e^{-min\left\\{\gamma_{1},\delta_{1}\right\\}t},$

(4.22) readily results from (4.27) with
$\eta_{1}=\frac{1}{2}\min\\{\gamma,\mu_{u}\gamma,\frac{1-\beta}{2},\frac{\delta_{z}}{2}\\}$
and some $C(p)>0$. ∎

At this position, due to the fact that the integrability exponent in (4.19)
does not exceed the considered spatial dimension $n=2$, the uniform decay of
$w$ is not achieved herein, however a somewhat optimal decay rate thereof with
respect to $L^{p}(\Omega)$ may be derived by the argument similar to that in
Lemma 4.7 instead of the simple interpolation. The desired result can be
stated below and the corresponding proof is omitted herein.

######  Lemma 4.8.

Let the condition in Lemma 4.1 hold. Then there exists $\varrho_{1}>0$ such
that for every $p\geq 2$,

$\|w(\cdot,t)\|_{L^{p}(\Omega)}\leq C(p)e^{-\varrho_{1}t}$ (4.28)

with some $C(p)>0$ for all $t>0$.

Next we proceed to establish the convergence properties in (1.6)–(1.7) stated
in Theorem 1.1, which are beyond that in Lemma 4.7 and Lemma 4.8. To this end,
thanks to Lemma 4.7, Lemma 4.8, Lemma 4.1, Lemma 4.2 and Lemma 4.3, we turn to
make sure that $\int_{\Omega}|\nabla v|^{4}$ decays exponentially, which
results from a series of testing procedures.

######  Lemma 4.9.

Let conditions in Theorem 1.1 hold. Then there exist $\eta_{2}>0$ and $C>0$
such that

$\int_{\Omega}|\nabla v|^{4}\leq Ce^{-\eta_{2}t}~{}~{}\hbox{for all $t>0$.}$
(4.29)

###### Proof.

Testing the identity

$a_{t}=D_{u}\triangle a+D_{u}\nabla v\cdot\nabla a+f(x,t),\quad
x\in\Omega,\quad t>0$

with $f(x,t)=\mu_{u}a(1-u)-\displaystyle\rho
az+\chi_{u}a(\alpha_{u}u+\alpha_{w}w)v$ by $-\triangle a$, and using Young’s
inequality, we get

$\begin{array}[]{rl}\displaystyle\frac{d}{dt}\int_{\Omega}|\nabla
a|^{2}+2D_{u}\displaystyle\int_{\Omega}|\triangle
a|^{2}=&-2D_{u}\chi_{u}\displaystyle\int_{\Omega}(\nabla a\cdot\nabla
v)\triangle a-2\int_{\Omega}f\triangle a\\\
\leq&D_{u}\displaystyle\int_{\Omega}|\triangle
a|^{2}+2D_{u}\chi_{u}^{2}\int_{\Omega}|\nabla a|^{2}|\nabla
v|^{2}+\frac{2}{D_{u}}\int_{\Omega}|f|^{2}.\end{array}$ (4.30)

Note that by the Gagliardo–Nirenberg type interpolation with standard elliptic
regularity theory and Poincaré’s inequality, one can find constants $c_{1}>0$
and $c_{2}>0$ such that for all $\varphi\in W^{2,2}(\Omega)$ with
$\frac{\partial\varphi}{\partial\nu}=0$ on $\partial\Omega$,

$\|\nabla\varphi\|^{4}_{L^{4}(\Omega)}\leq
c_{1}\|\Delta\varphi\|^{2}_{L^{2}(\Omega)}\|\varphi\|^{2}_{L^{\infty}(\Omega)}$

and

$\|\nabla\varphi\|^{2}_{L^{2}(\Omega)}\leq
c_{2}\|\Delta\varphi\|^{2}_{L^{2}(\Omega)}$

(see Lemma A.1 and A.3 in [7]). Hence thanks to Lemma 3.5, we can pick
$c_{3}>0$ such that

$a(x,t)\leq c_{3},\quad b(x,t)\leq c_{3}$

for all $x\in\Omega$ and $t>0$, and thereby have

$\|\nabla a\|^{4}_{L^{4}(\Omega)}\leq c_{1}c_{3}^{2}\|\Delta
a\|^{2}_{L^{2}(\Omega)},\quad\|\nabla a\|^{2}_{L^{2}(\Omega)}\leq
c_{2}\|\Delta a\|^{2}_{L^{2}(\Omega)}$ (4.31)

as well as

$\|\nabla b\|^{4}_{L^{4}(\Omega)}\leq c_{1}c_{3}^{2}\|\Delta
b\|^{2}_{L^{2}(\Omega)},\quad\|\nabla b\|^{2}_{L^{2}(\Omega)}\leq
c_{2}\|\Delta b\|^{2}_{L^{2}(\Omega)}$ (4.32)

Combining (4.31) with (4.30), the Young inequality shows that for every
$\varepsilon>0$

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\frac{D_{u}}{2c_{2}}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+D_{u}\displaystyle\|\triangle
a\|^{2}_{L^{2}(\Omega)}\\\\[8.53581pt] \leq&\varepsilon\|\nabla
a\|^{4}_{L^{4}(\Omega)}+\displaystyle\frac{D_{u}}{2c_{2}}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{u}^{2}\chi_{u}^{4}}{\varepsilon}\|\nabla
v\|^{4}_{L^{4}(\Omega)}+\frac{2}{D_{u}}\|f\|_{L^{2}(\Omega)}^{2}\\\
\leq&\varepsilon c_{1}c_{3}^{2}\displaystyle\|\Delta
a\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{u}}{2}\|\triangle
a\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{u}^{2}\chi_{u}^{4}}{\varepsilon}\|\nabla
v\|^{4}_{L^{4}(\Omega)}+\frac{2}{D_{u}}\|f\|_{L^{2}(\Omega)}^{2},\end{array}$
(4.33)

which along with the choice of $\varepsilon=\frac{D_{u}}{4c_{1}c_{3}^{2}}$
implies that

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\frac{D_{u}}{2c_{2}}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{u}}{4}\|\triangle
a\|^{2}_{L^{2}(\Omega)}\\\\[8.53581pt] \leq&\displaystyle
4D_{u}\chi_{u}^{4}c_{1}c_{3}^{2}\|\nabla
v\|^{4}_{L^{4}(\Omega)}+\frac{2}{D_{u}}\|f\|_{L^{2}(\Omega)}^{2}.\end{array}$
(4.34)

Likely, we can get

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\|\nabla
b\|^{2}_{L^{2}(\Omega)}+\frac{D_{w}}{2c_{2}}\|\nabla
b\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{w}}{4}\|\triangle
b\|^{2}_{L^{2}(\Omega)}\\\\[8.53581pt] \leq&\displaystyle
4D_{w}\chi_{w}^{4}c_{1}c_{3}^{2}\|\nabla
v\|^{4}_{L^{4}(\Omega)}+\frac{2}{D_{w}}\|g\|_{L^{2}(\Omega)}^{2}.\end{array}$
(4.35)

with $g(x,t)=-b+\displaystyle\rho
ue^{-\chi_{w}v}z+\chi_{w}b(\alpha_{u}u+\alpha_{w}w)v$.

Now in order to appropriately compensate the first summand on right-hand side
of (4.34) and (4.35), we use the third equation in (2.2) to see that

$\begin{array}[]{rl}&\displaystyle\frac{1}{4}\frac{d}{dt}\int_{\Omega}|\nabla
v|^{4}\\\ =&-\displaystyle\int_{\Omega}|\nabla v|^{2}\nabla v\cdot\nabla
v_{t}\\\
=&-\alpha_{u}\displaystyle\int_{\Omega}a(v+\chi_{u})e^{\chi_{u}v}|\nabla
v|^{4}-\alpha_{u}\displaystyle\int_{\Omega}ve^{\chi_{u}v}|\nabla v|^{2}\nabla
v\cdot\nabla a\\\
&-\alpha_{w}\displaystyle\int_{\Omega}b(v+\chi_{w})e^{\chi_{w}v}|\nabla
v|^{4}-\alpha_{w}\displaystyle\int_{\Omega}ve^{\chi_{w}v}|\nabla v|^{2}\nabla
v\cdot\nabla b.\end{array}$ (4.36)

Here recalling the uniform positivity of $a$ stated in Lemma 4.2, we can pick
$c_{4}>0$ fulfilling

$\alpha_{u}\displaystyle\int_{\Omega}a(v+\chi_{u})e^{\chi_{u}v}|\nabla
v|^{4}\geq\alpha_{u}\chi_{u}c_{4}\displaystyle\int_{\Omega}|\nabla v|^{4}$

and thus infer by the Young inequality and Lemma 2.2 that for all $t>t_{0}$

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\int_{\Omega}|\nabla
v|^{4}+c_{5}\int_{\Omega}|\nabla v|^{4}\\\
\leq&c_{6}\|v(\cdot,t_{0})\|_{L^{\infty}(\Omega)}\displaystyle\int_{\Omega}(|\nabla
a|^{4}+|\nabla b|^{4})\end{array}$ (4.37)

with constants $c_{5}>0$, $c_{6}>0$.

Now if we write
$d_{1}:=\frac{8c_{1}c_{3}^{2}(D_{u}\chi_{u}^{4}+D_{w}\chi_{w}^{4})}{c_{5}}$,
combining (4.37), (4.34) with (4.35) yields

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\left(\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\|\nabla b\|^{2}_{L^{2}(\Omega)}+d_{1}\|\nabla
v\|_{L^{4}(\Omega)}^{4}\right)\\\\[8.53581pt]
&+\displaystyle\frac{D_{u}}{2c_{2}}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{u}}{4}\|\triangle
a\|^{2}_{L^{2}(\Omega)}+\frac{D_{w}}{2c_{2}}\|\nabla
b\|^{2}_{L^{2}(\Omega)}+\displaystyle\frac{D_{w}}{4}\|\triangle
b\|^{2}_{L^{2}(\Omega)}+\frac{c_{5}d_{1}}{2}\|\nabla
v\|_{L^{4}(\Omega)}^{4}\\\\[8.53581pt]
\leq&c_{6}d_{1}\|v(\cdot,t_{0})\|_{L^{\infty}(\Omega)}(\|\nabla
a\|_{L^{4}(\Omega)}^{4}+\|\nabla
b\|_{L^{4}(\Omega)}^{4})+\displaystyle\frac{2}{D_{u}}\|f\|_{L^{2}(\Omega)}^{2}+\frac{2}{D_{w}}\|g\|_{L^{2}(\Omega)}^{2},\end{array}$
(4.38)

which together with (4.31), (4.32) and Lemma 4.3 entails that there exists
$t_{1}>1$ suitably large such that for all $t>t_{1}$,

$\begin{array}[]{rl}&\displaystyle\frac{d}{dt}\left(\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\|\nabla b\|^{2}_{L^{2}(\Omega)}+d_{1}\|\nabla
v\|_{L^{4}(\Omega)}^{4}\right)\\\\[8.53581pt]
&+\displaystyle\frac{D_{u}}{2c_{2}}\|\nabla
a\|^{2}_{L^{2}(\Omega)}+\frac{D_{w}}{2c_{2}}\|\nabla
b\|^{2}_{L^{2}(\Omega)}+\frac{c_{5}d_{1}}{2}\|\nabla
v\|_{L^{4}(\Omega)}^{4}\\\\[8.53581pt]
\leq&\displaystyle\frac{2}{D_{u}}\|f\|_{L^{2}(\Omega)}^{2}+\frac{2}{D_{w}}\|g\|_{L^{2}(\Omega)}^{2}.\end{array}$
(4.39)

Due to Lemma 3.5, we see that

$|f(x,t)|^{2}+|g(x,t)|^{2}\leq
c_{7}(|u(x,t)-1|^{2}+|w(x,t)|^{2}+|z(x,t)|^{2}+|v(x,t)|^{2})$

with some $c_{7}>0$, and thus there exist $\eta_{3}>0$ and $c_{8}>0$ such that

$\int_{\Omega}|f(x,t)|^{2}+|g(x,t)|^{2}\leq c_{8}e^{-\eta_{3}t}~{}~{}\hbox{for
all $t>t_{1},$}$ (4.40)

thanks to Lemma 4.7, Lemma 4.8, Lemma 4.1 and Lemma 4.3. Therefore (4.29)
readily results from (4.40) and (4.39). ∎

At this position, as an application of known smoothing estimates for the
Neumann heat semigroup, the latter readily turns to the exponential decay
property of $u-1$ as well as $w$ with respect to $L^{\infty}(\Omega)$-norm.

######  Lemma 4.10.

Let the conditions in Theorem 1.1 hold. Then there exist $\eta,\varrho>0$ and
$C>0$ fulfilling

$\|u(\cdot,t)-1\|_{L^{\infty}(\Omega)}\leq Ce^{-\eta t}$ (4.41)

as well as

$\|w(\cdot,t)\|_{L^{\infty}(\Omega)}\leq Ce^{-\varrho t}$ (4.42)

for all $t>0$.

###### Proof.

Since the proof is similar to that of Lemma 4.7, we only give a short proof of
(4.41). In view to known smoothing estimates for the Neumann heat semigroup on
$\Omega\subset\mathbb{R}^{2}$ ([31]), there exist $c_{1}>0$, $c_{2}>0$
fulfilling

$\left\|e^{\sigma D_{u}\Delta}\varphi\right\|_{L^{\infty}(\Omega)}\leq
c_{1}\sigma^{-\frac{1}{2}}\|\varphi\|_{L^{2}(\Omega)}$ (4.43)

for each $\varphi\in C^{0}(\Omega)$, and for all
$\varphi\in\left(L^{4}\left(\Omega\right)\right)^{2}$,

$\left\|e^{\sigma
D_{u}\Delta}\nabla\cdot\varphi\right\|_{L^{\infty}\left(\Omega\right)}\leq
c_{2}(1+\sigma^{-\frac{3}{4}})e^{-\lambda_{1}\sigma}\left\|\varphi\right\|_{L^{4}\left(\Omega\right)}$
(4.44)

with $\lambda_{1}>0$ the first nonzero eigenvalue of $-\Delta$ in $\Omega$
under the Neumann boundary condition.

According to the variation-of-constants representation of $u$ related to the
the first equation in (1.3), we utilize (4.43) and (4.44) to infer that

$\begin{array}[]{rl}&\|(u-1)(\cdot,t)\|_{L^{\infty}(\Omega)}\\\
\leq&\|e^{t(D_{u}\Delta-1)}(u_{0}-1)\|_{L^{\infty}(\Omega)}+\xi_{u}\displaystyle\int^{t}_{0}\|e^{(t-s)(D_{u}\Delta-1)}\nabla\cdot(u\nabla
v)\|_{L^{\infty}(\Omega)}ds\\\
&+\displaystyle\int^{t}_{0}\|e^{(t-s)(D_{u}\Delta-1)}((\mu_{u}u-1)(1-u)-\rho
uz)\|_{L^{\infty}(\Omega)}ds\\\ \leq&\displaystyle
e^{-t}\|u_{0}-1\|_{L^{\infty}(\Omega)}+c_{3}\int^{t}_{0}(1+(t-s)^{-\frac{3}{4}})e^{-(1+\lambda_{1})(t-s)}\|\nabla
v(\cdot,s)\|_{L^{4}(\Omega)}ds\\\
&+c_{3}\displaystyle\int^{t}_{0}(1+(t-s)^{-\frac{1}{2}})e^{-(t-s)}(\|(u-1)(\cdot,s)\|_{L^{2}(\Omega)}+\|z(\cdot,s)\|_{L^{2}(\Omega)})ds\end{array}$
(4.45)

with some $c_{3}>0$. This readily establish (4.41) with appropriate $\eta>0$
in view of (4.29), (4.22) and (4.1). ∎

Thereby our main result has essentially been proved already.

Proof of Theorem 1.1. The statement on global boundedness of classical
solutions has been asserted by Lemma 3.5. The convergence properties in
(1.6)–(1.9) are precisely established by Lemma 4.10, Lemma 4.1 and Lemma 4.3,
respectively.

##  5 Acknowledgments

This work is partially supported by NSFC (No.12071030).

## References

  * [1] T. Alzahrani, R. Eftimie, D. Trucu,  _Multiscale modelling of cancer response to oncolytic viral therapy_ , Math. Bioci. 310(2019), 76–95.
  * [2] A. R. Anderson, M. A. J. Chaplain, E. L. Newman, R. J. C. Steele, A. M. Thompson,  _Mathematical modelling of tumour invasion and metastasis_ , J. Theor. Med. 2(2000), 129–154.
  * [3] X. Cao, _Boundedness in a three-dimensional chemotaxis–haptotaxis model_ , Z. Angew. Math. Phys. 67(2016), 67:11.
  * [4] Z. Chen,  _Dampening effect of logistic source in a two-dimensional haptotaxis system with nonlinear zero-order interaction_ , J. Math. Anal. Appl. 492(2020), 124435.
  * [5] M.A.J. Chaplain, G. Lolas,  _Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system_ , Math. Mod. Meth. Appl. Sci. 18(2005), 1685–1734.
  * [6] M.A. Fontelos, A. Friedman, B. Hu, _Mathematical analysis of a model for the initiation of angiogenesis_ , SIAM J. Math. Anal. 33(2002), 1330–1355.
  * [7] M. Fuest, _Global solutions near homogeneous steady states in a multi-dimensional population model with both predator-and prey-taxis_ , SIAM J. Math. Anal. 52(2020), 5863–5891.
  * [8] H. Fukuhara, Y. Ino, T. Todo, _Oncolytic virus therapy: A new era of cancer treatment at dawn_ , Cancer Sci. 107(2016), 1373–1379.
  * [9] S. Gujar, J. G. Pol, Y. Kim, P. W. Lee, G. Kroemer, _Antitumor benefits of antiviral immunity: An underappreciated aspect of oncolytic virotherapies_ , Trends Immunol. 39(2018), 209–221.
  * [10] I. Ganly, D. Kirn, et al.,  _A phase I study of Onyx-015, an E1B-attenuated adenovirus, administered intratumorally to patients with with recurrent head and neck cancer_ , Clinical Cancer Res. 6(2000), 798–806.
  * [11] C. Jin,  _Global classical solutions and convergence to a mathematical model for cancer cells invasion and metastatic spread_ , J. Diff. Equations 269(2020), 3987–4021.
  * [12] H.Y. Jin, T. Xiang,  _Negligibility of haptotaxis effect in a chemotaxis–haptotaxis model_ , Math. Mod. Meth. Appl. Sci. 31(7)(2021), 1373–1417.
  * [13] N. L. Komarova,  _Viral reproductive strategies: how can lytic viruses be evolutionarily competitive?_ , J. Theor. Biol. 249(2007), 766–784.
  * [14] Y. Li, J. Lankeit, _Boundedness in a chemotaxis–haptotaxis model with nonlinear diffusion_ , Nonlinearity 29(2016), 1564–1595.
  * [15] J. Li, Y. Wang, _Boundedness in a haptotactic cross–diffusion system modeling oncolytic virotherapy_. J. Diff. Equation 270(2021), 94–113.
  * [16] J. Nemunaitis, I. Ganly, et al.,  _Selective replication and oncolysis in p53 mutant tumors with ONYX-015, an E1B-55kD gene-deleted adenovirus, in patients with advanced head and neck cancer: a phase II trial_ , Cancer Res. 60(2000), 6359–6366.
  * [17] P. Y. H. Pang, Y. Wang,  _Global boundedness of solutions to a chemotaxis–haptotaxis model with tissue remodeling_ , Math. Models Methods Appl. Sci. 28(2018), 2211–2235.
  * [18] P. Y. H. Pang, Y. Wang,  _Asymptotic behavior of solutions to a tumor angiogenesis model with chemotaxis–haptotaxis_ , Math. Models Methods Appl. Sci. 29(2019), 1387–1412.
  * [19] J. Prüss, R. Zacher, R. Schnaubelt,  _Global asymptotic stability of equilibria in models for virus dynamics_ , Math. Model. Nat. Phenom. 3(2008), 126–142.
  * [20] G. Ren, B. Liu,  _Global classical solvability in a three-dimensional haptotaxis system modelong oncolytic virotherapy_ , Math. Methods Appl. Sci. 44(2021), 9275–9291.
  * [21] C. Stinner, C. Surulescu, M. Winkler,  _Global weak solutions in a PDE-ODE system modeling multiscale cancer cell invasion_ , SIAM J. Math. Anal. 46(2014), 1969–2007.
  * [22] Y. Tao, M. Winkler,  _Global classical solutions to a doubly haptotactic cross–diffusion system modeling oncolytic virotherapy_ , J. Diff. Equations 268(2020), 4973–4997.
  * [23] Y.Tao, M. Winkler,  _Asymptotic stability of spatial homogeneity in a haptotaxis model for oncolytic virotherapy_ , Proc. Roy. Soc. Edinburgh Sect. A (2021), 1–21. DOI: https://doi.org/10.1017/prm.2020.97
  * [24] Y.Tao, M. Winkler,  _Critical mass for infinite-time blow-up in a haptotaxis system with nonlinear zero-order interaction_ , Discrete Contin. Dyn. Syst. A 41(2021), 439–454.
  * [25] Y.Tao, M.Winkler,  _A critical virus production rate for efficiency of oncolytic virotherapy_ , European J. Appl. Math. 32(2021), 301–316.
  * [26] Y.Tao, M. Winkler,  _A critical virus production rate for blow-up suppression in a haptotaxis model for oncolytic virotherapy_ , Nonlinear Anal. 198(2020), 111870.
  * [27] Y. Tao, M. Winkler,  _Large time behavior in a multidimensional chemotaxis–haptotaxis model with slow signal diffusion_ , SIAM J. Math. Anal. 47(2015), 4229–4250.
  * [28] Y. Tao, M. Winkler,  _Energy-type estimates and global solvability in a two-dimensional chemotaxis-haptotaxis model with remodeling of non-diffusible attractant_ , J. Diff. Equations 257(2014), 784–815.
  * [29] C. Walker, G.F. Webb, _Global existence of classical solutions for a haptotaxis model_ , SIAM J. Math. Anal. 38(2007), 1694–1713.
  * [30] Y. Wang, _Boundedness in the higher-dimensional chemotaxis–haptotaxis model with nonlinear diffusion_ , J. Diff. Equations 260(2016), 1975–1989.
  * [31] M. Winkler,  _Aggregation vs. Global diffusive behavior in the higher-dimensional Keller–Segel model_ , J. Diff. Equations 12(2010), 2889–2905.
  * [32] J. Zheng, Y. Ke, _Large time behavior of solutions to a fully parabolic chemotaxis–haptotaxis model in N dimensions_ , J. Diff. Equations 266(2019), 1969–2018.
  * [33] A. Zhigun, C. Surulescu, A. Uatay, _Global existence for a degenerate haptotaxis model of cancer invasion_ , Z. Angew. Math. Phys. 6(2016), 67:146.