An analysis of a mathematical model describing the growth of a tumor treated with chemotherapy

TL;DR

This paper analyzes a mathematical model combining ODEs and PDEs to describe tumor and normal cell growth under localized chemotherapy, providing theoretical results and numerical simulations.

Contribution

It introduces a novel mixed ODE-PDE model for tumor growth with localized drug delivery, including existence, uniqueness, and simulation results.

Findings

Model exhibits different behaviors depending on parameters.

Existence and uniqueness of solutions are established.

Numerical simulations illustrate tumor response to therapy.

Abstract

We present a mathematical analysis of a mixed ODE-PDE model describing the spatial distribution and temporal evolution of tumor and normal cells within a tissue subject to the effects of a chemotherapeutic drug. The model assumes that the influx of chemotherapy is restricted to a limited region of the tissue, mimicking a blood vessel passing transversely. We provide results on the existence and uniqueness of the model solution and numerical simulations illustrating different model behaviors.

Figures

Tables

Table 1. Table 1. Set-up of different simulations an their biological outcomes. Each row indicates the numerical values used for the chemotherapeutic parameters α A subscript 𝛼 𝐴 \alpha_{A} (cytotoxicity), σ 𝜎 \sigma (diffusion coefficient), μ 𝜇 \mu (infusion rate), and the position of ω ⊂ Ω ⊂ ℝ 2 𝜔 Ω superscript ℝ 2 \omega\subset\Omega\subset\mathbb{R}^{2} . Simulation 1 was performed in a two-dimensional domain Ω = [ 0 , 1 ] × [ 0 , 1 ] Ω 0 1 0 1 \Omega=[0,1]\times[0,1] , while simulations 2-5 were performed in a one-dimensional domain Ω = [ 0 , 1 ] Ω 0 1 \Omega=[0,1] .

Simulation	Figure	Outcome	${\alpha }_A$	$\mu$	$\sigma$	$\omega$
1	1	tumor persistence	5	3	0.1	$[0.45,0.55]\times [0.45,0.55]$
2	2	tumor persistence	10	3	0.1	$[0,0.1]$
3	3	tumor extinction	10	6	0.1	$[0,0.1]$
4	4	tumor extinction	10	3	0.2	$[0,0.1]$
5	5	tumor extinction	20	3	0.1	$[0,0.1]$

Equations

$$\left\{\begin{array}[]{lcl}\displaystyle\frac{\partial N}{\partial t}=r_{N}-\mu_{N}N-\beta_{1}NA-\alpha_{N}\gamma_{N}DN,&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial A}{\partial t}=r_{A}A\left(1-\frac{A}{k_{A}}\right)-(\mu_{A}+\epsilon_{A})A-\alpha_{A}\gamma_{A}DA,&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial D}{\partial t}=\sigma\Delta D+\mu\chi_{\omega}-\gamma_{A}DA-\gamma_{N}DN-\tau D,&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial D}{\partial\eta}=0,&\textup{on}&\Gamma,\vspace{0.2cm}\\ \displaystyle N(\cdot,0)=N_{0}(\cdot),A(\cdot,0)=A_{0}(\cdot),D(\cdot,0)=D_{0}(\cdot),&\textup{in}&\Omega.\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}\displaystyle\frac{dN}{dt}=r_{N}-\mu_{N}N-\beta_{1}NA-\alpha_{N}\gamma_{N}DN,\vspace{0.2cm}\\ \displaystyle\frac{dA}{dt}=r_{A}A\left(1-\frac{A}{k_{A}}\right)-(\mu_{A}+\epsilon_{A})A-\beta_{3}NA-\alpha_{A}\gamma_{A}DA,\vspace{0.2cm}\\ \displaystyle\frac{dD}{dt}=\mu-\gamma_{A}DA-\gamma_{N}DN-\tau D,\end{array}\right.$$

$$\begin{gathered} W_p^k(\Omega )=\left\{f \\ \in L^p(\Omega ):D^{\alpha }f\in L^p(\Omega ),|\alpha |\le k\right\}; \end{gathered}$$

$$W_q^{2,1}(Q)$$

$$W$$

$$L^p(0,T;B)$$

$$N\le ||N_0|{|}_{L^{\infty }(Q)}+r_{N}T,A\le C_{\lambda }||A_0|{|}_{L^{\infty }(\Omega )}a.e.inQ$$

$$||N|{|}_W+||A|{|}_W+||D|{|}_{W_4^{2,1}(Q)}\le C,$$

$$\Vert u{\Vert }_{L^q(Q)}\le C\Vert u{\Vert }_{W_p^{2,1}(Q)},$$

$$\left\{\begin{array}[]{lcl}\displaystyle\frac{\partial u}{\partial t}-\sum\limits_{i,j=1}^{n}a_{ij}(x,t)\frac{\partial u^{2}}{\partial x_{i}x_{j}}+\sum\limits_{j=1}^{n}a_{i}(x,t)\frac{\partial u}{\partial x_{j}}+a(x,t)u=f&\textup{in}&Q,\\ \displaystyle\sum\limits_{i=1}^{n}b_{i}(x,t)\frac{\partial u}{\partial x_{i}}+b(x,t)u=0&\textup{on}&\Gamma,\\ \displaystyle u(\cdot,0)=u_{0}(\cdot)&\textup{in}&\Omega.\end{array}\right.$$

$$\Vert u{\Vert }_{W_p^{2,1}(Q)}\le C_p\left(\Vert f{\Vert }_{L^p(Q)}+\Vert u_0{\Vert }_{W_p^{2-\frac{2}{p}}(\Omega )}\right).$$

$$\left\{\begin{array}[]{lcl}\displaystyle\frac{\partial\hat{N}}{\partial t}=r_{N}-\mu_{N}\hat{N}-\beta_{1}\hat{N}\hat{A}-\alpha_{N}\gamma_{N}|\hat{D}|\hat{N},&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial\hat{A}}{\partial t}=r_{A}\hat{A}\left(1-\frac{\hat{A}}{k_{A}}\right)-(\mu_{A}+\epsilon_{A})\hat{A}-\alpha_{A}\gamma_{A}|\hat{D}|\hat{A},&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial\hat{D}}{\partial t}=\sigma\Delta\hat{D}+\mu\chi_{\omega}-\gamma\hat{D}\hat{A}-\gamma_{N}\hat{D}\hat{N}-\tau\hat{D},&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial\hat{D}}{\partial\eta}(\cdot)=0,&\textup{on}&\Gamma,\vspace{0.2cm}\\ \displaystyle\hat{N}(\cdot,0)=N_{0}(\cdot),\hat{A}(\cdot,0)=A_{0}(\cdot),\hat{D}(\cdot,0)=D_{0}(\cdot),&\textup{in}&\Omega.\end{array}\right.$$

