Mathematical Modeling of the Tumor–Immune System with Time Delay and Diffusion

: This paper proposes a partial differential equation model based on the model introduced by V. A. Kuznetsov and M. A. Taylor, which explains the dynamics of a tumor–immune interaction system, where the immune reactions are described by a Michaelis–Menten function. In this work, time delay and diffusion process are considered in order to make the studied model closer to reality. Firstly, we analyze the local stability of equilibria and the existence of Hopf bifurcation by using the delay as a bifurcation parameter. Secondly, we use the normal form theory and the center manifold reduction to determine the normal form of Hopf bifurcation for the studied model. Finally, some numerical simulations are provided to illustrate the analytic results. We show how diffusion has a signiﬁcant effect on the dynamics of the delayed interaction tumor–immune system.


Introduction
Some local tissue cells in the biological system lose their normal regulation or are infected by certain viruses in their growth and then become tumor cells. A primary cause of human death is malignant tumor cells. In order to minimize the possibility of a patient dying, it is necessary to use some specific treatments or a combination of therapies such as surgery, chemotherapy, radiotherapy, etc. However, the immune function of the body can inhibit the attack of normal cells by invaders and the transformation into potential cancer cells, which is a natural defense of biological organs against tumor growth [1]. For the early stages of cancer, the use of an immune response is effective to treat tumors [2,3]. Thus, immunotherapy is becoming an essential treatment for tumors [4]. Furthermore, the study of the dynamics of the tumor-immune interaction system by mathematical models is very important and is one of the biomedical models that has attracted considerable attention over the past few decades.
In order to understand the complex dynamics of the tumor-immune system, mathematical models, such as ordinary differential equations (ODEs) [5][6][7] and partial differential equations (PDEs) [8][9][10][11], are largely used to explore the dynamics of the immune system and tumor growth. Various valuable results have been achieved, and these results have an important theoretical and clinical meaning in tumor research. For example, in ref. [12], the authors present an analytic approach to describe and establish solutions to a porous medium system of equations, showing applications in invasive-invaded biological dynamics. Furthermore, interesting analytical results apropos invasive-invaded systems with nonlinear diffusion and advection are given in [13]. New results on the multidimensional stability of V-shaped traveling fronts for a reaction-diffusion equation with a nonlinear convection term in R n (n ≥ 3) are presented in [14]. It is certain that constructing mathematical models supplies us with a new view to understanding the tumor-immune system.
Scientists have found that the immunity of an organism itself, which can help the immune cells fight against tumor cells, may be delayed [15]. Furthermore, resting cells can hunt and kill malignant tumor cells more than they can degenerate and transform into hunting cells themselves. Meanwhile, such a process can be achieved within a certain time interval, sometimes even a short time interval. Moreover, numerous scholars researching tumor-immune system interaction models have used only current information for predicting the future state [16,17]; however, future development not only depends on present information but also on past phenomena. For improving this issue, numerous scholars have considered the time delay in tumor-immune systems [18,19]. Ref. [20] considered a system of three-dimensional nonlinear delay differential equations, with the authors including time delay as a representation of the time lag for the conversion of resting cells to hunting cells. Ref. [21] studied the dynamics of tumor-immune system competition with time delay, and showed that the dynamical behavior of the system depends on the time delay parameter. Ref. [22] constructed a system of a number of immune cells that cooperate to protect the organs from tumor cells; the model included time delay to illustrate the activation process of immune cells. Ref. [23] proposed a Lotka-Volterra competition model with a step-wise constant delay to illustrate the interaction between tumor cells and T lymphocytes, and investigated the necessary and sufficient conditions for local and global stability.
The motivation for this paper is as follows. First, due to the processes of degeneration, proliferation, and transformation of tumor cells that can be achieved within a certain time interval, the existence of a time delay in the tumor-immune interaction system is necessary. Moreover, many scholars have shown that a time delay induces a change in the dynamical behavior of models, for example in [24,25], which can be biologically interpreted. Furthermore, the models governed by DDEs use past phenomena for predicting future information, which is more realistic. Second, the interaction between tumor cells and the immune system not only depends on time, but also on the space they live in. Moreover, tumor cells could spread in organs and trigger some factors, initiating the diffusion of tumor cells. Therefore, it is important to include time delay and diffusion into the tumorimmune interaction system, although investigating the time delay in specific mathematical models, such as models including spatial factors, is significant and more reasonable when we consider tumor-immune interaction problems. Third, Galach [26] analyzed a simplified version of Kuznetsov and Talyor's model [27], where the immune reaction is described by a linear term, and then included the time delay in the interaction term between tumor cells and the immune system in the equation that describes the growth of tumor cells. However, the authors did not investigate the effect of time delay in Kuznetsov and Taylor's model with the Michaelis-Menten function. Moreover, the diffusion of tumor cells and immune cells was not taken into account. Both immune and tumor reactions take a certain amount of time to identify alien cells and then react after a certain time delay.
Kuznetsov and Taylor's model explains the response of effector cells to the growth of cancer cells; this model is different from other tumor-immune interaction system models because the authors consider the penetration of tumor cells by effector cells. Moreover, they use the Michaelis-Menten function to illustrate the production of tumor-specific immune cells. However, tumor invasion and metastatic spread are two essential and fundamental spatial processes that must be considered in such an important model. To establish the spatial effect issue, in this work, we consider Kuznetsov and Taylor's model governed by PDEs, considering the time delay in the interaction term. Our model is the first to consider the time delay and diffusion process together in the original Kuznetsov and Taylor model (with the Michaelis-Menten function). In this paper, we analyze the stability equilibria to obtain the parameters that affect tumor growth and immune response the most. Moreover, we provide a detailed study to display the effect of diffusion parameters and time delay on the stability of equilibria and the occurrence of Hopf bifurcation. This paper is organized as follows: We present a partial differential equation model with time delay and diffusion in Section 2. In Section 3, we carry out a qualitative analysis of the model, investigate the local stability of equilibria and detect the existence of Hopf bifurcation. In Section 4, by the use of the center manifold reduction method, we determine the normal form of Hopf bifurcation. We present some numerical simulations in Section 5. Finally, we present the conclusion of this study.

