www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Optimal control analysis of combined chemotherapy-immunotherapy treatment regimens in a PKPD cancer evolution model Anuraag Bukkuri University of Minnesota bukku001@umn.edu Received: 21 September 2019, accepted: 13 February 2020, published: 2 March 2020 Abstract—The author constructs a mathematical model capturing tumor-immune dynamics, incor- porating the evolution of drug resistance, PKPD (pharmacokinetics and pharmacodynamics) of ad- ministered drugs, and immunotherapy possibilities. Numerical simulations are performed to analyze the model under a variety of treatment possibilities. A sensitivity analysis is performed to determine the parameters contributing the most to the variance in effector cell, resistant, and sensitive tumor cell pop- ulations. Then, a detailed optimal control analysis is performed, along with a numerical simulation of optimal treatment profiles for a hypothetical patient. Keywords-Optimal Control, Oncology, Chemother- apy, Immunotherapy, Pharmacokinetics, Pharmaco- dynamics MSC: 92C45, 92C40, 92C50, 34H05 I. INTRODUCTION Mathematical modeling has the potential to sig- nificantly impact the treatment of cancer through optimization of drug administration schedules. Currently, clinical treatments are administered in an ad hoc way, modifying the treatment along the way, after observing how the patient is responding to it. Since it is not possible to clinically test all possibilities for drug dosage, combinations, sequence, timing, and duration for a patient, math- ematical modeling is able to help determine the optimal treatment profile for each patient, leading to truly personalized cancer treatment [Rockne et al. 2019]. In recent years, immunotherapy has become a viable and promising treatment option for cancer patients. This biological treatment focuses on us- ing the immune system to fight cancerous cells. Though the immune system can naturally pre- vent/slow cancer growth, cancer cells have evolved many ways to bypass the body’s immune system such as having genetic mutations which make them harder to detect by the immune system, having surface proteins which deactivate immune cells, and changing the healthy cells surrounding the tumor so that they interfere with the immune response to cancer. Immunotherapies elicit or am- Copyright: c© 2020 Bukkuri. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment regimens in a PKPD cancer evolution model, Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 1 of 12 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... plify the immune response in these cases to fight against the cancer more effectively in many ways. For example, immune checkpoint inhibitors block immune checkpoints, leading to a stronger immune response, adoptive cell therapy boosts the natural ability of T cells to fight against cancerous ones, and monoclonal antibodies help the immune cells recognize cancerous cells more accurately [Kruger et al. 2019]. Chemotherapy is another common cancer treat- ment plan. Like immunotherapy, there are many drugs with different mechanisms chemotherapy subsumes including doxorubicin, novantrone, and epirubicin. Each of these drugs targets and at- tempts to kill cells which are growing and di- viding rapidly, like those characteristic of cancer. However, a marked side effect of this is that the drug also kills many healthy cells in the process, especially those fast-growing cells in the bone marrow, hair, and skin [Schirrmacher 2019]. Current medical literature is focused on attempt- ing to combine these two treatments. Though much work has been done in this area, the question still remains: what is the best way to combine these treatments to ensure tumor remission while also attempting to minimize side effects? This paper attempts to address this question in a few ways. In section 2, a mathematical model capturing the key components of our system is presented. In section 3, parameter value estimates will be provided and simulations with varying amounts of continuous infusion chemotherapy will be presented. In section 4, a sensitivity analy- sis will be performed, allowing us to understand which parameter values impact the variance in the effector cell and tumor cell populations the most. In section 5, the author analyzes the existence and characterization of the optimal control. In section 6, numerical simulations depicting the optimal control profile under different treatment side effect combinations for a hypothetical patient will be provided. II. MODEL CREATION The model presented below can be thought of to consist of two parts: the first three equations de- scribe tumor-immune interactions. Separate equa- tions are given for sensitive and resistant (to the chemotherapy drug) tumor cells, and a term is pro- vided which converts the sensitive cells to resistant cells. For greatest generality, immunotherapy is included by adding to the effector cell population. The last two equations form the second part of the model, capturing the pharmacokinetics and pharmacodynamics of chemotherapy drug admin- istration. This component is linked to the first set of equations via the effect of the chemotherapy drug on the sensitive tumor cell population. Below is our model: dS dt =aS ( 1− TT Cmax ) −qES−S ( C γ 2 C γ 2 +wIC50 ) −lS −dTS (1) dR dt =aR(1− TT Cmax )−qER+lS−dTR (2) dE dt =s+pE TT g + TT −mETT−µEC1−dEE +s2v(t) (3) dC1 dt =−(k1+k2)C1+ s1u(t) V1 (4) dC2 dt =k12 V1 V2 C1−k2C2 (5) The terms u(t) and v(t) are termed controls for our treatment: u(t) refers to the chemotherapy treatment, and v(t) refers to the immunotherapy treatment. In this model, S and R represent the drug- sensitive and drug-resistant tumor cells, TT rep- resents the total number of tumor cells (S + R), E represents the effector cells, and C1 and C2 represent the drug concentrations in the plasma and at the tumor site, respectively. Below, we will explain the functional forms of each of the equations above: dS/dt: the first term represents the normal dynamics of these cells, taking into account the carrying capacity (captured by Cmax). The second term captures the killing rate of tumor cells by effector cells. The third term is the killing of sensitive tumor cells by the drug; IC50 is the Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 2 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... median inhibitory concentration, and w scales between the measurements in vitro to in vivo (we can assume this is 1); γ is provided to help scale the hyperbolic curve for PKPD data fitting purposes. The fourth term is used to incorporate drug resistance–l represents the mutation rate from sensitive to resistant tumor cells. Finally, the last term captures the natural death rate of the cancer cells. dR/dt: The first, second, and last terms serve the same purpose as in the dS/dt equation. The lS captures the transferal of sensitive to resistant cells via mutation. Note that we assume identical base growth dynamics for the sensitive and resistant cells, as well as identical interactions with effector cells. dE/dt: The first term is the rate of influx of effector cells into the region of tumor localization and the second term captures the tumor activation of effector cells while ensuring a maximum rate at which effector cells are produced and dE is their death rate. The third term is the death of effector cells due to interaction with the tumor cells. The fourth term is added to incorporate drug toxicity. Finally, the last term, s2v(t), captures the administration of the immunotherapy treatment, in which v(t) is the input function ranging from 0 to 1 and s2 is a scaling term. Since each immunotherapy treatment mechanistically works differently, the generalization made here is that the immunotherapy effectively increases the number of effector cells at the tumor site. dC1/dt: Here, k1 is the elimination from the plasma, and k2 is the elimination from the effect compartment. u is the input function (ranging from 0 to 1) for administration protocol of the chemotherapy treatment (with s1 serving as a scaling factor) and V1 is the volume of distribution of the compartment. dC2/dt: k12 is the link process between the plasma component and the tumor, V1 and k2 are as above, V2 is the effect component for C2(t). Together, these last two equations incorporate the pharmacokinetic components of drug administra- tion. III. IMPACTS OF CONTINUOUS INFUSION OF CHEMOTHERAPY The parameter values used in model simulations are given in Table I. The most standard clinical treatment for cancer is chemotherapy administration. As such, in this section, numerical simulations are performed for various chemotherapy treatment situations. Initial conditions of 300 effector cells and 10,000 sensi- tive cancer cells and an s1 value of 100 were used for all simulations below. Three simulations were performed for varying intensities of chemotherapy treatment. The results can be seen in Figure 1. Fig. 1. Chemotherapy Numerical Simulations Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 3 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... Parameter Meaning Estimated Value Source a Natural Growth Rate of Cancer Cells 0.616/day dePillis Cmax Carrying Capacity of Cancer Cells 9.804 ∗ 108cells dePillis q Killing Rate of Cancer Cells by Effector Cells 1.101 ∗ 10−7/cell ∗ day Lai estimation w Scaling Between In Vitro and In Vivio 1 Assumption IC50 Median Drug Inhibitory Concentration 10 Depends on drug γ Hyperbolic Curve Modification 1 N/A l Conversion Rate of Sensitive to Resistant Cancer Cells 10−9 Birkhead dT Natural Death Rate of Cancer Cells 0.17/day Liao; Lai dE Natural Death Rate of Effector Cells 0.0412/day Kuznetsov s Effector Cells at Tumor Site 1.3 ∗ 105cells/day Kuznetsov p Tumor Activation of Effector Cells 0.1245/day Kuznetsov g Maximum Rate of Effector Cell Production 2.019 ∗ 107cells Kuznetsov m Death of Effector Cells due to Cancer Cells 3.422 ∗ 10−10/cell ∗ day Kuznetsov µ Side Effect of Chemotherapy on Effector Cells 2 ∗ 10−2 liters/g*day Estimation k1 Elimination from Plasma 1.6/day Iliadis k2 Link Process between Plasma and Tumor 0.8/day Iliadis u Input function for Chemotherapy Protocol Depends on control N/A v Input function for Immunotherapy Protocol Depends on control N/A V1 Volume of the Plasma Compartment 25 Iliadis V2 Volume of the Tumor Compartment 15 Iliadis k12 Link Process between Plasma and Tumor Compartments 0.4/day Iliadis TABLE I PARAMETERS USED IN NUMERICAL SIMULATIONS In this figure, the top panel corresponds to a u value of 0.2, or 20% of the maximal chemotherapy administration possible. The middle panel repre- sents 50% and the bottom panel is 100%. There are a few things to notice in these figures. First, consider the key similarity: in all three treatment scenarios, the effector cell population reached the same equilibrium of 4∗106 cells at similar times. Note that, in these simulations, a low side effect of the chemotherapy drug on the immune system was assumed. The critical differences lie in the populations of sensitive and resistant cells. We see a common trend in that the higher the dose of chemotherapy, the faster the sensitive cancer cell population goes extinct, but the larger the final equilibrium of resistant cancer cells is. Since, in this model, we are dealing with a rapidly evolving cancer cell population, regardless of the intensity of continuous infusion chemother- apy, resistance is ensured to occur, without com- plete remission of the cancer. Thus, to more rig- orously determine which aspects of the cancer to attack, a sensitivity analysis was performed. IV. SENSITIVITY ANALYSIS A sensitivity analysis was performed on the above ODE model using the Sobol-Martinez method. Parameters which cannot be modified with treatments were given a range 10% higher and 10% lower than the original value used in the original numerical simulations. The Sobol- Martinez method is based on variance decomposi- tion techniques to provide a rigorous measure of the contributions of the input to output variance. The algorithm outlined by Zhang et al. was imple- mented for the sensitivity analysis given in Figures 2, 3 and 4. The images on the left are the first order Sobol indices: the contribution to the output variance by the single model input alone. The second image is the total Sobol index, one which measures the contribution to the output variance caused by a model input, including the first-order effects as well as all higher-order interactions. As one can see, in the case of the effector cells, the Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 4 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... Fig. 2. Sensitivity Analysis: Effector Cells Fig. 3. Sensitivity Analysis: Resistant Cancer Cells most sensitive parameter value for the first 50 or so time steps is s, a quantity that is increased by immunotherapy treatments such as dendritic cell therapy. However, after that time, until 400 time steps, the µ parameter (side effects of the drug) is just as sensitive, if not more sensitive. Since the side effects of immunotherapy do not directly impact the effector cells usually, this can be interpreted as the side effects of chemotherapy impacting the effector cells even more than the po- tential benefits of the immunotherapy. Medically, from the effector cell side, it would seem that the best way to progress would be to give a high dose of dendritic cell therapy, for example, at the beginning, then work to limit the side effects of other chemotherapy drugs that are given. In the case of the resistant cells, one can see Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 5 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... Fig. 4. Sensitivity Analysis: Sensitive Cancer Cells that all three parameter values shown are quite sensitive. Considering the total Sobol index, it seems that, for most of the time, a (growth rate of the cancer) and q (killing rate of tumor by effector cells) switch off between being the most and second most sensitive parameter values, while the mutation rate (l) becomes more important, relative to the other parameters, as time progresses. This implies that, for the resistant cells, it might be most effective to alternate some treatment that directly reduces the cancer’s ability to grow and immunotherapy for the first 200 time steps, and then just give immunotherapy for the last 200 time steps. However, it’s hard to directly reduce natural tumor growth rate using a drug. There is, though, a burgeoning area of medical research focused on the elimination of the Tudor-SN pro- tein from cancer cells using CRISPR-Cas9 gene editing techniques. It is known that Tudor-SN is the main influencer of cancer cell growth and thus, silencing this protein would suppress the rapid cellular proliferation characteristic of cancers. A recent study showed that when Tudor-SN was removed from human cells, the levels of several microRNAs increase, putting brakes on genes that encourage cell growth. This significantly slowed down cell progression from the preparatory to cell division phase [Elbarbary 2017]. As suggested by the sensitivity analysis of the resistant cancer cell output values, using such a treatment based on the suppression of Tudor-SN in conjunction with immunotherapy would be the most effective treat- ment to reduce the resistant tumor cell population. For the sensitive cells, it’s clear that, at the beginning 50 time steps, a is the most sensitive parameter. This then quickly switches to q being the most important parameter, with the IC50 value of the drug becoming more sensitive as time pro- gresses. Medically, this implies that it’s best to give immunotherapy throughout the treatment, while introducing increasing amounts of chemotherapy over time. Thus, overall, it seems that the parameters a, q, s, and l are the parameters that impact the dynamics the most. Medically, it then makes sense to target these parameter values. Though it is often not possible to target all of these parameters in a treatment plan, our analysis sug- gests the plausibility of a combined chemotherapy- immunotherapy treatment plan. So, to deter- mine what the best combination and timing of chemotherapy-immunotherapy administered proto- cols are, we perform an optimal control analysis. Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 6 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... V. OPTIMAL CONTROL Here, we formulate the problem of determining the most effective treatment regimen as an optimal control one. Specifically, we desire to minimize the amount of drug given and the tumor size over a fixed therapy interval [0,T] while maximizing the number of effector cells, thus maintaining an effective tumor-immune balance as done in [Bukkuri 2019]. Drug toxicity is accounted for by the u(t) term in the integral, denoting the total drug dosage over the given treatment period as a penalty term. We choose as our control class piecewise con- tinuous functions defined for all t such that 0 ≤ u(t) ≤ 1 where u(t)= 1 represents max- imal chemotherapy and u(t)= 0 represents no chemotherapy. Thus, we depict the class of ad- missible controls as U1(t) = u(t) piecewise continuous s.t. 0 ≤ u(t) ≤ 1,∀t ∈ [0,T] U2(t) = u(t),v(t) piecewise continuous s.t. 0 ≤ u(t),v(t) ≤ 1,∀t ∈ [0,T] Now, we define our objective functional and optimal control problem. We specifically consider two cases: one treatment (J1) with just chemother- apy and another treatment (J2) with a combined chemotherapy-immunotherapy (e.g. dendritic cell therapy) regimen. For a fixed therapy horizon [0,T], maximize the objective functional J1(u)= ∫ T 0 αE(t)−β1S(t))−β2R(t)−Bu(t))dt J2(u,v)= ∫ T 0 αE(t)−β1S(t))−β2R(t) −B1u(t)−B2v(t))dt over all Lebesgue-measurable functions u : [0,T] → [0,umax] subject to the above ODE dynamics and initial conditions. Theorem 1: Consider the objective functional J1(u) subject to the state system in section 2. Assume that: • there exists an admissible pair ( ~K,u(t)) • N( ~K,U1,t) is convex in U1 for each ( ~K,t) • U1 is closed and bounded • ∃ a number θ s.t. ||~K|| ≤ θ∀t ∈ [t0, t1] and all admissible pairs ( ~K,u(t)) Then there exists an optimal control pair ( ~K∗,v∗) that maximizes J1(u). Moreover, for suf- ficiently small T, the optimal control system has a unique solution. Proof. An admissible pair ( ~K,v(t)) is needed for the existence of optimal control. Since the system in equations (1) – (5) (our original set of equations) has bounded coefficients and any solu- tions are bounded on the finite time interval, using a result from Lukes [Lukes 1982], we obtain the existence of the solution of the system described by equations (1) – (5). Now, we need N( ~K,U1,t) to be convex in U1 for each ( ~K,t). We define: w1 =(αE −β1S −β2R−Bu1 + γ1, aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2+wIC50 ) −lS−dTS, aR ( 1− TT Cmax ) −nER + lS −dTR, s + pE TT g + TT −mETT −µEC1 −dEE, − (k1 + k2)C1 + s1u1 V1 , k12 V1 V2 C1 −k2C2) for some γ1 ≤ 0 and u1 ∈ U1, w2 =(αE −β1S −β2R−Bu2 + γ2, aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2+wIC50 ) −lS−dTS, aR ( 1− TT Cmax ) −nER + lS −dTR, s + pE TT g + TT −mETT −µEC1 −dEE, − (k1 + k2)C1 + s1u2 V1 , k12 V1 V2 C1 −k2C2) for some γ2 ≤ 0 and u2 ∈ U1. We must now prove Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 7 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... that for every λ ∈ [0,1] we have w3 = λw1 + (1−λ)w2 ∈ N( ~K,U1, t). To do this, let z1 =λ(αE −β1S −β2R−Bu1 + γ1) + (1−λ)(αE −β1S −β2R−Bu2 + γ2) = αE −β1S −β2R−B((1−λ)u2 + λu1) + λγ1 + (1−λ)γ2 and define γ3 = z1 − (αE −β1S −β2R) + Bu3 for u3 = (1−λ)u2 + λu1. Then, γ3 = λγ1 + (1−λ)γ2 ≤ 0, since γ1,γ2 ≤ 0 and λ ∈ [0,1]. Then, we see that z2 =λ ( aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2 + wIC50 ) − lS −dTS ) + (1−λ) ( aS ( 1− TT Cmax ) −nES −S ( C γ 2 C γ 2 + wIC50 ) − lS −dTS ) =aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2+wIC50 ) −lS−dTS z3 =λ(aR(1− TT Cmax )−nER + lS −dTR) + (1−λ) ( aR ( 1− TT Cmax ) −nER+lS−dTR ) = aR ( 1− TT Cmax ) −nER + lS −dTR z4 =λ ( s+pE TT g+TT −mETT−µEC1−dEE ) +(1−λ) ( s+pE TT g+TT −mETT−µEC1−dEE ) = s + pE TT g + TT −mETT −µEC1 −dEE z5 =λ ( −(k1 + k2)C1 + s1u1 V1 ) + (1−λ) ( −(k1 + k2)C1 + s1u2 V1 ) = −(k1 + k2)C1 + s1u3 V1 z6 =λ ( k12 V1 V2 C1 −k2C2 ) +(1−λ) ( k12 V1 V2 C1 −k2C2 )) = k12 V1 V2 C1 −k2C2 Combining this information, we find a v3 ∈ [0,1] and γ3 ≤ 0 such that λw1 + (1−λ)w2 =( αE −β1S −β2R−Bu3 + γ3, aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2+wIC50 ) −lS−dTS, aR ( 1− TT Cmax ) −nER + lS −dTR, s + pE TT g + TT −mETT −µEC1 −dEE, −(k1 + k2)C1 + s1u3 V1 , k12 V1 V2 C1 −k2C2 ) Therefore, λw1 + (1 − λ)w2 ∈ N( ~K,U1, t). Thus, N( ~K,U1, t) is convex in U1. The next requirement for the existence of opti- mal control is that U is closed and bounded. This is true, by definition. Finally, there exists a number θ s.t. ||~K|| ≤ θ∀t ∈ [t0, t1] and all admissible pairs ( ~K,u(t)). The boundedness argument is analogous to those previously performed [Fister et al. 1998, Burden et al. 2004, Fister et al. 2005]. A similar argument can be applied for the existence of an optimal control for J2(u,v). For the analysis, we shall continue with the J2(u,v) objective so that the proofs are more Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 8 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... insightful. Note that J1 is identical to the J2 objective, but with a zero v term. Since we’ve proven the existence of optimal controls to maximize the J1 and J2 functionals, first order necessary conditions for optimality can be determined by a version of the Pontryagin maximum principle. For the characterization of the optimal control, we first define the Hamiltonian associated with J2(u,v) and the system of ODEs as follows: H =αE −β1S −β2R−B1u +λ1 ( aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2 +wIC50 ) − lS−dTS ) +λ2 ( aR ( 1− TT Cmax ) −nER+lS−dTR ) +λ3 ( s+pE TT g+TT −mETT−µEC1−dEE +s2v(t) ) + λ4 ( −(k1 + k2)C1 + s1u(t) V1 ) + λ5 ( k12 V1 V2 C1 −k2C2 ) The existence of an optimal control for the state system given in section 2 associated with the ob- jective J1(u) can be determined from the Filippov- Cesari theorem. For the theorem, the following notation shall be used: ~K =   S R E C1 C2   and N( ~K,U1, t)=( αE(t)−β1S(t)−β2R(t)−B1u(t)−B2v(t)+φ, aS ( 1− TT Cmax ) −nES−S ( C γ 2 C γ 2 +wIC50 ) − lS −dTS, aR(1− TT Cmax )−nER + lS −dTR, s + pE TT g + TT −mETT −µEC1 −dEE, − (k1 + k2)C1 + s1u(t) V1 , k12 V1 V2 C1 −k2C2 ) where φ ≤ 0 and u ∈ U1. Regarding the uniqueness of the controls, simi- lar proofs are given in [Burden et al. 2004, Fister et al. 1998]. The additional constraint that T must be small is due to the fact that the state system is moving forward in time while the adjoint system is moving backwards. Theorem 2: Given optimal controls u∗ and v∗ and solutions of the corresponding state system, there exist adjoint variables λi for i = 1,2,3,4,5 satisfying the following: dλ1 dt = − ∂H ∂S dλ2 dt = − ∂H ∂R dλ3 dt = − ∂H ∂E dλ4 dt = − ∂H ∂C1 dλ5 dt = − ∂H ∂C2 where λi(T) = 0 for i=1,2,3,4,5 by the PMP transversality condition. Furthermore, from the optimality condition, u∗ is given by:{ 0, if s1λ4 V1 −B1 < 0 1, if s1λ4 V1 −B1 > 0 while v∗ is similarly given by:{ 0, if s2λ3 −B2 < 0 1, if s2λ3 −B2 > 0 Proof. From the Hamiltonian, the derivatives of the adjoints were calculated, and the following is seen. u∗ is given by:  0, if s1λ4 V1 −B1 < 0 1, if s1λ4 V1 −B1 > 0 singular, if s1λ4 V1 −B1 = 0 Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 9 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... while v∗ is similarly given by:  0, if s2λ3 −B2 < 0 1, if s2λ3 −B2 > 0 singular, if s2λ3 −B2 = 0 To determine the representation of the control, we first prove that both controls are bang-bang via proof by elimination. If both controls are singular on [0,T] and take time derivatives, then λ3(t) = 0 and λ4(t) = 0 on [0,T] due to the continuity of λ3(t) and λ4(t). This is a contradiction of the fact that λ3 = B2 s2 and λ4 = B1V1 s1 since we’re assuming the presence of a tumor and some side effects of the drug. Thus, it’s not possible for both controls to be singular. Now, let’s analyze the possibility of one of the controls being singular. For simplicity, we will work with the v∗ control, but the same argument applies to the u∗ control. If the v∗ control is sin- gular, λ3 = B2 s2 on [0,T]. Taking a time derivative of s2λ3 −B2 = 0, we get s2λ′3(t) = 0, implying that λ′3 = 0 and that λ3(t) is constant. However, since λ3(t) is continuous on [0,T] and λ3(T) = 0 by the transversality condition, then λ3(t) = 0 for any subset on [0,T]; however, again, λ3 = B2 s2 , a nonzero quantity for any meaningful therapy. Thus, we conclude that both controls must be bang-bang. Note that this agrees with current medical knowledge that chemotherapy is best given in bang-bang controls and hence why chemotherapy is given in cycles: a dose of one or more drugs fol- lowed by several days or weeks without treatment. This ensures that the normal cells have enough time to recover from the drug side effects. A similar argument can be given for immunotherapy, though the side effects are less direct than those for chemotherapy Intuitively, this makes sense since it’s stating that no treatment will be used when the V1 is high (i.e. essentially that the tumor/normal com- partment ratio is high) and the maximum effect of the drug is low compared to the patient’s treatment tolerance level. On the other hand, if the patient can tolerate great side effects, the tumor size is large, and there exists a high drug efficacy, maximum treatment should be given. VI. SIMULATED PATIENTS We can think of our optimality system as a two- point boundary value problem, which we solve using a fourth-order iterative Runge-Kutta scheme, as done in [Jung et al. 2002]. In this scheme, we perform a forward sweep of the state equations with initial guesses for u and v, before performing a backward calculation using the adjoint equation and an update of the controls. This method is done iteratively until convergence is obtained. Below are simulations of a few hypothetical patients. We specifically analyze the cases of low and high side effects for chemotherapy and im- munotherapy. The weighting terms in the objective functional were fixed at α = 0.03, b1 = 0.5, and b2 = 0.6. s1 was fixed at 100 and s2 was fixed at 100,000. The table below shows the values of B1 and B2 used to represent high and low side effects for chemotherapy and immunotherapy: Side Effect Chemotherapy (B1) Immunotherapy (B2) Low 1000 10,000 High 1,000,000 100,000,000 TABLE II B1 AND B2 VALUES USED In Figure 5 depictions of all four possible com- binations of side effects for chemotherapy and immunotherapy over 700 days are presented. Note that, in accordance with clinical oncol- ogy practices, the optimal controls in all side effect scenarios include bang-bang controls for both chemotherapy and immunotherapy. The key differences in optimal control protocols occur in the first 400 days of treatment. In the case of low side effects for both chemotherapy and immunotherapy, it is advised that chemotherapy be delivered for the first 150 days and immunother- apy be delivered for the last 380 days. This is in stark contrast to the other three scenarios, in which chemotherapy and immunotherapy treatments do Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 10 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... Fig. 5. Numerical Simulations of Optimal Control Profiles for a Hypothetical Patient not overlap for the first 400 days. For example, consider low side effects for chemotherapy and high side effects for immunotherapy. In this case, chemotherapy should be prescribed for the first 150 days and immunotherapy should be given for the next 250 days. When chemotherapy has high side effects for the patient, it is effectively removed from the optimal treatment plan (for the first 400 days). Instead immunotherapy is given for the last 380 days (if it has a low side effect) or the last 250 days (if it has a high side effect). Past day 400, optimal treatment protocols in all four cases are to give 70 days of chemotherapy followed by 70 days of a drug vacation and to also give 70 days of immunotherapy followed by 70 days of a drug vacation–note that these two treatments overlap and are often given together. VII. CONCLUSION In this paper, a model was created capturing dynamics among tumor cells sensitive and resistant to chemotherapy, effector cells, and chemotherapy PKPD. Evolution of resistance to chemotherapy was taken into count. Parameter values were esti- mated from medical literature and numerical simu- lations were performed displaying dynamics under control and treatment cases. Then, Sobol-Martinez sensitivity analyses were performed which found the growth rate of cancer cells, killing of cancer cells by the immune system, effector cells in the tumor compartment, and the resistance emergence in cancer cells to be the parameters which impact the dynamics of effector cells, sensitive tumor cells, and resistant tumor cells the most. Next, a characterization of the optimal treatment pro- file was given analytically, along with proofs for the existence and uniqueness of such a protocol. Numerical optimal control was then performed to illustrate the optimal treatment schedule for a hypothetical cancer patient. This showed a profile of bang-bang chemotherapy and immunotherapy controls, often with overlapping treatment regi- mens. The author hopes that the model and anal- ysis will help inspire further medical research to be conducted in creating drugs which reduce the natural growth rate of cancer cells, as well Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 11 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Anuraag Bukkuri, Optimal control analysis of combined chemotherapy-immunotherapy treatment ... as the clinical use of combined chemotherapy- immunotherapy treatments in accordance with the optimal protocols presented. REFERENCES [1] Birkhead BG, Rankin EM, Gallivan S, Dones L, Rubens RD (1987). A Mathematical Model of the Development of Drug Resistance to Cancer Chemotherapy. Eur. J. Cancer Clin. Oncology, 23(9), 1421-1427. [2] Bukkuri A (2019). Optimal control analysis of combined anti-angiogenic and tumor immunotherapy. Open Jour- nal of Mathematical Sciences, 3, 349-357. [3] Burden T, Ernstberger J, Fister R (2004). Optimal Con- trol Applied to Immunotherapy. Discrete and Continu- ous Dynamical Systems, 4(1), 135–146. [4] de Pillis LG, Gu W, Radunskaya AE (2006). Mixed immunotherapy and chemotherapy of tu- mors: Modeling, applications and biological interpreta- tions. Journal of Theoretical Biology, 238(4):841–862. https://doi.org/10.1016/j.jtbi.2005.06.037 [5] Elbarbary RA, Miyoshi K, Myers JR, Du P, Ashton JM, Tian B, Maquat LE (2017). Tudor-SN-mediated endonucleolytic decay of human cell microRNAs pro- motes G1/S phase transition.; Science (New York, N.Y.); Vol 356(6340). [6] Fister R, Donnelly, J (2005). Immunotherapy: An Op- timal Control Theory Approach, Math. Biosci. and En- grg., 2(3), 499-510. [7] Fister R, Lenhart S, McNally J (1998). Optimizing Chemotherapy in an HIV , Elec. J.DE., 32, 1–12. [8] Iliadis A, Barbolosi D (2000). Optimizing Drug Regi- mens in Cancer Chemotherapy by an Efficacy–Toxicity Mathematical Model [9] Jung E, Lenhart S, Feng Z (2002). Optimal control of treatments in a two-strain tuberculosis model, Discrete Contin.Dyn. Syst. Ser., 2, 473-482. [10] Kruger S, Ilmer M, Kobold S. et al. (2019). Advances in cancer immunotherapy 2019 – latest trends. J Exp Clin Cancer Res 38, 268 [11] Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS (1994). Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bull. Math. Biol., 56(2): 295—321. [12] Lai X, Friedman A (2017). Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitors: A mathematical model. PLoS ONE 12(5): e0178479. https://doi.org/10.1371/journal.pone.0178479 [13] Liao KL, Bai XF, Friedman A (2014). Mathematical modeling of interleukin-27 induction of anti-tumor T cells response. PLoS ONE 9(3):e91844 [14] Lukes D (1982). Differential Equations: Classical to Controlled, Math. Sci. Engrg. 162. [15] Rockne RC, Hawkins-Daarud A, Swanson KR, Sluka JP, Glazier JA, Macklin P, Hormuth D, Jarrett AM, da Fonseca Lima EAB, Oden J. et al. (2019). The 2019 mathematical oncology roadmap, Physical Biology 16(4), 041005. [16] Schirrmacher V. (2019). From chemotherapy to biolog- ical therapy: A review of novel concepts to reduce the side effects of systemic cancer treatment. Int. J. Oncol. 54(2), 407-419. [17] Zhang XY, Trame MN, Lesko LJ, Schmidt, S (2015). Sobol Sensitivity Analysis: A Tool to Guide the De- velopment and Evaluation of Systems Pharmacology Models. CPT Pharmacometrics Syst. Pharmacol. 4, 69- 79. Biomath 9 (2020), 2002137, http://dx.doi.org/10.11145/j.biomath.2020.02.137 Page 12 of 12 http://dx.doi.org/10.11145/j.biomath.2020.02.137 Introduction Model Creation Impacts of Continuous Infusion of Chemotherapy Sensitivity Analysis Optimal Control Simulated Patients Conclusion References