MATHEMATICAL AND NUMERICAL ANALYSIS OF A MATHEMATICAL MODEL OF MIXED IMMUNOTHERAPY AND CHEMOTHERAPY OF CANCER

. In this study, a previously published mathematical model of mixed immunotherapy and chemotherapy of tumors is considered. The stability anal- ysis of the tumor-free equilibrium obtained in the previous study of this model is ﬂawed. In this paper, a suitable analysis is performed to correct this error, and the parameter conditions for the stability of the tumor-free equilibrium are obtained. The stability condition gives an indicator of the host’s ability to ﬁght a cancer. The parameter conditions are examined using experimental data from clinical studies to show that the immune system is able to control a small tumor, and the host’s ability to ﬁght a cancer depends on individual variation. A numerical method based on the continuation technique is employed for one-parameter bifurcation analysis of the mathematical model with periodically pulsed therapies. The unstable ﬁxed point curve provides a good approximation of the maximum tumor burden as a function of the dosage. Chemotherapy-induced lymphocyte damage, which may cause treat- ment failure, is observed in the numerical simulation. The numerical method also produces a set of combined chemotherapy and immunotherapy dosages from which an eﬃcient and safe combination of dosages can be determined.


1.
Introduction. Multi-pronged approaches have become exciting new strategies for treating tumors. Chemotherapy and immunotherapy are among the most common options for cancer treatment. Chemotherapy is the use of drugs to kill a tumor, and immunotherapy is a treatment that stimulates the immune system in fighting cancer. Unfortunately, chemotherapy destroys not only the cancer cells but also the surrounding healthy cells; additionally, more effective drugs are likely to be more toxic. Immunotherapy is ineffective in patients with solid tumors or advanced cancer [4,18,33]. Several research studies have suggested that an appropriate combination of chemotherapy and immunotherapy produces a synergistic effect [25,27,29]. Other clinical studies have used a combination of chemotherapy and immunotherapy to treat patients with different types of cancers [13,32,33], and these research studies have found that the success of combined therapy strongly depends on the choice of the drugs used and the administration schedule.
Mathematical modeling of tumor growth is a powerful tool to understand, predict, and improve the outcome of a treatment regimen. Many mathematical models of tumor growth use nonlinear ordinary differential equations (ODEs), for example [5,6,7,9,10,11,21,22,24,28,30,31,37]. This work is motivated by a mathematical model developed by de Pillis et al. [9] which is more realistic, compared with the other models above. The mathematical model is a system of ODEs governing the tumor growth on a cell population level with a combination of immunotherapy and chemotherapy treatments. The development of this mathematical model was based on the experiments conducted by Diefenbach et al. [12] using EL4 (a thymoma), B16-BL6 (a melanoma), and RMA (a T-cell lymphoma) tumor cells.
The tumor-free case is the basic configuration of the mathematical model because it determines whether the immune system is able to control the growth of a small tumor without treatment and whether the tumor will relapse after attempts to remove the main tumor mass [30]. Evidence suggests that many of us have microscopic tumor cells that remain asymptomatic and never develop and progress to become large and lethal [1,16,17]. At this stage, a balance between the immune system and the tumor cells could result in long-term persistence of tumor dormancy. This implies that the tumor-free equilibrium may be locally stable.
In the work by de Pillis, the tumor-free equilibrium and its stability were carried out using local stability analysis. Several bifurcation diagrams were presented to study the relationship between the tumor growth and model parameters. Unfortunately, the results obtained by de Pillis et al. [9] are not correct because the mathematical model is not differentiable at the tumor-free equilibrium, and thus the linearization at the tumor-free equilibrium does not exist. In this paper, we propose other appropriate analysis techniques to prove the stability property of the tumor-free equilibrium.
The dosage and administration schedule for cancer therapy are important determinants for its potential success. Overdose may cause damage to the patient while underdose may result in residual tumor cells. An optimal dosage for a treatment is to maximize the effect and minimize the damage to the patient. Numerical simulations using data from mouse experiments [12] and data from patients with metastatic melanoma [13] were performed by de Pillis [9], and some significant results have been drawn. Bistability existed, meaning that the long-term tumor state depended on its initial size. Immunotherapy combined with chemotherapy was more effective than single-modality treatments. The combination therapy with vaccination altered system parameter values and enhanced the ability to eliminate a tumor. However, these numerical simulations were carried out using ODE integration with selected dosages and initial tumor sizes. It is not clear whether these results apply to other sets of parameter values.
Bifurcation analysis is a powerful tool to study the changes in dynamics as parameter values vary. Recently, numerical methods for bifurcation analysis have been constructed for ODE systems with periodically pulsed therapies [38,39,40]. In this paper, we propose the numerical method developed in [38] for one-parameter bifurcation analysis, which is an efficient technique to approximate the initial tumor size as a function of the dosage and to predict the outcome of a treatment over a range of parameter values.
To provide details, the mathematical model and mathematical analysis of the stability of the tumor-free equilibrium are given in Sec. 2. The basic mathematical results of the system with periodically pulsed therapy are presented in Sec. 3. Some numerical examples and a discussion are given in Sec. 4. Section 5 provides a brief conclusion.
2. Mathematical model and mathematical analysis.
2.1. Mathematical model. The dimensional form of the mathematical model proposed by de Pillis et al. [9] is as follows: The state variables T (t), N (t), L(t), C(t), M (t), and I(t) represent the tumor cell population, NK cell population, CD8 + T cell population, number of circulating lymphocytes, chemotherapy drug concentration in the bloodstream, and immunotherapy drug concentration in the bloodstream, respectively. The physical meanings of the parameters and the terms in Eqs. (1)-(7) were explained in detail by de Pillis et al. [9,10]. Eqs. (1)-(7) are non-dimensionalized using the following transformations: T * = bT, N * = a 2 N/(αe), L * = bL, C * = aC/α, M * = M, I * = I, t * = at, D * = D/a, c * = cαe/a 3 , d * = d/a, e * = e, f * = f /a, g * = g/a, h * = hb 2 , j * = j/a, k * = kb 2 /a 2 , l * = l, m * = m/a, p * = p/(ab) q * = q/(ab), r * 1 = r 1 αe/a 3 , r * 2 = r 2 α/a 2 , s * = s, u * = uαe/(a 3 b), Dropping the stars for notational clarity, the nondimensional form becomes the following:Ṫ The interactions between tumor cells and effector cells are often modeled using the linear response for its simplicity, for example [6,11,21,24,30]. Eq. (1) employs the linear response for tumor lysis by the NK cells. For the term involving tumor lysis by CD8 + T cells, the experimental data from [12] indicates that the percentage of lysis is a function of the ratio of CD8 + T cells to tumor cells and the percentage does not exceed a maximum [10]. Thus, the function D in Eq. (14) is of the rational form involving L/T . Note that D is not defined for T = 0. Here, we define DT = 0 at T = 0 so that DT becomes continuous at T = 0, and this continuity property is enough for the analysis, which will be presented later in Sec. 3, of the stability property of the tumor free equilibrium. Other interaction terms in Eqs. (2) and (3) have also been justified by experimental data [9].
Due to different choices of the terms in mathematical models of tumor growth, the tumor free equilibrium is never stable in some studies, for example [5,6,11,22]. This means the body is not able to control a tumor of any size. However, the study by Folkman and Kalluri [17] suggests that the tumor free equilibrium be stable under some parameter conditions. Although other choices of the terms in the mathematical models result in a stable tumor free equilibrium under some parameter conditions [10,21,24,30], these conditions have not been justified by experimental data.
In the work by de Pillis et al. [9], the numerical simulation using Eqs. (1)-(7) and human data [13] has shown that the patients were able to control a small tumor in the absence of treatment. Furthermore, the stability condition was obtained in their study using the linearization at the tumor-free equilibrium. However, the parameter values from the above human data [13] did not satisfy the stability condition. The results from the local stability analysis and numerical simulation did not agree because the linearization at the tumor-free equilibrium did not exist. Thus, we are motivated to analyze the stability property of the tumor-free equilibrium using an appropriate mathematical approach.