$$\Lambda (\varphi )(x,t)=\frac{A_0(x)k_A{e}^{\lambda t}{e}^{-{\alpha }_A{\gamma }_A{\int }_0^t|\varphi (\xi ,x)|d\xi }}{k_A+A_0(x)r_A{\int }_0^t{e}^{\lambda s}{e}^{-{\alpha }_A{\gamma }_A{\int }_0^s|\varphi (\xi ,x)|d\xi }ds}$$

$$\Theta (\varphi )(x,t)=\frac{N_0(x)+r_N{\int }_0^t{e}^{{\mu }_{N}s}{e}^{{\alpha }_N{\gamma }_N{\int }_0^s|\varphi (x,\xi )|d\xi }{e}^{{\beta }_1{\int }_0^s\Lambda (\varphi )(x,\xi )d\xi }ds}{e^{{\mu }_{N}t}{e}^{{\alpha }_N{\gamma }_N{\int }_0^t|\varphi (x,\xi )|d\xi }{e}^{{\beta }_1{\int }_0^t\Lambda (\varphi )(x,\xi )d\xi }},$$

$$\left\{\begin{array}[]{lcl}\displaystyle\frac{\partial\hat{D}}{\partial t}=\sigma\Delta\hat{D}+\mu\chi_{\omega}-\gamma\hat{D}\Lambda(\hat{D})-\gamma_{N}\hat{D}\Theta(\hat{D})-\tau\hat{D},&\textup{in}&Q,\vspace{0.2cm}\\ \displaystyle\frac{\partial\hat{D}}{\partial\eta}(\cdot)=0,&\textup{on}&\Gamma,\vspace{0.2cm}\\ \displaystyle\hat{D}(\cdot,0)=D_{0}(\cdot),&\textup{in}&\Omega.\end{array}\right.$$

$$||\hat{D}|{|}_{W_4^{2,1}(Q)}\le C,$$

$$g^{\prime }(t)$$

$$\begin{array}[]{ll}\mbox{\bf(i)}&0\leq\Theta(\phi)(x,t)\leq||N_{0}||_{L^{\infty}(\Omega)}+r_{N}T;\vspace{0.2cm}\\ \mbox{\bf(ii)}&0\leq\Lambda(\phi)(x,t)\leq C_{\lambda}||A_{0}||_{L^{\infty}(\Omega)};\vspace{0.2cm}\\ \mbox{\bf(iii)}&\|\Lambda(\phi_{1})-\Lambda(\phi_{2})\|_{L^{\infty}{(Q)}}\leq C_{1}\|\phi_{1}-\phi_{2}\|_{L^{\infty}(Q)},\\ &where\ C_{1}\ is\ a\ constant\ depending\ on\ r_{A},k_{A},\alpha_{A},\gamma_{A},C_{\lambda},T\ and\ ||A_{0}||_{L^{\infty}(\Omega)};\vspace{0.2cm}\\ \mbox{\bf(iv)}&\|\Theta(\phi_{1})-\Theta(\phi_{2})\|_{L^{\infty}{(Q)}}\leq C_{2}\|\phi_{1}-\phi_{2}\|_{L^{\infty}(Q)},\\ &where\ C_{2}\ is\ a\ constant\ depending\ on\ r_{N},\mu_{N},\beta_{1},\alpha_{N},\gamma_{N},C_{\lambda},C_{1},T,\\ &||\phi_{1}||_{L^{\infty}(Q)},||\phi_{2}||_{L^{\infty}(Q)},||N_{0}||_{L^{\infty}(\Omega)}\ and\ ||A_{0}||_{L^{\infty}(\Omega)}.\end{array}$$

$$\begin{array}[]{rcl}\displaystyle\Theta(\phi)(x,t)=\frac{N_{0}(x)+r_{N}\int_{0}^{t}e^{\mu_{N}s}e^{\alpha_{N}\gamma_{N}\int_{0}^{s}|\phi(x,\xi)|d\xi}e^{\beta_{1}\int_{0}^{s}\Lambda(\phi)(x,\xi)d\xi}ds}{e^{\mu_{N}t}e^{\alpha_{N}\gamma_{N}\int_{0}^{t}|\phi(x,\xi)|d\xi}e^{\beta_{1}\int_{0}^{t}\Lambda(\phi)(x,\xi)d\xi}}\\ \displaystyle\hfil\\ \leq N_{0}(x)+r_{N}\frac{\int_{0}^{t}e^{\mu_{N}s}e^{\alpha_{N}\gamma_{N}\int_{0}^{s}|\phi(x,\xi)|d\xi}e^{\beta_{1}\int_{0}^{s}\Lambda(\phi)(x,\xi)d\xi}ds}{e^{\mu_{N}t}e^{\alpha_{N}\gamma_{N}\int_{0}^{t}|\phi(x,\xi)|d\xi}e^{\beta_{1}\int_{0}^{t}\Lambda(\phi)(x,\xi)d\xi}}.\end{array}$$

$$g(x,t)=\frac{{\int }_0^t{e}^{{\mu }_{N}s}{e}^{{\alpha }_N{\gamma }_N{\int }_0^s|\varphi (x,\xi )|d\xi }{e}^{{\beta }_1{\int }_0^s\Lambda (\varphi )(x,\xi )d\xi }ds}{e^{{\mu }_{N}t}{e}^{{\alpha }_N{\gamma }_N{\int }_0^t|\varphi (x,\xi )|d\xi }{e}^{{\beta }_1{\int }_0^t\Lambda (\varphi )(x,\xi )d\xi }},$$

$$\Theta (\varphi ,\phi )(x,t)\le N_0(x)+r_{N}T$$

$$\le ||N_0|{|}_{L^{\infty }(\Omega )}+r_{N}T.$$

$$\begin{array}[]{rcl}\displaystyle\Lambda(\phi)(x,t)=\frac{A_{0}(x)k_{A}e^{\lambda t}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi(\xi,x)|d\xi}}{k_{A}+A_{0}(x)r_{A}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi(\xi,x)|d\xi}ds}\\ \displaystyle\hfil\\ \leq\frac{1}{k_{A}}A_{0}(x)k_{A}e^{\lambda t}e^{-\alpha_{A}\gamma_{A}\int^{t}_{0}|\phi(x,\xi)|d\xi}\\ \\ \displaystyle\leq A_{0}(x)e^{\lambda t}\leq C_{\lambda}A_{0}(x)\leq C_{\lambda}||A_{0}||_{L^{\infty}(\Omega)}.\end{array}$$

