An Optimal Control Strategy for Mathematically Modeling Cancer Combination Therapy
AA N O PTIMAL C ONTROL S TRATEGY FOR M ATHEMATICALLY M ODELING C ANCER C OMBINATION T HERAPY
Jayanth S. Pratap
Dulles Math and Science AcademyPresented at the 3 rd Systems Approaches to Cancer Biology ConferenceNovember 11 th - 13 th , 2020. A BSTRACT
While the use of combination therapy is increasing in prevalence for cancer treatment, it isoften difficult to predict the exact interactions between different treatment forms, and theirsynergistic/antagonistic effects on patient health and therapy outcome. In this research, a systemof ordinary differential equations is constructed to model nonlinear dynamics between tumor cells,immune cells, and three forms of therapy: chemotherapy, immunotherapy, and radiotherapy. Thismodel is then used to generate optimized combination therapy plans using optimal control theory.In-silico experiments are conducted to simulate the response of the patient model to varioustreatment plans. This is the first mathematical model in current literature to introduce radiotherapyas an option alongside immuno- and chemotherapy, permitting more flexible and effective treatmentplans that reflect modern therapeutic approaches.I. I
NTRODUCTION
Combination therapy in cancer treatment is currentlya highly active area of research. Immunotherapy,chemotherapy, and radiotherapy are often used togetherin some combination to maximize the likelihoodof tumor regression, but in these cases it is ofparamount importance to consider both the synergisticand antagonistic effects of these drugs when combined.Chemotherapy and radiotherapy often have toxic effectson the patient’s immune system, and immunotherapy canbe used to counteract these effects by stimulating immunesystem activation.Quantitative modeling provides a promising avenuefor analyzing and optimizing the synergistic effectsof different therapies. In this project, the firstphase consists of constructing a model of ordinarydifferential equations to capture interactions betweentumor cells, immune effector cells, and three formsof therapeutic intervention—immuno-, chemo-, andradiotherapy. Biological parameters fitted to murinetrials and patient data are used to simulate the model.The second phase is setting up different optimal controlscenarios in fixed final-time form to minimize differentquantities over the course of treatment, such as averagetumor burden, final tumor burden, and therapy delivered,as well as introducing a term to maximize patient healthin the form of average immune cell population. II. ODE M
ODEL F ORMULATION
The model is comprised of a two-dimensional system,and is modeled around the Lotka-Volterra predator-preysystem, with form ˙ x = αx − βxy ˙ y = δxy − γy, where x represents the prey species and y represents thepredator species.In the system of tumor-immune dynamics, the predatorspecies y is considered as the effector cell population(immune cells with the ability to kill tumor cells), and theprey species x is represented as the tumor cell population. A. ODE Model: General Form of Equations
There are five populations whose dynamics are consideredin the ODE model:• N ( t ) , tumor cell population• E ( t ) , effector cell population• I ( t ) , immunotherapy dose level• C ( t ) , chemotherapy dose level• R ( t ) , radiotherapy dose levelFurthermore, there are several forms of terms that mustbe included in the equations:1 a r X i v : . [ m a t h . O C ] J a n PREPRINT - J
ANUARY
29, 2021• Net growth term ( G N , G E , G I , G C , G R )• Tumor-immune agonism/antagonism ( A E → N , A N → E )• Chemo-therapeutic effects ( C N , C E )• Radio-therapeutic effects ( R N , R E )• Immuno-therapeutic effects ( I E )• Therapy intervention ( u I , u C , u R ) A.i. Growth and death terms
In the absence of therapy, tumor cells ( N ) grow insome fashion limited by in vivo carrying capacity. Thelimited-growth model used in describing this is logisticgrowth, based on murine data. This gives us the form G N = αN (1 − βN ) . Immune cells are assumed to be generated at a constantrate, and to have a natural lifespan until death. This givesus the form G E = ζ − λE. For all three forms of therapy analyzed, we assumean exponential decay function adequately models theelimination of the dose from the body, thus G I = − τ E,G C = − ωC, and G R = − χR. A.ii. Tumor-immune agonism/antagonism
Tumor cells are assumed to be deactivated by circulatingeffector cells in a manner that scales with the probabilityof an effector/tumor cell encounter, thus giving the form A E → N = − γN E. There is an immunostimulatory effect as effector cellscome into contact with tumor cells. This is due tocytokine-induced activation, where chemokines (immunesignaling molecules) and inflammatory markers inducethe increased production and activation of effector cells.This term is modeled with Michaelis-Menten kineticsbecause immunostimulation asymptotically approaches amaximum rate due to limitation in the rate of immune cellcirculation and cell signaling. Thus we use form A N → E = ηENθ + N where η represents maximum rate of activation and θ represents the half-saturation constant, equivalent to thevalue of tumor cells for which the activation rate is half ofits maximum. A.iii. Therapeutic effects
The form of immunotherapy analyzed in this project isIL-2 (interleukin-2) immunotherapy; IL-2 is a naturallyoccurring cytokine in the body, and its effect on immunesystem activation is thus described in the same form as M N → E : I E = νEIρ + I .
Chemotherapy kill terms reflect fractional cell kill:due to chemotherapeutic agents such as doxorubicin,cisplatin, and bleomycin being mostly or exclusivelyeffective towards tumor cells undergoing specific phasesof the cell cycle, there is always only a particularfraction of cells that are killed with each dosage,proportional to the dosage itself. Additionally, thedose-response curves of chemotherapeutic drugs oftendemonstrate bounded efficacy of chemotherapy, and thusa saturation term is used to give near-linear kill rate forlow drug concentrations, with plateauing kill rate forhigher concentrations. With these considerations, we useterms I N = − δN (1 − e − C ) and C E = − µE (1 − e − C ) . Radiotherapy kill terms reflect the linear-quadraticmodel, used in radiation biology to reflect thephenomenon that cell survival fractions after irradiationoften take the form of an exponential function with alinear and quadratic term. We use terms R N = − N ( (cid:15)R + κR ) and R E = − E ( σR + φR ) . The ratios (cid:15)κ and σφ , as per clinical trial data in [4], areset as (cid:15)κ = σφ = 10 . A.iv. Therapy intervention
We allow functions u I ( t ) , u C ( t ) , u R ( t ) to model theinfusion rate of immuno-, chemo-, and radio-therapyrespectively, defined in terms of time. These are intendedto serve as controls to optimize treatment, which will beexamined later in the project. B. ODE Model: Specific Forms of Equations
Considering the developments in the previous section, theODE system has been written as follows: ˙ N ( t ) = αN (1 − βN ) − γN E − δN (1 − e − C ) − N ( (cid:15) ˆ R + κ ˆ R ) (1a) ˙ E ( t ) = ζ − λE + ηENθ + N + νEIρ + I − µE (1 − e − C ) − E ( σR + φR ) (1b)2 PREPRINT - J
ANUARY
29, 2021 ˙ I ( t ) = u I ( t ) − τ I (1c) ˙ C ( t ) = u C ( t ) − ωC (1d) ˙ R ( t ) = u R ( t ) − χR (1e) C. Model Parameters
The parameters used in the project are sourced from various past models which have analyzed data from murineexperiments and human clinical trials for curve-fitting. See below for a table of parameters, their descriptions, and thevalues used. Table 1: Estimated Model Parameter Values
Parameter Description Units Estimated Value Source α Maximum growth rate of tumor cells day − × − [8] β Reciprocal of carrying capacity of tumor cells cells − × − [8] γ Effector-cell-induced tumor death rate cells − day − × − [8] δ Chemotherapy (Bleomycin) kill rate coefficient for tumor cells day − × − [10] (cid:15) Radiotherapy linear kill coefficient for tumor cells Gy − × − [2] κ Radiotherapy quadratic kill coefficient for tumor cells Gy − × − [2] ζ Constant effector cell production rate cells day − × [8] η Maximum cytokine activation rate of effector cells day − × − [8] θ Half-saturation constant for cytokine activation of effector cells cells 2.019 × [8] λ Death rate of effector cells day − × − [8] ν Maximum immunotherapy (IL-2) activation rate of effector cells day − × − [7] ρ Half-saturation constant for immunotherapy (IL-2) activation of effector cells cells 2.00 × [7] µ Chemotherapy (Bleomycin) kill rate coefficient for effector cells day − × − [10] σ Radiotherapy linear kill coefficient for effector cells Gy − × − Assumed φ Radiotherapy quadratic kill coefficient for effector cells Gy − × − Assumed τ Decay rate of immunotherapy (IL-2) day − × [7] ω Decay rate of chemotherapy (Bleomycin) day − × − [10] χ Decay rate of radiotherapy day − × − [9] I max Maximum tolerable dose of immunotherapy (IL-2) IU 7.20 × [3] C max Maximum tolerable dose of chemotherapy (Bleomycin) IU 3.00 × [1] R max Maximum tolerable dose of radiotherapy Gy 4.50 × [4]
1, 3
Assumed based on ratio from [4]. Assumed as equal to tumor cell kill coefficient as per [9].
C.i. Non-dimensionalization of Parameters
We define the nondimensionalized state variables asfollows: ˆ N = N/N , ˆ E = E/E , ˆ I = I/I , ˆ C = C/C , ˆ R = R/R , where the scaling factors are chosen such that N = 10 , E = 10 ,I = I max , C = C max ,R = R max , to rescale the populations by a reasonable order ofmagnitude. Time is rescaled with respect to tumor cell deactivationsuch that t (cid:48) = t/t , where t = ( γN ) − , the time scaling factor. This eliminates the parameter γ .The parameter values are non-dimensionalized with thefollowing transformations: α ∗ = α · t , β ∗ = β · N , γ ∗ = γ · N · t , δ ∗ = δ · t ,(cid:15) ∗ = (cid:15) · t , κ ∗ = κ · t , ζ ∗ = ζ · t E , λ ∗ = λ · t ,η ∗ = ηE , θ ∗ = θ · t , ν ∗ = ν · t , ρ ∗ = ρI ,µ ∗ = µ · t , σ ∗ = σ · R , φ ∗ = φ · R ,τ ∗ = τ · t , ω ∗ = ω · t , χ ∗ = χ · t . Dropping the stars for notational clarity, the finalnon-dimensionalized system is now given by:3
PREPRINT - J
ANUARY
29, 2021 ˙ˆ N ( t (cid:48) ) = α ˆ N (1 − β ˆ N ) − ˆ N ˆ E − δ ˆ N (1 − e − ˆ C ) − ˆ N ( (cid:15) ˆ R + κ ˆ R ) (2a) ˙ˆ E ( t (cid:48) ) = ζ − λ ˆ E + η ˆ E ˆ Nθ + ˆ N + ν ˆ E ˆ Iρ + ˆ I − µ ˆ E (1 − e − ˆ C ) − ˆ E ( σ ˆ R + φ ˆ R ) (2b) ˙ˆ I ( t (cid:48) ) = u ˆ I ( t (cid:48) ) − τ ˆ I (2c) ˙ˆ C ( t (cid:48) ) = u ˆ C ( t (cid:48) ) − ω ˆ C (2d) ˙ˆ R ( t (cid:48) ) = u ˆ R ( t (cid:48) ) − χ ˆ R (2e) D. Treatment-Free Equilibria and Stability
In this section we consider the equilibrium behavior ofthe system in the absence of treatment. For convenienceof notation, let x ≡ ˆ N and x ≡ ˆ E . First we must findthe nullclines of this system, or the curves along which d x d t = d x d t = 0 . The intersection of these nullclines givesthe fixed points of the system.Setting d x d t = 0 gives two nullclines of the system, oneof which represents a tumor-free equilibrium and the othera non-zero tumor population at equilibrium: x = 0 (3a) x = α (1 − βx ) ≡ h ( x ) (3b)Setting d x d t = 0 gives the last nullcline: x = ζλ − ηx θ + x ≡ j ( x ) (4)Using these nullclines, we can now analyze the variousequilibrium states of the system. D.i. Tumor-Free Equilibrium
The intersection of h ( x ) with x = 0 gives a fixed pointwith coordinates ( x , x ) = (0 , ζλ ) ≡ ξ .To permit stability analysis, we linearize about point ξ .Let f ( x , x ) represent d x d t , and f ( x , x ) represent d x d t for ease of notation. The Jacobian matrix for thissystem is of the form J ( x , x ) = ∂f ∂x ∂f ∂x ∂f ∂x ∂f ∂x Evaluating at ξ gives a Jacobian of J ( ξ ) = α − ζλ ηζθλ − λ By the Hartman-Grobman theorem, the real componentsof both eigenvalues of the Jacobian matrix must benegative for the given equilibrium state to be stable.Eigenvalues are λ = α − ζλλ = − λ < Thus, for the tumor-free fixed point ξ to be a stableequilibrium point for the system, it holds that α < ζλ This allows us to see that to eliminate the tumor cellpopulation at the stable equilibrium of the system, themaximum rate of tumor proliferation must be less thanthe equilibrium tumor cell population. If this condition isnot met, the tumor-free equilibrium will be unstable.
D.ii. Nonzero Equilibria
There may be zero, one, or two additional equilibriabesides ξ , depending on the relationship between h ( x ) and j ( x ) .Intersecting h ( x ) and j ( x ) gives a quadratic of thefollowing form, whose solutions are values of x for allother equilibria of the system. Ax + Bx + C = 0 , (5)where A = β ( λ − η ) , B = ζα + βλη + η − λ, and C = θ (cid:16) ζα − λ (cid:17) . The number of solutions to this quadratic depends onvalues of A , B , and C , thus depending on the parametervalues being used. For the chosen parameter values,the quadratic formula yields two positive values for x at equilibrium. However, one of these values yields a4 PREPRINT - J
ANUARY
29, 2021negative x and is thus eliminated. Therefore, in thiscase our analysis will only include one additional nonzerofixed point ξ . E. Limit Cycles
The Dulac-Bendixson criterion can be used to provethe absence of closed orbits in this system. TheDulac-Bendixson criterion states that if there are noclosed orbits, there exists a function φ such that ∂∂x ( φ ˙ x ) + ∂∂x ( φ ˙ x ) (cid:54) = 0 (6)We choose φ = x x . ∂∂x ( φ ˙ x ) + ∂∂x ( φ ˙ x ) (7) = − ( αβx + ζx x ) (cid:54) = 0 . (8) Therefore, the Dulac-Bendixson criterion is satisfiedand it holds that there are no limit cycles or homoclinicconnections in the system. Evaluating the criterion forboth treatment-free and treatment-present cases gives thisconclusion. This proves that the system will alwayscollapse to an equilibrium point, with or without anytherapeutic intervention. F. Numerical Solution
We numerically solve the ODE system using the PythonPyomo library. As mentioned before, the estimatedparameter values are such that there are two steady statesof the system: ξ being an unstable tumor-free state and ξ being a stable steady state with a non-zero, dormanttumor population. Graphs of dynamics are shown belowwith various initial conditions.Figure 1: Time dynamics (left) and phase portrait (right) of tumor cells and immune cells without treatment. Initialconditions ( N, E ) = (6 × , × ) , (5 × , . × ) , (1 . × , × ) , (6 × , × ) III. O
PTIMAL C ONTROL
A general optimal control problem is stated as thefollowing: a system of ODEs is given that models thedynamics of a system in terms of state variables andcontrol inputs, such as ˙ x = f ( x , u , t ) , (9)where f is a vector-valued function and x , u are vectorsthat describe the state of the system and values of controlinputs, respectively. The goal is to drive the systemto a final state given initial states of the system, whileminimizing an objective function J through the timeduration [ t , t f ] . There may also be additional constraintsadded that must be satisfied for this duration, called path constraints. In the following sections, we aim todefine and solve an optimal control problem to optimizetreatment. A. Defining Objective Function
In this section, we explore a variant of the problem inwhich the final time t f is fixed at a particular value, knownas a fixed-final time problem in optimal control theory.The objective function is the function that must beminimized through the course of treatment. In thisproblem, there are several possible metrics that can beused to represent the goals for treatment duration andoutcome.5 PREPRINT - J
ANUARY
29, 2021
A.i. Incorporating Physician Preference
There must be some margin provided in the model toallow for physician preference to be incorporated intooptimization. We introduce weighting coefficients w and w to allow weighted preference for two different parts ofthe objective function. Thus, the objective function is inthe form J = w K + w L , (10)where K and L represent chosen expressions involvingdifferent quantities in the model that are to be minimized.The rest of this section will explore various choicesfor these expressions, and the impact they have onoptimization results. Since there are two terms in theobjective function, the results will be Pareto-optimal. A.ii. Maximizing Average Immune Cell Count
By selecting weighting coefficients such that w < , onecan choose to maximize a quantity throughout treatmentduration. Ideally, patient health should be maximizedthroughout treatment, and a good representative markerfor patient health is immune cell count. Thus, we select L to be average immune cell count, in the form L = 1 t f (cid:90) t f ˆ E ( t ) dt. (11)This is the L used for all optimization scenarios exploredin the rest of the project. Note that with varying values of w < , a physician can select an optimization setup thatreflects the preferred approach for treatment, choosing amore cautious route with higher values of | w | , and a moreaggressive treatment plan with lower values of | w | . A.iii. Minimizing Final Tumor Cell Count
One possible aim is to drive the tumor population to aslow a value as possible at the final time of treatment. Inthis case, the value of K would be defined as the tumorcell count at the end of treatment. K = ˆ N ( t f ) . (12) A.iv. Minimizing Average Tumor Cell Count
Another possible aim is to keep the tumor population aslow as possible for not just the end, but the entire durationof treatment, to prevent development of malignancy. Assuch, K would be defined as the average tumor populationover the treatment time. K = 1 t f (cid:90) t f ˆ N ( t ) dt. (13) A.v. Minimizing Treatment
There is a limited set of initial conditions and parametervalues where a solution is possible such that N ( t f ) = 0 .In these cases, it is possible to explore the solution spaceand optimize for the solution wherein the least treatmentis used. Minimizing the provided treatment serves to alsominimize toxicity to the patient and cost of the therapy.In this case, K would be defined as the integral of alltreatment options used over the duration of therapy: K = (cid:90) t f u I ( t ) + u C ( t ) + u R ( t ) dt. (14)When using this objective function, there must also be aconstraint placed on the optimization setup, such that ˆ N ( t f ) = 0 . (15)This restricts the use of this objective function to onlythe cases where it is possible to eliminate the tumor– otherwise, the constraint is infeasible and one of theprevious two objective functions must be used instead. A.vi. Pontryagin’s Maximum Principle
Recall that u I , u C , and u R represent the infusion rateof immuno-, chemo-, and radiotherapy, respectively, at agiven time.For future notation, we define u to be the vector [ u I u c u R ] T , and x to be [ x x x x x ] T ,where x , x , x represent ˆ I, ˆ C, ˆ R respectively, and asshown earlier, x , x represent ˆ N , ˆ E respectively. Thisallows us to capture the dynamics in Equations 2a-2e as ˙ x = f ( x , u , t ) .Now that we have the correct form for the dynamics, wecan use Pontryagin’s Maximum Principle. Pontryagin’sMaximum Principle defines the first-order conditionsnecessary for a particular solution to the optimal controlproblem to be optimal.The Hamiltonian of the system is defined as H = J ( u , t ) + λ T ( t ) f ( x , u , t ) . (16) λ ( t ) is the vector of co-state (adjoint) variables, whichvary with time. The Maximum Principle states that to findan optimal control, the Hamiltonian in this case must beminimized. The necessary conditions for optimality areas follows:1. ˙ x = ∂H∂λ = f ( x , u , t ) (State Equation)2. ˙ λ = − ∂H∂ x (Costate Equation)3. ∂H∂ u = 0 (Optimal Control)4. ≤ u ≤ (Constraint on Control Inputs)5. x ( t ) = x (Initial Conditions)6 PREPRINT - J
ANUARY
29, 2021
A.vii. Numerical Solutions
The Python Pyomo library [6] provides a robustframework for numerical optimization. The co-state variables are calculated to minimize the Hamiltonian andmeet the aforementioned constraints. Below are shownthe results after optimizing with the chosen constraintsand objective function. (a)
Minimize final tumor population. (b)
Minimize average tumor population. (c)
Minimize treatment administered, paired with constraint that final tumor population is zero.
Figure 2:
Time dynamics (left) and phase portrait (middle) of tumor cells and immune cells with treatment. Treatmentregimen (right) showing the calculated control inputs for the duration of treatment. Initial conditions ( N, E ) =(6 × , × ) . Weighting coefficients w = 1 , w = 0 . . PREPRINT - J
ANUARY
29, 2021
B. Results and Discussion
As shown in the computational experiments above, theoptimal treatment regimen differs depending on thechosen objective function. When minimizing final tumorpopulation, the optimal treatment plan involves peaks ofbleomycin and radiotherapy near the last 10-day period,paired with high-dose IL-2 after the first month oftreatment. When minimizing average tumor population,the regimen involves aggressive early peaks of bleomycinand radiotherapy for the first 2-week period and high-doseIL-2 for the entire duration. Finally, since eliminating thetumor is feasible for these parameters, the most optimaltreatment regimen – to eliminate the tumor with minimumimpact to patient health – is one where a single 10-daytreatment cycle of radiotherapy is used. The dynamicswill vary for different parameter values, due to thepresence of additional equilibria states. Regardless, themodel with the given parameters is an effective case studyof generating an optimal treatment plan from pre-definedpatient data values.IV. C
ONCLUSION
This project builds upon past contributions tocomputational oncology by extending the use of optimalcontrol theory to combination therapy, as opposed tosingle-therapy treatment plans. This is a stepping stone tothe possibility of personalized, optimized combinationtherapy regimens built around patient data. Futurework will include optimization of treatment time andincorporating more biological phenomena into the model.R
EFERENCES [1] A
STON , W. J., H
OPE , D. E., N
OWAK , A. K.,R
OBINSON , B. W., L
AKE , R. A.,
AND L ESTERHUIS , W. J. A systematic investigationof the maximum tolerated dose of cytotoxicchemotherapy with and without supportive care inmice.
BMC Cancer 17 , 1 (2017).[2] C
HANGRAN , G., P
AGANETTI , H.,
AND G RASSBERGER , C. Prediction of treatmentresponse for combined chemo- and radiationtherapy for non-small cell lung cancer patients usinga bio-mathematical model.
Scientific Reports 7 (122017).[3] D
UTCHER , J. P., S
CHWARTZENTRUBER , D. J.,K
AUFMAN , H. L., A
GARWALA , S. S., T
ARHINI , A. A., L
OWDER , J. N.,
AND A TKINS , M. B. Highdose interleukin-2 (aldesleukin) - expert consensuson best management practices-2014.
Journal forImmunoTherapy of Cancer 2 , 1 (2014).[4] F
OWLER , J. F., T OM ´ E , W. A., F ENWICK , J. D.,
AND M EHTA , M. P. A challenge to traditionalradiation oncology.
International Journal ofRadiation Oncology*Biology*Physics 60 , 4 (2004),1241–1256.[5] G
LUZMAN , M., S
COTT , J. G.,
AND V LADIMIRSKY , A. Optimizing adaptive cancertherapy: dynamic programming and evolutionarygame theory.
Proceedings of the Royal Society B:Biological Sciences 287 , 1925 (2020), 20192454.[6] H
ART , W. E., W
ATSON , J.-P.,
AND W OODRUFF ,D. L. Pyomo: modeling and solving mathematicalprograms in python.
Mathematical ProgrammingComputation 3 , 3 (2011), 219–260.[7] K
IRSCHNER , D.,
AND P ANETTA , J. C. Modelingimmunotherapy of the tumor - immune interaction.
Journal of Mathematical Biology 37 , 3 (1998),235–252.[8] K
UZNETSOV , V., M
AKALKIN , I., T
AYLOR , M.,
AND P ERELSON , A. Nonlinear dynamics ofimmunogenic tumors: Parameter estimation andglobal bifurcation analysis.
Bulletin of mathematicalbiology 56 (04 1994), 295–321.[9] M
KANGO , S. B., S
HABAN , N., M
UREITHI ,E.,
AND N GOMA , T. Dynamics of breastcancer under different rates of chemoradiotherapy.
Computational and Mathematical Methods inMedicine 2019 (2019), 1–12.[10] P
ILLIS , L. D., G U , W., AND R ADUNSKAYA ,A. Mixed immunotherapy and chemotherapyof tumors: modeling, applications and biologicalinterpretations.
Journal of Theoretical Biology 238 ,4 (2006), 841–862.[11] P
ILLIS , L. D.,
AND R ADUNSKAYA , A. Thedynamics of an optimally controlled tumor model: Acase study.
Mathematical and Computer Modelling37 , 11 (2003), 1221–1244.[12] V
ILLASANA , M.,
AND R ADUNSKAYA , A. Adelay differential equation model for tumor growth.