2.2.
Stability of the tumor-free equilibrium. In this subsection, we consider Eqs. (8)- (11) with Eq. (14), and let v L (t) = 0, M (t) = 0, and I(t) = 0 for all t > 0. The tumor-free equilibrium is given by (T, N, L, C) = (0, 1/f β, 0, 1/β). Since C is independent of T , N , and L, it may be assumed that C(t) is kept at its equilibrium value 1/β. Thus, the ODE system considered in this subsection becomes as follows: where C = 1/β. The tumor-free equilibrium of Eqs. (15)- (17) is denoted by x 0 = (0, 1/f β, 0). In the following two lemmas, we will give the condition under which a tumor will never grow to a size twice as large as its initial size. The proofs of Lemmas 2.1 and 2.2 are given in Appendix A.
Lemma 2.2. Let the notations and assumptions in Lemma 2.1 be adopted in this lemma. If T (0) = , then T (t) ≤ 2 for all t > 0.
Using Lemmas 2.1 and 2.2, we will prove that the tumor-free equilibrium is stable.
Theorem 2.3. Let the notations and assumptions in Lemma 2.1 be adopted in this theorem. Then, the tumor-free equilibrium is locally stable.
The proof of Theorem 2.3 is given in Appendix A.  Table 1 shows the parameter values estimated by de Pillis et al. [9] using two sets of human data in [13], which is referred to as patients 9 and 10. These parameter values satisfy the stability condition,kt 0 > 2m, in Theorem 2.3, indicating that the tumor-free equilibrium for both patients is locally stable. Remark 1. The quantitykt 0 /(2m) in the stability condition not only shows whether the body is able to fight a small tumor without treatment but also indicates the strength of the body itself in fighting a cancer.