$$|e^{z_2}-e^{z_1}|\le |z_2-z_1|,\forall z_1,z_2\le 0.$$

$$\begin{array}[]{rcl}\big{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}-e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\big{|}&\leq&\big{|}-\alpha_{A}\gamma_{A}\int_{0}^{t}(|\phi_{1}(x,\xi)|-|\phi_{2}(x,\xi)|)d\xi\big{|}\\ \\ &\leq&\alpha_{A}\gamma_{A}T||\phi_{1}-\phi_{2}||_{L^{\infty}(Q)}.\end{array}$$

$$\displaystyle\big{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}ds$$

$$\displaystyle\displaystyle e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}ds\big{|}$$

$$\displaystyle\displaystyle\big{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}-e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\big{|}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}ds$$

$$\displaystyle\displaystyle e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}\big{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{1}(x,\xi)|d\xi}-e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}\big{|}ds.$$

$$\begin{array}[]{rcl}\displaystyle\big{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}ds&-&\\ \\ \displaystyle e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{1}(x,\xi)|d\xi}ds\big{|}&\leq&\\ \\ \displaystyle 2\alpha_{A}\gamma_{A}C_{\lambda}T^{2}||\phi_{1}-\phi_{2}||_{L^{\infty}(Q)}.\end{array}$$

$$|\Lambda ({\varphi }_1)(x,t)-\Lambda ({\varphi }_2)(x,t)|$$

$$\displaystyle A_{0}(x)e^{\lambda t}\bigg{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}-e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{2}(x,\xi)|d\xi}\bigg{|}$$

$$\displaystyle\frac{1}{k_{A}}{A_{0}(x)}^{2}r_{A}e^{\lambda t}\bigg{|}e^{-\alpha_{A}\gamma_{A}\int_{0}^{t}|\phi_{1}(x,\xi)|d\xi}\int_{0}^{t}e^{\lambda s}e^{-\alpha_{A}\gamma_{A}\int_{0}^{s}|\phi_{2}(x,\xi)|d\xi}ds$$

Full text

1. Introduction

The main objective of this work is to perform a rigorous mathematical analysis of a system of nonlinear partial differential equations corresponding to a generalization of a mathematical model describing the growth of a tumor proposed in [6].

To describe the model, let $\Omega\subset\mathrm{I\!R\!}^{2}$, be an open and bounded set; let also $0<T<\infty$ be a given final time of interest and denote $t$ the times between $[0,T]$ and $Q=\Omega\times(0,T)$, the space-time cylinder and $\bar{\Gamma}=\partial\Omega\times(0,T)$, the space-time boundary. Then, the system of equations we are considering is the following:

[TABLE]

In [6], Fassoni studied an ODE system corresponding to system (1.1) in a spatially homogeneous setting. Such model describes the growth of a tumor and its effect on the normal tissue, the tissue response to the tumor and the application of chemotherapeutic treatments, without spatial heterogeneity. The aim of the authors was to understand the phenomena of cancer onset and treatment as transitions between different basins of attraction of the underlying ODE system. The equations of the model that were studied in [6] are

[TABLE]

where $N$ represents the number of normal cells in a given tissue of the human body, $A$ represents the number of tumor cells in the tissue and $D$ represents the concentration of a chemotherapeutic drug used to treat such a tumor.

Parameter $r_{N}$ represents a constant influx of new normal cells produced by the tissue stem cells and $\mu_{N}$ presents the natural mortality of normal cells. A constant influx is considered because the imperative dynamics within a formed tissue is the maintenance of a homeostatic state through the natural replenishment of old and dead cells, see [14].

On the other hand, tumor cells maintain their own growth program [7]. Thus, a density dependent growth is considered for tumor cells. The logistic growth is chosen due to its simplicity. Parameter $\mu_{A}$ represents the natural mortality of tumor cells, and $\epsilon_{A}$ represents an extra mortality rate due to apoptosis [4].

Parameters $\beta_{1}$ and $\beta_{3}$ encompass the many negative interactions exerted by tumor cells on normal cells and vice-versa, such as competition for nutrients and oxygen. Besides competition, parameter $\beta_{3}$ encompasses also the effects on normal cells of anti-growth and death signals released by normal cells. In the same way, the parameter $\beta_{1}$ encompasses also mechanisms developed by tumor cells that damage normal tissue, such as increased local acidity, growth suppression, and release of death signals [9].

The third equation of (1.2) describes the dynamics of chemotherapeutic drug concentration according the following assumptions. The drug has a constant infusion rate $\mu$ and a clearance rate $\tau$. Such constant infusion rate mimics a metronomic dosage, i.e., a near continuous and long-term administration of the drug. The absorption and deactivation of the drug by normal and cancerous cells are described in terms of the law of mass action with rates $\gamma_{N}$ and $\gamma_{A}$. Following the log-linear hypothesis [3], it is assumed that the amounts of drug absorbed by normal ($\gamma_{N}ND$) and cancerous cells ($\gamma_{A}AD$) kill such cells with rates $\alpha_{N}$ and $\alpha_{A}$, respectively. Although many models of cancer treatment do not consider drug absorption explicitly, in [6], the authors believe that it is an important fact to be considered, since, this phenomenon contributes to decrease the concentration of drug as time passes.

System (1.2) is similar to the classical Lotka-Volterra competition model, frequently used in models for tumor growth and population dynamics. The fundamental difference here is the use of a constant flux for normal cells instead of a logistic growth. Such constant flux, also used in other well-known models of cancer [5], removes the symmetry observed in the Lotka-Volterra equations, so that there is no steady state with $N=0$. Thus, it is impossible to observe the extinction of one of the populations (the normal cells in this case), as opposed to the Lotka-Volterra models. The authors of [6] claim that this is a realistic result since, roughly speaking, cancer ”does not win” by killing all the cells in the tissue, but by reaching a dangerous size that disrupts the proper functioning of the tissue and threatens the health of the individual.

In this work, we are not interested in analyzing the dynamics (stability, asymptotic behavior) of the model, as such study has already been made in [6]. Our objective is to study the existence and uniqueness of the solution of system (1.1). Such system extends the ODE model (1.2) to a more realistic situation by considering spatial variation of normal and cancer cells and the diffusion of the chemotherapeutic drug through the tissue, with diffusion coefficient $\sigma$ [2]. Further, it is also assumed that the drug influx is restricted to a limited region of the tissue, corresponding to a blood vessel passing transversely in such region. This is mathematically described in the model by the expression $\mu\chi_{\omega}$, where $\chi_{\omega}$ is the characteristic function of the subset $\omega\subset\Omega$. Finally, due to mathematical necessity to simplify the model, we set $\beta_{3}=0$. This corresponds to a situation where normal cells do not exert negative effects on tumor cells, and is a plausible biological assumption, since there are many tumors that develop resistance to the normal tissue’ mechanisms which suppress tumor growth [9].

