Mathematical Modelling for the Role of CD4+T Cells in Tumor-Immune Interactions

Mathematical modelling has been used to study tumor-immune cell interaction. Some models were proposed to examine the effect of circulating lymphocytes, natural killer cells, and CD8+T cells, but they neglected the role of CD4+T cells. Other models were constructed to study the role of CD4+T cells but did not consider the role of other immune cells. In this study, we propose a mathematical model, in the form of a system of nonlinear ordinary differential equations, that predicts the interaction between tumor cells and natural killer cells, CD4+T cells, CD8+T cells, and circulating lymphocytes with or without immunotherapy and/or chemotherapy. This system is stiff, and the Runge–Kutta method failed to solve it. Consequently, the “Adams predictor-corrector” method is used. The results reveal that the patient's immune system can overcome small tumors; however, if the tumor is large, adoptive therapy with CD4+T cells can be an alternative to both CD8+T cell therapy and cytokines in some cases. Moreover, CD4+T cell therapy could replace chemotherapy depending upon tumor size. Even if a combination of chemotherapy and immunotherapy is necessary, using CD4+T cell therapy can better reduce the dose of the associated chemotherapy compared to using combined CD8+T cells and cytokine therapy. Stability analysis is performed for the studied patients. It has been found that all equilibrium points are unstable, and a condition for preventing tumor recurrence after treatment has been deduced. Finally, a bifurcation analysis is performed to study the effect of varying system parameters on the stability, and bifurcation points are specified. New equilibrium points are created or demolished at some bifurcation points, and stability is changed at some others. Hence, for systems turning to be stable, tumors can be eradicated without the possibility of recurrence. The proposed mathematical model provides a valuable tool for designing patients' treatment intervention strategies.