Rewriting the indicatork
we have found that the strength of the body itself in fighting a cancer depends on the following factors. First, the parameterm = (s/(d − 1)) 1/l indicates the effectiveness of tumor lysis by CD8 + T cells. An individual with a smallerm value is able to kill tumor cells more effectively. Second, r 2 /β is related to the production of the CD8 + T cells where r 2 is the production rate. Recall that the number of circulating lymphocytes equals 1/β at the tumor-free steady state. This implies that an individual with a smaller β value has a stronger immune system and is able to produce more CD8 + T cells to fight against the tumor cells. Third, the death rate m is also a factor that affects the number of active CD8 + T cells. A large m value causes the CD8 + T cells to die out quickly and may result in an insufficient number of CD8 + T cells. Both patient 9 and patient 10 have similarm and r 2 /β values. This means that they have similar production of CD8 + T cells and similar effectiveness of tumor lysis. However, Patient 10 has a very large m value compared with patient 9. The CD8 + T cells in patient 10 die out very quickly. We havekt 0 /(2m) = 21.6 for patient 9 andkt 0 /(2m) = 1.9 for patient 10. It is expected that patient 9 is stronger than patient 10 in fighting a cancer. This agrees with the numerical simulation performed by de Pillis et al. [9] in which patient 9 is able to kill a tumor of 10 6 tumor cells while patient 10 is able to kill a tumor of 10 5 tumor cells.
3. The τ -shift map and one-parameter bifurcation analysis. The drug administration is usually modeled by applying an external input using the Dirac delta function or the step function. The functions v M , v L , and v I in Eqs. (8)-(14) take the following forms: where δ is the Dirac delta function, r is the number of doses within a period τ , and d • is the dosage given at nτ where 0 ≤ t 1 <t 1 < t 2 <t 2 < · · · < t r <t r < τ . Figs. 1(a) and (b) show examples of Eqs. (20) and (21), respectively, for two periods. We define the time-τ map, F , by where (T (t), N (t), L(t), C(t), M (t), I(t)) is the solution to Eqs. (8)- (14) with the initial condition The boundedness of F is proven in the next theorem.
and v I (t) are defined in the forms of Eqs. (20) or (21).
The proof of Theorem 3.1 is given in Appendix A.
One-parameter bifurcation analysis is the trace of a smooth fixed-point curve. Using the dosage of a treatment regimen as a bifurcation parameter, a fixed-point curve can be used to approximate the maximum tumor burden that a treatment regimen is able to eliminate [38,40]. In [38], one-parameter bifurcation analysis was performed for a mathematical model of superficial bladder cancer with weekly Bacillus Calmette Guerin therapy. The ODE system is turned into a discrete dynamical system using a τ -shift map. Then, the numerical method for one-parameter bifurcation analysis of discrete dynamical systems constructed in [23] can be directly applied to the τ -shift map. We refer to [38] for the detailed procedure and implementation. Note that because the drug administration schedule in Eqs. (8)-(14) is much more complicated than that of the mathematical model studied in [38], the Jacobian matrix of the τ -shift map needs to be calculated with care. 4. Numerical simulations and discussion. In the following, we are most interested in whether a mixed therapy can significantly reduce the dosage of chemotherapy drugs. Therefore, bifurcation diagrams for treatments using only chemotherapy and treatments using mixed therapy will be computed. A chemotherapy treatment is often given with courses repeated every three or four weeks [8,20,36]. For a mixed therapy, many clinical studies use a regimen in which chemotherapy starts a week before immunotherapy, and it is followed by an intravenous infusion of TIL and a bolus dose of IL2 once every 8 hours [3,13,14,15,34,35]. In our simulation, a cycle of three weeks, which corresponds to τ = 9.051, will be employed. Chemotherapy starts from day 1 to day 7, and it is then followed by an intravenous infusion of TIL on day 8 and a bolus dose of IL2 once every 8 hours from day 8 to day 10. The parameter values are shown in Table 2 where t i = 3.017 + 0.144i, i = 0, 1, . . . , 8. Fig. 2 shows an example of the treatment schedule for one period.
Several d M values less than 14.86 are selected to study the maximum tumor burden that a given cancer regimen is able to eliminate. At each of the selected d M values, the numerical integration of Eqs. (8)- (13) is performed using the initial condition x s0 = (T 0 , 119, 10 −7 , 34.48, 0, 0) with different T 0 values. The critical valueT 0 is identified so that an initial condition x s0 with T 0 <T 0 is attracted to the tumor-free fixed point and an initial condition x s0 with T 0 >T 0 is attracted to the high-tumor fixed point. TheT 0 values at these selected d M values are plotted with circles in Fig. 3(a). Then, the numerical integration is repeated at the selected d M values for a weakened immune system with the initial condition In both examples, the circles are located above the diamonds indicating that a healthy immune system is able to eliminate a larger tumor than a weakened immune system. A patient with a healthy immune system should not employ chemotherapy if the initial tumor size is smaller than T 1 . When the initial tumor size is greater than T 1 , the chemotherapy dosage obtained from the unstable fixed-point curve is sufficient for treating the tumor. A patient with a weakened immune system should employ a chemotherapy dosage d M = d 1 if the initial tumor size is less than T 2 . A dosage slightly greater than the dosage obtained from the unstable fixed-point curve should be used if the initial tumor size is greater than T 2 .
In the following example, the same treatment schedule is applied to patient 10. Using the same treatment schedule allows us to study the effect of interindividual variability. One-parameter bifurcation diagrams at d L = 0 and d L = 100 are shown in Figs. 4(a) and (b), respectively. The parameter values for patient 10 are shown in Table 2. Patient 10 has a lower ability to eliminate tumor cells than patient 9, and thus treating patient 10 is more difficult than treating patient 9. This is also shown in Figs. 3(a) and (b) and Figs. 4(a) and (b). Fig. 4(a) shows that chemotherapy alone is not able to control a small tumor unless high dosages of drugs, d M > 19.8, are administrated. The unstable fixed-point curve is plotted using a log scale for the tumor size as shown in Fig. 4(c). Recall that the right hand side of Eqs. (1)- (7) is not differentiable at T = 0. The numerical method employed in this paper uses the continuation technique and the differentiation of the right hand side of Eqs. (1)- (7) is required [23,38]. Therefore, the numerical results are not reliable when T is close to zero. Here, we show only the part of the curve with T > 10 −8 in Fig. 4(c). The unstable fixed-point curve decreases  very quickly near d M = 0.0, and chemotherapy actually becomes a detriment to treatment. Clinical studies reported that chemotherapy-induced tumor remission is often combined with chemotherapy-induced lymphocyte damage, which influences the duration of chemotherapy-obtained remission and causes tumor recurrence [2,19,26]. Recall that patient 10 is able to eliminate a tumor burden of 10 5 cells without treatment. However, chemotherapy-induced damage to the immune system may, on the contrary, help the tumor escape.
To show the possible outcome of the treatment for the part where the unstable fixed-point curve does not exist, we rely on the numerical integration of Eqs. (8)- (13). In Fig. 5, we set v M = 1 and use initial condition (T (0), N (0), L(0), C(0)) = (10 −9 , 178, 10 −7 , 51.72) which represents a healthy immune system for patient 10 with one tumor cell. Fig. 5 shows that the use of chemotherapy drugs reduces the number of circulating lymphocytes. This causes the decrease in the recruitment of CD8 + T cells and the ratio L(t)/T (t), and thus the tumor grows. Although L(t) starts to increase as T (t) increases, L(t)/T (t) remains a small value with time and the tumor continues to grow.
From Figs. 4(a) and (c), a treatment with chemotherapy alone is not recommended for patient 10. Strategies to boost the immune system subsequent to chemotherapy are needed to sustain remission. Recall that an intravenous infusion of TIL on is given on day eight after the treatment of chemotherapy. Fig.  4(b) shows that the treatment becomes effective if d L = 100. Fig. 4(d) shows the unstable fixed-point curve plotted using a log scale for the tumor size. From Fig.  4(d), we propose that a patient with a healthy immune system uses immunotherapy alone with d L = 100 for a tumor size T < 10 −3 and a combination therapy with the chemotherapy dosage obtained from the unstable fixed-point curve for a tumor size T > 10 −3 . For a patient with a weakened immune system, the part of the unstable fixed-point curve with d M > 1.4 can be used to estimate the chemotherapy dosage. The above numerical simulations show that combination therapy allows the usage of lower dosages of each therapy drug and is more effective. Finally, we will investigate the best range of combined dosages for treating a tumor. Fig. 6 shows combined dosages, d M and d L , that are sufficient for treating a tumor burden of 10 8 cells, T 0 = 0.1, using the data of patient 9. A large dosage, d M > 14.9, is needed if d L = 0, and a treatment is unlikely to be successful if d M < 1.2. High dosages sometimes cause undesirable side-effects, and low dosages result in ineffective treatment. Fig. 6 also shows the trade-off effect between both dosages. From Fig. 6, we are able to determine the combined dosages that are successful and able to prevent the patient from harm due to therapy. The combined dosages (d M , d L ) ∈ [1.25, 2] × [0, 400], located on the curve of Fig. 6, are preferable for their efficacy and safety because they contain a small dosage of each cancer therapy.