The paper is organized as follows. In Section 2 we present the technical hypothesis and state our main result. In Section 3 we study an auxiliary problem. Using its solution, we prove our main result in Section 4. In Section 5 we present numerical simulations illustrating model behavior.

2. Technical hypotheses and main result

Let $\Omega\subset\mathrm{I\!R\!}^{2}$ be a domain with boundary $\partial\Omega$, $0\leq T<\infty$, and denote $Q=\Omega\times(0,T)$ and $\Gamma=\partial\Omega\times(0,T)$. We will use standard notations for Sobolev spaces, i.e., given $1\leavevmode\nobreak\ \leq\leavevmode\nobreak\ p\leavevmode\nobreak\ \leq\leavevmode\nobreak\ +\infty$ and $k\in\mathbb{N}$, we denote

[TABLE]

when $p=2$, as usual we denote $W_{2}^{k}(\Omega)=H^{k}(\Omega)$; properties of these spaces can be found for instance in Adams [1, Theorem 5.4, p. 97]. Problem (1.1) will be studied in the standard functional spaces denoted by

[TABLE]

[TABLE]

and

[TABLE]

where $B$ is suitable Banach space, and the norm is given by $\|f(t)\|_{L^{p}(0,T;B)}=\|\ \|f(t)\|_{B}\ \|_{L^{p}((0,T))}$. We remark that $L^{p}(Q)=L^{p}((0,T);L^{p}(\Omega))$. Results concerning these spaces can be found for instance in Ladyzhenskaya [10] and Mikhaylov [15].

Next, we state some hypotheses that will be assumed throughout this article.

2.1. Technical Hypotheses:

• (i)

$\Omega\subset\mathbb{R}^{2}$ is a bounded $C^{2}$-domain;

• (ii)

$0<T<\infty$, and $Q=\Omega\times(0,T)$;

• (iii)

$N_{0},A_{0}\in L^{\infty}(\Omega)$ and $D_{0}\in W^{\frac{3}{2}}_{4}(\Omega)$, satisfying $\frac{\partial D_{0}}{\partial\eta}(\cdot)=0,\textup{ on }\partial\Omega$;

• (iv)

$0\leq D_{0}\leq\frac{\mu}{\tau}$ and $N_{0},A_{0}\geq 0$ a.e. on $\Omega$.

Remark 2.1**.**

The constraints imposed in (iv) on the initial conditions are natural biological requirements.

2.2. Main result:

Theorem 2.2**.**

Assume that the Technical Hypotheses 2.1 hold; then, there exists a unique nonnegative solution $(N,A,D)\in W\times W\times W^{2,1}_{4}(Q)$ of Problem (1.1). Moreover, $N,A$ and $D$ are functions satisfying

[TABLE]

and

[TABLE]

where $C$ is a constant depending on $r_{N}$, $\mu_{N}$, $\beta_{1}$, $\alpha_{N}$, $\gamma_{N}$, $C_{\lambda}$, $r_{A}$, $k_{A}$, $\mu_{A}$, $\epsilon_{A}$, $\alpha_{A}$, $\gamma_{A}$, $\mu$, $\tau$, $T$, $\omega$, $||N_{0}||_{L^{\infty}(\Omega)}$, $||A_{0}||_{L^{\infty}(\Omega)}$ and $||D_{0}||_{W^{\frac{3}{2}}_{4}(\Omega)}$.

Remark 2.3**.**

The explicit knowledge on how the constant $C$ appearing in the above estimates depends on the given data is important for applications in related control problems.

2.3. Known technical results:

To ease the references, we also state some technical results to be used in this paper. The first one is sometimes called the Lions-Peetre embedding theorem (see Lions [11], pp.15); it is also a particular case of Lemma 3.3, pp.80, in Ladyzhenskaya [10]: (obtained by taking $l=1$ and $r=s=0$).

Lemma 2.4**.**

Let $\Omega$ be a domain of $\mathrm{I\!R\!}^{n}$ with boundary $\partial\Omega$ satisfying the cone property. Then, the functional space $W^{2,1}_{p}(Q)$ is continuously embedded in $u\in L^{q}(Q)$ for $q$ satisfying: (i) $1\leq q\leq\frac{p(n+2)}{n+2-2p}$, if $p<\frac{n+2}{2}$; (ii) $1\leq q<\infty$, if $p=\frac{n+2}{2}$ and (iii) $q=\infty$, if $p>\frac{n+2}{2}$. In particular, for such $q$ and any function $u\in W^{2,1}_{p}(Q)$ we have that

[TABLE]

with a constant $C$ depending only on $\Omega$, $T$, $p$, $q$, $n$.

In the cases (ii), (iii) or in (i) when $\displaystyle 1\leq q<\frac{p(n+2)}{n+2-2p}$, the referred embedding is compact.

Next, we consider the following simple parabolic initial-boundary value problem:

[TABLE]

Existence and uniqueness of solutions for this problem is a particular case of Theorem $9.1$, pp.$341$, in Ladyzenskaya [10] for the case of Neumann boundary condition, according to the remarks at the end Chapter IV, section 9, p. 351 in [10]. In the following, we state this particular result, stressing the dependencies certain norms of the coefficients, that will be important in our future arguments.

Proposition 2.5**.**

Let $\Omega$ be a bounded domain in $\mathbb{R}^{n}$, with a $C^{2}$ boundary $\partial\Omega$, $a_{ij}$ be bounded continuous functions in $Q$, and $q>1$. Assume that

1. (1)

$a_{ij}\in C(\bar{Q})$, $i,j=1,\ldots,n$; $[a_{ij}]_{n\times n}$ is a real positive matrix such that for some positive constant $\beta$ we have $\sum\limits_{i,j=1}^{n}a_{ij}(x,t)\xi_{i}\xi_{j}\geq\beta|\xi|^{2}$ for all $(x,t)\in Q$ and all $\xi\in R^{n}$, ; 2. (2)

$\displaystyle f\in L^{p}(Q)$; 3. (3)

