Mathematical modeling on helper T cells in a tumor immune system

Activation of CD$8^+$ cytotoxic T lymphocytes (CTLs) is naturally regarded as a major antitumor mechanism of the immune system. 
In contrast, CD$4^+$ T cells are commonly classified as helper T cells (HTCs) on the basis of their roles in providing help to the generation and maintenance of effective CD$8^+$ cytotoxic and memory T cells. 
In order to get a better insight on the role of HTCs in a tumor immune system, 
we incorporate the third population of HTCs into a previous two dimensional ordinary differential equations (ODEs) model. 
Further we introduce the adoptive cellular immunotherapy (ACI) as the treatment to boost the immune system to fight against tumors. 
Compared tumor cells (TCs) and effector cells (ECs), the recruitment of HTCs changes the dynamics of the system substantially, 
by the effects through particular parameters, i.e., the activation rate of ECs by HTCs, $p$ (scaled as $\rho$), 
and the HTCs stimulation rate by the presence of identified tumor antigens, $k_2$ (scaled as $\omega_2$). 
We describe the stability regions of the interior equilibria $E^*$ (no treatment case) 
and $E^+$ (treatment case) in the scaled $(\rho,\omega_2)$ parameter space respectively. 
Both $\rho$ and $\omega_2$ can destabilize $E^*$ and $E^+$ and cause Hopf bifurcations. 
Our results show that HTCs might play a crucial role in the long term periodic oscillation behaviors 
of tumor immune system interactions. 
They also show that TCs may be eradicated from the patient's body under the ACI treatment.


