MATHEMATICAL ANALYSIS AND NUMERICAL SIMULATION OF A TUMOR-HOST MODEL WITH CHEMOTHERAPY APPLICATION

. In this paper, a system of non-linear quasi-parabolic partial differential system, modeling the chemotherapy application of spatial tumor-host interaction is considered. At some certain parameters, we derive the steady state of the anti-angiogenic therapy, baseline therapy and anti-cytotoxic therapy models as well as their local stability condition. We use the method of upper and lower solutions to show that the steady states are globally stable. Since the system of non-linear quasi-parabolic partial differential cannot be solved analytically, we formulate a robust numerical scheme based on the semi-ﬁtted ﬁnite difference operator. Analysis of the basic properties of the method shows that it is consistent, stable and convergent. Our numerical results are in agreement with our theoretical ﬁndings.


Introduction
A tumor is scientifically referred to as a neoplasm, or an abnormal tissue area that is either fluid-filled or solid in appearance.A tumor does not necessarily mean cancer, it can be classified into the benign type (which means non-progressive or cannot metastasize, for example, uterine fibroids and moles), pre-malignant (or pre-cancerous growth which include: actinic keratosis, dysplasia of the cervix, metaplasia of the lung and leukoplakia), or malignant (cancerous), and the malignant tumor which is a cancerous case.Record have shown that there are various types of tumors, which are made up of specific types of cancer cells; these include the carcinoma, sarcoma, lymphoma or leukemia, blastoma and germ cell tumor [13].
In practice, it requires an expert to differentiate between an ordinary tumor and a cancerous type.Its treatment can be classified into stages I-IV, the first two stages is called a low-grade tumor which can be treated by watchful monitoring or surgery.The high-grade class (III, IV) are malignant and can spread quickly from the affecting area to other surrounding tissues.Its treatment include the use of radiation therapy, chemotherapy, targeted therapy and tumor treating fields.According to statistics conducted on patients with the tumor related issues, chemotherapy approach forms the integral part of tumor treatment, because it requires administering drugs that can destroy the affected cancer cells by impeding their growth and reproduction.Chemotherapy delivery methods include: injection, orally (by mouth as a pill or liquid), intravenously (through infusion into a vein), topically (as a cream on the skin), direct placement (through a lumbar puncture or device placed under the scalp) [7].
In this paper, we consider the models derived by Hinow et al.,in [6,.The models consider cancer as a complex multiscale disease.Thus, to simplify their work toward finding an efficient healing of cancer, Hinow et al. derived a model for scenario of tumor growth only limited by intrinsic constraints such as oxygen supply and the surrounding stromal tissue, which they referred to as the baseline model.In order to cure such a complex multi-scale disease, Hinow et al., modified the baseline model to a cytostatic model in the sense that the drugs used are not toxic to the cells, but instead inhibit some mechanism essential for cell division or a specific function, which they referred to as an anti-angiogenic chemotherapy model and a drug that affects cells in the proliferative state which they referred to as cytotoxic chemotherapy model.Their findings are that when no treatment is applied to the model, their model reproduces a typical dynamics of early tumor growth and initially avascular tumor reaches a diffusion limited size of the order of millimeters which initiates angiogenesis through the release of vascular endothelial growth factor (VEGF) secreted by hypoxic cells in the core of the tumor.They also stated that the VEGF stimulates endothelial cells to migrate towards the tumor and establishes a nutrient supply sufficient for sustained invasion.However, when they applied cytostatic treatment in the form of a VEGF-inhibitor, they found out that the treatment has the capability to reduce tumor mass.Moreover, they were able to determine that inhibition of endothelial cell proliferation as the more important of the two cellular functions targeted by the drug.On the other hand, considering the application of a cytotoxic drug as a diffusible substance entering the tissue from the blood vessels, they concluded that the drug can either reduce the tumor mass significantly or in fact accelerate the growth rate of the tumor.All these findings are based on experiments without mathematical analysis.Thus, to this end, we believe that mathematical analysis of their model can elaborate more clearly to the audience in this field.Thus, in the next paragraphs, we highlight their models.
Let t, x to denote time, space, w, n, h, a, m, f , g, c denote oxygen concentration, normoxic cells, hypoxic cells, apoptic cells, endothelial cells, extracellular matrix, angiogenic factor, concentration of a drug, and the density of all types of cells and matrix combined is For convenient of explaining the models, we group the dependent variables and diffusion terms into vector u = [w, n, h, a, m, f , g, c] T and the diffusion coefficients D = [D w , D n max{n−v c , 0}+ D m , 0, 0, D m , 0, D g , D c ] T respectively, then the cytotoxic chemotherapy (drug that affects cells in the proliferative state) and anti-angiogenic chemotherapy (also considered as cytostatic in the sense that the drugs used are not toxic to the cells, but instead inhibit some mechanism essential for cell division or a specific function) models in [6, page 525-527] was derived through the baseline model (scenario when tumor growth only limited by intrinsic constraints such as oxygen supply and the surrounding stromal tissue) which we rewrite in a vector form as as the initial conditions, F(u) = F(w, n, h, a, m, f , g, c) and the boundary conditions for the reaction-diffusion system in equation ( 1) are the no flux-boundary conditions.Thus, for u 1 in equation ( 1) we have where, D w , α w , w max , β w , γ w denote the coefficient of oxygen diffusion, rate of diffusability of oxygen, maximum oxygen density, uptake rate of oxygen by normoxic, hypoxic, endothelial cells and loss rate of oxygen.The reason for the choice of the source term α w m(w max − w) is that, at high environmental levels of oxygen, less oxygen is released through the vessel walls.
For u 2 in equation ( 1) we have where, v c , D n max{n−v c , 0}n x , D m , χ n , α n , v max , w h , α h , H denote a threshold which adds to the dispersion of these cells through crowding-driven motion, nonlinear diffusion, random motility, haptotactic movement, rate of logistic growth, maximal density, critical value, rate of gain or loss, Heaviside function which is 1 for positive arguments and zero otherwise.We also see that the growth of normoxic cells levels off in regions where the sum of all cells and matrix approaches the maximal density v max .In the regions where the concentration of oxygen drops below a certain critical value w h , normoxic cells enter the hypoxic class at a rate α h and this transition process is reversible, and the reverse transition is denoted by the the reduced rate α h /10.
For u 3 in equation ( 1) we have The first and the second terms are in equation ( 14) are dictated by conservation of mass and correspond to terms in equation (13).The third term in equation ( 14) denotes the transition of hypoxic cells to apoptotic cells at rate β h as the level of oxygen falls below a second threshold, which is w a < w h .In [6, page 526] is reported that hypoxic cells are less active in general due to reduced availability of oxygen and other nutrients and thus, they assumed that lack of energy causes them to be immobile.
For u 4 in equation ( 1) we have The first term denotes the transition from hypoxic class to apoptic class.
For u 5 in equation ( 1) we have where, D m denotes the random motility for endothelial cells and endothelial cells respond via chemotaxis to gradients of angiogenic factor and this require the presence of angiogenic factor for proliferation.Proliferation is capped by the total density of cells.The proliferation constant for endothelial cells is α m .
For u 6 in equation (1) we have where the tissue matrix is degraded by the tumor cells according to a mass-action law with rate constant β f .For u 7 in equation (1) we have where, D g denotes diffusion coefficient, also see that angiogenic factor is produced by hypoxic cells (alone) at rate α g and consumed by endothelial cells with a mass-action coefficient of β g .
To ensure that they model is universal, Hinow et al., dimensionelise the baseline model (see in [6, page 538-539]) and the dimensionless baseline model is When in equation ( 9), the baseline model becomes the anti-angiogenic chemotherapy model.
However, when a drug that affects cells in the proliferative state is introduced then the normoxic cells are driven into apoptosis at a rate proportional to the drug concentration, then the baseline becomes the cytotoxic chemotherapy model.This implies that the reaction function equation in (13) becomes with rate constant γ > 0, whereas, in order to balance the loss of normoxic cells, the reaction function for apoptotic cells in equation ( 15) is amended accordingly to and whenever, the concentration of the drug denoted by c ∈ [0, 1] which is assumed to be delivered through blood infusion and thus enters the tissue from the blood stream.The production rate of the drug is therefore proportional to the density of endothelial cells, which implies that regions with higher vascular density will experience a higher concentration of the drug.Hence it is described by a reaction-diffusion equation similar to that of oxygen as The drug diffuses at a constant rate D c in the tissue and decays at a constant rate c , it also affects the normoxic cells in such a way that it is consumed at a rate kγ n when it kills normoxic cells.
Instead of looking at a constant drug supply Hinow et al.,in [6, page 532] included a scheduling of the drug by letting the production rate α c (t) > 0 be a time-dependent function, corresponding to the drug being delivered to the tissue.
It is well known fact that many research work have been done with regard to solve dynamic models either analytically or numerically.For instance Cristini et al., [4] developed a thermodynamically consistent mixture model for avascular solid tumor growth which takes into account the effects of cell-to-cell adhesion, and taxis inducing chemical and molecular species.They showed that their model is well-posed and the governing equations are of Cahn-Hilliard type.
To solve the governing equations, they developed a numerical algorithm based on an adaptive Cartesian block-structured mesh refinement scheme and found out that the their numerical solutions presented that taxis may play a role in tumor invasion, because when nutrient plays the role of a chemoattractant, they found out that diffusional instability is exacerbated by nutrient gradients.Accordingly, they also believed that their model is capable of describing complex invasive patterns observed in experiments.
However, based on biological and ecological fields, Tuan et al.,in [20] studied the identification of the population density of a logistic equation backwards in time associated with non-local diffusion and nonlinear reaction.Since their original model is ill-posed, they applied the quasireversibility method to convert their model to a stable approximation model.Their numerical results presented the convergence and stability of the regularized solution for the backward parabolic problem with non-local diffusion and nonlinear reaction term.
However, some of the solution to a model is considered as a cure to a tumor disease, because they are compared to experimental data in most cases, see in [3,23].Due to limited availability of the literature, we are only able to mention a few.
Since, Hinow et al., [6] presented their results without presenting their stability analysis of the models in their work in [6].Thus, we believe that they were not able to explain their results by referring to direct mathematical analysis.This, we believe has undermined the qualitative features of the underlying mathematical models at hand.Therefore in this paper, our aim is to present our mathematical analysis of the models on their merits and establish the stability condition for each model.We then use stability conditions to present our results.Due to the non-linearities of terms in the models, we construct a reliable numerical method capable of capturing the qualitative features of these models.
Thus, it is a well known fact that solving the problems like the one in equation ( 1) with the standard finite difference methods (SFDMs) has limitations.That is, explicit methods finite difference methods (EFDMs), can solve such differential equations with low computational cost, with very small stability regions, which in turn implies severe restrictions on the meshes sizes, which are required in order to achieve the desired results.On the other hand, the implicit finite difference methods (IFDMs) do solve such differential equations with wider stability regions as compared to the EFDMS.However, their associated computational complexity can achieve only one order as compared to EFDMs that use the same number of stages [2].
The rest of the paper is organized as follows.Mathematical analysis of the main model is presented in Section 2. A robust numerical scheme based on the semi-fitted finite difference technique is formulated in Section 3, analysis of the basic properties of this scheme is also examined for convergence.To justify the effectiveness of the proposed schemes, we present some numerical results in Section 4. Section 5 concludes the paper.

Mathematical analysis of the models
In this work, our major attention is on the baseline, anti-angiogenic and cytotoxic chemotherapy models.We shall adapt linear stability analysis method to discuss the general dynamics of each model.At the steady state, the cytotoxic chemotherapy model in equation ( 1), becomes where we see that the steady states of the baseline and anti-angiogenic models are the same.
Solving the system of nonlinear equations in equation ( 13), we obtain We see the oxygen and endothelial cells staedy states (w * and m * ) are positive for the baseline, anti-angiogenic and cytotoxic chemotherapy models.Consequently, the hypoxic, apoptic cells' , angiogenic factor (VEGF) and the extra-cellular matrix steady states are zero.On the other hand, when a host is tumor free, we find that which implies that the extra-cellular matrix ( f ) in all three models is well and healthy.We also see that the oxygen cells' steady state is zero, rendering the tumor free case unworthy to investigate.
Therefore, for the baseline, anti-angiogenic and cytotoxic chemotherapy models, we have the uniform steady states as Stability of the steady states in equation ( 16) can be investigated from the corresponding homogeneous baseline, anti-angiogenic and cytotoxic chemotherapy models.Thus, extracting the Jacobian matrices for all three models, from the system in equation ( 1), we obtain the non-zero entries of the Jacobian matrices for all the three models evaluated at the critical points as which clearly reduces to the characteristic equation of the form where The steady states in equation ( 16) are asymptotically stable if and only if the trace and the determinant in equation ( 19) is negative and postive respectively, ( [8]).This implies that tr(J 2×2 ) < 0 gives and det(J 2×2 ) > 0 gives Thus, for the baseline model, we have from equations (19)(20) that which implies that the product of the logistic growth of the normoxic cells with the maximal density is bounded by the rates of the transitions of normoxic cells into hypoxic cells (and vice versa), hypoxic cells into apoptic cells (vice-versa).
Hence, the baseline model stability makes the anti-angiogenic chemotherapy model uncoditional stable.Consequently, the cytotoxic chemotherapy model asymptotic stability is which implies that the product of the logistic growth of the normoxic cells with the maximal density is bounded by the sum of death rate of the normoxic cells (due to the presence of the drug), rates of the transitions of normoxic cells into hypoxic cells (and vice versa), hypoxic cells into apoptic cells (vice-versa).Remark 2.1 The stability conditions present the following facts: (a) heavyweight functions should not be all zero for the baseline model to exists.
(b) the stability conditions mainly only affects maximal density, normoxic, hypoxic cells and the cytotoxic drug.
(c) the apoptic cells and agiogenic factor are not affected by the stability conditions.

Global stability of the steady states
In this segment, we mainly prove that the equilibrium points in equation ( 16) are globally asymptotically stable with the upper and lower solution method in [14,15].Denoting the reaction functions in equation ( 1) by l j (w, n, h, a, m, f , g, c) for j = 1, 2, equation (1) we let and let U ⊂ R 8 + such that U = {u ∈ R 8 + : u ≤ 0 ≤ ū} and K j be any positive constant satisfying then we have the following results (see [22]).
Proof.From the maximum principle of parabolic equations, it is known that for any initial value (w 0 (t, x), n 0 (t, x), h 0 (t, x), a 0 (t, x), m 0 (t, x), f 0 (t, x), g 0 (t, x), c 0 (t, x)) > (0, 0, 0, 0, 0) the corresponding non-negative solution (w(t, x), n(t, x), h(t, x), a(t, x), m(t, x), f (t, x), g(t, x), c(t, x)) is strictly positive for t > 0. Since the equilibrium points in equation ( 16) are non-negative, then let ε 0 ∈ (0, 1).Then according to Lemma (2.1) and the comparison principle of parabolic equations, there exists t 1 > t 0 > 0 such that, for any t > t 1 , and Thus, for t > t 0 , it is possible to obtain Since l j (w, n, h, a, m, f , g, c) in equation ( 26) is a C 1 function of w, n, h, a, m, f , g, c, where l 1 is quasi-monotone non-increasing in n, h, a, m, f , g, c, l 2 is mixed quasi-monotone in w, h, a, m, f , g, c, l 3 is quasi-monotone non-decreasing in w, n, a, m, f , g, c, l 4 is mixed quasi-monotone in w, n, h, m, f , g, c, h 5 is mixed quasi-monotone in w, n, h, a, f , g, c, l 6 is quasi-monotone non-decreasing in w, n, h, a, m, g, c, l 7 is mixed quasi-monotone in w, n, h, a, m, f , c and l 8 is mixed quasi-monotone in w, n, h, a, m, f , g, then by the method of upper and lower solutions we know that the system in (1) has a unique global non-negative solution w, n, h, a, m, f , g, c, [14].Thus, Therefore, ( w, n, h, ā, m, f , ḡ, c) and (w, n, h, a, m, f , g, g), are a pair of coupled upper and lower solutions of system (1), [22], respectively.

A construction and analysis of the numerical method
In this section, we describe the derivation of the fitted numerical method for solving the system in equation (1).We determine an approximation to the derivatives of the functions w(t, x), n(x,t), h(x,t), a(x,t), m(x,t), f (x,t), g(x,t), c(x,t), with respect to the spatial variable x.
Then we approximate the second order spatial derivative in the system in (1) by where and the denominator functions where φ w → ∆x, φ n → ∆x, φ m → ∆x, φ c → ∆x, as ∆x → 0. Let t f be a positive integer and j = 1/t f where 0 < t < t f .Discretizing the time interval [0,t f ] through the points We approximate the time derivative at t i by where where ψ w → ∆t, ψ c → ∆t as ∆t → 0. The denominator functions in equations ( 36) and (37) are used explicitly to remove the inherent stiffness in the central finite derivatives parts and can be derived by using the theory of nonstandard finite difference methods, see, e.g., [10,16,17] and references therein.
Combining the equation (36) for the spatial with equation (37) for time derivatives, we obtain and denoting the piecewise approximations of w(x j ,t i ), n(x j ,t i ), h(x j ,t i ), a(x j ,t i ), m(x j ,t i ), f (x j ,t i ), g(x j ,t i ), c(x j ,t i ) in the interval [x j ,t i ] with w i j , n i j , h i j , a i j , m i j , f i j , g i j , c i j , then we let Thus, in view of the FOFDM, we see that the local truncation errors (ς w , ς n , ς h , ς a , ς m , ς f , ς g , ς c ) are given by    eter values are as given in Table 1.  1.  1.

TABLE 1 .
Parameter values