5.
Conclusion. This paper, which is based on a mathematical model first developed by de Pillis et al. [9], studies mixed immunotherapy and chemotherapy of cancer. Suitable analytical and numerical techniques have been proposed to study their mathematical model. From this study, we have found that some phenomena, which were not originally included in their study, match clinical observations.
In this paper, we first prove that the tumor-free equilibrium is stable if the quantitykt 0 /(2m) is greater than one. This quantity can be used to indicate the ability of a host to kill a tumor without treatment and the degree of response to a treatment. An individual with a largerkt 0 /(2m) value is able to eliminate a larger tumor. The value of the indicator mainly depends on the effectiveness of tumor lysis by the CD8 + T cells, the number of the active CD8 + T cells, and the strength of the immune system. This stability condition is demonstrated with two patients using clinical data [13]. Both patients are able to fight a small tumor without treatment.
Periodically pulsed immunotherapy and chemotherapy are considered for the treatment terms in the mathematical model. The boundedness of solutions with periodically pulsed therapies is proven by mathematical analysis. Bifurcation diagrams show that the unstable tumor fixed-point curve is a good approximation of the maximum tumor burden as a function of the dosage. The bifurcation diagrams also show that (i) a treatment is more likely to be successful for an individual with a largerkt 0 /(2m) value; (ii) the use of chemotherapy drugs may cause treatment failure because it damages the immune system and consequently results in a low production of CD8 + T cells; and (iii) a combination therapy significantly reduces the dosage of the chemotherapy drug and consequently improves the safety of the treatment. The numerical results agree with clinical studies that have suggested that the combination therapy using chemotherapy and immunotherapy may overcome the chemotherapy-induced lymphocyte damage and also improve the outcomes of treatments [2,19,26].
Finally, a curve representing the combined dosages of chemotherapy and immunotherapy drugs for treating a tumor of a given size is computed. From this curve, we are able to determine the desired dosage of each therapy so that the treatment is successful, safe, and efficient. It is expected that these results could be applied to clinical practice. Further studies including clinical tests are needed to verify these results.