(Communicated by Gail Wolkowicz)
Abstract. Activation of CD8 + cytotoxic T lymphocytes (CTLs) is naturally regarded as a major antitumor mechanism of the immune system. In contrast, CD4 + T cells are commonly classified as helper T cells (HTCs) on the basis of their roles in providing help to the generation and maintenance of effective CD8 + cytotoxic and memory T cells. In order to get a better insight on the role of HTCs in a tumor immune system, we incorporate the third population of HTCs into a previous two dimensional ordinary differential equations (ODEs) model. Further we introduce the adoptive cellular immunotherapy (ACI) as the treatment to boost the immune system to fight against tumors. Compared tumor cells (TCs) and effector cells (ECs), the recruitment of HTCs changes the dynamics of the system substantially, by the effects through particular parameters, i.e., the activation rate of ECs by HTCs, p (scaled as ρ), and the HTCs stimulation rate by the presence of identified tumor antigens, k 2 (scaled as ω 2 ). We describe the stability regions of the interior equilibria E * (no treatment case) and E + (treatment case) in the scaled (ρ, ω 2 ) parameter space respectively. Both ρ and ω 2 can destabilize E * and E + and cause Hopf bifurcations. Our results show that HTCs might play a crucial role in the long term periodic oscillation behaviors of tumor immune system interactions. They also show that TCs may be eradicated from the patient's body under the ACI treatment.
1. Introduction. A tumor, also known as a neoplasm, is an abnormal mass of tissue which may be solid or fluid filled, but it does not mean cancer. Tumors can be benign (not cancerous), pre-malignant (pre-cancerous), or malignant (cancerous). The immune system can monitor and protect the host from tumor onset, and immune deficiencies are accompanied by an increased risk of cancer. However, the mechanisms of tumor evasion and the immune response are still not completely understood. Cell-mediated immunity with cytotoxic T lymphocytes (CTLs) [3] and natural killer cells (NKs) [35], generally called effector cells (ECs) that are cytotoxic to the tumor cells (TCs), plays an essential role in immune responses against tumors [7,28]. Besides the cellular immunity, the humoral antibody response in the control of tumor growth has been reported [31]. Moreover, efficient antitumor immunity requires the concerted action of both CD4 + T cells and CD8 + T cells. CD4 + T cells are crucial for the induction, generation and maintenance of effective CD8 + cytotoxic and memory T cells, a phenomenon known as CD4 + T cell help [18]. They are commonly classified as helper T cells (HTCs). When CD4 + T cells are stimulated through antigen presentation by macrophages or dendritic cells, they can secrete the cytokine interleukin-2 (IL-2 known as cytokines which are information signaling molecules) to stimulate proliferation of T cells. CD4 + T cells orchestrate the antitumor CD8 + cytotoxic T cell responses through direct cellto-cell interaction and interleukin-2 stimulation. They can directly help CD8 + T cell activation via CD40-CD154 interaction [4,23]. Activated CD8 + T cells differentiate into CTLs and can directly kill tumor cells. Recently, it has been reported that CD4 + T cells may provide greater antitumor effect than CD8 + T cells [14,26]. Emerging evidence suggests that CD4 + T cells can even develop cytotoxic activity to tumors [25,27]. In this paper, we mainly focus on the helper role of HTCs to ECs, rather than the cytotoxic role of HTCs to TCs. The CD4 + T cell help is well known to play a vital role in the dynamics of tumor immune system interactions, which is one motivation for this work.
The study how the immune system interacts with cancer cells continues to lead to new progress in the treatment of cancer, including the development and application of many different immunotherapy strategies such as cytokine-based therapies, antibody-based therapies, cellular therapies, therapeutic vaccines and so on. One of the most promising advances is a therapeutic class called adoptive cellular immunotherapy (ACI). It uses a patient's own cells to stimulate the immune system and is usually performed in conjunction with large amounts of IL-2. There are two types of immune cells that are cultured for this purpose: lymphokine activated killer cells (LAK): these cells are taken from patients and in vitro cultured with high concentrations of IL-2 of peripheral blood leukocytes, then injected back to cancer site; tumor infiltrating lymphocytes (TIL): these cells are taken from lymphocytes recovered from the patient, and incubated with high concentrations of IL-2 before injected back to the patient [15,33]. The use of ACI may provide patients with a promising alternative treatment in the fight against cancer, which is another motivation for this work.
There is an extensive body of work which develops models to help understand the interactions between the immune system and a growing tumor [2,5,6,10,13,15,17,20,21,22]. Tumor immune models have been around since the early 1990s. A good summary of early works of tumor immune dynamics can be found in [1]. Recently Eftimie et al. (2011) [10] gave a detailed review of non-spatial models described by ODEs of interactions between the immune system and cancer. They reviewed models starting from the very simplest (involving a single equation) and built up to much more complex models that include greater biological details. The simplest mathematical model one can develop for tumor immune system interactions consists of two compartments: tumor cells and immune cells. A simple but classical mathematical model of a cell mediated response to a growing tumor cell population was proposed and analyzed by Kuznetsov et al. [17] in 1994, which differs from most others because it takes into account the penetration of TCs by ECs as well as the inactivation of ECs. Kuznetsov-Taylor model [17] can be applied to describe two different mechanisms of the tumor: tumor dormancy and sneaking through. Galach [13] firstly simplified Kuznetsov-Taylor model by replacing the Michaelis-Menten form with a Lotka-Volterra form for the immune reactions. Then she introduced the time delay into the simplified model and observed a state of the "returning" tumor. Commonly, the time delay can induce oscillatory behaviors. In 1998, Kirschner and Panetta [15] generalized Kuznetsov-Taylor model and illustrated the dynamics between TCs, ECs and IL-2. They firstly introduced ACI treatments into their models which can explain both short-term tumor oscillations in tumor sizes as well as long-term tumor relapse. In 2009, Kirschner and Tsygvintsev [16] initiated and performed a global analysis of Kirschner-Panetta model [15] by using the generalized Lyapunov method. Further, Starkov and Coria [30] continued studies of global dynamics of Kirschner-Panetta model with the help of applying the localization method of compact invariant sets. Arciero et al. [2] extended Kirschner-Panetta model by including a suppressive cytokine known as TGF-β. de Pillis and colleagues [5,6] presented several mathematical models to investigate the effect of NKs and CTLs in tumor surveillance, with the goal of understanding the dynamics of immune-mediated tumor rejection. However, there are few mathematical models that investigate on the role of helper T cells (HTCs) in tumor immune system interactions. Szymanska [32] proposed two mathematical models of the immune response to the identification of cancer antigens, which included the population of HTCs to compare two types of immunotherapy. Eftimie et al. [9] investigated the interactions between tumors and the two T helper (Th) cell types, Th1 and Th2 cells, which are defined by the cytokines produced by the T cells upon stimulation. In this paper, we construct the model at the cellular level without cytokines which were well investigated in [15]. But we use the novel treatment approach proposed in [15] to include ACI, which represents the introduction of immune cells into cancer patients that have been stimulated to have specific antitumor activity. Besides we follow the line in [32] to consider HTCs as a whole population to activate more ECs, unlike in [9] which divided HTCs into Th1 and Th2 cells to extend model complexity.
The purpose of the present paper is to formulate mathematical models to describe the tumor immune dynamics under the effects of HTCs and treatments. We provide theoretical analysis and implement numerical simulations to investigate the interactions between TCs and ECs and HTCs. Compared to two dimensional Kuznetsov-Taylor model [17] and Galach model without the time delay [13], the recruitment of HTCs changes the dynamics of the system substantially, by the effects through particular parameters (i.e., the activation rate of ECs by HTCs and the HTCs stimulation rate by the presence of identified tumor antigens). The numerical results show that an unstable periodic solution can bifurcate from the interior equilibrium E * (no treatment case), however a stable periodic solution can bifurcate from the interior equilibrium E + (treatment case), which implies that all the trajectories in the neighborhood of this equilibrium spiral towards the bifurcated limit cycle as time increases. This oscillatory phenomenon of tumor levels is known as Jeff's Phenomenon, and it has been observed clinically [12]. Periodic oscillations of tumor levels also have been observed in previous models of the interaction between tumors and the immune system (see, e.g., [13,15,20,21,22,34]). These results may be beneficial for the further study in mathematical modeling of the tumor immune system, especially for the better understanding of the role of HTCs in the tumor immune system interactions. Our results show that HTCs might play a crucial role in the long term periodic oscillation behaviors of tumor immune system interactions. They also show that TCs may be eradicated from the patient's body under the ACI treatment.
The paper is organized as follows. In the next section, we develop a mathematical model to investigate the interactions between TCs and ECs and HTCs. And we give the local stability analysis of three equilibria as well as local Hopf bifurcation at the interior equilibrium E * . In section 3, we incorporate a treatment term into our previous three dimensional model. And we give the local stability analysis of two equilibria as well as local Hopf bifurcation at the interior equilibrium E + . In section 4, numerical simulations are given to support the analytical results. In section 5, we provide a brief summary and discussions.
2. Mathematical modeling. The idea to use the qualitative theory of ODEs in mathematical biology goes back to 1920s when Lotka and Volterra developed a simple mathematical model in population dynamics [33]. In 1994, Kuznetsov et al. [17] applied Lotka-Volterra ideas to cancer modeling. The model was expressed as follows: where T (t) represents the population of TCs (prey) and E(t) represents the population of ECs (predator). Here the logistic form aT (1 − bT ) describes the growth of TCs, where a is the maximal growth rate and b −1 is the maximal carrying capacity for TCs. The mass action form −nET describes the loss of TCs due to the presence of ECs. s is the normal (i.e., not increased by the presence of the tumor) rate of the flow of adult ECs into the tumor site, and d is natural death rate of ECs. The Michaelis-Menten form p ET g+T describes the growth of ECs in response to TCs. The mass action form −mET describes the decay of ECs due to the interaction with TCs. All parameters a, b, n, s, d, p, g, m are positive.
Inspired from [17], Galach simplified system (1) by replacing the Michaelis-Menten form p ET g+T with Lotka-Voterra form θET . Therefore, she obtained the following model: where k = θ − m. All parameters a, b, n, s, d except k are positive. The sign of k depends on the relation between θ and m. If the stimulation coefficient of the immune system exceeds the neutralization coefficient of ECs in the process of the formation of EC-TC complexes, then k > 0. Further she introduced the time delay τ into system (2) by using the term kE(t − τ )T (t − τ ), because the immune system needs some time to develop a suitable response after the recognition of non-self cells. The Dulac-Bendixon criterion was used to show that there are no nonnegative periodic solutions to (1) and (2), which means that no Hopf bifurcations giving rise to limit cycles occur.
The idea in this paper is to investigate the role of HTCs in the tumor immune system. Therefore we incorporate the third population of HTCs into system (1). For the mathematical simplicity, bilinear term also has been used to describe the growth of the immune response to tumors, which is similar to system (2). As well-known, ECs can provide direct protective immunity by attacking and destroying the malignant TCs. In contrast, HTCs play major roles in providing help to promote or dampen cellular and humoral immune responses. HTCs can release various cytokines which stimulate ECs so that ECs can hunt and kill more and more malignant TCs. In our model we assume that HTCs can help ECs by the direct cell-to-cell interaction according to researches [4,18,23]. Several other assumptions are made such as the immune response to the appearance of TCs is positive [13].
Here we adopt the constant influx of HTCs [32], but we neglect the constant influx of effector ECs. Since ECs are cytotoxic T-cells produced as "naive" cells, that is, as not being able to produce any response to TCs unless they are activated by antigen presenting cells [19]. Then we shall establish a three dimensional ODEs model described as below: where T (t) and E(t) are the same variables as defined in system (1), and H(t) represents the population of HTCs. The first equation describes the rate change of the TCs population. Here the logistic growth term aT (1 − bT ) is chosen, where a is the maximal growth rate of the TCs population and b −1 is the carrying capacity of the biological environment for TCs. n represents the loss rate of TCs by ECs interaction. The second equation describes the rate change of the ECs population.
ECs have a natural lifespan of an average 1/d 1 days. k 1 is the ECs simulation rate by ECs-lysed TCs debris. p is the activation rate of ECs by HTCs. The third equation gives the rate change of the HTCs population. s 2 is the birth rate of HTCs produced in the bone marrow. HTCs have a natural lifespan of an average 1/d 2 days. k 2 is the HTCs stimulation rate by the presence of identified tumor antigens. We nondimensionalize the model (3) by taking the following scaling: As suggested in [17], we choose the scaling E 0 = T 0 = H 0 = 10 6 cells to improve the performance of numerical simulations. Then, after the substitution of the scaling into (3) and replacing τ by t, we obtain the following scaled model: where x, y and z denotes the dimensionless density of the TCs population, the ECs population and the HTCs population respectively.
Without the HTCs population, system (4) becomes the dimensionless form: There always exists a trivial equilibrium P 0 (0, 0), which is unstable. Besides there always exists a boundary equilibrium P 1 (1/β, 0), which is locally asymptotically stable when ω 1 ≤ βδ 1 . When ω 1 > βδ 1 , there exists an interior equilibrium , which is locally asymptotically stable. We call P 0 , P 1 and P 2 the tumor-free equilibrium, the tumor-dominant equilibrium and the immune-control equilibrium respectively.
In order to investigate the local stability of the above three equilibria E 0 , E 1 and E * of system (4), we linearize the system and obtain the characteristic equation evaluated atĒ(x,ȳ,z), given by the following determinant: Theorem 2.1. System (4) always has one tumor-free equilibrium E 0 , which is unstable.
By the Routh-Hurwitz criterion, the roots of equation (9) have negative real parts if and only if a 1 > 0, a 3 > 0 and a 1 a 2 − a 3 > 0.
Proof. We follow the analysis of a third-order characteristic equation by Ruan and Wolkowicz in [29] to determine the occurrence of Hopf bifurcation and choose ρ as a bifurcation parameter.
The abscissas x of intersection points of the two graphs y = f (x) and y = g(x) are the roots of equation (17). Thus there exists only one positive equilibrium E + if α > σ1δ2 δ1δ2−ρσ2 (g(0) < f (0)), that is ρ < δ1δ2 σ2 − σ1δ2 ασ2 ≡ ρ 1 , and x + < 1/β. Next we show the stability of the above two equilibria E − and E + of system (14). We linearize the system (14) and obtain the characteristic equation whose expression is the same as the equation (8). (14) has one tumor-free equilibrium E − when ρ < ρ 0 , which is locally asymptotically stable if ρ > ρ 1 .