Introduction
Cancer is one of the leading causes of death worldwide. According to the World Health Organization (WHO), there were 8.8 million deaths in 2015 due to cancer [1]. e global cancer burden is expected to rise to nearly 21.4 million cases and 13.5 million deaths by 2030 [2]. Cancer treatment includes surgery, radiation [3], hormonal therapy [4], virotherapy [5], chemotherapy, and recently immunotherapy [6]. Chemotherapy is a well-known method for treatment of cancer. It is based on the administration of drugs that can kill tumor cells.
ese drugs do not only kill tumor cells but can also kill normal cells. us, many patients suffer from side effects of treatment as well as resistance to therapy and recurrence [7].
New approaches to treatment have been investigated, and immunotherapy has been recently approved for the treatment of many types of cancers [8,9]. Immunotherapy is based on enhancing the effectiveness of the immune system to identify and kill tumor cells using two strategies: (1) passive immunotherapy, where effector components of the immune system are used to directly attack tumor cells. is strategy includes antibody-targeted therapy and genetically engineered T cells (e.g., chimeric antigen receptor [CAR]-T).
(2) Active immunotherapy, where the activity of the immune system is enhanced.
is includes cancer vaccines, cytokines, and adoptive cell therapy [10]. Cancer vaccines enhance cytotoxic T lymphocytes response to specific antigens produced by the tumor cells. Cytokines are proteins important for cell signalling and are produced by many cells, such as macrophages, B lymphocytes, and T lymphocytes. However, not all cytokines are approved for the treatment of cancer. e Food and Drug Administration (FDA) in USA approved only two cytokines, interleukin-2 (IL-2) and interferon alpha (IFN-α), since they demonstrated moderate clinical benefit [11]. IL-2 is most commonly used to trigger the immune system by activating T-cells and natural killer cells that identify and kill tumor cells. INF-α has a similar response rate to IL-2, but it does not achieve long-term patient survival compared to IL-2 [2,10,12,13]. Adoptive cell therapy (ACT) involves ex vivo stimulation and expansion of tumor infiltrated T cells, then infusing the cells back into cancer patients [10]. e majority of ACT clinical trials and animal studies comprise CD4 + T helper 1 cells and CD8 + cytotoxic T cells [14].
Mathematical modelling provides a powerful tool to describe and analyse many engineering and physical problems. It has also been used to describe some biological processes such as heart beats [15], diffusion of drugs [15], hepatitis C virus [16], management of HIV infection [17], treatment of diabetes [18], clearing the antibiotic resistant infection [19], aortic aneurysm formation [20], and tumor growth and cancer treatment. Studies in the field of cancer biology focus on developing suitable mathematical models involving quantitative approaches to understand various aspects of tumor growth and the response of cancer cells to clinical interventions. Jones et al. proposed a model describing tumor growth due to its internal pressure [15]. e model was in the form of a system of partial differential equations and included tumor size and internal pressure and nutrients' concentration without any treatment. Similar models, in the form of systems of partial differential equations, were also proposed by Tao et al. [21], Wei and Cui [22], Wise et al. [23], Frieboes et al. [24], Lee et al. [25], Zhang et al. [26], Knopoff et al. [27], de Pillis et al. [28], Villasana and Radunskaya [29], and Xu and Bai [30]. e last two models involved time delayed equations as well. A statistical model for cancer mortality has been studied by Ghosh and Samanta [31].
An important division of the mathematical models about cancer intends to investigate the effect of varying the concentrations of immunotherapy and chemotherapy and to explore the influence of various parameters. Valuable and comprehensive reviews can be found in [32][33][34][35]. Most of these models are in the form of a system of coupled nonlinear ordinary differential equations (see, for example, references [12,[36][37][38][39]). Along this line, de Pillis and colleagues developed a useful model for the interaction between the tumor and immune system [13,[40][41][42][43]. ey obtained a numerical solution and performed both stability analysis and sensitivity analysis. ey took into account the following assumptions: (a) a tumor grows logistically in absence of an immune response, (b) both natural killer and CD8 + T cells are capable of killing tumor cells, (c) both natural killers and CD8 + T cells respond to tumor cells by expanding and increasing cytolytic activity, (d) natural killers are normally present in the body, even when no tumor cells are present, (e) active tumor specific CD8 + Tcells are only present in large numbers when tumor cells are present, and (f ) both CD8 + T cells and natural killers become inactive after some number of encounters with tumor cells [40,42]. Arabameri et al. [8] proposed a model for immunotherapy using mature dendritic cells. Mamat et al. [2] proposed a modified model that includes the effect of both IL-2 and INF-α. A model was proposed by Eftimie et al. [44] and was modified by Anderson et al. [45] and by Hu and Jang [46] to study the effect of CD4 + T cells on the system. CD4 + T cells are helper cells that activate CD8 + T cells to achieve a regulated immune response and protect against tumors. However, CD4 + T cells can play a direct role in killing tumor cells via secretion of cytokines in addition to its traditional helper role [45][46][47][48][49].
Intriguingly, previous studies modelled the effect of natural killer cells, circulating lymphocytes and CD8 + T cells, and neglected the role of CD4 + T cells, or studied the role of CD4 + T cells, neglecting the role of other immune cells. In this study, a mathematical model is proposed for active immunotherapy that includes the effects of natural killer cells, circulating lymphocytes, CD8 + T cells, CD4 + T cells, and chemotherapy on the tumor. Both roles of CD4 + T cells are considered in the model. e work of de Pillis et al. [13,[40][41][42], Eftimie et al. [44], Anderson et al. [45], and Hu and Jang [46] are the main references for the proposed model. e model is solved numerically using the "Adams predictor-corrector" method with Adams-Bashforth predictor and Adams-Moulton corrector. A stability analysis is performed for patients' data to determine if the system is stable. For a nonzero tumor case, the tumor is totally eradicated by the treatment intervention and the possibility of recurrence is to be studied. At the zero tumor case, the system must be stable to avoid recurrence and a condition for stability is deduced. It is recommended to work on modifying the patient's body parameters, using cancer vaccines, for example, before applying the treatment intervention in order to satisfy the condition needed for stability [42]. Finally, bifurcation analysis is performed to study the effect of changing patient's parameters on stability.

Mathematical Model.
To illustrate the interaction between tumor cells T(t) and immune cells present in the microenvironment of the tumor without external effects, such as treatment, the following model is proposed: 2 Computational and Mathematical Methods in Medicine where N(t) is the natural killer cell population, L(t) is the CD8 + T cell population, Y(t) is the CD4 + T cell population, C(t) is the circulating lymphocyte cell population not including natural killer cells and active CD8 + Tand CD4 + Tcells [50,51], I(t) is the concentration of the IL-2 cytokine in the body, and D is defined as and the parameters a, a 1 , b, c, c 1 , d, e, f, g, g i , h, j, k, l, m, p, p i , q, r 1 , r 2 , s, u, α, α 1 , α 2 , β, β 1 , β 2 , μ 1 , μ i , and δ 2 are given in Table 1. e role of CD4 + Tcells is included in this model, and it depends mainly upon the work of Eftimie et al. [44], Anderson et al. [45], and Hu and Jang [46]. In this model, CD4 + T is included in equation (4), and its effect is included in the terms (c 1 T/(a 1 + T)) I and (β 2 T/(α 2 + T)) Y in equations (1) and (6) [46,54]. In equation (1), a T (1 − b T) represents tumor growth and is assumed to be logistic [40,42,46]. e term − cNT represents the natural killer-induced tumor death [40,42]. e next term − D T is the tumor lysis by CD8 + T cells [40,42]. Unlike CD8 + T cells, CD4 + T cells cannot kill tumor cells directly but through the cytokines that they produce. e role of cytokine in killing the tumor is represented by the term -(c 1 T/(a 1 + T)) I in equation (1) [46,54]. In equation (2), eC represents the production of natural killer cells from circulating lymphocytes [40,42]. is growth term is tied to the overall immune health levels as measured by the population of circulating lymphocytes [54][55][56]. − fN represents natural killers' turnover [51]. e term (g T 2 /(h + T 2 )) N represents the recruitment of natural killer cells induced due to the existence of the tumor [42,54,57]. e next term − p N T is natural killer death by exhaustion of tumor-killing resources [40,42,51]. In equation (3), − m L represents the inactivation of CD8 + T cells consisting only of natural death rate since no CD8 + T cells are assumed to be present in the absence of tumor cells [40,42]. e term (j D 2 T 2 /(k + D 2 T 2 ))L represents the recruitment in CD8 + T cell population due to the presence of the tumor [40,42,58]. e next term − qLT is the CD8 + T cell death by exhaustion of tumor-killing resources [40,42,51]. e term r 1 NT represents CD8 + Tcell stimulation by natural killerlysed tumor cell debris, and the term r 2 CT represents the activation of native CD8 + T cells in the general lymphocyte population [42,51]. e term − u N L 2 is the natural killer cell regulation of CD8 + T cells. It describes regulation and suppression of CD8 + T cell activity, which occurs when there are very high levels of activated CD8 + T cells without responsiveness to cytokines present in the system [40,42,59,60]. e activating effect of cytokines on CD8 + T cells is represented in equation (3) by the term (p i I/(g i + I))L [40,42,61]. CD4 + T cell population is represented by equation (4) which has a constant decaying rate μ 1 . e other decaying term − δ 2 T Y exists due to the interaction with the tumor. Activating effect of cytokines on CD4 + T cells due to existence of the tumor is represented by the term (β 1 T/(α 1 + T)) I in equation (4) [46,54]. Equation (5) represents the circulating lymphocytes cell population which are assumed to be generated at a constant rate α and death rate β [40,42,[54][55][56]. Equation (6) represents the concentration of cytokines which is decaying with a constant rate μ i . e activating role of the CD4 + T cells on cytokines due to the existence of the tumor is represented by the term ((β 2 T)/(α 2 + T)) Y in equation (6) [46,54]. e previous model represents the interaction of patient's body with the tumor only without any treatment. It can be generalized to introduce the effects of chemotherapy and/or immunotherapy; refer, for example, de Pillis and Radunskaya [40] and Hu and Jang [46]. e general model takes the following form: and v Y (t) are the concentrations of chemotherapy, IL-2 cytokine therapy, CD8 + T, and CD4 + T adoptive immunotherapy, respectively. e effect of chemotherapy on immune cells is represented by the terms e three terms illustrate the reduction in number of each type of immune cell populations, while tumor cell killing due to chemotherapy is represented by − K T (1 − e − M ) T and the concentration of the chemotherapeutic agent is represented by equation (14) [40,42]. e parameters K T , K N , K L , K C , and c are defined in Table 1.

Method of Solution.
In numerical methods for solving differential equations, a problem may occur when small Computational and Mathematical Methods in Medicine changes in the independent variable lead to large changes in the dependent variable. If this problem exists in a system of differential equations, the system is called a stiff system [62,63]. e treatment doses v L (t), v Y (t), v I (t), and v M (t) change suddenly from zero to very large values within small intervals of t, and systems (8)-(14) become a stiff system. We tried to solve the system using the Runge-Kutta method, but the solution was divergent in most cases due to the stiffness.
e Adams-Bashforth method [62][63][64][65] is a powerful method for solving stiff problems and was originally proposed in 1883. For the following nonlinear initial value problem, e approximate solution at t n is given by where Δt is the step size and k is the order of the method. e coefficients β 1 , β 2 , . . . , β k can be deduced using Taylor's expansion of y(t) near the point t � t 0 with error of order O((Δt) k+1 ) [65]. is means that, at any point, we can find the solution numerically using any number of the previous points. e Adams-Moulton method [62,63,65,66] is a modification of the Adams-Bashforth method, which was Death rate of natural killer cells [40,42].
Rate of which CD8 + T cells are stimulated to be produced; as a result, tumor cells killed by natural killer cells [40,42].
Computational and Mathematical Methods in Medicine proposed in 1926 by adding an implicit term to equation (16) to become Coefficients can be found in the same way of Adams-Bashforth method.
Finally, a combined two-staged method is proposed by first computing the predictor stage using the Adams-Bashforth method [62,65,67]: and then the corrector stage using the Adams-Moulton formula: is combined method is called the "Adams predictorcorrector method." In this paper, we use it to solve the systems of nonlinear ordinary differential equations presented by equations (1)-(6) and (8)-(14).

Results and Discussion
e model proposed in this study is solved numerically using the Adams predictor-corrector method up to the 12 th order, and the solutions are interpolated. e method has nonzero stability ordinates for this order [68]. Values of the parameters a, a 1 , b, c, c 1 , d, e, f, g, g i , h, j, k, K T , K N , K L , K C , l, m, p, p i , q, r 1 , r 2 , s, u, α, α 1 , α 2 , β, β 1 , β 2 , c, μ 1 , μ i , and δ 2 are given in Table 1 for humans. e solutions are obtained for three cases: without any treatment, with continuous treatment, and with pulsed treatment as detailed below.

First: A Case of No Treatment.
In this case, we illustrate the interaction between tumor cells and the immune system without any treatment. In Figures 1(a) and 1(b), simulation was performed for first patient's data where the initial natural killer cells N(0) � 10 3 , initial CD8 + T cells L(0) � 10, initial CD4 + T cells Y(0) � 10 6 , and initial circulating lymphocytes C(0) � 6 × 10 8 . e immune system could overcome tumor of size T(0) � 10 5 , as shown in Figure 1(a), but cannot overcome a tumor of size T(0) � 10 7 , as shown in Figure 1(b). e simulation is performed again with tumor size T(0) � 10 7 but with a stronger immune system such that N(0) � 10 3 , L(0) � 10 2 , Y(0) � 10 6 , and C(0) � 6 × 10 10 , which means that initial CD8 + T cells and initial circulating lymphocytes are increased. Figure 1(c) shows that tumor cells were destroyed in such case.
In this section, we demonstrate the role of different immune cells in abolishing the tumor. e results suggest that patient's body can overcome cancer cells up to a certain size. In case of larger tumors, however, a stronger immune system is necessary to overcome cancer cells.

Second: A Case of Continuous Treatment.
In this case, we simulated the model with continuous immunotherapy. is is performed using the first patient's data, with initial natural killer cells N(0) � 10 3 , initial CD8 + T cells L(0) � 10, initial CD4 + T cells Y(0) � 10 6 , initial circulating lymphocytes C(0) � 6 × 10 8 , and initial tumor size T(0) � 10 6 . First, a continuous adoptive therapy is introduced using CD8 + T cells with concentration v L � 10 4 and cytokine IL-2 with concentration v I � 10 4 . Figure 2(a) demonstrates that the tumor cannot be eliminated using this treatment. Second, a continuous adoptive therapy is introduced using CD4 + T cells with concentration v Y � 10 6 instead of both CD8 + T cells and IL-2 cytokine therapies, and the tumor could be eliminated, as shown in Figure 2 Treatment is applied using the second patient's data with initial conditions: N(0) � 10 3 , L(0) � 10 2 , Y(0) � 10 6 , C(0) � 6 × 10 8 , and tumor size T(0) � 10 6 . Figure 3(a) illustrates the results of using continuous CD8 + T cell therapy v L � 10 6 together with cytokine IL-2 treatment v I � 10 6 , and Figure 3(b) shows the result when using only CD4 + T cell therapy v Y � 10 6 . In both cases, the tumor is not eliminated, but when using the three therapies together, v L � 10 6 , v I � 10 6 , and v Y � 10 6 , the tumor is eliminated in about 10 days, as shown in Figure 3(c).

ird: A Case of Pulsed Treatment.
In this case, the patient is administered therapy on certain days only. In Figures 4(b), 5(c), and 5(d), a tumor of initial size T(0) � 8 × 10 5 is assumed to exist in the first patient with initial conditions: N(0) � 10 3 , L(0) � 10 2 , Y(0) � 10 6 , and C(0) � 6 × 10 8 . In each figure, we use a certain regimen of treatment. In Figure 4(a), five pulses of chemotherapy of intensity v M � 5 starting from the 6 th day with 5 days period can eliminate the tumor, as shown in Figure 4(b). An alternative to chemotherapy, 7 pulses of CD8 + T cell treatment v L � 10 5 combined with IL-2 cytokine v I � 10 5 with 2 days period starting from day 6 were proposed, as illustrated in Figures 5(a) and 5(b). Figure 5(c) demonstrates that this combination is not enough to eradicate the tumor. Oscillations of CD4 + T cell population and IL-2 cytokine concentration are experienced due to fast consumption of the treatment by the resisting tumor cells. ese large changes in small time intervals explain the reason of the failure of the Runge-Kutta method in obtaining the numerical solution. When we use the same treatment schedule with additional 5 pulses of CD4 + T cell treatment v Y � 10 7 with 2 days period starting from day 6, the tumor is eliminated, as shown in Figure 5(d).
ese simulations suggest that CD4 + T cells could replace chemotherapy and thus could spare patients the agony of chemotherapy-induced toxicities and side effects.
Pulsed treatment is applied to the first patient with initial conditions: N(0) � 10 3 , L(0) � 10, Y(0) � 10 6 , and C(0) � 6 × 10 8 but with a larger tumor size T(0) � 3 × 10 6 . In Figure 5           is applied starting from the 6 th day with 5 days period. Chemotherapy alone is insufficient to overcome tumor growth. In Figure 5(f ), the same chemotherapy drug is applied combined with 5 pulses of CD4 + T cell therapy v Y � 6 × 10 6 with 2 days period starting from day 6. Combined treatment could abolish the tumor, while chemotherapy alone was not able to do so. e results suggest that, depending on tumor size, CD4 + T cell treatment can be an alternative to chemotherapy or is needed in combination with therapy.
De Pillis and Radunskaya [40] studied the case of a larger tumor of initial size T(0) � 10 7 . Pulsed treatment is applied to the first patient with initial conditions: N(0) � 10 3 , L(0) � 10, and C(0) � 6 × 10 8 . ey applied 9 pulses of chemotherapy of intensity v M � 5 starting from the 1 st day with 10 days period combined with 6 pulses of the IL-2 cytokine therapy v I � 5 × 10 5 from day 8 to day 11 and CD8 + T cells treatment v L � 10 9 during days 7 and 8. Combined treatment could abolish the tumor in this case, as shown in Figure 5(g). As an alternative to this therapy, with the same initial conditions, however, we added the effect of CD4 + Tcells (Y(0) � 10 6 ). We suggest 5 pulses of CD4 + T cell therapy v Y � 4 × 10 7 with 2 days period starting from day 6 combined with three pulses of chemotherapy of intensity v M � 5 applied starting from the 1 st day with 10 days period. Figure 5(h) illustrates that, in this case, the tumor can be also abolished. It is observed that we used a lower dose of the chemotherapeutic agent. Also, we used CD4 + T cell therapy as an alternative to both CD8 + T cell therapy and cytokine therapy.

Stability Analysis.
In this section, stability analysis is performed for the system of equations (1) to (6). e main advantage of studying stability is to investigate the possibility of eradicating the tumor without the occurrence of recurrence. If an equilibrium point is stable, it will attract the system towards it, in its domain of attraction. However, if the equilibrium point is unstable, any small perturbation from equilibrium will cause the system to move away from this equilibrium point; i.e., upon terminating treatment, the tumor will escape immune surveillance unless every single tumor cell is killed. To avoid recurrence of tumor after treatment, the zero tumor case has to be stable. In this section, a condition for stability of the zero tumor case is deduced. e first step is equating the right-hand sides of these unforced autonomous equations to zero [69][70][71][72]. From equation (1), we can deduce that either T � 0 or a (1 − b T) − c N − D − (c 1 /(a 1 + T)) I � 0 is at equilibrium. In the following part, the two equilibrium conditions are considered; first, the zero tumor equilibrium, T eq � 0, and then the nonzero tumor equilibrium, a (1 − b T eq ) − c N eq − D eq − (c 1 /(a 1 + T eq )) I eq � 0 and T eq ≠ 0.

Equilibrium and Stability for Nonzero
Tumor. In this case, then, 8 Computational and Mathematical Methods in Medicine 80 80

Computational and Mathematical Methods in Medicine
From equations (2) and (5) at equilibrium, From equation (7) at equilibrium, From equation (3) at equilibrium, From equations (4) and (6) at equilibrium, en, from equation (28), the nonzero tumor equilibrium will be divided into two cases: the zero cytokine case at which I eq � 0 and nonzero cytokine case at which ((β 1 T eq / (α 1 + T eq )) − (μ i (α 2 + T eq )(μ 1 + δ 2 T eq )/(β 2 T eq ))) � 0 and I eq ≠ 0. In the nonzero tumor equilibrium with zero cytokine, substitution in equation (4) at equilibrium gives Y eq � 0. Equilibrium points in this case can be found by solving equations (26) and (27) for T eq and L eq and then calculating other variables. Data for the two patients are used for numerical solution. Since equations (26) and (27) are highly nonlinear, the solution is obtained by plotting both of them in the T − L plane and determining the intersections. In the nonzero tumor equilibrium with nonzero cytokine, Rearranging, e equilibrium points in this case can be found by solving equation (30) for T eq . For each value of T eq , equations (26) and (27) are solved for I eq and L eq and then substituted for other variables. Due to the high nonlinearity of the two equations, the solution is obtained by plotting both of them in the I − L plane and determining the intersections. For each equilibrium point, a linearization method is used to write the Jacobian and to determine the eigenvalues. Consequently, the stability of the system is specified. Table 2 illustrates the results for the parameter values given in Table 1. For more illustration, streamlines are plotted. We use the values of N eq , L eq , C eq ,   and I eq of each equilibrium point and allow Y and T to vary with time according to equations (1) and (4). For the first patient, the streamlines are plotted for equilibrium point 1, T � 0, and Y � 0, in Figure 6(a). One can notice that, for any state of the system in the vicinity of this unstable equilibrium point, the system moves away from that point. e streamlines are plotted for equilibrium point 2, T � 1.90609 × 10 7 and Y � 0, in Figure 6(b). Also, for this unstable equilibrium point, any state of the system in the vicinity of this point, the system moves away from it. In the same manner, the streamlines for the second patient are obtained; however, to be concise, they are not included.

Results of Stability Analysis for Both the Zero Tumor and Nonzero Tumor.
In this analysis, it is shown that all equilibrium points are unstable for both patients. is means that, for zero tumor case, the tumor can recur, and in nonzero tumor case, the size of the tumor is uncontrollable; i.e., it could increase without a limit. For the zero tumor, an analytical expression for the stability of the equilibrium point is obtained. Simply, if the condition a − (ceα)/(fβ) < 0 is satisfied, then the zero tumor case is a stable state and the tumor cannot recur. If the condition is not satisfied, then the tumor-free state is an unstable state and the tumor is recurrent. ough for the nonzero tumor the stability analysis is analysed numerically, changing the system parameters could lead to stable equilibrium. is is tackled in the next section. Hence, it is recommended for any successful treatment to work on modifying the patient's body parameters in order to satisfy the stability condition.

Bifurcation Analysis.
A bifurcation analysis is performed to study the effect of changing patient's parameters on the system stability [15,40,42,73]. is change could lead to turning unstable point(s) into stable; hence, for nonzero tumor, the size is limited and can be reduced by surgery or radiation and then immunotherapy and/or chemotherapy to reach a stable zero tumor case. In the zero tumor case, it is required to change patient's parameters, using cancer vaccines, for example, to satisfy the condition of stability stated in the previous section and prevent recurrence. Bifurcation diagram is plotted for parameters a and c for the first patient and the second patient. In the same manner, bifurcation analysis can be studied for other system parameters.
First, the value of parameter a is changed from 0 to 10. Equilibrium points are monitored, and then bifurcation diagrams are plotted. In Figure 7(a), the bifurcation diagram for parameter a with first patient's data is illustrated. For a < 2.02257 × 10 − 5 , there is only one equilibrium point at   Figure 6: Streamlines for first patient's data: (a) at equilibrium point 1, N eq � 315534, L eq � 0, C eq � 6.25 × 10 10 , and I eq � 0; (b) at equilibrium point 2, N eq � 199.335, L eq � 2.82461 × 10 6 , C eq � 6.25 × 10 10 , and I eq � 0.
which T � 0 and the system is stable.
e value a � 2.02257 × 10 − 5 is a bifurcation point because the system is turned into unstable after it. For 0.01 ≤ a ≤ 2.33, there is another equilibrium point due to the nonzero tumor case with zero cytokine. is point moves on curve 2 according to the value of a, and it is always unstable. e boundaries of curve 2, a � 0.01 and a � 2.33, are bifurcation points because this equilibrium point is created at a � 0.01 and vanished at a � 2.33. For a > 2.33, there is only one equilibrium point at the zero tumor state and it is unstable.
In Figure 7(b), the bifurcation diagram for the parameter a with second patient's data is plotted. It shows that, for a < 2.02257 × 10 − 5 , there is only one equilibrium point at the zero tumor case and the system is stable. e value a � 2.02257 × 10 − 5 is a bifurcation point because the system is turned into unstable after it. For 0.01 ≤ a ≤ 0.59, there is another equilibrium point due to the nonzero tumor case with zero cytokine. is point moves on curve 2 according to the value of a and is always unstable. e points a � 0.01 and a � 0.59 are bifurcation points. For the nonzero tumor with nonzero cytokine, there are two equilibrium points, illustrated by curves 3 and 4. e two equilibrium points created at the bifurcation points a � 0.3 and a � 0.6, respectively. Curve 3 is always unstable and exists for a ≥ 0.3. Curve 4 appears for a ≥ 0.6, and it represents the case with small size tumor and is always stable. A healthy state might be one which maintains the system at this stable low-tumor size level. Parameter a can be altered through treatments such as vaccination, and the tumor can be abolished/reduced by surgery or radiation and then immunotherapy and/or chemotherapy to reach the attraction zone of the tumor-free equilibrium state in the stable region of curve 1.
Bifurcation analysis is then performed for parameter c. It is changed from 0 to 4 × 10 − 6 . Figure 7(c) shows that, for first patient's data, there are two equilibrium points for all values of c. e equilibrium point moving on curve 2 for the nonzero tumor is always unstable. However, the equilibrium point moving on curve 1 for the zero tumor is unstable for  c < 1.365 × 10 − 6 and stable for c > 1.365 × 10 − 6 . e point c � 1.365 × 10 − 6 is a bifurcation point. Figure 7(d) shows that, for second patient's data, there are two equilibrium points for c < 0.4 × 10 − 6 and three equilibrium points for c ≥ 0.4 × 10 − 6 . e equilibrium points moving on curves 2 and 3 are always unstable; however, the equilibrium point moving on curve 1 is unstable for c < 1.356 × 10 − 6 and stable for c > 1.356 × 10 − 6 . e point c � 1.356 × 10 − 6 is a bifurcation point. Point c � 0.4 × 10 − 6 at which curve 3 is created is a bifurcation point also. It can be deduced that changing the value of parameter c leads to changing the stability of the zero tumor case only. For other equilibrium points, changing c is not practically effective because other points are unstable. So, surgery is slight effective.
Bifurcation analysis performed in this section shows that stability of the tumor can be controlled by changing patient's parameters. It is important to perform this analysis before deciding the treatment strategy and thus avoiding tumor recurrence.

Conclusions
In this study, a new mathematical model has been proposed in the form of a system of nonlinear ordinary differential equations to study the interaction between tumor cells and some immune cells. Effects of natural killer cells, circulating lymphocytes, CD8 + T cells, and CD4 + T cells on cancer cells have been investigated. Traditionally, CD4 + T cells are considered helper cells that activate CD8 + T cells. However, in the current study, the indirect role of CD4 + T cells in killing tumor cells by secreting cytokines has been considered. Due to stiffness of the system of nonlinear ordinary differential equations, the model has been solved numerically using the method of Adams predictor-corrector with Adams-Bashforth predictor and Adams-Moulton corrector. e Adams method is usually used to solve a single nonlinear ordinary differential equation; however, in the current study, it is used to solve a system of nonlinear ordinary differential equations.
ree cases have been studied: first, the no-treatment case where interaction between immune cells and tumors with different sizes has been discussed. e results reveal that the ability of the patient's body to overcome tumor cells depending upon the patient's parameters, initial tumor size, and initial immune cells populations. For large tumors, the immune system alone fails to overcome the tumor and needs the intervention of treatment. Second, the case of continuous immunotherapy is considered where the role of CD4 + T cells has been found to be very important and significant. In some cases, adoptive CD4 + T cell therapy can be an alternative to combined CD8 + T cells and cytokine therapy. In other cases, combined CD8 + T cells and cytokine therapy is not enough to eliminate the tumor and introducing CD4 + T cells is necessary.
ird, the case of pulsed chemotherapy and immunotherapy is considered where CD4 + T cell therapy has shown to be an alternative to chemotherapy in some cases. In other cases, specific pulses of chemotherapy alone could not be enough for treatment and immunotherapy could be introduced in order to completely eradicate the tumor. For large tumors, previous models abolished the cancer cells using intensive combined chemotherapy, CD8 + T cells, and cytokine. In the current study, the same tumor has been abolished with a lower dose of chemotherapy combined with CD4 + T cell therapy. Stability analysis has been performed for the mathematical model. A condition for the stability of the zero tumor point has been deduced, and it depends on some immunological parameters of the patient. e stability analysis has answered the question of whether the tumor can recur after reaching size zero or not. If a − (ceα)/(fβ) < 0, then the zero tumor state will be a stable state and the tumor cannot recur.
e results reveal that all equilibrium points are unstable for both the first and the second patients. is means that the size of the tumor is uncontrollable; i.e., it could increase without a limit. Also, after reaching the zero tumor equilibrium point using treatment, recurrence of the tumor can occur unless every single tumor cell is killed. We would like here to emphasize that the previously mentioned conclusions about the three studied cases would not be disturbed by the fact that the equilibrium points are unstable. For instance, it is still valid that CD4 + T cells can replace chemotherapy or at least reduce its amount and, hence, spare patients the complications of chemotherapy-induced toxicities and side effects.
Finally, bifurcation analysis has been performed to show the effect of changing patient's parameters on the stability. Bifurcation analysis has been considered for parameters a and c with the data for the first patient and the second patient. It has specified the values of a and c for the system to be stable, and hence, tumor cannot recur after reaching the zero tumor case. e analysis has also shown that, for the second patient, there is a nonzero small amount of tumor case which is stable. Bifurcation analysis provides a powerful tool to control the stability of the tumor by changing patient's parameters. e proposed model offers a powerful method for predicting the body response to immunotherapy and/or chemotherapy. us, the model provides a valuable tool for planning patients' treatment intervention strategies.

Data Availability
Previously reported human data were used to support this study and are available in Table 1. ese prior studies are cited at relevant places within Table 1.