EFFICIENT NUMERICAL METHOD FOR A MODEL ARISING IN BIOLOGICAL STOICHIOMETRY OF TUMOUR DYNAMICS

. In this paper, we extend a system of coupled ﬁrst order non-linear system of delay diﬀerential equations (DDEs) arising in modeling of stoichiom- etry of tumour dynamics, to a system of diﬀusion-reaction system of partial delay diﬀerential equations (PDDEs). Since tumor cells are further modiﬁed by blood supply through the vascularization process, we determine the local uniform steady states of the homogeneous tumour growth model with respect to the vascularization process. We show that the steady states are globally stable, determine the existence of Hopf bifurcation of the homogeneous tumour growth model with respect to the vascularization process. We derive, analyse and implement a ﬁtted operator ﬁnite diﬀerence method (FOFDM) to solve the extended model. This FOFDM is analyzed for convergence and we observe seen that it has second-order accuracy. Some numerical results conﬁrming the- oretical observations are also presented. These results are comparable with those obtained in the literature.


1.
Introduction. When a body develops a disease, the situation can be thought of as that of a predator and prey ecology. However some diseases are very much complex and the manner in which their formation and the manner they reinforce their presence in a body as compare to a situation in a physical environment. Such diseases are for example HIV and cancer. In this paper, we focus on a tumour. Therefore when a host is affected by a tumour, the system requires interaction with its environment. Therefore, it is not strange for the ecological system of cancer cells to interacts with surrounding cells, both healthy and malignant. This means cancerogenous and healthy cells start competing for resources, such as oxygen, nutrient and space. The cancerogenous cells not only compete against the healthy cells but also among themselves for the same resources. It is this competition that necessitates the consideration of the effects of biological stoichiometry to supplement the first principles of malignant tumours' cells [7]. Similar considerations of biological stoichiometry can be found in [4,5,23]. In this paper, we consider the biological stoichiometry derived in [7].
In [7], it is established that the qualitative dynamics of the models are essentially unchanged whether phosphorus limits blood vessel construction or not. Thus, Elser et al. [7] claimed that varying the supply of phosphorus continuously seems to create no new dynamical behavior. The aforementioned arguments prompted us to investigate the extended homogeneous tumour growth model directly.
Hence, our first aim in this paper is to investigate the qualitative features of the model with regard to blood vessel construction and determine the possibility of the time delay τ on the dynamics of the homogeneous tumour growth model. By applying the Poincaré normal form and center manifold theorem [8,24], we determine conditions on the functions and derive formulas which determine the properties of Hopf bifurcation [19] such as the direction of bifurcation, the period of periodic solutions and the stability of solutions. More specifically, we show that the positive equilibrium point losses its stability and the system exhibits Hopf bifurcation under certain conditions.
Our second aim is to solve the extended model. Thus, we develop an efficient numerical method for solving the extended model with respect to the qualitative features of the extended model. To this end we highlight our motivation for our numerical method. The deficiencies of the standard finite difference methods in solving the problems like the one in equation (2) are well-known. While explicit methods can solve such differential equations with low computational cost, they have the drawback that their stability regions are very small. This implies that severe restrictions on the time and space step-sizes will be required in order to achieve satisfactorily converging results. On the other hand, implicit schemes do have wider stability regions but the associated computational complexity is very high and they cannot achieve more than one order as compared to explicit methods that use the same number of stages [3].
When a single solid tumour is growing within an organ Elser et al. [7] mentioned that the initial mass starts near some genetically determined carrying capacity (k h ) and its vascularization process takes place at approximately 0.01kg. Since the parenchyma cells may contain distinct cell types that differ in their nutrient use and growth rates, then Elser et al. [7] developed the heterogeneous tumor model with dietary regulation as dx dt = x a min 1, Pe where x, y i (i = 1, 2 in this case), z, f , β, P e , P , L, n and m i denote mass of healthy cells, tumour mass contributed by the i th parenchyma cell type, mass of tumour micro-vessels, fraction of the total fluid within an organ, therapeutic intervention, extracellular phosphorus within the organ, the homeostatic regulation of the total amount of phosphorus, maximum proliferation rate of tumour cells, the mean amount of phosphorus in healthy cells and mean amount of phosphorus in parenchyma cells in that order. The cells proliferation and death at maximum per capita rates are denoted by a, b i and d x , d i , respectively, whereas α and g denote mass of cancer cells of which one unit of blood vessels can barely maintain and measurements of sensitivity of tumour tissue due to lack blood.
When a tumour has only one parenchyma cell type then the system of first order delay differential equations in equations (DDEs) (1) is known as homogeneous tumour growth model and heterogeneous tumour growth model when a tumour has more than one parenchyma cell type.
We observe that, all of the above models did not take spatial effects into account. In fact, cells can move around subject to many factors including diffusion. Thus, instead of depicting the models with purely time dependent ordinary differential equations (ODEs) with delay, it is more realistic to introduce the diffusion of the cells into the system, and the simplest way to reach this goal is to use the concept of reaction-diffusion equations [10,11]. Elser et al., [7] showed that at a steady state, tumor growth is no longer limited by its blood vessel infrastructure. For that reason, they also found out that when the tumor is viewed as a single entity, the homogeneous and heterogeneous models essentially generate the same dynamics. Thus, in this paper we consider the homogeneous tumor growth model. Therefore, incorporating spatial effects into equation (1), the homogeneous tumour growth model in (1) becomes ∂X ∂t − Dx∆X = X a min 1, Pe where ∆ denotes the Laplace operator, , Ω ∈ R 3 denotes a bounded domain with smooth boundary ∂Ω and ν denotes the outward unit normal on ∂Ω. The initial function η j (x, t) is Holder continuous on [−τ, 0] [18]. We imposed the no flux-boundary conditions in order to ensure that we exclude the external effects. The rest of the paper is as follows, Section 2 we analyse the steady state and existence of Hopf bifurcation analysis for two possible blood limiting cases. We derive, analyse our numerical method in Section 3 and present our numerical results in Section 4 and conclude the paper with Section 5.
2. Mathematical analysis of the homogeneous tumour growth model. To proceed, we recall from Elser et al. in [7], that in a phosphorus-rich environment, to denote such possibilities of limitations. Since maximum proliferation rate of tumour cells is dictated by the value of min 1, g Z−αY1 Y1 [7], then, at the steady state, the homogeneous tumour growth model in equation (2) becomes which gives the trivial equilibrium point (0, 0, 0), and the system if min 1, g Z−αY1 Y1 = 1. Thus, we consider both cases in the next two sections.
2.1. Local stability for min 1, g Z−αY1 Y1 = 1. When min 1, g z−αY1 and the determinant of the matrix in equation (7) is provided that a > d x , aϕ Px > d x and b 1 = d 1 . Thus, where the determinant of the numerator in equation (9) is Therefore, equation (9) becomes . For Y * 1 , we then have from equations in (7-8) that Similarly for Z * , from equations in (7)(8) we have which requires b 1 > d 1 and b 1 ϕ Py 1 > d 1 , as anticipated. Therefore, we have two positive equilibrium points and they are This gives us the following results.
Theorem 2.1. When min 1, g Z−αY1 Y1 = 1 and the following conditions simultaneously hold the homogeneous tumour growth model in equation (2) Therefore, when the maximum proliferation rate of tumor cells is greater than unity, then the steady states are positive as long as the genetically determined, carrying capacity for healthy cells is bigger than that of the tumour cells.
then the homogeneous tumor growth model has a unique positive constant solution Similarly, when the maximum proliferation rate of tumor cells drops below unity, then the steady states are positive as long as the genetically determined, carrying capacity for healthy cells is bigger than that of the tumour cells.