Numerical simulations.
In this section we demonstrate through numerical simulations that stability change may occur at the branch of the interior equilibria E * and E + , leading to Hopf bifurcations. First, we expect to find stability regions of E * and E + by simulations with the following set of parameter values: α = 1.636, β = 0.002, δ 1 = 0.3743, ω 1 = 0.04, σ 2 = 0.38, δ 2 = 0.055.
In the following we choose one point in Fig. 2a to show the stability of E * . And then we choose four points in Fig. 2b to show the stability of E + and the existence of stable periodic solutions. The existence of periodic solutions is common in mathematical modeling of tumor immune system interactions [13,15,20,21,22].
SoẼ * is stable (see Fig. 3). SoẼ + 2 is unstable. Fig. 4c shows that HTCs have a periodic solution with a small amplitude (value of amplitude is about 3.5). And a limit cycle is formed around the fixed pointẼ + 2 through the Hopf bifurcation, thus resulting in a stable periodic solution (see Fig. 4e).
In Fig. 2b, if we choose ρ = 0.042 and ω 2 = 0.015, there exists one interior equilibrium,Ẽ   SoẼ + 4 is unstable. Fig. 4d shows that HTCs have a periodic solution with a relatively large amplitude (value of amplitude is about 125). And a limit cycle is formed around the fixed pointẼ + 4 through the Hopf bifurcation, thus resulting in a stable periodic solution (see Fig. 4f).
Next, we use bifurcation diagrams to find new dynamics as the parameters vary. The variable parameter is ρ, which is used as a bifurcation parameter. For no treatment system (4), we let ρ vary between 0 and ρ 0 while fixing ω 2 = 0.002 and investigate possible bifurcations. The bifurcation diagram in ρ − x plane is given in Fig. 5a, which exhibits one Hopf bifurcation at ρ 4 = 0.02945. It is interesting to observe that for changing ρ, system (4) bifurcates from a stable equilibrium to an unstable periodic oscillation. However for the treatment system (14), we let ρ vary between 0 and ρ 1 while fixing ω 2 = 0.015 and investigate possible bifurcations. For E + to become unstable, characteristic roots have to cross the imaginary axis to the right when ρ increases. The bifurcation diagram in ρ − x plane is given in Fig. 5b, which exhibits two Hopf bifurcations at ρ 2 = 0.04018 and ρ 3 = 0.01323. We note that the branch of interior equilibrium E + can change from stable to unstable as parameter ρ increases. It is interesting to observe that for changing ρ, system (14) bifurcates from a stable equilibrium to a stable periodic oscillation with large amplitude of malignant cell density. A stable limit cycle is observed between ρ 2 and ρ 3 . This means that tumors and immune cells can coexist for a quite long time and their development exhibits stable periodic oscillation phenomenon. Bifurcation curve of period of cycles versus ρ is depicted in Fig. 7. For no treatment system (4), the tumor-free equilibrium E 0 is always unstable, which means that tumors can not be totally eliminated from patient's body. But for the treatment system (14), we note that the tumor-free equilibrium E − exists when 0 < ρ < ρ 0 and is stable when ρ > ρ 1 , and the immune-control equilibrium E + exists when 0 < ρ < ρ 1 . From Fig. 5b we can see that TCs level x reduces gradually with the increase of the parameter ρ. And x can reach at an extremely low level when ρ is very close to ρ 1 . As ρ increases across ρ 1 , E + disappears and turns to E − , which means that tumors have been eliminated and patient is definitely cured (Fig. 6a). However when the activation rate of ECs by HTCs is too high, that is ρ > ρ 0 , there has no any positive equilibrium. As a consequence the population of ECs will grow uncontrollable (Fig. 6b). It can be explained that too strong immune response to tumors is also harmful for normal cells.

