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

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.


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].
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 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 µ 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 (γ N N D) and cancerous cells (γ A AD) 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.
Next, we state some hypotheses that will be assumed throughout this article.

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 IR n with boundary ∂Ω satisfying the cone property. Then, the functional space W 2,1 p (Q) is continuously embedded in u ∈ L q (Q) for q satisfying: 2 . In particular, for such q and any function u ∈ W 2,1 p (Q) we have that u L q (Q) ≤ C u W 2,1 p (Q) , with a constant C depending only on Ω, T , p, q, n.
Next, we consider the following simple parabolic initial-boundary value problem: in Ω.
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.
. . , n, and the coefficients b i (x, t) satisfy the condition (Ω) with p = 3 and satisfying the compatibility condition Then, there exists a unique solution u ∈ W 2,1 p (Q) of Problem (2.1); moreover, there is a positive constant C p such that the solution satisfies . Such constant C p depends only on Ω, T , p, r, s, β, δ and on the norms b i C 2 (Γ) , b C 2 (Γ) , a ij C(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 ∂Ω 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: in Ω.
Now we observe that, since the equation forN in this last problem is, for each x ∈ Ω, an ordinary differential equation which is linear inN , we can find an explicit expression for it in terms ofÂ and |D|. However,Â is, for each x ∈ Ω, a nonlinear differential equation inÂ, and we can determine its explicit expression in terms of |D| using Bernoulli's method. Using these observations and setting andD satisfies the following integro-differential system: in Ω.
where C is a constant depending on µ, T , ω and ||D 0 || .
Proof: Since f is continuous in (0, T ), it follows that As f (t) > 0 we have g(t) ≥ 0 and using the fact that f (t) ≥ 0 we obtain 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:

Finally, the expression in (3.3) suggests
. and using the estimates obtained in (3.8), (3.9) and (3.10) and making the possible simplifications, we obtain for almost everything (x, t) ∈ Q, i.e.,

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 Ψ defined as follows: where D is the unique solution of (3.12) 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 Θ(φ) − τ ∈ L 4 (Q), because by Lemma 3.5, Λ(φ), Θ(φ) ∈ L ∞ (Q). Thus, we conclude that there is a unique solution D ∈ W 2,1 4 (Q) of problem 3.12. Moreover, D satisfies the following estimate: .
Using the Proposition 2.5 and the fact that L ∞ (Q) → L 4 (Q) and D 1 ≤ µ τ , we get Then, by Lemmas 3.5 and 2.4, we finally have where C depends onC p , C 1 , C 2 , γ, γ N , µ, τ and the immersion constant.

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) = 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 (ν = 0). The expressions for N 2 and A 2 are: 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/(r N /µ N ) and A ← A/K A . Therefore, the population sizes range from 0 to 1. The re-scaled parameter values used in the model simulations were fixed to These values were chosen to describe: normal cells that reach the equilibrium N = r N /µ N = 1 at absence of tumor cells; a tumor with the same carrying capacity of normal cells (K A = r N /µ 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.
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).
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  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.
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.  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.  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.   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.  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.