2.3.
Global stability of the uniform steady states. In this section we show that the uniform positive steady states are globally uniform. Since the homogeneous tumour growth model is of quasi-Lotka-Vorterra type [6], we let V : upon multiplying the last equation with Z. When min 1, g Z−αY1 Y1 ≡ 1, then equation in (17) becomes after dropping the bar signs. Similarly, min 1, g Z−αY1 Thus, we have the following results.
Thus, in this section we determine the positive solutions of both cases and found out that similar conditions hold for the survival of tumour cells which play an integral part in this stoichometric dynamics.
In the next sections we consider the existence of Hopf bifurcation.

2.4.
Stability of uniform equilibrium points and the existence of Hopf bifurcation for min 1, g Z−αY1 Y1 = 1 and min 1, g Z−αY1 Y1 = 1. In this section, we concentrate on the dynamical behavior of equation (2). Our goal is to investigate the stability of the equilibrium points of (2) and also the existence of Hopf bifurcation. This is achieved by taking the delay time τ as a bifurcation parameter. Thus, we study effects of the time delay on the dynamics of (2) as follows.
2.4.1. Stability of positive equilibrium points and the existence of Hopf bifurcation for min 1, g Z−αY1 Y1 = 1. Let (X * , Y * 1 , Z * ) be an equilibrium point for the system in equation (2) and ). Linearizing the system in equation (2) around (X * , Y * 1 , Z * ), and dropping bars again, we obtain Let 0 = µ 0 < µ 1 < · · · be the eigenvalues of the operator ∆ on Ω with the homogeneous Neumann boundary condition, then the characteristic equation for equation in (22) is given by where Negative real parts for the first equation in (23) requires that whereas in the second equation in (23), when τ = 0, Routh-Hurwitz criteria requires that   (2) is locally asymptotically stable when τ = 0.
Next, we examine when equation in (23) has pure imaginary roots λ = ±iω with ω real number and ω > 0. This is given by the following lemma. Proof. If λ = iω is a root of the characteristic equation (23) where ω > 0, then we have Separating real and imaginary parts, we have the following two equations Equations in (27) give possible values of τ and ω for which the characteristic equation in (23) can have pure imaginary roots. To see it we square each equation and we obtain Adding the two equations, we obtain which implies that Hence, Thus, in view of the first equation in (27) and the first equation in (28), we obtain as the critical values of τ , for j = 0, 1, 2, . . . , and this complete the proves. This gives us the following lemma.
EFFICIENT NUMERICAL METHOD 601 We have verified that the hypotheses for Hopf bifurcation occurs at τ * = τ * 0 except for the transversality condition. Differentiating the second equation in (23) with respect to τ , we have Re Thus, we obtain the following results.
Lemma 2.7. The transversality condition Re is satisfied.