Model Formulation
The immune system contains many mechanisms for the identification and attack of tumor cells. The immune cells and tumor cells cannot recognize each other immediately, and there exists a time delay, τ, for responding. This response time may be small; however, it cannot be neglected. Furthermore, under the conditions of bounded resources and space, tumor cells and immune cells diffuse in the organs. Although the spatial environment affects the growth and interaction between tumor cells and immune cells, it does not only depend on time. Hence, it is also necessary to provide the complex interaction caused by spatial characteristics. Considering the above aspect, we construct the following system of non-linear equations governed by the partial differential equations in the presence of a time delay to study the dynamics of the tumor-immune interaction system: where Ω is a bounded domain in R m (m = 1, 2, 3) with smooth boundaries, ∂Ω. For simplicity, we chose Ω = [0, π], but all calculations can be extended for higher dimensions. u and v denote the density of cancer cells and immune cells at time t and location x. The descriptions of the parameters are given in Table 1.  The discrete time delay in interconnection terms

Equilibria
In this section, we will find the positive steady states of model 1 which are biologically feasible. The steady state of model 1 is the solution of the system: which are:

1.
Free equilibrium (i.e., the absence of cancer cells and immune cells) The axial equilibrium (i.e., the cancer cells exist), The interior equilibrium (i.e., the coexistence of cancer cells and immune cells), If (H 2 ) β + η > ζ holds, then E 0 and E 1 are the only equilibria points. If (H 1 ) β + η < ζ holds, then there is one interior equilibrium point, E 2 = (u * , v * ).

Stability and Bifurcation Analysis
To understand the dynamics of the suggested model and detect the effect of time delay and diffusion in the proposed model, we study the stability of equilibria and then we develop the bifurcation theory to describe how small parameters can change the qualitative behavior of the model.

Stability of E 0
The linearized system for model 1 about steady state E 0 has the form: We can see that the linearized system is independent of the time delay; therefore, the stability of the equilibrium E 1 will depend only on the diffusion parameters.
We can investigate the stability of E 0 by substituting: into the linearized system of model 1 to obtain the characteristic equation associated with E 0 as follows: Then, the stability result of the equilibrium E 0 is resumed in the following lemma.
Proof. The characteristic equation associated with equilibrium E 0 has two roots: λ 1 = −c 2 n 2 − β, which is always negative, and λ 2 = −c 1 n 2 + α, which is negative if c 1 > α n 2 = c * 1 , (n = 0) and positive if c 1 < c * 1 for all τ ≥ 0. If n = 0 (without diffusion), then the equilibrium E 0 is unstable for all τ ≥ 0. Remark 1. The above result shows the necessity of including diffusion in the model; however, the diffusion parameters are responsible for the stability of the equilibrium E 0 .

Stability of E 1
We pass now to the stability analysis of the equilibrium E 1 to show the parameters that have an impact on the dynamics of the system.
The linearized system about the steady state E 1 = (1, 0) has the form: The characteristic equation for τ = 0 is given by: The characteristic equation above has two roots: Under hypothesis (H 2 ), λ 2 is negative. Then, we have the stability results as follows.

Stability and Bifurcation Analysis of the Interior Equilibrium
When (H 1 ) holds, there exists an interior equilibrium. In what follows, we present the stability results for E 2 = (u * , v * ).
The characteristic equation around the equilibrium E 2 is given as follows: where A n , B n , C n and E n are given as follows: For τ = 0, the characteristic Equation (2) becomes: where Note that under the assumption (H 3 ) α(2u * − 1) + cv * > 0, φ 1 is positive. We are now able to formulate the next result.
When τ > 0, we suppose that ±iω (ω > 0) are a pair of purely imaginary roots of Equation (2), and by separating the real and imaginary parts we obtain: Which implies that Let ξ = ω 2 , then the above equation can be rewritten in the following form: The equation G(ξ) = 0: 1.