$\displaystyle a_{i}\in L^{r}(Q)$* with either r=\max\big{(}p,n+2\big{)} if $p\neq n+2$ or $r=n+2+\varepsilon$, for any $\varepsilon>0$, if $p=n+2$;* 4. (4)

$\displaystyle a\in L^{s}(Q)$* with either s=\max\big{(}p,(n+2)/2\big{)} if $\displaystyle p\neq(n+2)/2$ or $s=(n+2)/2+\varepsilon$, for any $\varepsilon>0$, if $\displaystyle p=(n+2)/2$.* 5. (5)

$b_{i},b\in C^{2}(\bar{\Gamma})$, $i=1,\ldots,n$, and the coefficients $b_{i}(x,t)$ satisfy the condition $\left|\sum\limits_{i=1}^{n}b_{i}(x,t)\eta_{i}(x)\right|\geq\delta>0$ for $a.e.$ in $\partial\Omega\times(0,T)$, where $\eta_{i}(x)$ is the $i^{th}$-component of the unitary outer normal vector to $\partial\Omega$ in $x\in\partial\Omega$; 6. (6)

$u_{0}\in\ W^{2-\frac{2}{p}}_{p}(\Omega)$* with $p\neq 3$ and satisfying the compatibility condition *

$\displaystyle\sum\limits_{i=1}^{n}b_{i}\frac{\partial u_{0}}{\partial x_{i}}+b\ u_{0}=0$* on $\partial\Omega$ when $p>3$.*

Then, there exists a unique solution $u\in W^{2,1}_{p}(Q)$ of Problem (2.1); moreover, there is a positive constant $C_{p}$ such that the solution satisfies

[TABLE]

Such constant $C_{p}$ depends only on $\Omega$, $T$, $p$, $r$, $s$, $\beta$, $\delta$ and on the norms $\|b_{i}\|_{C^{2}(\bar{\Gamma})}$, $\|b\|_{C^{2}(\bar{\Gamma})}$, $\|a_{ij}\|_{C(\bar{Q})}$, $\|a_{i}\|_{L^{r}(Q)}$ and $\|a\|_{L^{s}(Q)}$. Moreover, we may assume that the dependencies of $C_{p}$ on stated the norms are non decreasing.

Remark 2.6**.**

The result set out in Proposition 2.5 can be formulated for the parabolic problem with Dirichlet conditions (see Ladyzenskaya [10, Theorem 9.1, pp.$341$]). In the problem with Dirichlet condition the compatibility condition in Proposition 2.5-($6$) can be replaced by $u_{0}=0$ on $\partial\Omega$ when $p>3/2$. This way, all the results in this paper holds if we replaced the Neumann conditions by Dirichlet conditions.

3. An auxiliary problem

In this section we will prove an auxiliary result to be used in the proof of Theorem 2.2. To cope with difficulties with the signs of certain terms during the derivation of the estimates, we firstly have to consider the following modified problem:

[TABLE]

Now we observe that, since the equation for $\hat{N}$ in this last problem is, for each $x\in\Omega$, an ordinary differential equation which is linear in $\hat{N}$, we can find an explicit expression for it in terms of $\hat{A}$ and $|\hat{D}|$. However, $\hat{A}$ is, for each $x\in\Omega$, a nonlinear differential equation in $\hat{A}$, and we can determine its explicit expression in terms of $|\hat{D}|$ using Bernoulli’s method. Using these observations and setting $\lambda=r_{A}-(\mu_{A}+\epsilon_{A})$, we introduce operators $\Lambda:L^{\infty}(Q)\to L^{\infty}(Q)$ and $\Theta:L^{\infty}(Q)\to L^{\infty}(Q)$, defined respectively by

[TABLE]

and

[TABLE]

where $0\leq s\leq t\leq T$.

Remark 3.1**.**

Thus, $(\hat{N},\hat{A},\hat{D})$ is a solution of (3.1) if, and only if, $\hat{N}=\Theta(\hat{D})$, $\hat{A}=\Lambda(\hat{D})$ and $\hat{D}$ satisfies the following integro-differential system:

[TABLE]

Remark 3.2**.**

Notice that, to guarantee that $(N,A,D)$, with $D=\hat{D}$, $N=\Theta(\hat{D})$ and $A=\Lambda(\hat{D})$ is also a solution of system (1.1), it is enough to prove that the solution $\hat{D}$ of Problem (3.4) is nonnegative.

For the Problem 3.4, we have the following existence result:

Proposition 3.3**.**

Assuming that the Technical Hypotheses 2.1 hold, there exists at least one nonnegative solution $\hat{D}\in W^{2,1}_{4}(Q)$ of Problem (3.4). Moreover, $\hat{D}\leq\frac{\mu}{\tau}$ a.e. in $Q$ and

[TABLE]

where $C$ is a constant depending on $\mu$, $T$, $\omega$ and $||D_{0}||_{W^{\frac{3}{2}}_{4}(\Omega)}$.

Lemma 3.4**.**

Let $f:(0,T)\to\mathbb{R}$ differentiable such that $f(t)>0$ and $f^{\prime}(t)\geq 0$. If $g(t)=\frac{\int_{0}^{t}f(x)dx}{f(t)}$, then $g(t)\leq T$, for all $t\in(0,T)$.

Proof: Since $f$ is continuous in $(0,T)$, it follows that

[TABLE]

As $f(t)>0$ we have $g(t)\geq 0$ and using the fact that $f^{\prime}(t)\geq 0$ we obtain $\frac{f^{\prime}(t)}{f(t)}g(t)\geq 0$. Therefore, $g^{\prime}(t)\leq 1$, which suggests $g(t)\leq t$, for all $t\in(0,T)$. Thus, $g(t)\leq T$, as intended. $\Box$

Since in the proof of existence of solutions of (3.4) the expression of $\Lambda$ and $\Theta$ will play important roles, we state some of their properties in the following:

Lemma 3.5**.**

If $N_{0},A_{0}\in L^{\infty}(\Omega)$ and $C_{\lambda}=\max\{1,e^{\lambda T}\}$, then for any $\phi,\phi_{1},\phi_{2}\in L^{\infty}(Q)$ and for almost every $(x,t)\in Q$, there holds

[TABLE]

Proof (i) and (ii): By the expressions (3.2) and (3.3) it is immediate that $\Lambda(\phi)(x,t),\Theta(\phi)(x,t)\geq 0$. To prove that $\Theta(\phi)(x,t)\leq||N_{0}||_{L^{\infty}(\Omega)}+r_{N}T$, we observe that

[TABLE]

Fixed $x\in\Omega$, we define

[TABLE]