5.
Conclusion. In this paper, we explore interactions between tumor cells and immune cells through a system of bilinear ordinary differential equations. Especially we pay attention to the helper role of CD4 + T cells in the tumor immune system. CD4 + T helper cells play an important role in orchestrating an immune response to tumors by stimulating the activation of the antitumor ECs which target and destroy the invading TCs. In this paper, we mainly consider the helper role of HTCs to ECs, rather than the cytotoxic role of HTCs to TCs. The key idea is the incorporation of the direct cell-to-cell interaction between ECs and HTCs [4,18,23]. Recently new experimental results show that CD4 + T cells may provide even greater antitumor effect than CD8 + T cells and NK cells [25,26,27]. It will be very interesting to consider directly cytotoxic role of the CD4 + T cells on tumors and use these cells to attack tumors [10].
The tumor immune response dynamics in vivo are very complex, and not well understood primarily because the measurements of the necessary parameters are difficult in vivo. In model (4), main roles are played by two parameters that are the activation rate of ECs by HTCs, p, (scaled as ρ), and the HTCs stimulation rate by the presence of identified tumor antigens, k 2 (scaled as ω 2 ). Besides we consider the model (14) under the ACI treatment, and still change ρ and ω 2 to show what happens. Our studies reveal certain thresholds for scaled parameters ρ and ω 2 , which are effective to control the unlimited growth of malignant tumor cells so as to control the oscillations in the system (see Fig. 2). Further, our results show that it is possible to reach the tumor-free steady state which represents that TCs can be eliminated by activating HTCs into ECs in some parameter regions(e.g. ρ 1 < ρ < ρ 0 ) under the treatment. In comparison to system (1) and (2), the recruitment of HTCs changes the dynamics of no treatment system (4) substantially, which undergoes subcritical Hopf bifurcation resulting in the creation of an unstable periodic orbit. Besides the recruitment of HTCs also changes the dynamics of the treatment system (14) substantially, which undergoes supercritical Hopf bifurcation resulting in the creation of a stable periodic orbit. When a stable periodic orbit exists, it can be seen as "safe" in the neighborhood of this closed orbit since the trajectories originating from there will overwhelmingly go towards it [22]. Therefore tumors and immune cells can coexist for a long term although the tumors are not eradicated. The oscillatory dynamics in the tumor immune system have been observed in some mathematical models [2,8,10,15,17,20,21,22]. Periodic oscillation behaviors can be explained biologically that while the immune system fights with tumors in the host, there exists a balance between them, since the periodic changes in the internal tissues and the external circumstances such that they can coexist in a bounded region [21,22]. Our results show that HTCs might play a crucial role in the long term periodic oscillation behaviors of tumor immune system interactions. They also show that by using the ACI treatment, the tumor cells may be eradicated from the patient's body under some conditions.