3.
Has two positive roots, ω 1,2 , if (H 6 ) : holds. If the assumption (H 6 ) holds, then the critical time delay is determined as follows: for l = 0, 1, 2 k = 0, 1, 2, . . . . (6) Now, the transversality conditions yield: This will signify that at least one eigenvalue exists with a non-negative real part for l . Furthermore, the conditions for the existence of Hopf bifurcation are necessary to prove that periodic solutions exist. (1), τ > 0. If the hypotheses (H 1 ) and (H 3 ) hold, and the stability conditions of Theorem 1 hold, then:

3.
If (H 6 ) is satisfied, then there exists m ∈ N such that: l , l = 1, 2, and k = 0, 1, 2, · · · . Remark 2. If (H 1 ) holds, then the stability analysis for the equilibrium E 1 and E 2 in the presence of time delay τ are similar.

Normal Form of Hopf Bifurcation
In this section, we try to derive the normal formal of Hopf bifurcation at the interior equilibrium E 2 = (u * , v * ) when τ = 0 for the system (1). We define u(x, t) = u(x, t) − u * and v(x, t) = v(x, t) − v * . Herein, we have omitted the bar for convenience. The system (1) is then given by: Let h = (h 11 , h 12 ) T be the eigenvector of the linear operator corresponding to eigenvalue iωτ, < h * , h >= h * T .h = 1, where h * = (h 21 , h 22 ) T is the normalized eigenvector of the adjoint operator of the linear operator corresponding to the eigenvalues −iωτ. We then obtain: be the solution to Equation (7), where Denote We then have We take perturbation τ = τ + + µ to deal with the delayed terms, we expand U(x, t − 1) and v(x, t − 1) at u(x, T 0 − 1, T 1 , T 2 , . . .) and v(x, T 0 − 1, T 1 , T 2 , . . .), respectively; that is, where . .) and j = 1, 2, 3, . . .. Substituting (9) and (13) into Equation (7), for the -order terms, we obtain: Since ±iωτ are the eigenvalues of the linear part of (7), the solution to (14) can be expressed in the following form: where rest means the conjugate of the preceding terms and h is given in (8) . For the 2 order terms, we obtain: We substitute Equation (16) into (15), we obtain the coefficient vector of term e iωτ T 0 , then we obtain: Suppose that the solution of Equation (16) is where with p k =< cos(nx) cos(nx), cos(kx) >= π 0 cos(nx) cos(nx) cos(kx)dx.
For the 3 order terms, we obtain: We substitute Equations (15) and (17) into Equation (18), we obtain the coefficient vector of term e iωτ T 0 , then we obtain: In conclusion, the normal form of Hopf bifurcation for the system (1) reduced on the center manifold is ∂F ∂T Let F → F ; thus, the above equation becomes: To obtain the normal form of Hopf bifurcation in polar coordinates, we let F = re iθ , we thus obtain: We state the following theorem. 1. If (M)µ < 0, the bifurcating periodic solutions are unstable.

Numerical Simulations
This part is devoted to the exposition of the results throughout numerical simulations. The values of parameters are chosen according to [27] and are stated in Table 2. In the absence of delay (i.e., τ = 0), the "built-in pdepe" function was used to plot the results. We provide both the spatial patterns of the solutions in addition to solution profiles for t ∈ [0, 100]. When the time delay was positive (i.e., τ > 0), numerical simulations were obtained using a forward-in-time and centered-in-space finite differences scheme. The spatial mesh interval was ∆x = 0.03 and the time step was ∆t = 0.01.
(a) (b)  Observing the figures of the numerical simulations with a time delay and in the absence of a time delay, we notice that when the time delay is inferior to the critical value, the interior equilibrium is locally asymptotically stable, which means that the immune system can keep tumor cell growth under control at this moment. Tumor cells may not proliferate extensively in the organs. It is established that the immune and tumor cells could coexist in the body. Although the efficiency of the immune system is shown (see Figures 2 and 4), as the time delay increases, periodic variation is observed (Figures 3 and 5). This demonstrates that the hyperplasia time for tumor cells is large, which diminishes the efficiency of the immune system.

Conclusions
In this work, we explore a mathematical model which explains the tumor-immune interaction system. The model which refers to Kuznetsov-Taylor's model with the Michaelis-Menten function is based on partial differential equations with a time delay. First, we investigate the stability of the equilibria in the absence of a time delay. Using the time delay as a bifurcation parameter, the necessary conditions for the stability of the interior equilibrium in the presence of a time delay are established. The critical value for which Hopf bifurcation occurs is obtained. We show that when the bifurcation parameter, τ, passes through the critical value, the stability of an interior equilibrium changes from stable to unstable. Applying the normal form together with the center manifold reduction, the normal form of Hopf bifurcation of the considered model is determined. In this paper, we consider spatial factors and the essential time for the reaction of the immune system and the tumor cells; then, by analysis of the results and numerical simulations, we observe that the critical value of Hopf bifurcation depends on diffusion. The stability of bifurcating periodic solutions can be also changed. Therefore, the dynamical behavior of the tumor-immune interaction system can be affected by spatial factors and time delay.