This work is licensed under the Creative Commons Attribution 4.0 International License.
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 Ω ⊂ IR2, be an open and bounded set; let also 0 < T < ∞ be a given final time of interest and denote t the times between [0, T] and Q = Ω × (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:
\left\{{\matrix{{{{\partial N} \over {\partial t}} = {r_N} - {\mu_N}N - {\beta_1}NA - {\alpha_N}{\gamma_N}DN,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial A} \over {\partial t}} = {r_A}A\left({1 - {A \over {{k_A}}}} \right) - ({\mu_A} + {\varepsilon_A})A - {\alpha_A}{\gamma_A}DA,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial D} \over {\partial t}} = \sigma \Delta D + \mu {\chi_\omega} - {\gamma_A}DA - {\gamma_N}DN - \tau D,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial D} \over {\partial \eta}} = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{N(\cdot,0) = {N_0}(\cdot),A(\cdot,0) = {A_0}(\cdot),D(\cdot,0) = {D_0}(\cdot),} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
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
\left\{{\matrix{{{{dN} \over {dt}} = {r_N} - {\mu_N}N - {\beta_1}NA - {\alpha_N}{\gamma_N}DN,} \hfill\cr{{{dA} \over {dt}} = {r_A}A\left({1 - {A \over {{k_A}}}} \right) - ({\mu_A} + {\varepsilon_A})A - {\beta_3}NA - {\alpha_A}{\gamma_A}DA,} \hfill\cr{{{dD} \over {dt}} = \mu- {\gamma_A}DA - {\gamma_N}DN - \tau D,} \hfill\cr}} \right.
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 rN represents a constant influx of new normal cells produced by the tissue stem cells and μ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 μA represents the natural mortality of tumor cells, and εA represents an extra mortality rate due to apoptosis [4].
Parameters β1 and β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 β3 encompasses also the effects on normal cells of anti-growth and death signals released by normal cells. In the same way, the parameter β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 μ and a clearance rate τ. 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 γN and γA. Following the log-linear hypothesis [3], it is assumed that the amounts of drug absorbed by normal (γNND) and cancerous cells (γAAD) kill such cells with rates αN and α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 σ [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 μχω, where χω is the characteristic function of the subset ω ⊂ Ω. Finally, due to mathematical necessity to simplify the model, we set β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.
Technical hypotheses and main result
Let Ω ⊂ IR2 be a domain with boundary ∂Ω, 0 ≤ T < ∞, and denote Q = Ω × (0, T) and Γ = ∂Ω × (0, T). We will use standard notations for Sobolev spaces, i.e., given 1 ≤ p ≤ +∞ and k ∈ ℕ, we denote
W_p^k(\Omega) = \left\{{f \in {L^p}(\Omega):{D^\alpha}f \in {L^p}(\Omega),|\alpha | \le k} \right\};
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
\matrix{{W_q^{2,1}(Q) = \left\{{f \in {L^q}(Q):{D^\alpha}f \in {L^q}(Q),{\kern 1pt} \forall 1 \le |\alpha | \le 2,{f_t} \in {L^q}(Q)} \right\},}\cr{W = \left\{{f \in {L^\infty}(Q):{f_t} \in {L^\infty}(Q)} \right\}}\cr}
and
{L^p}(0,T;B) = \left\{{f:(0,T) \to B:\parallel f(t){\parallel_{{L^p}(0,T;B)}} <+ \infty} \right\},
where B is suitable Banach space, and the norm is given by ||f(t)||Lp (0, T; B) = ||||f(t)||B ||Lp ((0, T)). We remark that Lp(Q) = Lp((0, T);Lp(Ω)). 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.
0 \le {D_0} \le {\mu\over \tau}
and N0, A0 ≥ 0 a.e. on Ω.
Remark 2.1.
The constraints imposed in(iv)on the initial conditions are natural biological requirements.
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_4^{2,1}(Q)of Problem (1.1). Moreover, N, A and D are functions satisfyingN \le ||{N_0}|{|_{{L^\infty}(Q)}} + {r_N}T,\;A \le {C_\lambda}||{A_0}|{|_{{L^\infty}(\Omega)}}\;a.e.\;in\;Qand||N|{|_W} + ||A|{|_W} + ||D|{|_{W_4^{2,1}(Q)}} \le C,where C is a constant depending on rN, μN, β1, αN, γN, Cλ, rA, kA, μA, εA, αA, γA, μ, τ, T, ω, ‖N0‖L∞ (Ω), ‖A0‖L∞ (Ω)and||{D_0}|{|_{W_4^{{3 \over 2}}(\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.
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 Ω be a domain of IRn with boundary ∂Ω satisfying the cone property. Then, the functional spaceW_p^{2,1}(Q)is continuously embedded in u ∈ Lq(Q) for q satisfying:(i)1 \le q \le {{p(n + 2)} \over {n + 2 - 2p}}
, ifp < {{n + 2} \over 2}
; (ii) 1 ≤ q < ∞, ifp = {{n + 2} \over 2}and(iii)q = ∞, ifp > {{n + 2} \over 2}
. In particular, for such q and any functionu \in W_p^{2,1}(Q)we have that\parallel u{\parallel_{{L^q}(Q)}} \le C\parallel u{\parallel_{W_p^{2,1}(Q)}},with a constant C depending only on Ω, T , p, q, n.
In the cases(ii), (iii)or in(i)when1 \le q < {{p(n + 2)} \over {n + 2 - 2p}}
, the referred embedding is compact.
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 Ω be a bounded domain in ℝn, with a C2boundary ∂Ω, aij be bounded continuous functions in Q, and p > 1. Assume that
{a_{ij}} \in C(\bar Q)
, i, j = 1,..., n; [aij]n×nis a real positive matrix such that for some positive constant β we have
\sum\limits_{i,j = 1}^n{a_{ij}}(x,t){\xi_i}{\xi_j} \ge \beta |\xi {|^2}for all (x, t) ∈ Q and all ξ ∈ Rn,;
f ∈ Lp(Q);
ai ∈ Lr(Q) with either r = max p, n + 2 if p ≠ n + 2 or r = n + 2 + ε, for any ε > 0, if p = n + 2;
a ∈ Ls(Q) with either s = max (p, (n + 2)/2) if p ≠ (n + 2)/2 or s = (n + 2)/2 + ε, for any ε > 0, if p = (n + 2)/2.
{b_i},b \in {C^2}(\bar \Gamma)
, i = 1,..., n, and the coefficients bi(x, t) satisfy the condition\left| {\sum\limits_{i = 1}^n{b_i}(x,t){\eta_i}(x)} \right| \ge \delta> 0for a.e. in ∂Ω × (0, T), where ηi(x) is the ith-component of the unitary outer normal vector to ∂Ω in x ∈ ∂ Ω;
{u_0} \in \;W_p^{2 - {2 \over p}}(\Omega)with p ≠ 3 and satisfying the compatibility condition\sum\limits_{i = 1}^n{b_i}{{\partial {u_0}} \over {\partial {x_i}}} + b\;{u_0} = 0on ∂ Ω when p > 3.
Then, there exists a unique solutionu \in W_p^{2,1}(Q)of Problem (2.1); moreover, there is a positive constant Cp such that the solution satisfies\parallel u{\parallel_{W_p^{2,1}(Q)}} \le {C_p}\left({\parallel f{\parallel_{{L^p}(Q)}} + \parallel {u_0}{\parallel_{W_p^{2 - {2 \over p}}(\Omega)}}} \right).
Such constant Cp depends only on Ω, T, p, r, s, β, δ and on the norms\parallel {b_i}{\parallel_{{C^2}(\bar \Gamma)}}
,
\parallel b{\parallel_{{C^2}(\bar \Gamma)}}
,
\parallel {a_{ij}}{\parallel_{C(\bar Q)}}
, ||ai||Lr (Q)and ||ai||Ls (Q). Moreover, we may assume that the dependencies of Cp 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 u0 = 0 on ∂Ω when p > 3/2. This way, all the results in this paper holds if we replaced the Neumann conditions by Dirichlet conditions.
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:
\left\{{\matrix{{{{\partial \hat N} \over {\partial t}} = {r_N} - {\mu_N}\hat N - {\beta_1}\hat N\hat A - {\alpha_N}{\gamma_N}|\hat D|\hat N,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \hat A} \over {\partial t}} = {r_A}\hat A\left({1 - {{\hat A} \over {{k_A}}}} \right) - ({\mu_A} + {\varepsilon_A})\hat A - {\alpha_A}{\gamma_A}|\hat D|\hat A,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \hat D} \over {\partial t}} = \sigma \Delta \hat D + \mu {\chi_\omega} - \gamma \hat D\hat A - {\gamma_N}\hat D\hat N - \tau \hat D,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \hat D} \over {\partial \eta}}(\cdot) = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{\hat N(\cdot,0) = {N_0}(\cdot),\hat A(\cdot,0) = {A_0}(\cdot),\hat D(\cdot,0) = {D_0}(\cdot),} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
Now we observe that, since the equation for
\hat N
in this last problem is, for each x ∈ Ω, an ordinary differential equation which is linear in
\hat N
, we can find an explicit expression for it in terms of  and
|\hat D|
. However, Â is, for each x ∈ Ω, a nonlinear differential equation inÂ, and we can determine its explicit expression in terms of
|\hat D|
using Bernoulli’s method. Using these observations and setting λ = rA − (μA + εA), we introduce operators Λ : L∞ (Q) → L∞ (Q) and Θ : L∞ (Q) → L∞ (Q), defined respectively by
\Lambda (\phi)(x,t) = {{{A_0}(x){k_A}{e^{\lambda t}}{e^{- {\alpha_A}{\gamma_A}\int_0^t|\phi (\xi,x)|d\xi}}} \over {{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}}
and
\Theta (\phi)(x,t) = {{{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} \over {{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}}}},
where 0 ≤ s ≤ t ≤ T.
Remark 3.1
Thus, (
\hat N
, Â,
\hat D
) is a solution of (3.1) if, and only if,\hat N = \Theta (\hat D)
,
\hat A = \Lambda (\hat D)and\hat Dsatisfies the following integro-differential system:\left\{{\matrix{{{{\partial \hat D} \over {\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,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \hat D} \over {\partial \eta}}(\cdot) = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{\hat D(\cdot,0) = {D_0}(\cdot),} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
Remark 3.2
Notice that, to guarantee that (N, A, D), withD = \hat D
,
N = \Theta (\hat D)andA = \Lambda (\hat D)is also a solution of system (1.1), it is enough to prove that the solution\hat Dof 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_4^{2,1}(Q)of Problem(3.4). Moreover,\hat D \le {\mu\over \tau}a.e. in Q and||\hat D|{|_{W_4^{2,1}(Q)}} \le C,where C is a constant depending on μ, T , ω and||{D_0}|{|_{W_4^{{3 \over 2}}(\Omega)}}
.
Lemma 3.4
Let f : (0, T) → ℝ differentiable such that f (t) > 0 and f′(t) ≥ 0. Ifg(t) = {{\int_0^tf(x)dx} \over {f(t)}}
, then g(t) ≤ T, for all t ∈ (0, T).
Proof
Since f is continuous in (0, T), it follows that
\matrix{{g'(t) = {{f{{(t)}^2} - f'(t)\int_0^tf(x)dx} \over {f{{(t)}^2}}}} \hfill\cr{\,\,\,\,\,\,\,\,\,\,\,\,\, = 1 - {{f'(t)} \over {f(t)}}g(t).} \hfill\cr}
As f (t) > 0 we have g(t) ≥ 0 and using the fact that f′(t) ≥ 0 we obtain
{{f'(t)} \over {f(t)}}g(t) \ge 0
. Therefore, g′(t) ≤ 1, which suggests g(t) ≤ t, for all t ∈ (0, T). Thus, g(t) ≤ T, as intended.
Since in the proof of existence of solutions of (3.4) the expression of Λ and Θ will play important roles, we state some of their properties in the following:
Lemma 3.5
If N0, A0 ∈ L∞ (Ω) and Cλ = max {1, eλT}, then for any ϕ, ϕ1, ϕ2 ∈ L∞ (Q) and for almost every (x, t) ∈ Q, there holds
0 ≤ Θ(ϕ)(x, t) ≤ ||N0||L∞ (Ω) + rNT;
0 ≤ Λ(ϕ)(x, t) ≤ Cλ ||A0||L∞ (Ω);
||Λ(ϕ1) − Λ(ϕ2)||L∞ (Q) ≤ C1||ϕ1 − ϕ2||L∞ (Q), where C1is a constant depending on rA, kA, αA, γA, Cλ, T and ||A0||L∞ (Ω);
We firstly need to observe that, due to the mean value inequality, given any z1, z2 ∈ IR, there is θ = θ (z1, z2) such that ez2 − ez1 = e(1−θ)z1 + θz2(z2 − z1); in particular, for any z1, z2 ≤ 0 we also have (1 − θ)z1 + θz2 ≤ 0 and thus
|{e^{{z_2}}} - {e^{{z_1}}}| \le |{z_2} - {z_1}|,\quad \forall {z_1},{z_2} \le 0.
Secondly, we note that by the inequality (3.5) and by ϕi ∈ L∞ (Q), i = 1, 2, we obtain
\matrix{{|{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}}|} \hfill &\le\hfill & {| - {\alpha_A}{\gamma_A}\int_0^t(|{\phi_1}(x,\xi)| - |{\phi_2}(x,\xi)|)d\xi |} \hfill\cr{} \hfill & {} \hfill & {} \hfill\cr{} \hfill &\le\hfill & {{\alpha_A}{\gamma_A}T||{\phi_1} - {\phi_2}{{||}_{{L^\infty}(Q)}}.} \hfill\cr}
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 Ψ defined as follows:
\matrix{\hfill {\Psi :} & \hfill {[0,1] \times {L^\infty}(Q)} & \hfill\to& \hfill {{L^\infty}(Q)} \cr\hfill {} & \hfill {(l,\phi)} & \hfill\mapsto& \hfill {D,} \cr}
where D is the unique solution of
\left\{{\matrix{{{{\partial D} \over {\partial t}} = \sigma \Delta D + \mu {\chi_\omega} - l\gamma D\Lambda (\phi) - l{\gamma_N}D\Theta (\phi) - \tau D,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial D} \over {\partial \eta}}(\cdot) = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{D(\cdot,0) = {D_0}(\cdot),} \hfill & {{\rm{in}}} \hfill & {\Omega,} \hfill\cr}} \right.
with Λ(ϕ) and Θ(ϕ) given by (3.2) and (3.3), respectively.
To apply such theorem we present next a sequence of lemmas:
Lemma 3.6
Suppose N0, A0 ∈ L∞ (Ω) and{D_0} \in W_4^{{3 \over 2}}(\Omega)
. Then the mapping Ψ : [0, 1] × L∞ (Q) → L∞ (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γΛ(ϕ) −lγNΘ(ϕ) − τ ∈ L4(Q), because by Lemma 3.5, Λ(ϕ), Θ(ϕ) ∈ L∞ (Q). Thus, we conclude that there is a unique solution
D \in W_4^{2,1}(Q)
of problem 3.12. Moreover, D satisfies the following estimate:
\matrix{{||D{{||}_{W_4^{2,1}(Q)}} \le {C_p}(||\mu {\chi_\omega}{{||}_{{L^4}(Q)}} + ||{D_0}{{||}_{W_4^{{3 \over 2}}(\Omega)}})} & {} & {} & {}\cr{\le {C_p}(\mu |\omega {|^{{1 \over 4}}}{T^{{1 \over 4}}} + ||{D_0}{{||}_{W_4^{{3 \over 2}}(\Omega)}}).} & {} & {} & {}\cr}
Finally, from Lemma 2.4, we have
W_4^{2,1}(Q) \hookrightarrow {L^\infty}(Q)
, and we conclude that the operator Ψ in well defined.
Lemma 3.7
Suppose D is a solution of (3.12) and0 \le {D_0} \le {\mu\over \tau}a.e. in Ω, then0 \le D \le {\mu\over \tau}a.e. in Q.
Proof
Multiplying the first equation in (3.12) by D− and integrating into Ω, we get
\matrix{\hfill {{1 \over 2}{d \over {dt}}\int_\Omega {{({D^ -})}^2}\;dx =- \sigma \int_\Omega |\nabla {D^ -}{|^2}\;dx - \mu \int_\omega {D^ -}\;dx} \cr\hfill {- l\gamma \int_\Omega \Lambda (\phi){{({D^ -})}^2}\;dx - l{\gamma_N}\int_\Omega \Theta (\phi){{({D^ -})}^2}\;dx - \tau \int_\Omega {{({D^ -})}^2}\;dx.} \cr}
Thus,
{d \over {dt}}\int_\Omega {({D^ -})^2}dx \le 0,
and using Gronwall’s inequality and the fact that D0 ≥ 0 a.e. in Ω, we obtain
\int_\Omega {({D^ -})^2}dx \le \int_\Omega {({D_0}^ -)^2}dx = 0,
that is, ||D− (·, t)||L2 (Ω) = 0 for all t ∈ (0, T), where we conclude that D− = 0 a.e. in Q and therefore D ≥ 0 a.e. in Q.
Now, we observe that the first equation in (3.12) can be rewritten as
{\partial\over {\partial t}}(D - {\mu\over \tau}) = \sigma \Delta (D - {\mu\over \tau}) - l\gamma \Lambda (\phi)D - l{\gamma_N}\Theta (\phi)D - \tau (D - {{\mu {\chi_\omega}} \over \tau}).
Thus, using Gronwall’s inequality and the fact that
{D_0} \le {\mu\over \tau}
a.e. in Ω, it follows that
\int_\Omega {((D - {\mu\over \tau}{)^ +})^2}\;dx \le \int_\Omega {(({D_0} - {\mu\over \tau}{)^ +})^2}\;dx = 0,
that is,
||(D(\cdot,t) - {\mu\over \tau}{)^ +}{||_{{L^2}(\Omega)}} = 0
for all t ∈ (0, T), and therefore
{(D - {\mu\over \tau})^ +} = 0
a.e. in Q, and we conclude that
D \le {\mu\over \tau}
a.e. in Q.
Lemma 3.8
For each fixed l ∈ [0, 1], the mapping Ψ(l, ·) : L∞ (Q) → L∞ (Q) is compact, i.e., it is continuous and maps bounded sets into relatively compacts sets.
Using the Proposition 2.5 and the fact that L∞ (Q) ↪ L4(Q) and
{D_1} \le {\mu\over \tau}
, we get
\matrix{\hfill {||\tilde D|{|_{W_4^{2,1}(Q)}} \le {C_p}|| - l\gamma {D_1}(\Lambda ({\phi_1}) - \Lambda ({\phi_2})) - l{\gamma_N}{D_1}(\Theta ({D_1}) - \Theta ({D_2}))|{|_{{L^4}(Q)}}} \cr\hfill {\le {{\bar C}_p}|| - l\gamma {D_1}(\Lambda ({\phi_1}) - \Lambda ({\phi_2})) - l{\gamma_N}{D_1}(\Theta ({\phi_1}) - \Theta ({\phi_2}))|{|_{{L^\infty}(Q)}}} \cr\hfill {\le {{\bar C}_p}\gamma {\mu\over \tau}||\Lambda ({\phi_1}) - \Lambda ({\phi_2})|{|_{{L^\infty}(Q)}} + {{\bar C}_p}{\gamma_N}{\mu\over \tau}||\Theta ({\phi_1}) - \Theta ({\phi_2})|{|_{{L^\infty}(Q)}}.} \cr}
Then, by Lemmas 3.5 and 2.4, we finally have
||\Psi (l,{\phi_1}) - \Psi (l,{\phi_2}{)||_{{L^\infty}(Q)}} \le C||{\phi_1} - {\phi_2}{||_{{L^\infty}(Q)}},
where C depends on
{\bar C_p}
, C1, C2, γ, γN, μ, τ and the immersion constant.
To show that Ψ(l, ·) is compact, we use the fact that the immersion
W_4^{2,1}(Q) \hookrightarrow {L^\infty}(Q)
is compact and that Ψ(l, ·) is the composition between the inclusion operator and the solution operator, i.e.,
\Psi (l, \cdot):{L^\infty}(Q) \to W_4^{2,1}(Q) \to {L^\infty}(Q)
.
Lemma 3.9
Given a bounded subset B ⊂ L∞ (Q), for each ϕ ∈ B, the mapping Ψ(·, ϕ) : [0, 1] → L∞ (Q) is uniformly continuous with respect to B.
Proof
Since B ∈ L∞ (Q) is bounded, there is rB ≥ 0 such that, for any ϕ ∈ B, we have ||ϕ||L∞ (Q) ≤ rB. Now, let us fix ϕ ∈ L∞ (Q) and consider l1, l2 ∈ [0, 1] and denote Ψ(l1, ϕ) = D1, Ψ(l2, ϕ) = D2 and
\tilde D = {D_1} - {D_2}
. Then,
\tilde D
satisfies
\left\{{\matrix{{{{\partial \tilde D} \over {\partial t}} - \sigma \Delta \tilde D + \gamma {l_2}\Lambda (\phi)\tilde D + {\gamma_N}{l_2}\Theta (\phi)\tilde D + \tau \tilde D =} \hfill & {} \hfill & {} \hfill\cr{\gamma \Lambda (\phi){D_1}({l_1} - {l_2}) - {\gamma_N}\Theta (\phi){D_1}({l_1} - {l_2}),} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \tilde D} \over {\partial \eta}} = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{\tilde D(\cdot,0) = {{\tilde D}_0}(\cdot) = 0,} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
Using the Proposition 2.5 and the fact that L∞ (Q) ↪ L4(Q),
{D_1} \le {\mu\over \tau}
, we get
\matrix{\hfill {||\tilde D{{||}_{W_4^{2,1}(Q)}} \le {C_p}||\gamma \Lambda (\phi){D_1}({l_1} - {l_2}) - {\gamma_N}\Theta (\phi){D_1}({l_1} - {l_2}{{)||}_{{L^4}(Q)}}} \cr\hfill {\le {{\bar C}_p}\gamma {\mu\over \tau}|{l_1} - {l_2}|||\Lambda (\phi {{)||}_{{L^\infty}(Q)}} + {{\bar C}_p}{\gamma_N}{\mu\over \tau}|{l_1} - {l_2}|||\Theta (\phi {{)||}_{{L^\infty}(Q)}}.} \cr}
Then, by Lemmas 3.5 and 2.4, we finally have
||\Psi ({l_1},\phi) - \Psi ({l_2},\phi {)||_{{L^\infty}(Q)}} \le C|{l_1} - {l_2}|,
where C depends on
{\bar C_p}
, γ, γN, μ, τ, rN, T, Cλ, ||N0||L∞ (Ω), ||A0||L∞ (Ω) and the immersion constant.
Lemma 3.10
Suppose{D_0} \le {\mu\over \tau}a.e. in Ω, then there exists a number ρ > 0 such that, for any l ∈ [0, 1] and any possible fixed point D ∈ L∞ (Q) of Ψ(l, ·), there holds ||D||L∞ (Q) < ρ.
Proof
Let D ∈ L∞ (Q) such that Ψ(l, D) = D. The analogous demonstration made in Proposition 3.7 guarantees us
||D{||_{{L^\infty}(Q)}} \le {\mu\over \tau}
. Therefore, just take
\rho= {\mu\over \tau} + 1
.
Lemma 3.11
The mapping Ψ(0, ·) : L∞ (Q) → L∞ (Q) has a unique fixed point.
Proof
Indeed, letting l = 0 in 3.12, D is a fixed point of Ψ(0, ·) if, and only if, D is the unique solution to the problem
\left\{{\matrix{{{{\partial D} \over {\partial t}} = \sigma \Delta D + \mu {\chi_\omega} - \tau D,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial D} \over {\partial \eta}}(\cdot) = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{D(\cdot,0) = {D_0}(\cdot),} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
But Proposition 2.5 guarantees the existence of a unique solution
D \in W_4^{2,1}(Q) \hookrightarrow {L^\infty}(Q)
of this last problem; therefore Ψ(0, ·) has a unique fixed point in L∞ (Q).
Proposition 3.12
There is a nonnegative solution\hat D \in W_4^{2,1}(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 Ψ : [0, 1] × L∞ (Q) → L∞ (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_4^{2,1}(Q)
is nonnegative and
\hat D
is the required solution of (3.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_4^{2,1}(Q)of the modified problem (3.1).
Proof
Just combine the Proposition 3.12, the Remark 3.1 and the Lemma 3.5.
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:\matrix{{|{{\partial \hat N} \over {\partial t}}| \le {r_N} + {\mu_N}(||{N_0}{{||}_{{L^\infty}(\Omega)}} + {r_N}T) + {\beta_1}(||{N_0}{{||}_{{L^\infty}(\Omega)}} + {r_N}T){C_\lambda}||{A_0}{{||}_{{L^\infty}(\Omega)}}}\cr{+ {\alpha_N}{\gamma_N}{\mu\over \tau}(||{N_0}{{||}_{{L^\infty}(\Omega)}} + {r_N}T),}\cr}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:\matrix{{|{{\partial \hat A} \over {\partial t}}| \le {r_A}{C_\lambda}||{A_0}{{||}_{{L^\infty}(\Omega)}} + {{{r_A}} \over {{k_A}}}{{({C_\lambda}||{A_0}{{||}_{{L^\infty}(\Omega)}})}^2} + ({\mu_A} + {\varepsilon_A}){C_\lambda}||{A_0}{{||}_{{L^\infty}(\Omega)}}}\cr{+ {\alpha_A}{\gamma_A}{\mu\over \tau}{C_\lambda}||{A_0}{{||}_{{L^\infty}(\Omega)}},}\cr}a.e. in Q, i.e.,Ât ∈ L∞ (Q).
Proposition 4.3
There is a nonnegative solution(N,A,D) \in W \times W \times W_4^{2,1}(Q)of problem (1.1).
Proof
Just combine the Proposition 4.1 and the Remarks 3.2 and 4.2.
Proposition 4.4
The solution (N, A, D) of the problem (1.1) is unique.
Proof
Let (N1, A1, D1) and (N2, A2, D2) be solutions to the problem (1.1); if Ñ = N1 − N2, Ã = A1 − A2 and
\tilde D = {D_1} - {D_2}
, then Ñ, Ã and
\tilde H
satisfy the following problems, respectively:
\left\{{\matrix{{{{\partial \tilde N} \over {\partial t}} =- {\mu_N}\tilde N - {\beta_1}{A_1}\tilde N - {\beta_1}{N_2}\tilde A - {\alpha_N}{\gamma_N}{N_1}\tilde D - {\alpha_N}{\gamma_N}{D_2}\tilde N,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{\tilde N(\cdot,0) = {{\tilde N}_0}(\cdot) = 0,} \hfill & {{\rm{in}}} \hfill & {\Omega,} \hfill\cr}} \right.\left\{{\matrix{{{{\partial \tilde A} \over {\partial t}} = {r_A}\tilde A - {{{r_A}} \over {{k_A}}}({A_1} + {A_2})\tilde A - ({\mu_A} + {\varepsilon_A})\tilde A - {\alpha_A}{\gamma_A}{A_1}\tilde D - {\alpha_A}{\gamma_A}{D_2}\tilde A,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{\tilde A(\cdot,0) = {{\tilde A}_0}(\cdot) = 0,} \hfill & {{\rm{in}}} \hfill & {\Omega,} \hfill\cr}} \right.\left\{{\matrix{{{{\partial \tilde D} \over {\partial t}} = \sigma \Delta \tilde D - \gamma {A_1}\tilde D - \gamma {D_2}\tilde A - {\gamma_N}{N_1}\tilde D - {\gamma_N}{D_2}\tilde N - \tau \tilde D,} \hfill & {{\rm{in}}} \hfill & {Q,} \hfill\cr{{{\partial \tilde D} \over {\partial \eta}}(\cdot) = 0,} \hfill & {{\rm{on}}} \hfill & {\Gamma,} \hfill\cr{\tilde D(\cdot,0) = {{\tilde D}_0}(\cdot) = 0,} \hfill & {{\rm{in}}} \hfill & {\Omega.} \hfill\cr}} \right.
Multiplying the first equation of (4.3) byÑ, integrating into Ω, using the fact that N1 ≤ ||N0||L∞ (Ω) + rNT and the inequality of Young, we have
\matrix{{{1 \over 2}{d \over {dt}}\int_\Omega {{\tilde N}^2}dx} \hfill & {=- {\mu_N}\int_\Omega {{\tilde N}^2}dx - {\beta_1}\int_\Omega {N_2}\tilde A\tilde Ndx - {\alpha_N}{\gamma_N}\int_\Omega {N_1}\tilde D\tilde Ndx} \hfill\cr{} \hfill & {- {\alpha_N}{\gamma_N}\int_\Omega {D_2}{{\tilde N}^2}dx} \hfill\cr{} \hfill & {\le (||{N_0}|{|_{{L^\infty}(\Omega)}} + {r_N}T)({\beta_1}\int_\Omega |\tilde A||\tilde N|dx + {\alpha_N}{\gamma_N}\int_\Omega |\tilde D||\tilde N|dx)} \hfill\cr{} \hfill & {\le C\int_\Omega ({{\tilde A}^2} + {{\tilde N}^2} + {{\tilde H}^2})dx,} \hfill\cr}
where C depends on β1, αN, γN, rN, T and ||N0||L∞ (Ω).
Now, multiplying the first equation of (4.4) byÃ, integrating into Ω, using the fact that A1 ≤ Cλ ||A0||L∞ (Ω) and the inequality of Young, we obtain
\matrix{{{1 \over 2}{d \over {dt}}\int_\Omega {{\tilde A}^2}dx} \hfill & {= {r_A}\int_\Omega {{\tilde A}^2}dx - {{{r_A}} \over {{k_A}}}\int_\Omega ({A_1} + {A_2}){{\tilde A}^2}dx - ({\mu_A} + {\varepsilon_A})\int_\Omega {{\tilde A}^2}dx} \hfill\cr{} \hfill & {- {\alpha_A}{\gamma_A}\int_\Omega {A_1}\tilde D\tilde Adx - {\alpha_A}{\gamma_A}\int_\Omega {D_2}{{\tilde A}^2}dx} \hfill\cr{} \hfill & {\le {r_A}\int_\Omega |\tilde A{|^2}dx + {\alpha_A}{\gamma_A}{C_\lambda}||{A_0}|{|_{{L^\infty}(\Omega)}}\int_\Omega |\tilde D||\tilde A|dx} \hfill\cr{} \hfill & {\le C\int_\Omega ({{\tilde A}^2} + {{\tilde N}^2} + {{\tilde H}^2})dx,} \hfill\cr}
where C depends on rA, αA, γA Cλ and ||A0||L∞ (Ω).
Lastly, multiplying the first equation of (4.5) by
\tilde H
, integrating into Ω, using the fact that
D \le {\mu\over \tau}
and the inequality of Young, we obtain
\matrix{{{1 \over 2}{d \over {dt}}\int_\Omega {{\tilde D}^2}dx} \hfill & {=- \sigma \int_\Omega |\nabla \tilde D{|^2}dx - \gamma \int_\Omega {A_1}{{\tilde D}^2}dx - \gamma \int_\Omega {D_2}\tilde A\tilde Ddx} \hfill\cr{} \hfill & {- {\gamma_N}\int_\Omega {N_1}{{\tilde D}^2}dx - {\gamma_N}\int_\Omega {D_2}\tilde N\tilde Ddx - \tau \int_\Omega {{\tilde D}^2}dx} \hfill\cr{} \hfill & {\le \gamma {\mu\over \tau}\int_\Omega |\tilde A||\tilde D|dx + {\gamma_N}{\mu\over \tau}\int_\Omega |\tilde N||\tilde D|dx} \hfill\cr{} \hfill & {\le C\int_\Omega ({{\tilde A}^2} + {{\tilde N}^2} + {{\tilde H}^2})dx,} \hfill\cr}
where C depends on γ, γN, μ and τ.
Thus,
{d \over {dt}}(\int_\Omega (|\tilde N{|^2} + |\tilde A{|^2} + |\tilde D{|^2})dx) \le C\int_\Omega (|\tilde N{|^2} + |\tilde A{|^2} + |\tilde D{|^2})dx,
and using the Gronwall’s inequality, we finally
\int_\Omega (|\tilde N{|^2} + |\tilde A{|^2} + |\tilde D{|^2})dx \le {e^{CT}}\int_\Omega (|{\tilde N_0}{|^2} + |{\tilde A_0}{|^2} + |{\tilde D_0}{|^2})dx = 0,
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 ∈ (0, T). Where we concludeÑ =Ã =
\tilde D
= 0 a.e. in Q and therefore N1 = N2, A1 = A2 and D1 = D2 a.e. in Q.
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 Ω = [0, L] × [0, L], with L = 1, discretized with n = 50 steps Δx = Δy = L/n = 0.02. The Laplacian Δ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) = N2, A(x, 0) = A2, D(x, 0) = 0, where (N2, A2, 0) is a globally asymptotically stable equilibrium point for the ODE system (1.2) without treatment (ν = 0). The expressions for N2 and A2 are:
{N_2} = {{{r_N}} \over {{\mu_N} + {\beta_1}{A_2}}},\;{\;_2} = {{{r_A} - {\mu_A} - {\varepsilon_A}} \over {{r_A}}}{K_A}.
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 ← N/(rN/μN) and A ← A/KA. Therefore, the population sizes range from 0 to 1. The re-scaled parameter values used in the model simulations were fixed to
\matrix{{{r_N} = 1,{\mu_N} = 1{,_A} = 1{,_A} = 1,{\beta_1} = 1.5,{\mu_A} = 0.05,{\varepsilon_A} = 0.05,}\cr{{\tau_H} = 0.9,{\gamma_N} = 0.1,{\alpha_N} = 1,{\gamma_A} = 1.}\cr}
These values were chosen to describe: normal cells that reach the equilibrium N = rN/μN = 1 at absence of tumor cells; a tumor with the same carrying capacity of normal cells (KA = rN/μN = 1) and a greater absorption of the chemotherapeutic drug by tumor cells in comparison with normal cells (γA > γ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 αA, the diffusion coefficient of the chemotherapeutic drug σ and the chemotherapy infusion rate μ. 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 ω, which is a mathematical description of a blood vessel crossing the tissue, from where the chemotherapy enters the tissue. The values for parameters αA, σ, μ and the position of σ used in each simulation are indicated in Table 1. We present the following results.
Set-up of different simulations an their biological outcomes. Each row indicates the numerical values used for the chemotherapeutic parameters αA (cytotoxicity), σ (diffusion coefficient), μ (infusion rate), and the position of ω ⊂ Ω ⊂ ℝ2. Simulation 1 was performed in a two-dimensional domain Ω = [0, 1] × [0, 1], while simulations 2–5 were performed in a one-dimensional domain Ω = [0, 1].
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., ω = [0.45, 0.55] × [0.45, 0.55]. We use the following parameter values: αA = 5, μ = 3, and σ = 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).
Fig. 1
Results of Simulation 1, for model (1.1) within a two-dimensional domain Ω = [0, L] × [0, L] = [0, 1] × [0, 1]. Plots of model solutions N(x, y, t) (normal cells, blue, top row), A(x, y, t) (cancer cells, red, middle row) and D(x, y, t) (chemotherapeutic drug concentration, green, bottom row) at time points t = 0, 1, 15 (columns 1,2 and 3, respectively). See Table 1 for parameter values used here. At time t = 0, the tumor is spread trough the tissue, and as chemotherapy is applied (t > 0), the tumor cells are reduced in the vicinity of the blood vessel, while the distant tumor cells persist along time (the shape of the solution at time t = 15 is stationary). Within the vicinity of the blood vessel, the removal of tumor cells allows the normal tissue to recover and grow.
In order to make easier to illustrate the model dynamics, we present the results of next simulations in a one-dimensional domain Ω = [0, 1]. In Simulation 2, we use the same parameters values used in Simulation 1 (see Table 1), but increase the chemotherapy toxicity αA and move the blood vessel to the left side of the tissue, ω = [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 μ (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 σ, 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 α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 (α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.
Fig. 2
Results of Simulation 2, with Ω = [0, L] = [0, 1]. Plots of model solutions A(x, t) (cancer cells, red), N(x, t) (normal cells, blue) and D(x, t) (chemotherapeutic drug concentration, green) at time points t = 0, 3, 6, 9, 12, 15. See Table 1 for parameter values used here. At time t = 0, the tumor is spread trough the tissue, and as chemotherapy is applied (t > 0), the tumor cells are reduced and extinct within a given distance from the blood vessel (x < 0.6), but not in the entire tissue (x > 0.6). Within the region of tumor extinction, the removal of tumor cells release the normal tissue to recover and grow.
Fig. 3
Results of Simulation 3, with Ω = [0, L] = [0, 1]. Plots of model solutions A(x, t) (cancer cells, red), N(x, t) (normal cells, blue) and D(x, t) (chemotherapeutic drug concentration, green) at time points t = 0, 3, 6, 9, 12, 15. See Table 1 for parameter values used here. At time t = 0, the tumor is spread trough the tissue, and as chemotherapy is applied (t > 0), the tumor cells are reduced and in the entire tissue. In comparison with Simulation 2, the tumor extinction is reached because the drug infusion rate μ is increased here. Within the entire tissue, the removal of tumor cells release the normal tissue to recover and grow.
Fig. 4
Results of Simulation 4, with Ω = [0, L] = [0, 1]. Plots of model solutions A(x, t) (cancer cells, red), N(x, t) (normal cells, blue) and D(x, t) (chemotherapeutic drug concentration, green) at time points t = 0, 3, 6, 9, 12, 25. See Table 1 for parameter values used here. At time t = 0, the tumor is spread trough the tissue, and as chemotherapy is applied (t > 0), the tumor cells are reduced and in the entire tissue. In comparison with Simulation 2, the tumor extinction is reached because the drug diffusion coefficient σ is increased here.
Fig. 5
Results of Simulation 5, with Ω = [0, L] = [0, 1]. Plots of model solutions A(x, t) (cancer cells, red), N(x, t) (normal cells, blue) and D(x, t) (chemotherapeutic drug concentration, green) at time points t = 0, 3, 6, 9, 12, 15. See Table 1 for parameter values used here. At time t = 0, the tumor is spread trough the tissue, and as chemotherapy is applied (t > 0), the tumor cells are reduced and in the entire tissue. In comparison with Simulation 2, the tumor extinction is reached because the chemotherapy toxicity against tumor cells, αA, was increased.