2.4.2.
Stability of positive equilibrium points and the existence of Hopf bifurcation for min 1, g Z−αY1 Y1 = 1. Similarly, after linearizing the system in equation (2) around (X * , Y * 1 , Z * ), we obtain the following characteristic equation where

Negative real parts for the first equation in (35) requires that
whereas in the second equation in (35), when τ = 0, Routh-Hurwitz criteria requires that Therefore, we have the following result.  (2) is locally asymptotically stable when τ = 0.
Next, we examine when equation in (35) has pure imaginary roots λ = ±iω with ω real number and ω > 0. This is given by the following lemma. Proof. If λ = iω be a root of the characteristic equation (35) where ω > 0, then we have Separating real and imaginary parts, we have the following two equations Equations in (39) give possible values of τ and ω for which the characteristic equation in (35) can have pure imaginary roots. To see it we square each equation and we obtain Adding the two equations, we obtain which implies that Hence, Thus, in view of the first equation in (40) and equation in (43), we conclude that are the critical values of τ , for j = 0, 1, 2, . . . , and this complete the proves. This gives us the following lemma. where τ * = τ * 0 . We have verified the hypotheses for Hopf bifurcation to occur at τ * = τ * 0 except for the transversality condition. Differentiating the second equation in (35) with respect to τ , we have Re Hence, the following results.
Lemma 2.11. The transversality condition Re is satisfied.
In the next section we derive our efficient numerical method.
3. Derivation of the numerical method. In this section, we describe the derivation of the fitted numerical method for solving the system in equation (2). We determine an approximation to the derivatives of the functions X(t, x), Y 1 (x, t) and Z(t, x) with respect to the spatial variable x.
Let N x be a positive integer. Discretize the interval [0, x f ] through the points Z j (t) denote the numerical approximations of X(t, j), Y 1 (t, j), Z(t, j), then we approximate the second order spatial derivatives by where It is obvious that all the φ → ∆x as ∆x → 0. Let N t be a positive integer and ∆t = T /N t where 0 < t < T . Discretizing the time interval [0, T ] through the points where t n+1 − t n = ∆t, n = 0, 1, . . . , (t Nt − 1).
We approximate the time derivatives at t n by where where we see that all the ψ → ∆t as ∆t → 0. The denominator functions in (47) and (49) are used explicitly to remove the inherent stiffness in the central finite derivatives parts and are derived by using the theory of nonstandard finite difference methods, see, e.g., [9,16,17] and references therein.
Combining the equation (47) for the spatial derivatives with equation (49) for time derivatives, we obtain where, the no-flux boundary conditions are discretised by means of the central finite difference [2], j = 1, 2, . . . , x Nx − 1, n = 0, 1, . . . , t Nt − 1 and is denoting the history functions corresponding to the equation in Y 1 .
The system in equation (51) can further be simplified as which can be written as a tridiagonal system given by where j = 1, . . . , x Nx − 1, n = 0, . . . , t Nt − 1 and On the interval [0, τ ] the delayed arguments t n −τ belong to [−τ, 0], and therefore the delayed variables in equation (51) are evaluated directly from the history functions and equation (54) becomes Let s be the largest integer such that τ s ≤ τ . By using the system equation (56) we can compute X n j , (Y 1 ) n j , Z n j for 1 ≤ n ≤ s. Up to this stage, we interpolate the data . . , (t s , (Y 1 ) s j ), using an interpolating cubic Hermite spline ι j (t). Then (Y 1 ) n j = ι Y1 (t n , x j ), for all n = 0, 1, . . . , s and j = 1, 2, . . . , x Nx − 1.
For n = s + 1, s + 2, . . . , t Nt − 1, when we move from level n to level n + 1 we extend the definitions of the cubic Hermite spline ι j (t) to the point (t n + k, (Y 1 ) n j . Then the history term (H Y1 ) n j can be approximated by the functions ι j (t n − τ ) for n ≥ s. This implies that, and equation (56) becomes Our FOFDM is then consists of equations (51)-(58). This method is analyzed for convergence in next section and the corresponding numerical results are presented in Section 4.
3.1. Analysis of convergence. The convergence for the proposed FOFDM is proved via consistency and stability.
3.1.2. Stability of the numerical method. Substituting the exact solutions X(t n , x j ), Y 1 (t n , x j ), Z(t n , x j ) instead of the approximations X n j , (Y 1 ) n j , Z n j in equation (58), we obtain Subtracting equation (66) from (58) and taking the absolute values on both two sides and using equation in (65), we have where e n+1 P X = P X (t n+1 , ·) − P X (t n+1 , ·). Similar notations are used in the second and third equation in (67).
This proves that the method is unconditionally stable.

4.
Numerical results and discussions. In this section, we present our numerical results with respect to the maximum proliferation rate of tumor cells and blood supply limitation indicator to help visualize when there are indeed limiting healthy and tumour cells growth. In order to investigate it, we implement our numerical method such that the stability conditions given in Theorems 2.1-2.2 are satisfied. Due to the unavailability of diffusion values, we set the diffusion constants D x = 10 −3 , D y1 = 20 −4 , D x = 30 −5 . Following Elser et al. [7], we present our numerical results as follows.  In Figure 1, we present the case when the birth rate of healthy cell are decreased than that of parenchyma cell (a < b 1 ) and dearth rate of healthy cell is increased than that of parenchyma cell (d x > d 1 ) for (a) no treatment blocking phosphorus uptake by tumor cells, (b) lowered phosphorus P by 20%, (c) increased the time delay τ from 7 to 11 days and in (d) blocked tumor cell uptake of phosphorus by half.     In Figure 2, we present the case when the birth and dearth rates of healthy and parenchyma cells are equal (a = b 1 and d x = d 1 ) for (a) no treatment blocking phosphorus uptake by tumor cells, (b) lowered phosphorus P by 20%, (c) increased the time delay τ from 7 to 11 days and in (d) blocked tumor cell uptake of phosphorus by half.
In Figure 3, we present the case when the birth of healthy cells are increased than that of a parenchyma cell (a > b 1 ) and increased death rate of healthy cells than that of a parencyma cell (d x < d 1 ) for (a) no treatment blocking phosphorus uptake by tumor cells, (b) lowered phosphorus P by 20%, (c) increased the time delay τ from 7 to 11 days and in (d) blocked tumor cells uptake of phosphorus by half.
The numerical solutions in Figure 1 present a slight growth of healthy cells and fast growth of tumour cell during the first 40 days from the infection date. We also see that after 40 days of infection, both cells converges to their steady states, with a slight decrease in their growth. This is due to the competition for resources among cells, presented by this case. After 40 days the process of vascularization starts to grow because the tumour cell has reached some genetically size. We also see that lowered phosphorus and blocked tumor cell uptake of phosphorus do not differ from the no treatment blocking phosphorus uptake by tumor cells, whereas the above features are postponed when we increased the time delay τ from 7 to 11 days.
The numerical solutions in Figure 2 present an increase growth of healthy cells and tumour cell during the first 40 days from the infection date. We also see that after 40 days of infection, both cells converges to their steady states, with a slight increase for healthy cells and a decrease growth of tumour cells. This is due to the competition for resources among cells presented by this case. After 40 days the process of vascularization starts to grow because the tumour cell has reached some genetically size. We note that lowered phosphorus and blocked tumor cell uptake of phosphorus do not differ from the no treatment blocking phosphorus uptake by tumor cells, whereas the above features are postponed when we increased the time delay τ from 7 to 11 days.
The numerical solutions in Figure 3 present an increase growth of healthy cells as compare to the previous cases and a drastic decrease for tumour cell during the first 40 days from the infection date. After 40 days of infection, both cells converges to their respective steady states, with an increase of healthy cells and a decrease growth of tumour cell. This is due to the competition for resources among cells presented by this case. After 40 days the process of vascularization starts to grow, thus causing tumour cell to grow gradually. We note that lowered phosphorus and blocked tumor cell uptake of phosphorus do not differ from the no treatment blocking phosphorus uptake by tumor cell, whereas the above features are postponed when we increased the time delay τ from 7 to 11 days.

5.
Conclusion. In this paper, we have considered the biological stoichiometry of tumour dynamics, vascularized by a single solid tumor growing within the confinement of an organ and the environment provided by an organ. We examined the presence of positive steady state solutions for all the possible limiting cases with respect to proliferation rate of tumour cell and determine the existence of Hopf bifurcation. Thus, our numerical solutions clearly present the fact that the biological stoichiometry of tumour dynamics is real and can contribute quite a great deal toward the development of therapeutically drugs which can contribute toward healing tumour and tumour related diseases. Thus, our approach in this work should serve as a first numerical attempt to incorporate the detailed effects of healthy and tumor cells competing for both space and essentials, but limiting nutrients within a host.