and using the Lemma 3.4 with $f(x,t)=e^{\mu_{N}t}e^{\alpha_{N}\gamma_{N}\int_{0}^{t}|\phi(x,\xi)|d\xi}e^{\beta_{1}\int_{0}^{t}\Lambda(\phi)(x,\xi)d\xi}$, it follows that

[TABLE]

To prove that $\Lambda(\phi)(x,t)\leq C_{\lambda}||A_{0}||_{L^{\infty}(\Omega)}$, note that

[TABLE]

Proof (iii): We firstly need to observe that, due to the mean value inequality, given any $z_{1},z_{2}\in\mathrm{I\!R\!}$, there is $\theta=\theta(z_{1},z_{2})$ such that $e^{z_{2}}-e^{z_{1}}=e^{(1-\theta)z_{1}+\theta z_{2}}(z_{2}-z_{1})$; in particular, for any $z_{1},z_{2}\leq 0$ we also have $(1-\theta)z_{1}+\theta z_{2}\leq 0$ and thus

[TABLE]

Secondly, we note that by the inequality (3.5) and by $\phi_{i}\in L^{\infty}(Q)$, $i=1,2$, we obtain

[TABLE]

Thirdly, we observe that

[TABLE]

How $e^{\lambda T}\leq C_{\lambda}$ and $e^{-\alpha_{A}\gamma_{A}||\phi_{2}||_{L^{\infty}(Q)}}\leq 1$, and using study analogous to that done in (3.6), we obtain that

[TABLE]

Finally, the expression in (3.3) suggests

[TABLE]

and using the estimates obtained in (3.6) and (3.7) and making the possible simplifications, we obtain

[TABLE]

for almost everything $(x,t)\in Q$, i.e.,

[TABLE]

Proof (iv): First, note that

[TABLE]

and by the inequality (3.5) and by $\Lambda(\phi_{i}),\phi_{i}\in L^{\infty}(Q)$, $i=1,2$, we obtain

[TABLE]

Since

[TABLE]

doing $||\phi||_{L^{\infty}(Q)}=\max\{||\phi_{1}||_{L^{\infty}(Q)},||\phi_{2}||_{L^{\infty}(Q)}\}$ and study analogous to that done in (3.9), guarantees us

[TABLE]

Finally, the expression in (3.3) suggests

[TABLE]

and using the estimates obtained in (3.8), (3.9) and (3.10) and making the possible simplifications, we obtain

[TABLE]

for almost everything $(x,t)\in Q$, i.e.,

[TABLE]

$\Box$

3.1. Proof of Proposition 3.3

To not overburden the notation, in this subsection we denote $D$ as a generic solution of the equations that follows.

To get a solution of problem (3.4), we will apply the Leray-Schauder fixed point theorem to the mapping $\Psi$ defined as follows:

[TABLE]

where $D$ is the unique solution of

[TABLE]

with $\Lambda(\phi)$ and $\Theta(\phi)$ given by (3.2) and (3.3), respectively.

To apply such theorem we present next a sequence of lemmas:

Lemma 3.6**.**

Suppose $N_{0},A_{0}\in L^{\infty}(\Omega)$ and $D_{0}\in W^{\frac{3}{2}}_{4}(\Omega)$. Then the mapping $\Psi:[0,1]\times L^{\infty}(Q)\rightarrow L^{\infty}(Q)$ is well defined.

Proof: We affirm that the coefficients of the Problem 3.12 satisfy the hypotheses of the Proposition 2.5. For example, it is immediate that $-l\gamma\Lambda(\phi)-l\gamma_{N}\Theta(\phi)-\tau\in L^{4}(Q)$, because by Lemma 3.5, $\Lambda(\phi),\Theta(\phi)\in L^{\infty}(Q)$. Thus, we conclude that there is a unique solution $D\in W^{2,1}_{4}(Q)$ of problem 3.12. Moreover, $D$ satisfies the following estimate:

[TABLE]

Finally, from Lemma 2.4, we have $W^{2,1}_{4}(Q)\hookrightarrow L^{\infty}(Q)$, and we conclude that the operator $\Psi$ in well defined. $\Box$

Lemma 3.7**.**

Suppose $D$ is a solution of (3.12) and $0\leq D_{0}\leq\frac{\mu}{\tau}$ a.e. in $\Omega$, then $0\leq D\leq\frac{\mu}{\tau}$ a.e. in $Q$.

Proof: Multiplying the first equation in $(\ref{P6})$ by $D^{-}$ and integrating into $\Omega$, we get

[TABLE]

Thus,

[TABLE]

and using Gronwall’s inequality and the fact that $D_{0}\geq 0$ a.e. in $\Omega$, we obtain

[TABLE]

that is, $||D^{-}(\cdot,t)||_{L^{2}(\Omega)}=0$ for all $t\in(0,T)$, where we conclude that $D^{-}=0$ a.e. in $Q$ and therefore $D\geq 0$ a.e. in $Q$.

Now, we observe that the first equation in (3.12) can be rewritten as

[TABLE]

Multiplying by $(D-\frac{\mu}{\tau})^{+}$ and integrating in $\Omega$, we obtain

[TABLE]

that is,

[TABLE]

Thus, using Gronwall’s inequality and the fact that $D_{0}\leq\frac{\mu}{\tau}$ a.e. in $\Omega$, it follows that

[TABLE]

that is, ||\big{(}D(\cdot,t)-\frac{\mu}{\tau}\big{)}^{+}||_{L^{2}(\Omega)}=0 for all $t\in(0,T)$, and therefore \big{(}D-\frac{\mu}{\tau}\big{)}^{+}=0 a.e. in $Q$, and we conclude that $D\leq\frac{\mu}{\tau}$ a.e. in $Q$.

$\Box$

Lemma 3.8**.**

For each fixed $l\in[0,1]$, the mapping $\Psi(l,\cdot):L^{\infty}(Q)\rightarrow L^{\infty}(Q)$ is compact, i.e., it is continuous and maps bounded sets into relatively compacts sets.

Proof: The functions $\Psi(l,\phi_{1})=D_{1}$ and $\Psi(l,\phi_{2})=D_{2}$ satisfy the system

[TABLE]

with $i=1,2$; letting $\tilde{D}=D_{1}-D_{2}$, we have

[TABLE]

Using the Proposition 2.5 and the fact that $L^{\infty}(Q)\hookrightarrow L^{4}(Q)$ and $D_{1}\leq\frac{\mu}{\tau}$, we get

[TABLE]

Then, by Lemmas 3.5 and 2.4, we finally have

[TABLE]

where $C$ depends on $\bar{C}_{p}$, $C_{1}$, $C_{2}$, $\gamma$, $\gamma_{N}$, $\mu$, $\tau$ and the immersion constant.

To show that $\Psi(l,\cdot)$ is compact, we use the fact that the immersion $W^{2,1}_{4}(Q)\hookrightarrow L^{\infty}(Q)$ is compact and that $\Psi(l,\cdot)$ is the composition between the inclusion operator and the solution operator, i.e., $\Psi(l,\cdot):L^{\infty}(Q)\rightarrow W^{2,1}_{4}(Q)\rightarrow L^{\infty}(Q)$. $\Box$

Lemma 3.9**.**

Given a bounded subset $B\subset L^{\infty}(Q)$, for each $\phi\in B$, the mapping $\Psi(\cdot,\phi):[0,1]\rightarrow L^{\infty}(Q)$ is uniformly continuous with respect to $B$.

Proof: Since $B\in L^{\infty}(Q)$ is bounded, there is $r_{B}\geq 0$ such that, for any $\phi\in B$, we have $||\phi||_{L^{\infty}(Q)}\leq r_{B}$. Now, let us fix $\phi\in L^{\infty}(Q)$ and consider $l_{1},l_{2}\in[0,1]$ and denote $\Psi(l_{1},\phi)=D_{1}$, $\Psi(l_{2},\phi)=D_{2}$ and $\tilde{D}=D_{1}-D_{2}$. Then, $\tilde{D}$ satisfies

[TABLE]

Using the Proposition 2.5 and the fact that $L^{\infty}(Q)\hookrightarrow L^{4}(Q)$, $D_{1}\leq\frac{\mu}{\tau}$, we get

[TABLE]

Then, by Lemmas 3.5 and 2.4, we finally have

[TABLE]

where $C$ depends on $\bar{C}_{p}$, $\gamma$, $\gamma_{N}$, $\mu$, $\tau$, $r_{N}$, $T$, $C_{\lambda}$, $||N_{0}||_{L^{\infty}(\Omega)}$, $||A_{0}||_{L^{\infty}(\Omega)}$ and the immersion constant. $\Box$

Lemma 3.10**.**

Suppose $D_{0}\leq\frac{\mu}{\tau}$ a.e. in $\Omega$, then there exists a number $\rho>0$ such that, for any $l\in[0,1]$ and any possible fixed point $D\in L^{\infty}(Q)$ of $\Psi(l,\cdot)$, there holds $\|D\|_{L^{\infty}(Q)}<\rho$.

Proof: Let $D\in L^{\infty}(Q)$ such that $\Psi(l,D)=D$. The analogous demonstration made in Proposition 3.7 guarantees us $||D||_{L^{\infty}(Q)}\leq\frac{\mu}{\tau}$. Therefore, just take $\rho=\frac{\mu}{\tau}+1$. $\Box$

Lemma 3.11**.**

The mapping $\Psi(0,\cdot):L^{\infty}(Q)\rightarrow L^{\infty}(Q)$ has a unique fixed point.

Proof: Indeed, letting $l=0$ in 3.12, $D$ is a fixed point of $\Psi(0,\cdot)$ if, and only if, $D$ is the unique solution to the problem

[TABLE]

But Proposition 2.5 guarantees the existence of a unique solution $D\in W^{2,1}_{4}(Q)\hookrightarrow L^{\infty}(Q)$ of this last problem; therefore $\Psi(0,\cdot)$ has a unique fixed point in $L^{\infty}(Q)$. $\Box$

Proposition 3.12**.**

There is a nonnegative solution $\hat{D}\in W^{2,1}_{4}(Q)$ of the problem (3.4).

Proof: From Lemmas 3.6, 3.8, 3.9, 3.10 and 3.11, we conclude that the mapping $\Psi:[0,1]\times L^{\infty}(Q)\rightarrow L^{\infty}(Q)$ satisfies the hypotheses of the Leray-Schauder’s fixed point theorem (see Friedman [8, pp. 189, Theorem 3]). Thus, there exists $\hat{D}\in L^{\infty}(Q)$ such that $\Psi(1,\hat{D})=\hat{D}$. Moreover, by Lemmas 3.6 and 3.7, $\hat{D}\in W^{2,1}_{4}(Q)$ is nonnegative and $\hat{D}$ is the required solution of (3.4). $\Box$

4. Proof of Theorem 2.2

Proposition 4.1**.**

There is a nonnegative solution $(\hat{N},\hat{A},\hat{D})\in L^{\infty}(Q)\times L^{\infty}(Q)\times W^{2,1}_{4}(Q)$ of the modified problem (3.1).

Proof: Just combine the Proposition 3.12, the Remark 3.1 and the Lemma 3.5. $\Box$

Remark 4.2**.**

We affirm that $\hat{N},\hat{A}\in W$. Indeed, by Lemma 3.5 we know that $\hat{N}=\Theta(\hat{D}),\hat{A}=\Lambda(\hat{D})\in L^{\infty}(Q)$. Moreover, returning to the first equation of (3.1), using the Lemmas 3.5 and 3.7, it follows that:

[TABLE]

a.e. in $Q$, i.e., $\hat{N}_{t}\in L^{\infty}(Q)$.

Moreover, returning to the second equation of (3.1) and using, again, the Lemmas 3.5 and 3.7, we get:

[TABLE]

a.e. in $Q$, i.e., $\hat{A}_{t}\in L^{\infty}(Q)$.

Proposition 4.3**.**

There is a nonnegative solution $(N,A,D)\in W\times W\times W^{2,1}_{4}(Q)$ of problem (1.1).

Proof: Just combine the Proposition 4.1 and the Remarks 3.2 and 4.2. $\Box$

Proposition 4.4**.**

The solution $(N,A,D)$ of the problem (1.1) is unique.

Proof: Let $(N_{1},A_{1},D_{1})$ and $(N_{2},A_{2},D_{2})$ be solutions to the problem (1.1); if $\tilde{N}=N_{1}-N_{2},\tilde{A}=A_{1}-A_{2}$ and $\tilde{D}=D_{1}-D_{2}$, then $\tilde{N}$, $\tilde{A}$ and $\tilde{H}$ satisfy the following problems, respectively:

[TABLE]

[TABLE]

[TABLE]

Multiplying the first equation of (4.3) by $\tilde{N}$, integrating into $\Omega$, using the fact that $N_{1}\leq||N_{0}||_{L^{\infty}(\Omega)}+r_{N}T$ and the inequality of Young, we have

[TABLE]

where $C$ depends on $\beta_{1}$, $\alpha_{N}$, $\gamma_{N}$, $r_{N}$, $T$ and $||N_{0}||_{L^{\infty}(\Omega)}$.

Now, multiplying the first equation of (4.4) by $\tilde{A}$, integrating into $\Omega$, using the fact that $A_{1}\leq C_{\lambda}||A_{0}||_{L^{\infty}(\Omega)}$ and the inequality of Young, we obtain

[TABLE]

where $C$ depends on $r_{A}$, $\alpha_{A}$, $\gamma_{A}$ $C_{\lambda}$ and $||A_{0}||_{L^{\infty}(\Omega)}$.

Lastly, multiplying the first equation of (4.5) by $\tilde{H}$, integrating into $\Omega$, using the fact that $D\leq\frac{\mu}{\tau}$ and the inequality of Young, we obtain

[TABLE]

where $C$ depends on $\gamma$, $\gamma_{N}$, $\mu$ and $\tau$.

Thus,

[TABLE]

and using the Gronwall’s inequality, we finally

[TABLE]

that is, $||\tilde{N}(\cdot,t)||_{L^{2}(\Omega)}^{2}+||\tilde{A}(\cdot,t)||_{L^{2}(\Omega)}^{2}+||\tilde{D}(\cdot,t)||_{L^{2}(\Omega)}^{2}=0$, for all $t\in(0,T)$. Where we conclude $\tilde{N}=\tilde{A}=\tilde{D}=0$ a.e. in $Q$ and therefore $N_{1}=N_{2},A_{1}=A_{2}$ and $D_{1}=D_{2}$ a.e. in $Q$. $\Box$

5. Numerical simulations

In this section, we provide numerical simulations illustrating different model behaviors. The settings and methods used to implement the simulations are the following. We consider the spatial domain as a square $\Omega=[0,L]\times[0,L]$, with $L=1$, discretized with $n=50$ steps $\Delta x=\Delta y=L/n=0.02$. The Laplacian $\Delta D$ is approximated by second order centered finite differences and the coupled ODE system arising from such discretization is solved with the method of lines in the software Mathematica. The simulations run from time $t=0$ until $t=25$ (which is enough to achieve stationary behavior in all simulations).

The initial conditions for numerical simulations are $N(x,0)=N_{2}$, $A(x,0)=A_{2}$, $D(x,0)=0$, where $(N_{2},A_{2},0)$ is a globally asymptotically stable equilibrium point for the ODE system (1.2) without treatment ($\nu=0$). The expressions for $N_{2}$ and $A_{2}$ are:

[TABLE]

Such equilibrium is allways globally asymptotically stable in system (1.2) (see details in [6]). From the biological point of view, these initial conditions correspond to the start of chemotherapy application when a tumor is already a formed, where the normal cells were not able to control tumor growth, and no chemotherapy was applied until the tumor reached a stationary state.

To avoid large numbers and numerical instabilities, we re-scale the populations with respect to their possible maximum values, setting $N\leftarrow N/(r_{N}/\mu_{N})$ and $A\leftarrow A/K_{A}$. Therefore, the population sizes range from [math] to $1$. The re-scaled parameter values used in the model simulations were fixed to

[TABLE]

[TABLE]

These values were chosen to describe: normal cells that reach the equilibrium $N=r_{N}/\mu_{N}=1$ at absence of tumor cells; a tumor with the same carrying capacity of normal cells ($K_{A}=r_{N}/\mu_{N}=1$) and a greater absorption of the chemotherapeutic drug by tumor cells in comparison with normal cells ($\gamma_{A}>\gamma_{N}$), due to the drug specificity.

In order to illustrate different biological outcomes in the model simulations, we allowed the following parameters to assume different values: the chemotherapeutic drug cytotoxicity against cancer cells $\alpha_{A}$, the diffusion coefficient of the chemotherapeutic drug $\sigma$ and the chemotherapy infusion rate $\mu$. We will show that these properties of the drug and the infusion rate are crucial for determining an effective treatment. We also simulated different positions for the subset $\omega$, which is a mathematical description of a blood vessel crossing the tissue, from where the chemotherapy enters the tissue. The values for parameters $\alpha_{A}$, $\sigma$, $\mu$ and the position of $\omega$ used in each simulation are indicated in Table 1. We present the following results.

In the first simulation of system (1.1), we confirm that our model and numerical methods are able to reproduce the expected biological behavior (Figure 1). The blood vessel crosses the tissue at its center, i.e., $\omega=[0.45,0.55]\times[0.45,0.55]$. We use the following parameter values: $\alpha_{A}=5$, $\mu=3$, and $\sigma=0.1$. With such values, the chemotherapy is not able to lead to tumor extinction. We observe that tumor cells that are near the blood vessel are eliminated but not extinct by the chemotherapeutic effect, and those which are distant from the blood vessel persist (Figure 1).

In order to make easier to illustrate the model dynamics, we present the results of next simulations in a one-dimensional domain $\Omega=[0,1]$. In Simulation 2, we use the same parameters values used in Simulation 1 (see Table 1), but increase the chemotherapy toxicity $\alpha_{A}$ and move the blood vessel to the left side of the tissue, $\omega=[0,0.1]$. Although the tumor cells in the vicinity of the blood vessel are extinct, the chemotherapy is still not able to eliminate the distant tumor cells (Figure 2). Thus, we observe tumor persistence in the long-term. In Simulation 3, we keep the parameters as in Simulation 2, but increase the chemotherapy infusion rate $\mu$ (mimicking a higher dose). We observe that the tumor cells are extinct in the entire tissue (Figure 3). In Simulation 4, we illustrate other mechanism to achieve tumor extinction: instead of increasing drug dose, we adopt the parameter values of Simulation 2, but increase the drug diffusion $\sigma$, so that it is capable to spread over the entire tissue and effectively eliminate all tumor cells (Figure 4). Finally, in Simulation 5, we also adopt the parameter values of Simulation 2, but increase the chemotherapy toxicity against tumor cells $\alpha_{A}$. This also leads to tumor extinction (Figure 5). An advantage of the strategies adopted in Simulations 4 and 5, in comparison with Simulation 3 (increasing dose), is that the former lead to less side effects. Simulation 3 describes the use of a drug which spreads faster, while Simulation 5 illustrates the use of a more potent and specific drug, which targets more tumor cells but not more normal cells ($\alpha_{N}$ was not changed). Taken together, these simulations and the different outcomes observed for different parameter values confirm the ability of the model to consistently describe tumor chemotherapy and illustrate the potential of mathematical models to provide testable hypothesis that could be studied together with clinicians in order to achieve better results in the treatment of cancer.