Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance
(USC Thesis Other)
Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DETERMINISTIC EVOLUTIONARY GAME THEORY MODELS FOR THE NONLINEAR DYNAMICS AND CONTROL OF CHEMOTHERAPEUTIC RESISTANCE by Yongqian Ma A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) August 2020 Copyright 2020 Yongqian Ma Dedication To my dear parents. ii Acknowledgments Firstly thanks a lot for my advisor Paul Newton and my lab mates. You provided me a lot of ideas about research and my future career. Without your help I won’t be able to complete all the work smoothly. Secondly I want to thank all my friends, my roommates, my classmates, etc. They kept me accompany throughout my PhD life in the past six years in the US or from China. It wouldn’t be as colorful as I had without you. After that I received a lot of encourage and care from my relatives, my uncles, aunts and cousins. During this special time of global pandemic, it’s so warm to hear from the other side of the earth. Then as one of my numbered interest I had a lot fun watching anime and playing games during my spare time. Thanks for all the wonderful products and I hope this year won’t be the terminal, but a turning point for the industry. Last but not least, there are also many online friends and groups of netizens. We may or may not know each other, but during my struggling time in the midnight, I really got relieved by joking, talking or even trolling with each other. Best regards for everyone in the world to survive the Covid-19 pandemic and to live a happy and peaceful life. iii Contents Dedication ii Acknowledgments iii List of Tables vi List of Figures vii Abstract xx 1 Introduction 1 2 CapitalizingonCompetition: AnEvolutionaryModelofCompetitiveRelease in Metastatic Castration Resistant Prostate Cancer Treatment 6 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 The likelihood of pre-existing resistance . . . . . . . . . . . . . . . . . 8 2.1.2 Leveraging the cost of pre-existing resistance . . . . . . . . . . . . . . 9 2.1.3 Retrospective analysis of metastatic castration-resistant prostate cancer 11 2.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 The replicator equation model . . . . . . . . . . . . . . . . . . . . . . 16 2.2.2 The linearized system and the cost of resistance . . . . . . . . . . . . 18 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3.1 Managing competitive release . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Nonlinear adaptive control of competitive release and chemotherapeutic resistance 30 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 A three-component replicator system . . . . . . . . . . . . . . . . . . . . . . 34 3.3 The prisoner’s dilemma as a cancer model . . . . . . . . . . . . . . . . . . . 37 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.1 Continuous therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.2 Time-dependent therapy . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4.3 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 iv 4 The role of synergy and antagonism in designing multi-drug adaptive chemotherapy schedules 56 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2 A three-component replicator system . . . . . . . . . . . . . . . . . . . . . . 60 4.2.1 Why the prisoner’s dilemma game? . . . . . . . . . . . . . . . . . . . 60 4.2.2 The three-component model . . . . . . . . . . . . . . . . . . . . . . . 63 4.3 Constant chemotherapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3.1 Additive interactions e = 0 . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3.2 Non-additive interactions . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 Designing adaptive schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5 Optimal control of evolutionary games: The two-by-two case 79 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2 model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.2.1 evolutionary games and replicator dynamics . . . . . . . . . . . . . . 81 5.2.2 Pontryagin Maximum Principle . . . . . . . . . . . . . . . . . . . . . 84 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4 discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.1 f 2 decrease . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.2 f 1 increases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.4.3 Time re-scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6 Optimal control of levitation by blowing 94 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.3 Time-dependent blowing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.4 Optimal blowing schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.4.1 Minimum time strategy . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.4.2 Minimum air volume strategy . . . . . . . . . . . . . . . . . . . . . . 105 6.4.3 Minimum energy strategy . . . . . . . . . . . . . . . . . . . . . . . . 106 6.5 Pontryagin Maximum Principle . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.5.1 Singular control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Reference List 119 A replicator dynamics of tumor growth 129 v List of Tables 4.1 Total dose (D), Time period (T), and Average Dose (D/T) associated with adaptive therapies with antagonistic, additive, and synergistic drug interactions. 74 6.1 Time, Energy and Air volume values for two trajectories shown in figure 6.3. The results are simulated using variables: b 1 = 6,α = 1m −1 ,g = 9.8m/s 2 ,m = 1kg. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 vi List of Figures 2.1 Schematic of competitive release in a tumor — (a) Prior to treatment, a tumor consists of a large population of sensitive cells (red) and a small population of less fit resistant cells (green) competing for resources with the surrounding healthy cells (blue); (b) Chemotherapy targets the sensitive pop- ulation (middle), selecting for the less fit resistant population that thrives in the absence of competition from the sensitive population; (c) Upon regrowth, the tumor composition has larger numbers of resistant cells, rendering the subsequent rounds of treatment less effective. . . . . . . . . . . . . . . . . . . 12 2.2 Clonal evolution of competitive release — A fishplot (sometimes known as a Müller plot), showing the tumor quantity (vertical axis) and composition (sensitive: red; resistant: green) over time (horizontal axis, left to right) with important events annotated. After first driver mutation (left), initial exponentialgrowthofsensitivepopulationoccursuntildiagnosis(dashedline). Continuous therapy targeting the chemo-sensitive population responds well withadecreaseintumorburden. Intheabsenceofsensitivecells, theresistant population (existing in small numbers before the start of therapy) grows to become the dominant subpopulation at relapse, albeit typically with lower exponential growth rate due to the cost of resistance. . . . . . . . . . . . . . 14 vii 2.3 Model fits of PSA measurement data for metastatic castration- resistant prostate cancer — Prostate-specific antigen (PSA) measurement data of four representative patients from three randomised clinical trials with metastatic castration-resistant prostate cancer. (A) A single patient data shown with model fit (eqn. 2.4) of growing resistant population (green), decayingsensitivepopulation(red), wherethesumrepresentsanoverallinitial good response and eventual relapse (q(t), dashed line). Left column: treat- ment with prednisone only [75]; middle column: treatment with mitoxantrone and prednisone [112]; right column: treatment with docetaxel and prednisone [89]. PSA data is normalized by q(t = 0) (black dots). The data is fit using the exponential model (eqn. 2.1; blue curve) and the replicator model (eqns 2.2, 2.3). Each patient is fit reasonably well with both models. Data was fit by parameter sweep of cost (α−β, eqn. 2.15, 2.16), initial fractional resistance f and selection pressure w. Parameters used are α−β = [0.02, 0.07, 0.20, 0.00, 0.20, 0.00, 0.20, 0.03, 0.20, 0.01, 0.15, 0.19]; f = [0.20, 0.02, 0.09, 0.05, 0.03, 0.07, 0.10, 0.10, 0.09, 0.10, 0.03, 0.03]; w = [0.10, 0.20, 0.20, 0.30, 0.20, 0.10, 0.30, 0.15, 0.10, 0.15, 0.20, 0.20] for B - M respectively. . . . . . . . . . 26 2.4 Fitness landscape before and during therapy — (A) Before therapy, a driver mutation leads to a fitness advantage of the cancer cell (red) and a subsequent resistant-conferring mutation comes at a fitness cost (green). Parameters determined by the prisoner’s dilemma payoff matrix reflect these relative fitness differences. (B) Rather than change the elements of the game payoff matrix directly, the selection parameter is manipulated for healthy (increased)andcancer(decreased)subpopulationssuchthatw 1 −w 2 = 2c. (C) Changes in w i results in a new relative fitness of each subpopulation during therapy. The fitness of the resistant population is unaffected by therapy’s selective pressure, but the healthy population is given an advantage over the chemo-sensitive population. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 viii 2.5 Dynamic phase portraits before and during chemotherapy — (A) Trilinear coordinate phase space representation; (B) Schematic of proposed adaptive therapy concept using the resistant nullclines to determine therapy “on" and “off" times in order to trap the tumor in the controllable region 2, and reach approximate cycle that repeats back on itself in red. The contin- uous therapy is also plotted in dashed blue, for comparison. Two nullclines divide the triangle into 3 regions; region 1: ˙ x R < 0 for both therapy on and off; region 2: ˙ x R > 0 for therapy off and ˙ x R < 0 for therapy on; region 3: ˙ x R > 0forboththerapyonandoff. (C)Beforechemotherapy, thehealthy(H), sensitive (S), and resistant (R) populations compete on a dynamical fitness landscape, with several solution trajectories shown (black) and the instanta- neous relative velocity indicated by background color gradient (red to blue). All internal trajectories lead to tumor growth and eventual saturation of the sensitive population (bottom left corner). Each population nullcline (line of zero growth: ˙ x i = 0) is plotted: healthy (dashed blue), sensitive (dashed red), and resistant (dashed green). The nullclines divide the triangle into 3 regions. Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; Region 2: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0; Region 3: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0; (D) Chemotherapy alters the selection pressure to the disadvantage of chemo-sensitive cancer population and advantage of the healthy population (shown for c = 0.6,α = 0.020, β = 0.018, w = 0.1). In this case, the nullclines divide the triangle into 6 regions; Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; Region 2: ˙ x H > 0 ˙ x S < 0 ˙ x R < 0; Region 3: ˙ x H > 0 ˙ x S < 0 ˙ x R > 0;Region 4: ˙ x H < 0 ˙ x S < 0 ˙ x R > 0; Region 5: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0; Region 6: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0; Solution trajectories (black) show initial trajectory toward healthy saturation (triangle top) but eventual relapse toward resistant population (bottom right of triangle) upon passing the resistant nullcline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 ix 2.6 The effect of dose on tumor relapse and progression free survival under continuous and adaptive therapy — (A) PSA data is from a sin- gle patient (see figure 2.3B) under continuous treatment (Mitoxantrone and Prednisone) is replotted (black dots) along with the best model fit (eqns. 2.2, 2.3) in blue. Using identical patient parameters, continuous treatment of two “new" drugs with higher effectiveness (i.e. increased effective dose,c) is shown in red and yellow. Time to relapse significantly increases with increasing dose while the progression free survival shows marginal, but decreasing, difference. An adaptive therapy (see figure 2.5B) is also simulated (solid black line), showing an increased control over the tumor; Similar control paradigms were simulated with error in administration times: slightly shorter time on therapy (1 day; green) and slightly longer (1 day; purple). (B) The same four therapies are shown in a fish plot. Continous therapies show relapse to initial tumor size is dominated by chemo-resistant population (green). The adaptive ther- apy successfully suppressed the growth of the resistant population (bottom). Parameters used: α−β = 0.02; w = 0.10; f = 0.2 . . . . . . . . . . . . . . 29 3.1 Schematic dose schedules.(a) Maximum tolerated dose (MTD) schedule; (b) Low dose metronomic (LDM) schedule; (c) Adaptive schedule; (d) General time-dependent schedule accuated by piecewise constant dose concentrations. 33 3.2 The competitive release mechanism. (a) The three-component phase space associated with competing populations of (H,S,R) cells. (b) With no therapy, the S corner of the triangle is a globally attracting fixed point, while the H and R corners are unstable. All initial conditions lead to a saturated tumor. Filled circles are stable, unfilled circles are unstable. (c) For continuous chemotherapy above a threshold level, theR corner of the triangle isagloballyattractingfixedpoint, whiletheH andS cornersareunstable. All initial conditions (except those lying on the separatrix connecting the interior balanced fitness state to the S−R side) lead to a resistant tumor. . . . . . . 37 x 3.3 Fixed point locations. Tracking the location of the fixed points as a func- tion of chemo-concentration parameterC. The fractional values indicate that the values are analytically obtained. . . . . . . . . . . . . . . . . . . . . . . . 41 3.4 Basins of attraction and nullclines. Panel showing the separatrices and nullclines (dx H /dt = 0,dx S /dt = 0,dx R /dt = 0) through the balanced fitness interior fixed point that determines the basin boundaries of attraction of the S state and the R states. The interior fixed point is the one in which each of the subpopulation fitness levels exactly matches the average fitness of the entire population. Filled circles are stable, unfilled circles are unstable. . . . 42 3.5 Basin of attraction areas. Areas of basin of attraction of S fixed point and R fixed point as a function of chemo-concentration C. . . . . . . . . . . . . . 43 3.6 Fitness landscapes as a function of time. With no therapy (C = 0), the fitness curves are contunuously decreasing functions as the tumor satu- rates with the sensitive cell population in a sigmoidal shaped growth curve. With continuous therapy (C = 0.6), the fitnesses initially increase indicating tumorregression, buteventuallydecrease. Thehealthysubpopulationinitially increases before the resistant population eventually saturates the tumor. Ini- tial condition (0.9,0.09,0.01) with w 0 = 0.1 . . . . . . . . . . . . . . . . . . . 44 3.7 Dynamical trajectories for constant C.(a) With no chemotherapy (C = 0), the sensitive corner S is a globally attracting fixed point. All initial con- ditions inside the triangle move toS along the sample trajectories shown. (b) Above the chemotherapy threshold C > 0.5, all initial conditions inside the triangle move to the resistant corner R. Shown are sample trajectories for C = 0.7. (c) Overlay of the solution trajectories forC = 0 andC = 0.7 create a curvilinear grid throughout the triangle. By switching between these two values of C, we construct a global trajectory made up of segments of the two families of trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 xi 3.8 Constructing a closed loop trajectory using C = 0 and C = 0.5, with C avg = 0.366.(a) By using segments of a trajectory for C = 0 and C = 0.5, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories startingfrompointA. BotheventuallymovetotheS corner,althoughinitially the MTD trajectory moves toward the H corner (tumor regression) before recurrance. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitnesshfi for the three different chemo-schedules. The adaptive schedule is able to maintain a higher average fitness throughout. . . 48 3.9 Constructing a closed loop trajectory using C = 0 and C = 0.6, with C avg = 0.375.(a) By using segments of a trajectory for C = 0 and C = 0.6, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories starting from point A. The MTD trajectory eventually moves to the S corner, while the LDM trajectory moves to the R corner. Both trajectories initially move toward the H corner (tumor regression) before recurrance. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitness hfi for the three different chemo-schedules. The adaptive schedule is able to maintain a higher average fitness throughout, although LDM initially achieves higher average fitness before declining. . . . . . . . . . . . . . . . . . . . . . 49 xii 3.10 Constructing a closed loop trajectory using C = 0.3 and C = 0.6, with C avg = 0.396.(a) By using segments of a trajectory for C = 0.3 and C = 0.6, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories starting from point A. Both trajectories move toward the R corner after initially moving toward H. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitnesshfi for the three different chemo- schedules. The adaptive schedule is able to maintain a higher average fitness throughout. Both MTD and LDM initially show higher average fitness before declining. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Constructing orbit transfers. (a) We take two arbitrary points labeled A and B and construct the incoming and outgoing solution trajectories from each, using two values C = 0 and C = 0.7. There must be a crossing point, which we label point O. The two segments AO followed by OB is the two- switch trajectory that takes us fromA toB. (b) Shown is the curvilinear grid constructed by families of solution curves for the three valuesC = 0,C = 0.3, C = 0.7. Any point of intersection can be used as a starting point and an ending point to construct a three-switch orbit. . . . . . . . . . . . . . . . . . 51 3.12 Robust family of nested orbits. (a) One example of a (continuous) family of nested closed orbits using bang-bang strategy toggling between C = 0 and C = 0.7, all with similar total dose (C∼ 0.4). A closed loop can be made arbitrarily small around a center point described in Figure 3.13. (b) Schedule associated with orbits 1 and 4; (c) Schedule associated with orbits 1 and 5; (d) Schedule associated with orbits 1 and 6; (e) Schedule associated with orbits 2 and 4; (f)Schedule associated with orbits 2 and 5; (g) Schedule associated with orbits 2 and 6; (h)Schedule associated with orbits 3 and 4; (i)Schedule associated with orbits 3 and 5; (j) Schedule associated with orbits 3 and 6. . 54 xiii 3.13 All closed orbit families are centered on dashed line defined by f 1 = f 3 . (a) As a family of closed orbits enclose smaller and smaller loops, they intersect at a point of tangency. The point of tangency must lie on the dashed line. (b) Examples of several closed loop orbits centered at different points down the line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1 (a) e> 0 showing synergistic profiles; (b) e = 0 showing additive profiles; (c) e< 0 showing antagonistic profiles. . . . . . . . . . . . . . . . . . . . . . . . 59 4.2 Depiction of evolutionary stable states (ESS) for different chemotherapy val- ues. (a) With no chemotherapy, the tumor saturates to the S corner regard- less of the initial make-up of the three subpopulations. (b) With C 1 = 0, and C 2 = 0.8, competitive release of the resistant population R 1 drives all trajectories to the R 1 corner. (c) With C 1 = 0.8 and C 2 = 0, competitive release of the resistant population R 2 drives all trajectories to the R 2 corner. (d) Trajectories associated with three different constant combinations of C 1 and C 2 , depicting the overlap of the trajectories at different times. . . . . . . 67 4.3 Tumor growth curves (log-plots) for untreated and constant therapies. Tumor recurrence (dashed line) occurs at t∼ 75 dimensionless time units. . . . . . . 68 4.4 Panel showing the full range of trajectories, ESS and their basins of attraction forconstantadditivetherapiesintherange 0≤C 1 ≤ 1, 0≤C 2 ≤ 1,C 1 +C 2 ≤ 1. Bifurcation values along top row are C 2 = 1/3, C 2 = 1/2. Bifurcation values along left column are C 1 = 7/18, C 1 = 1/2. . . . . . . . . . . . . . . . 69 4.5 Theeffectofsynergy(e> 0)vs. antagonism(e< 0)ontheESSandthebasins of attraction for a fixed choice of constant combination therapies C 1 = 0.35, C 2 = 0.27. (a) e =−0.4; (b) e =−0.3; (c) e =−0.2; (d) e =−0.1; (e) e = 0; (f) e = 0.1; (g) e = 0.2; (h) e = 0.3; (i) e = 0.4. . . . . . . . . . . . . . . . . 70 xiv 4.6 The opening of the basin of attraction for R 1 via a transcritical bifurcation at e =−310/1701 =−0.1822, C 1 = 0.35, C 2 = 0.27. At the bifurcation point, the areas of the basins of attraction of the respective regions are S = 76.2%,R 2 = 13.8%. (a) e =−0.21; (b) e =−0.19; (c) e =−0.175; (d) e =−0.15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.7 Four ways of depicting the antagonistic transcritical bifurcation. (a) Eigen- value collision that takes place as the two fixed points collide at theR 1 corner. The other eigenvalue remains negative for both fixed points; (b) Transcritical bifurcation diagram; (c) Pre-bifurcation phase plane; (d) Bifurcation phase plane; (e) Post-bifurcation phase plane. . . . . . . . . . . . . . . . . . . . . . 73 4.8 The closing of the basin of attraction for S via a transcritical bifurcation at e = 30/189 = 0.1587, C 1 = 0.35, C 2 = 0.27. At the bifurcation point, the areasofthebasinsofattractionoftherespectiveregionsareR 1 = 72.3%,R 2 = 17.7%. (a) e = 0; (b) e = 0.1; (c) e = 0.2; (d) e = 0.3. . . . . . . . . . . . . . 74 4.9 Four ways of depicting the synergistic transcritical bifurcation. (a) Eigen- value collision that takes place as the two fixed points collide at theS corner. The other eigenvalue remains negative for both fixed points; (b) Transcritical bifurcation diagram; (c) Pre-bifurcation phase plane; (d) Bifurcation phase plane; (e) Post-bifurcation phase plane. . . . . . . . . . . . . . . . . . . . . . 75 4.10 Areas of the three basins of attraction as a function of the parameter e. Note the sensitivity of theS andR 1 regions in the antagonistic regime−0.2≤e≤ −0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.11 (a) Average fitness curves for e = 0.3, 0.15, 0,−0.15,−0.3; (b) Tumor growth curves (log-plot) for e = 0.3, 0.15, 0,−0.15,−0.3. . . . . . . . . . . . . . . . . 76 xv 4.12 Closed loop adaptive schedules OABO using three cycles for e = 0, 0.3,−0.3. (a)e = 0, OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (b)Correspondingadaptiveschedule; (c)e = 0.3,OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (d) Corresponding adaptive schedule; (e) e =−0.3, OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (f) Corresponding adaptive schedule. . . . . . . . . . . . 77 4.13 Tumor growth curves for the adaptive schedules from figure 4.12 as compared with untreated growth (exponential) and constant schedule which eventually leads to tumor recurrence. (a) e = 0; (b) e = 0.3; (c) e =−0.3. . . . . . . . . 77 4.14 Rates of adaptation for e = 0, 0.3,−0.3: (a) R 1 rate of adaptation; (b) R 2 rate of adaptation; (c) x 1 rate of adaptation; (d) x 2 rate of adaptation; . . . 78 5.1 Schema for the Maximum tolerated Dose (MTD) and Low Dose Metronomic (LMD) schedule. The integral of these two schedule are the same. . . . . . . 88 5.2 Simulation results with D = 0.3,T = 1,x(0) = 0.1 and payoff matrix in eq.5.43 (a) the corresponding trajectory for MTD and LDM simulated by standard replicator equation.(b)the corresponding trajectory for MTD and LDM simulated by Adjusted replicator equation. . . . . . . . . . . . . . . . . 89 5.3 Bang-bang optimal control for both standard replicator equation and adjusted replicator equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.4 Simulation results with D = 0.3,T = 1,x(0) = 0.1 and payoff matrix in eq.5.44 with Adjusted replicator equation. . . . . . . . . . . . . . . . . . . . 91 5.5 (a)optimal controlu ∗ (t) for adjusted replicator equation (b) the corresponding trajectory x ∗ (t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.1 Ping-pong ball levitated by blowing through a straw. The streaming air blows upwards along the positive x direction against the gravity. x ∗ marks the equilibrium position where the forces balance. For full video, see https: //www.youtube.com/watch?v=bCRjPFhlSYk . . . . . . . . . . . . . . . . . . 96 xvi 6.2 (a) Phase curves in the (x,v) plane for 5 different values of the energy E b . x max marks the maximum height that the ball can achieve on the energy curve E b =E 0 b , whiletheellipticpointx ∗ marksthelevitationpointwheretheforces balance E b = E ∗ b . (b) Energy curves E b = E 0 b going through the origin for four different values of b. (c) Phase curves of free fall (b = 0) with different values of energy. (d)x max marks the intersection of the two curves defined by the left and right hand sides of eqn (6.8). . . . . . . . . . . . . . . . . . . . . 113 6.3 Diagram showing two different paths from point O to point x ∗ . Orbits cor- responding to energy E b = E 0 b for two different values, b = b 1 (blue) and b =b 2 <b 1 (red). Heightx ∗ represents the elliptic fixed point associated with b = b 1 . The value of b 2 is chosen so that x max for the red orbit intersects x ∗ . The pink orbit corresponds to a free-fall trajectory b = 0 with an initial upward launch velocity v ∗ so that the trajectory intersects x ∗ . Point A, with coordinates (x A ,v A ), marks the intersection of the free-fall trajectory with the blue orbit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 xvii 6.4 Examples of piecewise constant blowing schedules that move the ball from the origin to the target point (x,v) = (x ∗ , 0). Parameter values are given by: b∈ [0, 6],α = 1m −1 ,g = 9.8m/s 2 ,x ∗ =log(2)meter,m = 1kg, and we choose T = 1.34s for time scale. (a) Two-step schedule: b(t) starts with b = αx ∗ e αx ∗ e αx ∗ −1 , then it switches to b = e αx ∗ . (b) The corresponding trajectories for the two- step schedule. Note that the ball will remain at equilibrium on second step, so the trajectory only has one segment. (c) Three-step schedule with blowing on the first segment, followed by free falling (b = 0) for the second segment. (d) The corresponding trajectory for (c). The trajectories will get closer and closer to vertical axis with the increase of initialb(t). (e) Examples of general three-step schedule, with b 1 for first segment and b 2 for the second. (f) The corresponding b(t) for the examples of general three step schedule. (g) A schematic of a four-step schedule, where t 1 t 2 and t 3 mark three switching points.(h) The corresponding trajectory for (g) . . . . . . . . . . . . . . . . 115 6.5 The schematic of the trajectory driven by a general blowing function b(t) as indicated in blue. For comparison, we show what we call the bang-bang strategy, starting with maximum blowing b max from t = 0 to t =t s , then free fallb = 0 from t =t s tot =t f . Any general trajectory associated with b(t) is trapped between the two segments of the other orbit. . . . . . . . . . . . . . 116 6.6 The bang-bang (on-off) strategy is optimal to attain the transfer in minimum timeandusingminimumairvolume. (a)Theactualtrajectoryforthetransfer; (b) The time-dependent blowing schedule. . . . . . . . . . . . . . . . . . . . 116 6.7 Minimum energy strategy.(a)phase space. Simulated with parameters b max = 20,α = 1m −1 ,x ∗ =ln(5)m,g = 9.8m/s 2 ,m = 1kg (b)corresponding schedule for (a) .(c)phase space simulated with parameters: b max = 5,α = 1m −1 ,x ∗ = ln(5)m,g = 9.8m/s 2 ,m = 1kg.(d) corresponding optimal control and it’s still bang-bang control. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 xviii 6.8 Constant levitation with dissipation.The numbers for the simulation are b = 5,α = 1m −1 ,g = 9.8m/s 2 ,x ∗ =log(5)m,m = 1kg,γ = 1kg/s . . . . . . . . . 118 xix Abstract In this thesis, we use a three by three replicator dynamical system with a payoff matrix of Prisoner’sDilemmatypetodevelopevolutionarygametheorymodelsofadaptivechemother- apy that avoid chemotherapeutic resistance in cancer. With the primary goal of controlling resistance and delaying tumor recurrence, we introduce and develop a trajectory design method from the orbital mechanics literature in order to construct families of closed loops, called ‘evolutionary cycles’, using adaptive time-dependent dose schedules. With time- dependent adaptive chemotherapy control, we attempt to alter the fitness landscape of the heterogeneous cell populations in order to delay tumor recurrence. We show that for a three-component population made up of healthy cells, sensitive cells, and resistant cells, this method outcompetes both maximum tolerated dose (MTD) schedules and low-dose metro- nomic (LDM) schedules, which are standard currently used protocols. We then consider a more complex system made up of two types of resistant cells and two drugs applied sequen- tially and in combination. In this case, the drugs can interact synergistically, additively, or antagonistically. We show how to develop multi-drug schedules in this context using a similar method of adaptive trajectory design. Generally speaking, we show that antagonistic drug interactions are more effective at delaying tumor recurrence than synergistic ones which effectively cause the two drugs to act as one, limiting the flexibility in schedule design. We then apply the Pontryagin Maximum Principle to two by two evolutionary games. We show examples of payoff matrices where bang-bang control is optimal, and others where singular control techniques need to be applied. As a final chapter, we use a Hamiltonian model along with the trajectory design method to show how to levitate a mass with specially designed time-dependent blowing schedules. xx Chapter 1 Introduction My Ph.D. thesis develops techniques of adaptive control applied to evolutionary game theory[108, 85, 52] models of tumor growth and chemotherapeutic resistance[58]. We use the Prisoner’s Dilemma (PD) game[66, 85] and the replicator dynamical system as a paradigm for the competition (cooperation or defection) among heterogeneous cell subpopulations in a tumor as the system undergoes evolution via natural selection. We use adaptive control[40] methods to shape the fitness landscape of the tumor and steer its evolutionary outcome. An overriding goal of our adaptive (time-dependent) chemotherapy schedules is to manage the development of chemo-resistance in a tumor. Although the models used in this thesis focus on the prevention of chemotherapeutic resistance in tumors, the techniques developed may also be useful in steering the evolution of microbial populations with the goal of mitigating and reversing antibiotic resistance, which has a large and specialized separate literature[58, 1, 17, 105, 39, 44, 27]. Because cancer is a rising threat to human health, it is important to understand the underlying biological and physical mechanisms of tumor growth and the complex pharma- cological interactions between cancer cells and the toxins used to kill them. From the view of evolution, the proliferation of cancer cells in the patient is widely recognized as a process of natural selection [28]. As cancer cells proliferate and die, they can mutate, changing their rate of reproduction and ability to survive chemotherapy. Macroscopically, a tumor consists of multiple subpopulations of cells with varying growth rates as they compete for resources. The most aggressive subgroups have enhanced ability to gain more resources, grow faster and 1 take over the majority of the population. The weaker subpopulations, on the other hand, are unable to reproduce as fast so their population size will remain low, reduced, or even be completely eliminated. The result, after months or years of growth, is a process of evolution by natural selection in a tumor. To quantitatively model this evolutionary process, researchers used evolutionary game theory. An important first step for a mathematical model is to quantify the ability to repro- duce, which is defined as fitness of a subpopulation of cells. Evolutionary game theory constructs the fitness landscape as the system evolves to an evolutionary stable state (ESS) [108, 52] based on the interaction of the subpopulation of cells. More specifically, the fitness function will have linear dependence [85] on the population size of surrounding cells, thus a frequency dependent selection. Compared with constant selection with fixed fitnesses, evolu- tionary game theory models the adaptation to the changing sizes of competing populations. There are many successful models that are based on either stochastic simulations with evo- lutionary games, such as the Moran process[69, 85], or deterministic simulations using the replicator equations[52] or Lotka–Volterra equations[72]. In this thesis, I will develop mod- els of adaptive chemotherapeutic response of a tumor based on the replicator equations for Chapter 2-5. We will use techniques based on adaptive, time-dependent control of nonlinear dynamical systems. An important mechanism we need to quantify is how tumors respond to chemother- apy. Tumor recurrence for a patient undergoing treatment is one of the main clini- cal challenges, and the primary reason of tumor regrowth is the development of chemo- therapeutic resistance[58, 1, 17]. As an evolutionary process, the drugs will alter the fitness landscapes[120] by reducing the fitness of drug-sensitive cells. Initially, due to the reduction of sensitive cells, the regression of the tumor is observed quickly. After the sensitive cells are eliminated, the drug-resistant cells (or resistant cells) are favored by the selection process and the tumor will start to regrow.This process from treatment to regrowth of the tumor is called competitive release phenomenon in the ecological literature [44, 104, 74], because the absence of drug-sensitive cells contributes to the growth of drug-resistance cells as they are 2 released from competition. In order to prevent the recurrence of the tumor, one approach is to find some adaptive chemotherapeutic schedule [40] for determining when to give the drug to the patients and how much drug should be administered. We use our evolutionary game theory model with the altered fitness from chemotherapy, to investigate the phase space of the system governing the interacting cells with different chemo-dose levels. In the phase space, we are able to construct the desired trajectory (steering and altering the balance of cells) from one point to another, and even design closed loops which keeps the cells in per- petual competition without any of the subpopulations saturating the tumor. The adaptive chemotherapy schedule corresponding to these closed loops can keep the tumor burden low and prevent or delay resistance from emerging. Compared with the classical Maximum- tolerated-Dose method (On-off piecewise constant function) or Low-Dose-Metronomic (con- stantvalue)method, whichwillsaturatewithlargetumorsizewithinfinitetime, theadaptive therapy outperforms both of them by trapping the tumor size below a certain level. Discus- sions and details for this model are provided in Chapter 2. This method is then compared with patient data of tumor relapse in Chapter 3. In Chapter 4, we extend our evolutionary game theory model to study multi-drug combi- nation therapies with two types of resistant populations. The goal is to achieve an adaptive schedule of two drugs to avoid competitive release of both subpopulations of resistant cells. The existence of more than one drug, however, can result in drug-drug interactions. These interactions can be additive, synergistic, or antagonistic. In addition, with a two-drug sched- ule, there are many new possible combinations that emerge. The discussion of multidrug interactions goes back to the 1930s[14]. These works help to determine the cocktails of multi- ple drugs and reach the desired outcome for treatments[50, 57, 11, 100, 3, 25, 71]. Synergistic effects are commonly assumed to be the best if the goal is to kill the maximum amount of cancer cells since synergistic interactions are analogous to adding extra drug dose. However for an evolving tumor, maximum killing of sensitive cells enhances selection for resistant cells and can actually encourage tumor resistance. Although tumor regression often occurs early in the treatment cycle, tumor recurrence inevitably occurs and our goal is to delay this event 3 to the maximum extent possible. After we design the adaptive control method based on the finding of closed loops in phase space, we show that antagonistic effects can provide better tumor regression for an evolving system as compared with the synergistic effects, mostly because the process of selection for the resistance population is slowed. After we design feasible schedules of controlling the tumor size as mentioned above, the next goal is to optimize the different schedules. For example, we may aim at producing a schedulethatresultsinthesmallestpossibletumor(volume)attheendofthetreatment. The question then becomes an optimal control[66] problem in the replicator dynamical system. Typical methods for this problem involve Pontryagin Maximum Principle (PMP) developed in 1956 by the Russian mathematician Lev Pontryagin and his students[91, 98, 21, 113]. It is common to assume that the optimal control for a replicator equation is a bang-bang control [98,21], whileforstochasticsimulationwithMoranprocess, itisreported[122]tohaveahigh entropy schedule (Low-Dose Metronomic schedule) to be optimal for minimal tumor volume. The continuous limit of the Moran process with infinite population size is proved [114, 115] to be the adjusted replicator equation, a rescaled form of replicator equation. In Chapter 4 we apply the maximum principle to a general 2-by-2 evolutionary game along with both the standard replicator equation and the adjusted replicator equation so we can compare with the stochastic model. We find that the behaviors of these two dynamic systems are quite different depending how the drug alters the fitness landscape. The optimal schedule for the standard replicator equation to minimize tumor size is bang-bang control, while for the adjusted replicator system, may not be piecewise constant if the drug favors one type of species instead of reducing the fitness of the other one. Details are presented in Chapter 5. Finally, in Chapter 6 we use our time-dependent control method for a different problem of levitation of a mass by non-constant blowing schedules using the optimal control techniques. This work is inspired by the time-dependent optimal drug schedules previously used to prevent tumor recurrence. We present a model of a mass (say a ping-pong ball) blown by streaming gas and floating in air. We show that constant blowing can never levitate the ball from its starting point to a floating point in mid-air. Our goal is to design a time-dependent 4 blowing schedule (analogous to our adaptive chemotherapy schedule) to levitate the ball to a target point that finds its rest point in minimal time. We prove that bang-bang control is optimal for minimizing time and streaming gas volume, while minimizing energy requires extrasegmentscorrespondingtothesingularcontrolfromoptimalcontroltheory. Theresults are identical to the solution from Pontryagin Maximum principle. For more details about how the problem is analyzed and discussed, see Chapter 6. 5 Chapter 2 Capitalizing on Competition: An Evolutionary Model of Competitive Release in Metastatic Castration Resistant Prostate Cancer Treatment The development of chemotherapeutic resistance resulting in tumor relapse is largely the consequence of the mechanism of competitive release of pre-existing resistant tumor cells selected for regrowth after chemotherapeutic agents attack the previously dominant chemo- sensitive population. We introduce a prisoner’s dilemma game theoretic mathematical model based on the replicator of three competing cell populations: healthy (cooperators), sensitive (defectors), and resistant (defectors) cells. The model is shown to recapitulate prostate- specificantigenmeasurementdatafromthreeclinicaltrialsformetastaticcastration-resistant prostate cancer patients treated with 1) prednisone, 2) mitoxantrone and prednisone and 3) docetaxel and prednisone. Continuous maximum tolerated dose schedules reduce the sensitive cell population, initially shrinking tumor burden, but subsequently “release" the resistantcellsfromcompetitiontore-populateandre-growthetumorinaresistantform. The evolutionary model allows us to quantify responses to conventional (continuous) therapeutic strategies as well as to design adaptive strategies. These novel adaptive strategies are robust to small perturbations in timing and extend simulated time to relapse from continuous therapy administration. 6 2.1 Introduction In his now classic 1961 study of competition for space between two species of barnacles in the intertidal zone off the Scottish coast, Joseph Connell [27] discovered something inter- esting. The blue barnacles Balanus normally occupied the intertidal zone, while the brown barnacles Chthamalus occupied the coast above high tide. Despite the commonly held belief that each occupied their own niche because of different adaptations to local micro-conditions, Connell hypothesized that the colonization of the intertidal zone by Balanus was actually preventing Chthamalus from inhabiting this region. To test this, he removed the blue bar- nacles from the intertidal zone and tracked the subsequent penetration of Chthamalus into this region. He concluded that Chthamalus had undergone relief from competition with Bal- anus which allowed it to flourish where previously it could not. The point, he emphasized, was there was nothing inherent about the micro-environment of the intertidal zone that was preventing Chthamalus from occupying this region. It was simply the competition against a more dominant species that was holding it back. Without the presence of that species, Chthamalus happily claimed both zones as fundamental niches. Thus, the important notion of competitive release was formulated (see Grant [45]). When two (or more) sub-species com- pete for the same resources, with one species dominating the other, if the dominant species is removed, this can provide the needed release from competition that can allow the less dom- inant species to flourish. The mirror image of competitive release is the related notion of character displacement developed by Brown and Wilson [20] in which competition can serve to displace one or more morphological, ecological, behavioral, or physiological characteristics of two closely related species that develop in close proximity. These concepts are now well established as part of the overall framework of co-evolutionary ecology theory and play an important role in the evolution of chemotherapeutic resistance in cancer [74, 6, 46, 32]. Since co-evolution among competing clones and subclones is now a well established [86] process in malignant tumors, the mechanism of competitive release should be expected to play a role and affect the chemotherapeutic strategies one might choose to eliminate or control tumor growth. Indeed, tumor relapse and the development of chemo-therapeutic 7 resistance is now thought largely to be a consequence of the evolutionary mechanism of competitive release of pre-existing resistant cells in the tumor which are selected for growth after chemotherapeutic agents attack the subpopulation of chemo-sensitive cells which had previously dominated the collection of competing subpopulations. This concept, perhaps most forcefully advocated by Gatenby and collaborators [38], is gaining acceptance by clin- icians. A recent (2012) systematic literature analysis of cancer relapse and therapeutic research showed that while evolutionary terms rarely appeared in papers studying therapeu- tic relapse before 1980 (< 1%), the language usage has steadily increased more recently, due to a huge potential benefit of studying therapeutic relapse from an evolutionary perspective [1]. Anticancer therapies strongly target sensitive cells in a tumor, selecting for resistance cell types and, if total eradication of all cancer cells is not accomplished, the tumor will recur as derived from resistant cells that survived initial therapy [88]. It is argued by Gatenby [38] that eradicating most disseminated cancers may be impossible, undermining the typical goal of cancer treatment of killing as many tumor cells as possible. The underlying assumption of this approach has been that a maximum cell-kill will either lead to a cure or, at worst, max- imum life extension. Taking cues from agriculturists who have long abandoned the goal of complete eradication of pests in favor of applying insecticides only when infestation exceeds a threshold in the name of “control" over “cure," there are those who advocate for a shift from the cure paradigm in cancer treatments to a control paradigm [38, 12]. 2.1.1 The likelihood of pre-existing resistance Pre-existing resistant cells should generally be present in all patients with late-stage metastatic disease (for single point mutations which confer resistance), a conclusion sup- ported by probabilistic models [17] and from tumor samples taken prior to treatment [61, 96] which have been reported for melanoma [56], prostate cancer [97], colorectal cancer [29, 62], ovarian cancer [103], and medulloblastoma [77]. According to this view, treatment failure 8 would not be due to evolution of resistance due to therapy, but rather the pre-existing pres- ence of resistant phenotypes that are relatively sheltered from the toxic effects of therapy [32]. The likelihood of pre-existing resistance has important therapeutic implications. If we assume no pre-existing resistance, then most models predict maximum dose-density therapy will reduce the probability of resistance largely because this treatment minimizes the number ofcell-divisions, therebyminimizingtheriskofamutationleadingtoacquiredresistance[32]. By contrast, in pre-existing resistance scenarios, the maximum dose-density therapy strat- egy lends itself to competitive release due to the evolutionary nature of tumor progression. Most pre-clinical efforts that aim to maximize the short-term effect of the drug on sensitive cells does not significantly affect the long-term control of cancer [17]. This is because the phenomenon of competitive release can occur via the harsh selective pressure imposed by the tumor microenvironment after cancer therapies diminish the presence of the dominant (i.e. the chemo-sensitive) subpopulation. Additionally, the process of metastasis may allow a resistant subpopulation in the primary tumor to emerge elsewhere [116]. 2.1.2 Leveraging the cost of pre-existing resistance Pre-existing mutations that are responsible for conferring resistance may be associated with a phenotypic cost, or a reduced fitness, compared to the average fitness of the sensitive cell population [40, 41]. Even factoring in this fitness cost, deleterious mutations are still expectedtobepresentinlate-stagemetastaticcancers[16]. Thiscostcancomeinmanyways, such as an increased rate of DNA repair, or an active pumping out of the toxic drug across cell membranes. All of these strategies use up a finite energy supply that would otherwise be availableforinvasionintonon-canceroustissuesorproliferation. Therapidremovalofchemo- sensitive cells during therapy releases the resistant population from unwanted competition and thereby permits unopposed proliferation of the resistant cell population. In contrast, the goal of an adaptive therapy is to maintain a stable tumor burden that permits a significant 9 population of chemo-sensitive cells for the purpose of suppressing the less fit but chemo- resistant populations, consistent with the philosophy that it takes an evolutionary strategy to combat an evolving tumor [116]. There is some evidence that adaptive therapies may be beneficial for a range of fitness cost (or lack thereof) and therapeutic sensitivity [36, 49]. A theoretical framework for these adaptive therapies first developed by Gatenby [40], leverages the notion that pre-existing resistance is typically present only in small population numbers due to a cost of resistance. This less fit phenotype is suppressed in the Darwinian environment of the untreated tumor but treatments that are designed to kill maximum num- bers of cells remove the competition for the resistant population and ultimately select for that population during tumor relapse 1 . While the goal of an adaptive therapy (to capitalize on the competition resistant subpopulations incur through maintaining a stable sensitive cell population) has gained some level of acceptance, the ideal adaptive therapy schedule in practice, fordifferenttumortypesandgrowthratesisfarfromsettled. Gatenby’spapermod- ulates dose density only, while stating that an ideal adaptive therapy would also modulate the drug, as well as the dose schedule (both dose and density) [40]. Some of these evolution- ary ideas were tested experimentally using mouse models with modulated dose strength and dose vacations designed to maintain a stable, controllable tumor volume [31, 105]. This two- phase adaptive therapy involved an initial high-dose phase to treat the exponential growth of the tumor and a second phase designed for stable tumor control using a variety of strate- gies (such as decreasing doses or skipping doses when stability is achieved). Several spatial, agent-based computational models have modulated dose strength with respect to a thresh- old value of tumor size (a fraction of the original tumor burden) [36, 8]. Findings suggest that adaptive therapies based on evolutionary treatment strategies that maintain a resid- ual population of chemo-sensitive cells may be clinically viable, and is currently extended to an on-going clinical trial (NCT02415621) which adaptively controls the on/off cycling of 1 Itisimportanttonotethatbothhigh-dose,maximumtolerateddoseschedulesandlow-dose,metronomic dose schedules have this cumulative goal of achieving maximum cell-kill over the course of many cycles of treatment. 10 abiraterone [134]. We refer the reader to similar approaches used in mitigating antibiotic resistance [4, 48]. With these advances in mind, the goal of this manuscript is to use the evolutionary framework introduced and advocated over the past 10 years [40, 38, 36, 8, 134] to mathe- matically model the important parameters of competitive release and use that framework to better understand therapeutic implications of the cost of developing resistance and to learn how to exploit competition between subpopulations. Here, we ignore clonal dynamics, mutations between subpopulations, and non-genetic (rapid, bidirectional) state transitions, to focus on the inter-subpopulation competition [106, 51, 90]. Specifically, we introduce a three-component replicator system with a prisoner’s dilemma payoff matrix [119] to model the three relevant subclonal populations: healthy cells (H), sensitive cells (S), and resistant cells (R). We consider healthy cells as the surrounding non-neoplastic cells that do not con- tribute to the overall growth of the tumor, but nevertheless compete with neoplastic cells for space and resources. Using the nullcline information in a triangular phase plane repre- sentation of the nonlinear dynamics of the system, we first show the essential ingredients that render competitive release possible. Then, using the parameters that control selection pressure (hence relative growth rates) on the three subclonal populations, we attempt to maintain the tumor burden at low levels so that the resistant population does not reach fixation. 2.1.3 Retrospective analysis of metastatic castration-resistant prostate cancer A recent retrospective analysis of tumor measurement data (PSA levels) from eight ran- domized clinical trials with metastatic castration-resistant prostate cancer (mCRPC) used a simple linear combination of exponentials model to estimate the growth and regression rates of disease burden over time [127]. In total, over 67% of patients were fit to models with a positive regrowth rate, indicating failure due to resistance. Prostate-specific antigen (PSA) measurement data for patients in each treatment silo (prednisone only, mitoxantrone 11 Healthy Sensitive Resistant 3 species model of competitive release therapy Tumor is dominated by sensitive cells which out- compete both healthy cells and resistant cells. Pre-treatment After chemotherapy the sensitive cell population is reduced, leaving resistant cells free to re-populate without competition. Resistant cells eventually repopulate a larger proportion of the tumor, rendering the subsequent rounds of chemotherapy less effective. re-growth Post-treatment Tumor re-growth Figure 2.1: Schematic of competitive release in a tumor — (a) Prior to treatment, a tumor consists of a large population of sensitive cells (red) and a small population of less fit resistant cells (green) competing for resources with the surrounding healthy cells (blue); (b) Chemotherapy targets the sensitive population (middle), selecting for the less fit resistant population that thrives in the absence of competition from the sensitive population; (c) Upon regrowth, the tumor composition has larger numbers of resistant cells, rendering the subsequent rounds of treatment less effective. and prednisone, docetaxel and prednisone) were obtained through the Project Data Sphere open data portal (https://www.projectdatasphere.org), and we show that this model is able to adequately fit data for each treatment type with the additional capability of allowing us to track responses to conventional therapeutic strategies and design new adaptive strategies as the tumor evolves. The model can be used to test previously proposed adaptive thera- pies, but we propose a novel schedule utilizing quantitative tools from nonlinear dynamical systems theory which use the current global state of the nonlinear replicator system with respect to the nullcline curves of the equations as well as parameters controlling relative fitness levels of the competing sub-populations. The simulated chemotherapeutic strategies that we implement, based on tracking the phase-space structure of of replicator system, are 12 ones that can adapt on the same timescale as the inherent timescale of evolution of the sub- populations comprising the tumor, i.e. are as dynamic as the tumor. While this approach cannot be preplanned by the oncologist at the beginning of therapy like classical strate- gies, we provide discussion to explain how an evolutionary game theory model describing the fitness landscape (described below) is useful to understand the underlying features of a dynamical fitness landscape associated with a cost to resistance: a three-way prisoner’s dilemma. Specifically, the model indicates a boundary over which an adaptive therapy will cease to be effective. Our model focuses on cost to resistance, as opposed to specific mecha- nisms of resistance (i.e. hormone resistance, for example). This has the advantage in some ways as being agnostic to resistance mechanisms, but on the other hand, in this manuscript, we can not distinguish between different mechanisms: the analysis explores the effect of a fitness competition between subpopulations on the competitive release phenomenon.Three common therapies used in mCRPC do show that a fit of our model to prostate data lead to theemergenceofacostofresistanceinamajorityofpatients, withoutpredeterminingavalue to the cost in the model a priori. Similarly, ref. [127] found that continuous administration of different therapies (i.e. different resistance mechanisms) is associated with slow tumor regrowth rates (i.e. a cost of resistance conferred to growth potential). We now describe the model. 2.2 Materials and Methods A schematic of a finite-population three component model (healthy, sensitive, resistant) of competitive release is shown in figure 2.1, where the tumor consisting of sensitive and resistant cells is competing with the surrounding healthy tissue. At diagnosis (see figure 2.1, left), the tumor is dominated by sensitive cells (red) which out compete the surrounding healthy population (blue) during unhindered tumor progression. A small portion of resistant cells (green) remains in small numbers, suppressed by the larger sensitive population. After several rounds of chemotherapy, the tumor shrinks, leaving the resistant population largely unaffected (figure 2.1, middle). Inevitably, the tumor relapses due to the small number of 13 Figure 0 total tumor burden sensitive subpopulation resistant subpopulation first driver mutation exponential growth of chemo-sensitive population first resistant mutation begin continuous chemotherapy treatment good initial response to chemotherapy competitive release of resistant population Chemo-sensitive Chemo-resistant time tumor burden Figure 2.2: Clonal evolution of competitive release — A fishplot (sometimes known as a Müller plot), showing the tumor quantity (vertical axis) and composition (sensitive: red; resistant: green) over time (horizontal axis, left to right) with important events annotated. After first driver mutation (left), initial exponential growth of sensitive population occurs until diagnosis (dashed line). Continuous therapy targeting the chemo-sensitive population responds well with a decrease in tumor burden. In the absence of sensitive cells, the resistant population (existing in small numbers before the start of therapy) grows to become the dominant subpopulation at relapse, albeit typically with lower exponential growth rate due to the cost of resistance. sensitive cancer cells remaining after therapy (figure 2.1, right). In the absence of competi- tion from the dominant sensitive population, the resistant cells grow unhindered, rendering subsequent rounds of chemotherapy less effective. Subsequent application of identical thera- pies will have a diminished effect. Figure 2.2 shows the process in a ‘Müller fishplot’, which we will use later to track the subclonal populations. This representation was first utilized in cancer to compare modes of clonal evolution in acute myeloid leukemia (see [30]) and is particularly useful in computational models where every cell type can be tracked in time. A fishplot shows the tumor burden (vertical axis) over time (horizontal axis). The first resis- tant cell is assumed to arise from a single sensitive cell, undetectable by current methods of measurement. The schematic depicts unhindered tumor growth after the first driver muta- tion (figure 2.2, left) where the tumor grows exponentially before diagnosis, during which time a resistant mutation occurs (figure 2.2, middle). After diagnosis (dashed line), a regi- men of continuous chemotherapy shows initial good response and tumor regression, but the 14 resistant population grows back (although at a slower growth rate) unhindered by competi- tion, leading to relapse (figure 2.2, right). Previously, a linear combination of exponentials model has been proposed to track the tumor quantity via prostate-specific antigen (PSA) measurement data [23]. The tumor quantity q(t) is a function of the exponential death rate of the sensitive cells, d, the exponential growth rate of the resistant cells, g, and the initial fraction of resistant cells, f [17]. In this equation, q(t) is the tumour burden at time t in days, normalized to the PSA level at t=0, (see figure 2.3, blue lines). The model can be written as follows: q(t) = (1−f)e −dt +fe gt . (2.1) This model, shown to be a reasonably good description of the changing tumor size dur- ing therapy for colorectal, prostate, and multiple myeloma cancers, identifies the important parameters in competitive release: initial fractional resistance (f), and birth/death rates (g,d) for the resistant and sensitive populations, respectively. The model is used to fit prostate-specific antigen (PSA) measurement data from retrospective analysis of three ran- domized clinical trials with metastatic castration-resistant prostate cancer to estimate the growth (g) and regression (d) rates of disease burden over time. Four representative patients are chosen from the control arms of each randomized trial and shown in figure 2.3: treatment with prednisone only [75] (left column: figures 2.3B,E,H,K); treatment with mitoxantrone and prednisone [112] (middle column: figures 2.3C,F,I,L) and treatment with docetaxel and prednisone [89] (right column: figures 2.3D,G,J,M). PSA data and model fits are normalized by q(t = 0) (black dots) and exponential fits are shown in blue. Despite the fact that the model (2.1) curve-fits data reasonably well (labeled “exp." in figure 2.3), it contains no evolutionary information or concepts, a keystone principle behind competitive release. We also include in the figure 2.3 fits of our model presented here (eqns. 2.2, 2.3) in red (labeled “rep." in figure 2.3, where ˜ q(t) is given by eqn. 2.4, below). Our evolutionary model is able to capture similar trends as the exponential model of equations (eqn. 2.1) but has the important property of allowing us to calculate cost of 15 resistance associated with model fit. While each drug (columns) may have very different resistance mechanisms, ten of the twelve patients shown have a cost of resistance (α−β, described below). In the last section we will describe the dynamical phase portrait and implement adaptive strategies to capitalize on competition between resistant and sensitive subpopulations in order to maintain tumor control. 2.2.1 The replicator equation model The dynamics of the fitness landscape of three competing cell types are described by the replicator equation (see [114]), which is a deterministic birth-death process in which birth and death rates are functions of cell fitness, and cell fitness is a function of prevalence in the population. Each ith cell type (i = 1, 2, 3) competes according to equation 2.2, where x 1 , x 2 , x 3 are the corresponding frequency of healthy (H), sensitive (S) and resistant (R) cells, respectively, such that P i x i = 1. ˙ x i = (f i −hfi)x i (2.2) f i = 1−w i +w i (A~ x) i (2.3) Here, ~ x is the vector ~ x = (x 1 ,x 2 ,x 3 ) T and (Ax) i is the ith element of vector Ax. The prevalence of each sub-population, x i , changes over time according to the chang- ing population fitness, f i , as compared to the average fitness of all three populations hfi = f 1 x 1 +f 2 x 2 +f 3 x 3 . If the fitness of the sub-population is greater than the aver- age fitness (f i −hfi > 0), that sub-population grows exponentially, whereas if it is less (f i −hfi< 0), it decays. In order to directly compare results to previously published models of tumor quantity, q(t), during competitive release (i.e. equation 2.1), the tumor quantity can be written as the sum of the resistant and sensitive tumor populations (x 2 ,x 3 respectively), ˜ q(t), 16 ˜ q(t) =x 2 (t) +x 3 (t), (2.4) normalized by it’s initial condition, ˜ q(0). This is shown in a schematic in figure 2.3A. These represent the same subpopulations shown in green (resistant) and red (sensitive) in figure 2.2. Before therapy, each subpopulation (healthy, chemo-sensitive, and resistant cells) the selection pressure is constant across all cell types (i.e. w i ≡ w, i = 1, 2, 3) at a level that represents the natural selection pressure the tumor environment imposes on the different subpopulations. These values discussed in the literature are typically small, in the range w i ≡w≈ 0.1− 0.3. This is a weighting parameter between neutral dynamics (w = 0) and fitness-dependent dynamics (w = 1). Typically in replicator systems, this parameter scales out with time. In our model, we treat each subpopulation selection parameter independently. We implement chemotherapy in our model by changing the selection pressure parameters on each of the subpopulations of cells. Therapy can be administered at different doses (i.e. values of the drug concentration: c; 0≤ c≤ 1). A higher value of c indicates a stronger dose of chemotherapy drug (described in more detail in [122]). This follows the schematic in figure2.4whichdepictsthechangeinthefitnesslandscapebeforeandaftertherapy. Infigure 2.3, dose concentration is assumed to be constant for a specific drug while patient-specific parameters are the selection pressure (w), cost of resistance (discussed below), and initial fraction of resistant cells (f). Values are altered as follows (see figure 2.4 for explanation of changing fitness landscape): w 1 = (1 +c)w (healthy) (2.5) w 2 = (1−c)w (sensitive) (2.6) w 3 = w (resistant) (2.7) 17 The fitness landscape (eqn. 2.3) is described in detail by the entries of the payoff matrix A, A = H S R H S R a b o h j k l m n (2.8) where each pairwise cell-cell interaction is described by the row and column values, which are parametersinthefitnessequation(2.3). The 3×3payoffmatrixisconstructedasastandard prisoner’s dilemma matrix where healthy cells are the cooperators and sensitive/resistant cells are the defectors. In the absence of resistant cells, we describe in more detail the dynamics of the system and its relevance to tumor growth in the appendix. The key features are the Gompertzian growth of the cancer cell population which saturates at a lower overall population fitness level than the initial all healthy cell population. This necessitates the following inequalities of the payoff matrix below (eqn. 2.8): h > a > j > b, l > a > n > o, and k > n > j > m. Note that these payoff entries remain constant before and during therapy: chemotherapy is viewed as affecting the selection balance among the three subpopulations. This is implemented in our model by a change in the the selection pressure parameters only (see eqns. 2.5, 2.6, 2.7). More discussion of why the prisoner’s dilemma matrix, which models the evolution of defection, is a useful paradigm for cancer can be found in [119, 118] and the appendix. 2.2.2 The linearized system and the cost of resistance An additional important feature of the payoff matrix is the notion of the cost of resistance which is highlighted in figure 2.4A. With no therapy, the sensitive cells exhibit fastest growth due to their higher fitness value relative to both the resistant population and the healthy population. The difference between the baseline fitness values of the sensitive cells and the resistant cells can be thought of as the ‘price paid’ by the resistant population to retain their resistance to toxins. This cost, in our model, is quantified as the difference in the (linearized) 18 growth rates of the two populations (type 2: sensitive; type 3: resistant). Linearizing eqn (2.2), (2.3), (which form a cubic nonlinear system if expanded out) gives rise to the sensitive- resistant uncoupled system: ˙ x 2 =αx 2 (2.9) ˙ x 3 =βx 3 , (2.10) with the growth parameters: α =w 1 (1−a) +w 2 (h− 1) (2.11) β =w 1 (1−a) +w 3 (l− 1). (2.12) Using eqns (2.5), (2.6), (2.7) gives: α =w(h−a)−cw(h +a− 2) (2.13) β =w(l−a) +cw(1−a). (2.14) With no therapy, c = 0, we have: α =w(h−a) (2.15) β =w(l−a). (2.16) We call the fitness cost of resistance the difference between these growth rates with no therapy, hence (α−β) = w(h−l), where we require h > l in the payoff matrix. However, the w parameter is patient-specific, allowing each patient a varied cost. 2.3 Results It is useful to view the nonlinear dynamical trajectories of the system using the trilinear coordinates shown in figure 2.5A, which gives a representation of the clonal phase space 19 for every possible value of ~ x [102]. The corners represent saturation of a single cell type (e.g. the top corner represents ~ x = [1, 0, 0], or all healthy cells). The healthy (H; top corner), sensitive (S; bottom left corner), and resistant (R; bottom right corner) populations compete according to equation (2.2) and follow trajectories shown (black) in figure 2.5C. Figure 2.5C shows the complete dynamical information for untreated patient dynamics (c = 0, α = 0.020, β = 0.018, w = 0.1). All internal trajectories (black lines) lead to tumor growth and eventual saturation of the sensitive population (bottom left corner, “S"). Instantaneous relative velocity is indicated by background color gradient (blue to red). Sub- population “nullclines," the curves for which ˙ x i = 0 are shown for healthy (dashed blue), resistant (dashed green), and sensitive (dashed red). On one side there is positive growth ( ˙ x i > 0); on the opposite side negative growth ( ˙ x i < 0). These nullclines delineate the phase space into three regions. • Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; • Region 2: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0; • Region 3: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0. Likewise, figure 2.5D shows the complete dynamical information for the patient treated with chemotherapy (c = 0.6,α = 0.020, β = 0.018, w = 0.1). Trajectories initially move away from the sensitive/resistant corners toward healthy, but subsequently relapse to the saturation of the resistant subpopulation. The phase space for treated dynamics is divided into six regions by the three nullclines. • Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; • Region 2: ˙ x H > 0 ˙ x S < 0 ˙ x R < 0; • Region 3: ˙ x H > 0 ˙ x S < 0 ˙ x R > 0; • Region 4: ˙ x H < 0 ˙ x S < 0 ˙ x R > 0; 20 • Region 5: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0; • Region 6: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0. An important detail emerges from the model: during treatment, the resistant nullcline is reached before the healthy nullcline (figure 2.5D). For example, a tumor with an initial diag- nosis in region 2 (see figure 2.5D) can be expected to respond to treatment. The trajectory will follow along the black line towards the healthy corner (i.e. decreasing tumor quantity) until the trajectory passes over the nullcline for the resistant subpopulation (dashed green line; ˙ x R = 0) reaching region 3. Subsequently, the trajectory passes the healthy nullcline (dashed blue line; ˙ x H = 0) and the tumor relapses ( ˙ x H < 0), this time toward saturation of the resistant subpopulation. If the desire is to maintain treatment until the point of positive progression, this occurs when the trajectory is far past the resistant nullcline. Appearances (based on tumor burden) can be deceiving – while the tumor may appear to be responding, the overall state may be well past the point of no return, and the resistant population is preparing to re-populate. 2.3.1 Managing competitive release Figure 2.5B shows a schematic of the two resistant nullclines from 2.5C and 2.5D (untreated: solid green line; treated: dashed green line) overlaid on the same phase por- trait. From the initial condition (purple dot) treatment begins as the trajectory (solid blue line) approaches the nullcline, positive resistant growth can be avoided by a well-timed drug holiday (solid red line). Here, we propose an adaptive therapy with the goal of trapping the trajectory between the treated resistant nullcline and the untreated resistant nullcline, cre- ating an orbit in a closed loop for a finite period of time to ‘steer’ tumor evolution in region 2 of figure 2.5B. We do this by altering the dose concentration parameter c (a parameter that can be accessed clinically) in eqns (2.5) in an off-on (bang-bang) fashion, (eqn. 2.6, 2.7) between the two resistant nullclines: the untreated nullcline from figure 2.5C and the treated nullcine from 2.5D. 21 While all details of the ‘tumor phase space’ may not yet be directly measurable in the clinic, we propose that all successful adaptive therapies will operate in regions 1 and 2 in figure 2.5D; otherwise they ultimately will not be successful. In this way, meaningful insight is gained into the dynamics behind the cost to resistance, regardless of the mechanism of that resistance. Previous adaptive therapy schedules (described in the introduction; see [36, 8, 31, 134]) have benefited from not crossing this “hidden" resistant nullcline. We now propose one example of an adaptive therapy schedule which actively captures and uses information from the dynamic phase space. The simple control paradigm proposed indirectly controls the resistant population by systematically choosing when to administer therapy and when to give drug holidays. These holidaysallowasufficientnumberofsensitivecellstoremaininordertosuppresstheresistant population. A continuous dose of therapy is administered until the treated nullcline ( ˙ x R = 0) is reached (see figure 2.5D, green dashed line). This is the starting point of positive growth for the resistant population (further therapy would result in ˙ x R > 0). A drug holiday is then imposed until the second nullcline is reached (see figure 2.5C, green dashed line). The sensitive population is allowed to regrow until it is large enough to suppress the resistant population once again (and when ˙ x R = 0). Therapy is administered to allow the tumor to cycle back and forth between the two nullclines. This bang-bang (on-off) strategy allows an extension of relapse times. We emphasize that the specific times we turn the therapy on and off in this bang-bang strategy cannot be pre-planned, but depend on the position of the trajectory in the tumor phase space as the disease evolves. Next, we measure the success of this control paradigm compared to continuous treatment viatwoimportantmeasuresoftherapyeffectiveness: progressionfreesurvival(PFS)andtime to relapse. PFS is the time to which the patient has a measured response, before relapse of the tumor: shown in figure 2.6A, bottom. A more useful metric is “time to relapse," where disease reaches initial state. Measuring the effectiveness of a chemotherapy schedule based on the killing rate or progression free survival alone are not sufficient predictive measures of long-term cancer control [17]. PSA data from a single patient who relapsed due to treatment 22 resistance (figure 2.3B) is replotted in figure 2.6A (black dots) along with the model best fit of continuous treatment (blue; therapy 1). Also shown are simulated new drugs with increased effectiveness of killing sensitive cells (red, yellow) simulated by increasing the effective dose concentration with identical patient-specific parameters (from the patient in figure 2.3B). Each increased dose corresponds to a slightly shorter PFS, but an increased time to relapse to the initial tumor burden. However, despite the increase in relapse times, none of these doses optimizes tumor control, as seen in the fishplots (figure 2.6, B,C,D). At the point of relapse to the initial tumor burden, the tumor is dominated by the presence of resistant subpopulations (green), rendering future treatments ineffective. Oftentimes, the effectiveness of a new chemotherapy drug is be determined by PFS times when drugs that have high killing rates of sensitive cells may have shorter times to progression and lower total tumor burden at all times (everything else equal). The figure clearly shows that all treatments have similar progression free times but with a greater range of relapse times (even though continuous treatment always eventually leads to relapse). Next, compare continuous treatments to the proposed control paradigm (solid black line) for identical initial conditions and identical drug dose. As seen in the fishplot (figure 2.6E), the resistant population (green) is suppressed during drug holidays, leading to an extended time without relapse. This adaptive technique is successful for two reasons. First, the drug holidays allow an adequate sensitive population size to suppress the growth of the lower-fitness resistant population. Second, the resistant population is never allowed to reach a positive growth under treatment ( ˙ x R > 0), therefore cannot take over the tumor cell population. In order to test the robustness of this method, similar control paradigms were simulated with error in administration times: slightly shorter time on therapy (1 day; figure 2.6A; green) and slightly longer (1 day; figure 2.6A, purple). All else equal, it is better to err on the side of longer treatment on periods without decreasing off treatment time (purple). 23 2.4 Discussion The chemotherapeutic scheduling strategies outlined in this paper cannot be pre-planned by the oncologist at the beginning of therapy like classical strategies [84], as they rely on significant decision making and continuous monitoring of the different subpopulations of cells that co-evolve as the tumor progresses. This means that the quality of the cell population monitoring system is crucial to the entire strategy, as has been pointed out in [33]. There can be no adaptive tumor control strategy without continuous monitoring of the subpopulations as it is not the tumor burden that is of primary interest, but the heterogeneous balance of the subpopulations comprising the tumor. In addition, the information gleaned from a detailed monitoring system cannot be acted upon unless the various administered drugs are sufficiently targeted to act efficiently and exclusively on specific subpopulations . These two systems must be in place (sensing and actuating) in order to successfully shape the fitness landscape and steer a growth trajectory in a desired direction. We also want to emphasize a separate point, which is that it is not enough to know in detail the current state of the system in order to steer it successfully. One must also have a description of all possible nearby states of the system, both under therapeutic pressure and without therapy. Better yet is to have a global picture of all possible states of the system, with nonlinear nullcline information, as one would obtain by analyzing the full phase space of the entire system. With this information, one would know where to steer the system to get to a desired state, even if one does not know how to achieve this (clinically). In current state-of-the-art medical practice, such sophisticated sensor-actuator capability is not yet sufficiently developed as it is in many engineering contexts where versions of adaptive control theory are routinely used. Many similar challenges, and the necessary steps towards their implementation, present themselves in the ecology and pest control communities, and we point to Gould’s article [44] for a nice early overview. More recently, connections between the approaches developed in the past by ecologists and possible future strategies for oncologists have been discussed by Gatenby and collaborators [39]. Other groups [64, 18, 65] have also developed highly mathematical approaches to tumor control from different points of view. Clearly not all 24 of the clinical steps are in place to effectively test and implement many of the strategies that have been explored theoretically. Yet it is still important to continue to develop the kinds of mathematical models and computer simulations that would serve to identify the many possible schemes, parameter ranges, and sensitivities that could be tested via clinical trials that focus on adaptive therapies with the goal of suppression of potential evolution of resistance. 25 0 50 100 150 200 250 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 exp. rep. 0 100 200 300 400 0 0.2 0.4 0.6 0.8 1 1.2 exp. rep. 0 100 200 300 400 500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 50 100 150 200 250 300 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 0.2 0.4 0.6 0.8 1 1.2 1.4 exp. rep. 0 100 200 300 400 500 600 0 0.2 0.4 0.6 0.8 1 1.2 exp. rep. 0 50 100 150 200 250 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 1.2 exp. rep. 0 50 100 150 200 250 300 350 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. 0 100 200 300 400 500 600 700 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 exp. rep. Prednisone Mitoxantrone and Prednisone Docetaxel and Prednisone B C D E F G H I J K L M A Normalized sum of growth and regression Normalized Normalized PSA data Sensitive decay, x 2 (t) Resistant growth, x 3 (t) ˜ q(t)/˜ q(0) <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="KooCSFWDasnL+KUxvM75EPzAAqk=">AAAB9HicbZDNSgMxFIXv1L86Vq1rN8EiuCozbnQpuHFZwf5AO5RMJtOGJpkhuSOUoSt3bn0Kd+Lb+BC+g+m0C2s9EDick3BvvjiXwmIQfHm1nd29/YP6oX/U8I9PTpuNns0Kw3iXZTIzg5haLoXmXRQo+SA3nKpY8n48u1/2/WdurMj0E85zHik60SIVjKKLOuNmK2gHlci2CdemBWuNm9+jJGOF4hqZpNYOwyDHqKQGBZN84Y8Ky3PKZnTCh85qqriNymrNBbl0SULSzLijkVTp7xclVdbOVexuKopT+7dbhv91wwLT26gUOi+Qa7YalBaSYEaWfyaJMJyhnDtDmRFuV8Km1FCGjszGFKZiIyZTXPiVHJ7wL4xt07tuh0E7fAygDudwAVcQwg3cwQN0oAsMEniFN+/Fe/c+Vhhr3prnGWzI+/wB8mOTuA==</latexit> <latexit sha1_base64="YtijgKxQECAD0Ef77P7BSGnWD4o=">AAACEnicbZDLSgMxGIX/8Vpr1dGduBksQrupM250KbhxWcFeoB1KJpNpQzMXk3+EMgw+iiu3+hTuxL0P4TuYTgva1gOBk3MS/uTzEsEV2vaXsba+sbm1Xdop71b29g/Mw0pbxamkrEVjEcuuRxQTPGIt5ChYN5GMhJ5gHW98M+07j0wqHkf3OEmYG5JhxANOCepoYB73kQufZQ95Devnvxu7PjCrdsMuZK0aZ26qMFdzYH73/ZimIYuQCqJUz7ETdDMikVPB8nI/VSwhdEyGrKdtREKm3Kz4Qm6d6cS3gljqFaFVpH9vZCRUahJ6+mRIcKSWu2n4X9dLMbhyMx4lKbKIzgYFqbAwtqY8LJ9LRlFMtCFUcv1Wi46IJBQ1tYUpNPQkH44wLxfSfJxlGqumfdFw7IZzZ0MJTuAUauDAJVzDLTShBRSe4AVe4c14Nt6NjxnJNWOO9AgWZHz+AMS6nts=</latexit> <latexit sha1_base64="YtijgKxQECAD0Ef77P7BSGnWD4o=">AAACEnicbZDLSgMxGIX/8Vpr1dGduBksQrupM250KbhxWcFeoB1KJpNpQzMXk3+EMgw+iiu3+hTuxL0P4TuYTgva1gOBk3MS/uTzEsEV2vaXsba+sbm1Xdop71b29g/Mw0pbxamkrEVjEcuuRxQTPGIt5ChYN5GMhJ5gHW98M+07j0wqHkf3OEmYG5JhxANOCepoYB73kQufZQ95Devnvxu7PjCrdsMuZK0aZ26qMFdzYH73/ZimIYuQCqJUz7ETdDMikVPB8nI/VSwhdEyGrKdtREKm3Kz4Qm6d6cS3gljqFaFVpH9vZCRUahJ6+mRIcKSWu2n4X9dLMbhyMx4lKbKIzgYFqbAwtqY8LJ9LRlFMtCFUcv1Wi46IJBQ1tYUpNPQkH44wLxfSfJxlGqumfdFw7IZzZ0MJTuAUauDAJVzDLTShBRSe4AVe4c14Nt6NjxnJNWOO9AgWZHz+AMS6nts=</latexit> <latexit sha1_base64="fHWoOWp1xLe1izqTRWjO6o0EvhU=">AAACHXicbVDLSsNAFJ3UV62vqCtxEyxCu6mJG10W3bisYB/QhjKZTNuhM0mcuRFKCH6KK7f6Fe7ErfgR/oPTNKBtPXDhzDn3cuceL+JMgW1/GYWV1bX1jeJmaWt7Z3fP3D9oqTCWhDZJyEPZ8bCinAW0CQw47USSYuFx2vbG11O//UClYmFwB5OIugIPAzZgBIOW+uZRDxj3aXKfVqB69vuwq32zbNfsDNYycXJSRjkaffO754ckFjQAwrFSXceOwE2wBEY4TUu9WNEIkzEe0q6mARZUuUl2QmqdasW3BqHUFYCVqX8nEiyUmghPdwoMI7XoTcX/vG4Mg0s3YUEUAw3IbNEg5haE1jQPy2eSEuATTTCRTP/VIiMsMQGd2twWIjzJhiNISxl0Ps5iGsukdV5z7Jpza5frV3lSRXSMTlAFOegC1dENaqAmIugRPaMX9Go8GW/Gu/Exay0Y+cwhmoPx+QOn4qBf</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> <latexit sha1_base64="M2sIkIItDRb47czJxAJgmkbwut4=">AAACHXicbVDLSsNAFJ34rPUVdSVugkVoNzURQZdFNy4r2Ae0oUwmk3boZBJnboQSgp/iyq1+hTtxK36E/+C0DWhbD1w4c8693LnHizlTYNtfxtLyyuraemGjuLm1vbNr7u03VZRIQhsk4pFse1hRzgRtAANO27GkOPQ4bXnD67HfeqBSsUjcwSimboj7ggWMYNBSzzzsAuM+Te+zMlROfx92pWeW7Ko9gbVInJyUUI56z/zu+hFJQiqAcKxUx7FjcFMsgRFOs2I3UTTGZIj7tKOpwCFVbjo5IbNOtOJbQSR1CbAm6t+JFIdKjUJPd4YYBmreG4v/eZ0Egks3ZSJOgAoyXRQk3ILIGudh+UxSAnykCSaS6b9aZIAlJqBTm9lCQk+y/gCy4gQ6H2c+jUXSPKs6dtW5PS/VrvKkCugIHaMyctAFqqEbVEcNRNAjekYv6NV4Mt6Md+Nj2rpk5DMHaAbG5w+pIqBj</latexit> ˜ q(t) <latexit sha1_base64="eA0zx410i6ePjnM+0DH9JiPM3K4=">AAACDnicbVDLSsNAFJ3UV62vqks3g0Wom5KIoMuiG5cV7APaUCaTSTt0JokzN0IJ+QdXbvUr3Ilbf8GP8B+cplnY1gMXDufcy7kcLxZcg21/W6W19Y3NrfJ2ZWd3b/+genjU0VGiKGvTSESq5xHNBA9ZGzgI1osVI9ITrOtNbmd+94kpzaPwAaYxcyUZhTzglICRBgPgwmfpY1aH82G1ZjfsHHiVOAWpoQKtYfVn4Ec0kSwEKojWfceOwU2JAk4FyyqDRLOY0AkZsb6hIZFMu2n+c4bPjOLjIFJmQsC5+vciJVLrqfTMpiQw1sveTPzP6ycQXLspD+MEWEjnQUEiMER4VgD2uWIUxNQQQhU3v2I6JopQMDUtpFDpKT4aQ1bJYfpxlttYJZ2LhmM3nPvLWvOmaKqMTtApqiMHXaEmukMt1EYUxegFvaI369l6tz6sz/lqySpujtECrK9fme2bKQ==</latexit> <latexit sha1_base64="eA0zx410i6ePjnM+0DH9JiPM3K4=">AAACDnicbVDLSsNAFJ3UV62vqks3g0Wom5KIoMuiG5cV7APaUCaTSTt0JokzN0IJ+QdXbvUr3Ilbf8GP8B+cplnY1gMXDufcy7kcLxZcg21/W6W19Y3NrfJ2ZWd3b/+genjU0VGiKGvTSESq5xHNBA9ZGzgI1osVI9ITrOtNbmd+94kpzaPwAaYxcyUZhTzglICRBgPgwmfpY1aH82G1ZjfsHHiVOAWpoQKtYfVn4Ec0kSwEKojWfceOwU2JAk4FyyqDRLOY0AkZsb6hIZFMu2n+c4bPjOLjIFJmQsC5+vciJVLrqfTMpiQw1sveTPzP6ycQXLspD+MEWEjnQUEiMER4VgD2uWIUxNQQQhU3v2I6JopQMDUtpFDpKT4aQ1bJYfpxlttYJZ2LhmM3nPvLWvOmaKqMTtApqiMHXaEmukMt1EYUxegFvaI369l6tz6sz/lqySpujtECrK9fme2bKQ==</latexit> <latexit sha1_base64="eA0zx410i6ePjnM+0DH9JiPM3K4=">AAACDnicbVDLSsNAFJ3UV62vqks3g0Wom5KIoMuiG5cV7APaUCaTSTt0JokzN0IJ+QdXbvUr3Ilbf8GP8B+cplnY1gMXDufcy7kcLxZcg21/W6W19Y3NrfJ2ZWd3b/+genjU0VGiKGvTSESq5xHNBA9ZGzgI1osVI9ITrOtNbmd+94kpzaPwAaYxcyUZhTzglICRBgPgwmfpY1aH82G1ZjfsHHiVOAWpoQKtYfVn4Ec0kSwEKojWfceOwU2JAk4FyyqDRLOY0AkZsb6hIZFMu2n+c4bPjOLjIFJmQsC5+vciJVLrqfTMpiQw1sveTPzP6ycQXLspD+MEWEjnQUEiMER4VgD2uWIUxNQQQhU3v2I6JopQMDUtpFDpKT4aQ1bJYfpxlttYJZ2LhmM3nPvLWvOmaKqMTtApqiMHXaEmukMt1EYUxegFvaI369l6tz6sz/lqySpujtECrK9fme2bKQ==</latexit> <latexit sha1_base64="eA0zx410i6ePjnM+0DH9JiPM3K4=">AAACDnicbVDLSsNAFJ3UV62vqks3g0Wom5KIoMuiG5cV7APaUCaTSTt0JokzN0IJ+QdXbvUr3Ilbf8GP8B+cplnY1gMXDufcy7kcLxZcg21/W6W19Y3NrfJ2ZWd3b/+genjU0VGiKGvTSESq5xHNBA9ZGzgI1osVI9ITrOtNbmd+94kpzaPwAaYxcyUZhTzglICRBgPgwmfpY1aH82G1ZjfsHHiVOAWpoQKtYfVn4Ec0kSwEKojWfceOwU2JAk4FyyqDRLOY0AkZsb6hIZFMu2n+c4bPjOLjIFJmQsC5+vciJVLrqfTMpiQw1sveTPzP6ycQXLspD+MEWEjnQUEiMER4VgD2uWIUxNQQQhU3v2I6JopQMDUtpFDpKT4aQ1bJYfpxlttYJZ2LhmM3nPvLWvOmaKqMTtApqiMHXaEmukMt1EYUxegFvaI369l6tz6sz/lqySpujtECrK9fme2bKQ==</latexit> Figure 2.3: Model fits of PSA measurement data for metastatic castration- resistant prostate cancer — Prostate-specific antigen (PSA) measurement data of four representative patients from three randomised clinical trials with metastatic castration- resistant prostate cancer. (A) A single patient data shown with model fit (eqn. 2.4) of growing resistant population (green), decaying sensitive population (red), where the sum represents an overall initial good response and eventual relapse (q(t), dashed line). Left column: treatment with prednisone only [75]; middle column: treatment with mitoxantrone and prednisone [112]; right column: treatment with docetaxel and prednisone [89]. PSA data is normalized by q(t = 0) (black dots). The data is fit using the exponential model (eqn. 2.1; blue curve) and the replicator model (eqns 2.2, 2.3). Each patient is fit reasonably well with both models. Data was fit by parameter sweep of cost (α−β, eqn. 2.15, 2.16), initial fractional resistance f and selection pressure w. Parameters used are α−β = [0.02, 0.07, 0.20, 0.00, 0.20, 0.00, 0.20, 0.03, 0.20, 0.01, 0.15, 0.19]; f = [0.20, 0.02, 0.09, 0.05, 0.03, 0.07, 0.10, 0.10, 0.09, 0.10, 0.03, 0.03];w = [0.10, 0.20, 0.20, 0.30, 0.20, 0.10, 0.30, 0.15, 0.10, 0.15, 0.20, 0.20] for B - M respectively. 26 During Therapy (w 2 < w 3 < w 1 ) Fitness Unaltered fitness of resistant subpopulation remains constant before, during therapy Low Fitness High Fitness Chemotherapy results in fitness gap between healthy, sensitive cells A B C w = 1 w = 0 w 1 = (1 + c)w w 2 = (1 – c)w w 3 = w Baseline value of selection pressure in unperturbed tumor w Strong selection No selection w 1 - Altered selection pressure on healthy cell population w 2 - Altered selection pressure on sensitive cancer cell population w 3 - Unaltered selection pressure on resistant cell population Selection difference due to therapy: w 1 - w 2 = 2c Fitness advantage of driver mutation (prisoner’s dilemma) No Therapy (w 1 = w 2 = w 3 ) Fitness Sensitive Cancer Cell Resistant Cancer Cell Healthy Cell Fitness cost of resistance Low Fitness High Fitness Figure 2.4: Fitness landscape before and during therapy — (A) Before therapy, a driver mutation leads to a fitness advantage of the cancer cell (red) and a subsequent resistant-conferring mutation comes at a fitness cost (green). Parameters determined by the prisoner’s dilemma payoff matrix reflect these relative fitness differences. (B) Rather than change the elements of the game payoff matrix directly, the selection parameter is manipu- lated for healthy (increased) and cancer (decreased) subpopulations such thatw 1 −w 2 = 2c. (C) Changes in w i results in a new relative fitness of each subpopulation during therapy. The fitness of the resistant population is unaffected by therapy’s selective pressure, but the healthy population is given an advantage over the chemo-sensitive population. 27 A B C D S H R S H R x R x H x S Trilinear Coordinates Schematic Adaptive Therapy Schematic Therapy ON Therapy OFF S H R S H R No therapy: c = 0 With therapy: c = 0.6 Final Version ˙ x R =0 ˙ x S =0 ˙ x H =0 ˙ x R =0 ˙ x S =0 ˙ x H =0 3 2 1 1 2 3 4 5 6 (Therapy ON) (therapy ON) (therapy OFF) ˙ x R =0 ˙ x R =0 Pre-therapy Initial condition Continuous therapy 2 1 3 Unstable equilibrium Stable equilibrium Relative velocity 0 1 Figure 2.5: Dynamic phase portraits before and during chemotherapy — (A) Tri- linear coordinate phase space representation; (B) Schematic of proposed adaptive therapy concept using the resistant nullclines to determine therapy “on" and “off" times in order to trap the tumor in the controllable region 2, and reach approximate cycle that repeats back on itself in red. The continuous therapy is also plotted in dashed blue, for comparison. Two nullclines divide the triangle into 3 regions; region 1: ˙ x R < 0 for both therapy on and off; region 2: ˙ x R > 0 for therapy off and ˙ x R < 0 for therapy on; region 3: ˙ x R > 0 for both therapy on and off. (C) Before chemotherapy, the healthy (H), sensitive (S), and resistant (R) pop- ulations compete on a dynamical fitness landscape, with several solution trajectories shown (black) and the instantaneous relative velocity indicated by background color gradient (red to blue). All internal trajectories lead to tumor growth and eventual saturation of the sensi- tive population (bottom left corner). Each population nullcline (line of zero growth: ˙ x i = 0) is plotted: healthy (dashed blue), sensitive (dashed red), and resistant (dashed green). The nullclines divide the triangle into 3 regions. Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; Region 2: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0; Region 3: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0; (D) Chemotherapy alters the selection pressure to the disadvantage of chemo-sensitive cancer population and advantage of the healthy population (shown for c = 0.6,α = 0.020, β = 0.018, w = 0.1). In this case, the nullclines divide the triangle into 6 regions; Region 1: ˙ x H > 0 ˙ x S > 0 ˙ x R < 0; Region 2: ˙ x H > 0 ˙ x S < 0 ˙ x R < 0; Region 3: ˙ x H > 0 ˙ x S < 0 ˙ x R > 0;Region 4: ˙ x H < 0 ˙ x S < 0 ˙ x R > 0; Region 5: ˙ x H < 0 ˙ x S > 0 ˙ x R > 0; Region 6: ˙ x H < 0 ˙ x S > 0 ˙ x R < 0; Solution trajectories (black) show initial trajectory toward healthy saturation (triangle top) but eventual relapse toward resistant population (bottom right of triangle) upon passing the resistant nullcline. 28 Figure 4: Time to relapse and progression free survival POST REVIEW 0 100 200 300 400 500 600 700 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time to relapse progression free survival (PFS) Therapy 1 (c = 0.5) Therapy 2 (c = 0.6) Therapy 3 (c = 0.7) Adaptive (c = 0.5) Adaptive, shorter t OFF Adaptive, longer t ON Therapy 1 Therapy 2 Therapy 3 Chemo-sensitive Chemo-resistant t = 0 100 200 300 400 500 600 700 Adaptive Therapy A B C D E Figure 2.6: The effect of dose on tumor relapse and progression free survival under continuous and adaptive therapy — (A) PSA data is from a single patient (see figure 2.3B) under continuous treatment (Mitoxantrone and Prednisone) is replotted (black dots) along with the best model fit (eqns. 2.2, 2.3) in blue. Using identical patient parameters, continuous treatment of two “new" drugs with higher effectiveness (i.e. increased effective dose, c) is shown in red and yellow. Time to relapse significantly increases with increasing dose while the progression free survival shows marginal, but decreasing, difference. An adaptive therapy (see figure 2.5B) is also simulated (solid black line), showing an increased controloverthetumor; Similarcontrolparadigmsweresimulatedwitherrorinadministration times: slightly shorter time on therapy (1 day; green) and slightly longer (1 day; purple). (B) The same four therapies are shown in a fish plot. Continous therapies show relapse to initial tumor size is dominated by chemo-resistant population (green). The adaptive therapy successfully suppressed the growth of the resistant population (bottom). Parameters used: α−β = 0.02; w = 0.10; f = 0.2 29 Chapter 3 Nonlinear adaptive control of competitive release and chemotherapeutic resistance We use a three-component replicator system with healthy cells, sensitive cells, and resis- tant cells, with a prisoner’s dilemma payoff matrix from evolutionary game theory to model and control the nonlinear dynamical system governing the ecological mechanism of com- petitive release by which tumors develop chemotherapeutic resistance. The control method we describe is based on nonlinear trajectory design and energy transfer methods first intro- duced in the orbital mechanics literature for Hamiltonian systems. For continuous ther- apy, the basin boundaries of attraction associated with the chemo-sensitive population and the chemo-resistant population for increasing values of chemo-concentrations have an inter- twined spiral structure with extreme sensitivity to changes in chemo-concentration level as well as sensitivity with respect to resistant mutations. For time-dependent therapies, we introduce an orbit transfer method to construct continuous families of periodic (closed) orbits by switching the chemo-dose at carefully chosen times and appropriate levels to design schedules that are superior to both maximum tolerated dose (MTD) schedules and low-dose metronomic (LDM) schedules, both of which ultimately lead to fixation of sensitive cells or resistant cells. By keeping the three subpopulations of cells in competition with each other indefinitely, we avoid fixation of the cancer cell population and re-growth of a resistant tumor. The method can be viewed as a way to dynamically shape the average population fitness landscape of a tumor to steer the chemotherapeutic response curve. We show that the 30 method is remarkably insensitive to initial conditions and small changes in chemo-dosages, an important criterion for turning the method into an actionable strategy. 3.1 Introduction It is widely appreciated that the development of chemotherapeutic resistance is the pri- mary reason for recurrance of cancer in patients undergoing treatment, and remains one of the primary challenges in the field of oncology [58, 1, 17]. As a tumor grows, and even as tumor cells spread throughout the system and metastasis ensues, standard pre- scheduled chemotherapeutic protocols such as maximum tolerated dose (MTD) and low-dose metronomic schedules (LDM) often show early success as the tumor regresses temporarily. Schematics for several kinds of chemotherapeutic schedules are shown in Figure 3.1. A typ- ical chemotherapeutic cycle might involve one strong dose every three weeks, or a dose for 5 consecutive days, followed by a 28 day rest period [88]. After months of fixed periodic cycles, the cancer often recurs and the tumor begins to regrow. Because of the genetic and cellular heterogeneity of a typical tumor [116], instead of killing all of the cancer cells and thereby eliminating the tumor, the chemotherapeutic regimen actually selects for a resistant pheno- type [105, 39, 44, 27]. The diversity of cells within a tumor effectively protects the tumor from single-line or pre-scheduled chemotherapeutic assaults by allowing for elimination of the chemo-sensitive population in order to accomplish the subsequent release of the chemo- resistant population. By reducing the relative fitness of the sensitive cells, chemotherapy acts as the primary mechanism of natural selection that selects specifically against rapidly dividing cells [68]. The characterization of a typical tumor as an adaptive landscape made up of competing cellsofvaryingdegreesoffitness, whichdeterminegrowthratesofthevarioussubpopulations is a more accurate characterization of a tumor and suggests an ecological or evolutionary approach [74, 6, 10, 9, 31, 107, 12, 118, 122, 119]. If one had access to time-resolved infor- mation [33] on the relative balance and growth rates of the subpopulations of cells making up the tumor, then one could use chemotherapy as a control device (accuator) to keep the 31 subpopulations in balance, competing with each other indefinitely, without any one of the cancerous subpopulations dominating the landscape [105]. Chemotherapy would then be regarded more as a maintenance mechanism than a cure [38, 54] and one would be imposing no more selection than is necessary as has been advocated by Read et al. and others [94]. We introduce an evolutionary game theory model of chemotherapeutic resistance [58] along with a method of adaptive control to design adventageous chemotherapeutic schedules that are able to overcome, or at least manage resistance. The mechanism of resistance that we model is based on the ecological notion of competitive release [44, 104, 74] of the resistant cellpopulationwhenthesensitivecellpopulationisreducedbelowacertainthreshold. Above the threshold, the sensitive cells are able to outcompete (on average) the resistant cells due to the inherent cost-of-resistance [5] which tilts the fitness landscape of the system in favor of the sensitive cell population allowing the tumor to grow. Under sufficient chemotherapeutic pressure, the sensitive cell population is reduced enough to allow the resistant population to begin to flourish and eventually regrow the tumor in a form that is much harder to treat. A quantitative understanding of this phenomena is necessary in order to develop chemotherapeutic strategies (i.e. adaptive therapies) to combat it, a point of view adopted and developed in [38, 40, 107, 17, 10, 9]. In this paper, we frame the problem as one in nonlinear dynamical systems theory and use trajectory transfer methods developed relatively recently in the orbital mechanics literature [60, 43, 47, 99, 93] to design chemotherapeutic schedules that have the potential to outperform more standard approaches [84]. The mathematical model we use is based on a three-component replicator dynamical sys- tem with a frequency dependent fitness function based on a prisoner’s dilemma (PD) payoff matrix [85, 52]. Cell interactions occur between three cell types that form our ecosystem: healthy cells (H); chemo-sensitive cells (S); and chemo-resistant cells (R). The healthy cell subpopulation can be thought of as non-neoplastic cells that have a lower fitness than the chemo-sensitive subpopulation. The independent variables represent relative frequencies in the total population, and the tumor alone would be made up of S +R cancer cell subpop- ulations. In a PD scenario, the healthy cells act as cooperators while the cancer cells act 32 (a) (b) (c) (d) Figure 3.1: Schematic dose schedules.(a) Maximum tolerated dose (MTD) schedule; (b) Low dose metronomic (LDM) schedule; (c) Adaptive schedule; (d) General time-dependent schedule accuated by piecewise constant dose concentrations. as defectors [7]. Unchecked, the defectors saturate the population as the average fitness decreases to a sub-optimal outcome. The goal of chemotherapeutics in this framework is to coax the defectors to cooperate, leading to a higher fitness Nash equilibrium [85, 52]. Growth rates are determined by cell fitness functionals which, in turn, depend on subpopu- lation sizes, i.e. they are frequency dependent. There is a time-dependent controller in our model which determines chemotherapeutic dose schedule and concentrations. Although the model system is nonlinear and the control parameter enters as a coefficient in most all of the terms of the three-component cubic nonlinear system (making classic control schemes [18] 33 non-applicable), we use nonlinear trajectory design techniques introduced and developed in an orbital mechanics context [60, 43, 47, 99, 93]. In that body of literature, time-dependent controllers are used to design orbit transfers in a Hamiltonian mechanics setting, piecing together partial orbits at different energy levels and switching energies at carefully chosen times, much like classic Homann transfers for satellite control [92]. While the replicator sys- tem we describe is not a Hamiltonian system in which orbits transfer from one energy level to another, the time-dependent chemotherapeutic parameter, C(t), can be used to design adventageous orbits in the replicator dynamics tri-linear phase space in a similar manner where piecewise constant dose concentrations are used,C i =const. for carefully chosen time intervals t i < t < t i+1 , (i = 0,...,n) with switching times t i chosen in such a way as to produce a periodic (closed), continuous, piecewise differentiable orbit that stays trapped in a desirable region of the phase space. Orbits designed this way are shown to maintain a higher average level of fitness for the full population and avoid tumor recurrance. The existence of such orbits for appropriately chosen chemotherapeutic schedules in our model system suggests the possibility that similar orbits may also exist in a more complex tumor environment with carefully designed adaptive schedules [134, 40]. The general technique of designing an orbit with a time-dependent controller should also work in other contexts such as microbial drug resistance [81] or pest management [44, 19], although with perhaps different numbers of competing subpopulations and different payoff matrices that determine the fitness landscapes. 3.2 A three-component replicator system The model we develop is based on a three-component replicator nonlinear dynamical system governing three competing subpopulations of cells: ~ x = (x 1 ,x 2 ,x 3 ) T = (x H ,x S ,x R ) T , where x 1 represents the proportion of healthy cells (H), x 2 represents the proportion of sensitivecells(S),andx 3 representstheproportionofresistantcells(R),withx 1 +x 2 +x 3 = 1. In this context, the model makes the well-mixed assumption on the cell population [85, 52], i.e. no spatial dependence is modeled in order to highlight the method in as clean a setting 34 as possible. For discussions and comparisons of well-mixed models vs. those with spatial dependence, see [132]. The equations which describe the subpopulation interactions are: ˙ x 1 = (f 1 −hfi)x 1 , (3.1) ˙ x 2 = (f 2 −hfi)x 2 , (3.2) ˙ x 3 = (f 3 −hfi)x 3 . (3.3) with f i representing the fitness of each of the subpopulations (i = 1, 2, 3) as their rela- tive populations change, andhfi representing the average fitness of the entire population. The exponential growth-decay rates of each of the subpopulations are then determined by (f i −hfi) which dictates whether the subpopulation fitness is above or below the average population fitness, hence whether the subpopulation decays or grows. The fitness of each of the three subpopulations are defined by the linear functionals: f 1 = 1−w 1 +w 1 (A~ x) 1 , (3.4) f 2 = 1−w 2 +w 2 (A~ x) 2 , (3.5) f 3 = 1−w 3 +w 3 (A~ x) 3 , (3.6) where 0≤ w i (t)≤ 1, (i = 1, 2, 3) are time-dependent selection parameters (that serve as our controllers) we use to shape the fitness landscape of the system. A is the payoff matrix associated with the cell-cell interactions, which we describe in §3.3. The time dependence in our model enters through a chemo-concentration parameter C(t): w 1 = w 0 , (3.7) w 2 = w 0 (1−C(t)), (3.8) w 3 = w 0 , (3.9) 35 wherew 0 scales time (we typically takew 0 = 1). Note that the chemotherapy parameter acts linearly on the sensitive cell population lowering its fitness, although the three populations are coupled nonlinearly through [3.1]-[3.3]. The average population fitness is defined by the quadratic functional: hfi =x 1 f 1 +x 2 f 2 +x 3 f 3 (3.10) in the usual way. It is straightforward to see (for fixed values of the chemo-concentration parameterC) that the fixed points of the system [3.1]-[3.3] are of three basic types: (i) There are the three fixed points at each of the corners of the triangular phase space diagram shown in Figure 3.2, when two of the subpopulation values are zero, and the third is saturated: (x 1 ,x 2 ,x 3 ) = (x H ,x S ,x R ) = (1, 0, 0); (0, 1, 0); (0, 0, 1); (ii) There are three possible fixed points on the triangle sides, which correspond to one of the subpopulation values equaling zero, with the other two having fitness values equal to the population average: x 1 = 0,f 2 = hfi,f 3 =hfi;x 2 = 0,f 1 =hfi,f 3 =hfi;x 3 = 0,f 1 =hfi,f 2 =hfi; (iii) There is the balanced fitness state, when none of the subpopulation values is zero, but each of the subpopulation fitness values equals the population average:f 1 =hfi ;f 2 =hfi ;f 3 =hfi. Which of these fixed points lies on or inside the triangular phase space, and their stability properties, depend in detail on the parametric values, which we describe in §3.4. By using a prisoner’s dilemma payoff matrix (see section 3.3) we ensure: (i) Gompertzian growthofthecancercells[83]; (ii)areductioninoverallfitnessofthepopulationasthetumor grows; and (iii) a fitness cost associated with resistance. We first study the details of the nonlinear dynamics associated with eqns [3.1]-[3.3] for constant values of the chemotherapy parameter 0≤C≤ 1 to demonstrate the mechanism of competitive release for threshold val- uesC≥ 1/3. Then we investigate piecewise constant time-dependent functionsC(t) to show how to avoid the evolution of resistance of the tumor. Figure 3.1 shows several examples of the chemotherapeutic schedules we consider. These include maximum tolerated dose (MTD) schedules (Figure 3.1(a)), low-dose metronomic schedules (LDM) (Figure 3.1(b)), adaptive 36 (a) (b) (c) Figure 3.2: The competitive release mechanism. (a) The three-component phase space associated with competing populations of (H,S,R) cells. (b) With no therapy, theS corner of the triangle is a globally attracting fixed point, while the H and R corners are unstable. All initial conditions lead to a saturated tumor. Filled circles are stable, unfilled circles are unstable. (c) For continuous chemotherapy above a threshold level, the R corner of the triangle is a globally attracting fixed point, while the H and S corners are unstable. All initial conditions (except those lying on the separatrix connecting the interior balanced fitness state to the S−R side) lead to a resistant tumor. schedules (Figure 3.1(c)), and more general time-dependent schedules (Figure 3.1(d)) which we break-up into piecewise constant doses as done in forming the Riemann sum approxima- tion to an area under a curve. In all cases, we compare outcomes of the different schedules holding the total dosage, D, fixed over time period τ, D = Z τ 0 C(t)dt =const. (3.11) 3.3 The prisoner’s dilemma as a cancer model We first describe the standard version of the PD payoff matrix [85] in a 2× 2 setting to make clear why the PD evolutionary game gives rise to a Gompertzian (i.e. ‘S-shaped’) growth curve when two sub-populations compete. To be specific, we first describe what happens when healthy cells compete with cancer cells using the 2× 2 payoff matrix: 37 A PD = a b c d = 3 0 5 1 , (3.12) withc>a>d>b. The first row and column correspond to the payoffs associated with the cooperator (C) in the PD evolutionary game, and the second row and column correspond to the payoffs associated with the defector (D). In the simplest tumor growth paradigm in which a population of healthy cells competes with a population of cancer cells, the healthy cells are the cooperators, while the cancer cells are the defectors. In any mixed population ~ x = (x C ,x D ) T , 0≤x C ≤ 1; 0≤x D ≤ 1; x C +x D = 1, the fitness functions, ~ f = (f C ,f D ) T , associated with the two subpopulations are: ~ f = A PD ~ x, (3.13) which in component form yields: f C = (A PD ~ x) 1 = 3·x C + 0·x D , (3.14) f D = (A PD ~ x) 2 = 5·x C + 1·x D , (3.15) while the average fitness of the total population is given by the quadratic form: hfi = ~ x T A PD ~ x = 3x 2 C + 5x C x D +x 2 D ≥ 1. (3.16) Note that the average fitness of the healthy state (x C ,x D ) = (1, 0) is given byhfi| (x C =1) = 3, whilethatofthecancerousstate (x C ,x D ) = (0, 1)isgivenbyhfi| (x D =1) = 1, whichminimizes 38 the average fitness. Tumor growth is then modeled as a 2× 2 evolutionary game governed by the replicator dynamical system: ˙ x C = (f C −hfi)x C , (3.17) ˙ x D = (f D −hfi)x D . (3.18) It is straightforward to show: ˙ x D = [(c−a)− (d−b)]x D (1−x D ) 1 1− d−b c−a −x D , with fixed points at x D = 0, 1, (c−a) (c−a)−(d−b) . From this, we can conclude that for any initial condition containing at least one cancer cell: 0<x D (0)≤ 1, we have: (i) x D → 1, x C → 0 as t→∞ (ii)hfi→ 1 as t→∞. Condition (i) guarantees that the cancer cell population will saturate, while condition (ii) guarantees that the saturated state is sub-optimal, sincehfi| (x D =1) <hfi| (x C =1) . For these two reasons, the prisoner’s dilemma evolutionary game serves as a simple model for tumor growth both in finite population models, as well as replicator system (infinite population) models [118, 119, 122, 114]. For our purposes, we now generalize to a three-component system where the fitness functions [3.4]-[3.6] are defined via a payoff matrix A which is a 3× 3 matrix defining the evolutionary game played by the cell population. For this, we take every 2× 2 sub-matrix to be a PD game, i.e. we take A to be the 3× 3 PD matrix: A = a b c d e f g h i = 3 1.5 1.5 4 2 2.8 3.9 −2 2.2 , (3.19) 39 with the PD inequalities [85]: g>a>i>c, (3.20) d>a>e>b, (3.21) f >i>e>h. (3.22) The numerical values in [3.19] are chosen for convenience and satisfy the constraints [3.20]- [3.22]. In each cell-cell interaction, the healthy cellsx 1 (healthy H) are cooperators, and the two-species of cancer cells, x 2 (sensitive S), and x 3 (resistant R), are the defectors. In any interaction between a chemo-sensitive cell (S) and a chemo-resistant cell (R), the sensitive cell is the defector, while the resistant cell is the cooperator. The payoff matrix [3.19] guarantees that for any interaction between two cells, the system retains features (i),(ii) detailed previously (i.e. tumor growth leading to sub-optimal fitness). In addition, the payoff matrix also imposes a cost to resistance if we add the extra constraint d > g which guarantees: f S = dx H +ex S +fx R (3.23) > f R =gx H +hx S +ix R It is worth pointing out that shaping the fitness landscape by adjusting the selection parameters (w 1 ,w 2 ,w 3 ) is equivalent to choosing f i = ( ˜ A~ x) i , (3.24) where: ˜ A = 1−w 1 +aw 1 1−w 1 +bw 1 1−w 1 +cw 1 1−w 2 +dw 2 1−w 2 +ew 2 1−w 2 +fw 2 1−w 3 +gw 3 1−w 3 +hw 3 1−w 3 +iw 3 . 40 To see this, use the fact that (1−w i )≡ (1−w i )(x 1 +x 2 +x 3 ) in the fitness equations [3.4]-[3.6]. Figure 3.3: Fixed point locations. Tracking the location of the fixed points as a function of chemo-concentration parameter C. The fractional values indicate that the values are analytically obtained. 3.4 Results 3.4.1 Continuous therapy For the purposes of understanding how to implement time-dependent therapies in our model, we first describe the phase space dynamics for constant values 0 ≤ C ≤ 1, so 41 (a) (b) (c) (d) (e) (f) Figure 3.4: Basins of attraction and nullclines. Panel showing the separatrices and nullclines (dx H /dt = 0, dx S /dt = 0, dx R /dt = 0) through the balanced fitness interior fixed point that determines the basin boundaries of attraction of the S state and the R states. The interior fixed point is the one in which each of the subpopulation fitness levels exactly matches the average fitness of the entire population. Filled circles are stable, unfilled circles are unstable. total dose delivered over time period t is D = C·t. Figure 3.2 shows the mechanism of chemotherapeutic resistance via competitive release determined by the system [3.1]-[3.3] (for constantC) as depicted in the triangular phase space diagram shown in Fig 3.2(a). With no therapy (Fig.3.2(b), C = 0), the sensitive corner S is globally attracting, while the H and R corners are unstable. When the continuous therapy parameter C = 0.6 (Figure 3.2(c)), the resistant corner is globally attracting, while the H and S corners are unstable. Filled corners are stable, unfilled corners are unstable. This is the basic mechanism of competitive release induced by sufficiently strong chemotherapeutic dose [105]. 42 Figure 3.5: Basin of attraction areas. Areas of basin of attraction of S fixed point and R fixed point as a function of chemo-concentration C. Intermediate values of C reveal a much more complex picture. We show in Fig 3.3 the location of the fixed points as a function ofC. For values 79/228 = 0.34649...<C < 0.7, the balanced fixed point state is an interior fixed point, which forms the central spiral associated with the basin boundaries between the asymptotically stable S corner and R corner, shown in the panel in Figure 3.4. Also shown in the figure are the nullclines defined by the curves dx H /dt = 0,dx S /dt = 0,dx R /dt = 0. Opposite sides of each of the nullclines mark a change in whether the particular subpopulation decays or grows. The mixed population state exists between values 1/3<C < 1/2 where the basin of attraction sizes (areas), shown in Figure 3.5, change sensitively (asC varies), one at the expense of the other, in an intertwined spiral 43 (a) (b) (c) (d) Figure 3.6: Fitness landscapes as a function of time. With no therapy (C = 0), the fitness curves are contunuously decreasing functions as the tumor saturates with the sensitive cell population in a sigmoidal shaped growth curve. With continuous therapy (C = 0.6), the fitnesses initially increase indicating tumor regression, but eventually decrease. The healthy subpopulation initially increases before the resistant population eventually saturates the tumor. Initial condition (0.9,0.09,0.01) with w 0 = 0.1 structure centered at the balanced fixed point state. The steep transition curves between the two states occuring for small changes in the chemo-concentration parameter highlights the sensitivity of the system to chemotherapeutic dosing levels. Also worth highlighting in Fig. 3.4(b),(c) is the sensitivity to resistant mutations of the final asymptotic state of the system, 44 even if no pre-existing mutations exist in the population. As the tumor grows, assuming no pre-existing resistant mutations, the dynamics would traverse down the left side of the triangle from the H corner to the S corner. As it does, the very thin passageway (on the H-S side of the triangle in Figure 3.4(b),(c)) associated with the basin of attraction of the resistant corner indicates that one single resistant mutation would push the trajectory off the side into the basin of attraction associated with the R corner instead of converging to the S corner. The fitness landscapes are shown in Figure 3.6 both for C = 0 and C = 0.6. With no therapy, the fitness curves are monotonically decreasing functions as the sensitive population saturates the tumor in a sigmoidal shaped [83, 118] growth curve. With continuous therapy above threshold, the fitness curves initially increase (i.e. tumor regression), but eventually decrease monotonically as the tumor relapses. The healthy cell population initially increases, but eventually the resistant population saturates the tumor. 3.4.2 Time-dependent therapy Figure 3.7 shows the key idea behind the method we use to design trajectories for eqns [3.1]-[3.3] for time-dependent chemotherapeutic schedulesC(t). With no chemotherapy,C = 0, since the sensitive corner S is a globally asymptotically attracting fixed point whose basin of attraction is the full region, all trajectories that start inside the triangle eventually get trapped in the left S corner. By contrast, with chemotherapeutic levels at C = 0.7, competitive release acts to create a basin of attraction for the right resistant corner R for all initial conditions inside the triangle. Using these families of solution trajectories, we overlay the solution curves in Figure 3.7(c) to show the underlying curvilinear grid that spans the full tri-linear phase plane. By switching between the two values C = 0 and C = 0.7 at times when two curves intersect, it is possible to transition from a trajectory associated with the C = 0 family to one associated with the C = 0.7 family. This creates multiple possibilities for designing complex orbits using piecewise constant values ofC, with switching at appropriately chosen times. 45 One such trajectory is shown in Figure 3.8(a) achieved by switching between the two values C = 0 and C = 0.5. The trajectory starts at point A ((H,S,R) = (0.5, 0.4, 0.1)) using a trajectory with C = 0.5. When the trajectory reaches point B, the chemotherapy is switched off C = 0, and the system then travels back down to point A in a closed loop. With switching times labeled T 1 and T 2 (which are computed by monitoring the relative balances of the three subpopulations) , this closed orbit can be maintained indefinitely. By contrast, Figure 3.8(b) shows the corresponding MTD and LDM trajectories with the same initial condition (point A) using the same total dosages over the same time-periods. The MTD trajectory eventually gets trapped in the left S corner, as does the LDM trajectory. Figure 3.8(c) shows the schedules for all three cases. In all cases, we first design the adaptive schedule, then we create the MTD and LDM schedules using the same total dosage D. Neither the MTD nor the LDM standard chemo-schedules is able to prevent the system from saturating with a full grown tumor, whereas the adaptive schedule keeps the system trapped indefinitely near the top H corner of the system. Figure 3.8(d) shows the average fitness of the system for the MTD, LDM, and adaptive schedules. It is clear that the adaptive schedule is able to maintain a higher average fitness throughout the full course of chemotherapy. Figure 3.9 shows the sensitivity of the system to the chemo-concentration levels chosen. Here we switch between C = 0 and C = 0.6 (higher average dose than in figure 3.8) to construct the closed loop (figure 3.9(a)), which looks very similar to that in Figure 3.8. But notice that for these values, the LDM schedule creates an orbit that saturates at the R corner, whereas the MTD schedule saturates at the S corner. The actual schedules are shown in Figure 3.9(c). Figure 3.9(d) compares the average population fitness of the three schedules, showing the adaptive schedule maintaines a higher average throughout. Notice that LDM initially achieves a higher average fitness before tumor regression occurs. In Figure 3.10 we show the result of toggeling between values C = 0.3 and C = 0.6 to maintain the periodic loop (Figure 3.10(a)). For this case, both the MTD and LDM schedules send the trajectory to the R corner (Figure 3.10(b)). Figure 3.10(d) shows the 46 initial benefit of the MTD and LDM schedules in terms of higher initial average fitness, but eventually the adaptive schedule shows its superiority over both. Figure 3.11(a) shows that we can actually design an orbit that starts at an arbitrary point A inside the triangle, and send it to an arbitrary point B. We accomplish this by constructing the incoming and outgoing orbits from point A for two different C values, and those associated with pointB for those same twoC values, showing that they must intersect at some point which we label O. By sending the orbit out from point A to point O, then switching values of C at point O until we arrive at point B, we complete the transfer. One can immediately see the potential richness in the possible design of different orbits one can construct by switching values ofC among two, three, or more values, at appropriately chosen times. Figure 3.11(b) shows the richness of the curvilinear grid that can be created with three values ofC = 0, 0.3, 0.7 and the multitude of possible paths from one point to another inside the triangle if one allows for switching among the three values shown. (a) (b) (c) Figure 3.7: Dynamical trajectories for constant C.(a) With no chemotherapy (C = 0), the sensitive corner S is a globally attracting fixed point. All initial conditions inside the triangle move to S along the sample trajectories shown. (b) Above the chemotherapy threshold C > 0.5, all initial conditions inside the triangle move to the resistant corner R. Shown are sample trajectories forC = 0.7. (c) Overlay of the solution trajectories forC = 0 and C = 0.7 create a curvilinear grid throughout the triangle. By switching between these two values ofC, we construct a global trajectory made up of segments of the two families of trajectories. 47 (a) (b) (c) (d) Figure 3.8: Constructing a closed loop trajectory using C = 0 and C = 0.5, with C avg = 0.366.(a) By using segments of a trajectory for C = 0 and C = 0.5, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories starting from pointA. Both eventually move to theS corner, although initially the MTD trajectory moves toward theH corner (tumor regression) before recurrance. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitnesshfi for the three different chemo-schedules. The adaptive schedule is able to maintain a higher average fitness throughout. 3.4.3 Robustness Tobeclinicallyactionable, attheveryleastitwouldbeimportantthatthemethodnotbe sensititve to specific initial conditions and parameter values. We describe several remarkably robust features of the strategy in this section. We show in Figure 3.12(a) an example of a 48 (a) (b) (c) (d) Figure 3.9: Constructing a closed loop trajectory using C = 0 and C = 0.6, with C avg = 0.375.(a) By using segments of a trajectory for C = 0 and C = 0.6, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories starting from pointA. The MTD trajectory eventually moves to the S corner, while the LDM trajectory moves to the R corner. Both trajectories initially move toward the H corner (tumor regression) before recurrance. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitnesshfi for the three different chemo-schedules. The adaptive schedule is able to maintain a higher average fitness throughout, although LDM initially achieves higher average fitness before declining. continuous family of closed orbits that are easily achievable by using a bang-bang strategy, with blue orbits (no therapy) and red orbits (C = 0.7) forming a continuous family of closed, and nearly closed loops. The associated dose schedules are shown in Figure 3.12(b)-(j), with average doses all very close to 0.4. Because of the continual interlacing of these orbits, and the robustness of their associated schedules, it is clear there is nothing particularly special about any of the values shown, or the average dose - any of the orbits and schedules would work. Interestingly, Figure 3.13 shows that the entire region on the right side of the triangle centered on the dashed curve can be populated by continous families of closed orbits each centered at any point on the line. Figure 3.13(a) depicts the defining feature of the line which requires that two orbits with different C values have a point of tangency lying on it. A necessary and sufficient condition for this is f 1 = f 3 , which defines a line, a condition obtained by requiring that tangent vectors associated with two different values of C point in opposite directions, as shown in Figure 3.13(a). Figure 3.13(b) shows examples of closed orbits obtained at 5 representative points on the dashed line. 49 (a) (b) (c) (d) Figure 3.10: Constructing a closed loop trajectory using C = 0.3 and C = 0.6, with C avg = 0.396.(a) By using segments of a trajectory forC = 0.3 andC = 0.6, switching values at points A and B, we construct a closed periodic orbit. (b) Using the same total dose, we show the MTD and LDM trajectories starting from point A. Both trajectories move toward the R corner after initially moving toward H. (c) The MTD, LDM, and adaptive schedules are depicted. (d) We plot the average fitnesshfi for the three different chemo-schedules. The adaptive schedule is able to maintain a higher average fitness throughout. Both MTD and LDM initially show higher average fitness before declining. 3.5 Discussion The model shows that if the chemo-dose exceeds a threshold value ofC > 0.5, the tumor may regress for a period of time, but eventually regrows to form a resistant tumor (Figure 3.2(c)) as long as there is at least one resistant cell in the population (i.e. a pre-existing 50 (a) (b) Figure 3.11: Constructing orbit transfers. (a) We take two arbitrary points labeled A and B and construct the incoming and outgoing solution trajectories from each, using two valuesC = 0 andC = 0.7. There must be a crossing point, which we label pointO. The two segments AO followed by OB is the two-switch trajectory that takes us from A to B. (b) Shown is the curvilinear grid constructed by families of solution curves for the three values C = 0, C = 0.3, C = 0.7. Any point of intersection can be used as a starting point and an ending point to construct a three-switch orbit. resistant mutation). This process of competitive release is very robust and occurs for all initial distributions of the three subpopulations. For chemo-doses that are not as large, the results are much more sensitive to small changes in C since the phase space diagram generally has two basins of attraction, one associated with a saturated sensitive state, the other with a saturated resistant state, and these basins of attraction are spirally intertwined. Small differences in the balances among the three types of cells comprising our model, and a single resistant mutation can determine whether the long-time dynamics converges to the S state or the R state. The relative areas of the basins of attraction also sharply change with small changes in chemo-dose values between 1/3 < C < 1/2 reflecting the sensitivity in choosing the optimal chemo-dose levels. There is no constant value of C for which the system converges to the H state, as this state is always unstable. By contrast, for a time-dependent chemotherapy parameter C(t), one can design a rich, robust, and potentially endless array of trajectories (i.e. tumor responses) that remain in 51 desirable regions of the triangular phase space indefinitely, depending on how many time- switches and different dose levels one is willing to accept. One can, in principle, steer the system trajectory along any path (in the model). Identifying a given trajectory in our model as well as the chemo-therapeutic schedule that produces it is fairly straightforward and amounts to superposing the response curves associated with different constant therapies, and piecing together a global trajectory made up of sections generated from each of the constant therapy curves. Of course, in a clinical setting this procedure will be much more complex but suggests the possibility of using actual chemo-therapeutic responses with differ- ent chemo-schedules used on different patients as a design template to combine in new ways to predict possible responses to various time-dependent switching strategies. By using one dosing schedule to create a time-limited response, then a different dosing schedule to create a different (time-limited) response, one could think of piecing together finite-time limited responses to create new outcomes, superior to what would have been achieved by sticking to one single chemo-schedule. Oneofthechallengesoftestinganddesigningnewschedulesviaclinicaltrialsisthatshort term gain in average fitness with LDM or MTD does not always result in more long-term increased averaged fitness levels, i.e. recurrance sets in if the schedule does not completely eradicate the tumor. This is clearly shown in our simulations. The possibility of designing complexorbitswithvariouspotentiallyadventageousfeaturesisvirtuallyendlessifoneallows for enough switches among many different levels of the chemo-concentration parameter C. In this manuscript, we have focused on comparisons of average fitness levels to evaluate the quality of the schedule. One could, of course, imagine designing orbits that maximize a (time-averaged) population fitness while minimizing total dose, or perhaps delay recurrance as long as possible while avoiding regions of the phase space (i.e. relative balances of the subpopulation levels) where the total tumor burden becomes unsustainable. The fact that closed loop trajectories can be designed in our model three-component replicator system by chemo-scheduling alone suggests the possibility that similar orbits could potentially be created in an actual tumor environment, microbial environment, or pest management setting 52 with the right accuation (chemo-schedule/antibiotic regimen/pest control schedule), where fixation of the sensitive population or fixation of the resistant population are both avoided indefinitely as the balance is managed with proper adaptive therapy. This, of course, all hingesonourabilitytocarefullyandfrequentlymonitorthetumorenvironment[33]. Clinical trialsthatmakeuseofadaptiveschedulingideasarecurrentlybeingcarriedoutattheMoffitt Cancer Center, and results show great promise [134]. More possibilities exist with the use of combination multi-drug therapies [10, 81, 130], or immunotherapies [95], but these are outside the scope of our model which focuses only on the use of dosing schedules to shape the fitness landscape of a tumor to steer its evolutionary response. 53 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) Figure 3.12: Robust family of nested orbits. (a) One example of a (continuous) family of nested closed orbits using bang-bang strategy toggling between C = 0 andC = 0.7, all with similar total dose (C∼ 0.4). A closed loop can be made arbitrarily small around a center point described in Figure 3.13. (b) Schedule associated with orbits 1 and 4; (c) Schedule associated with orbits 1 and 5; (d) Schedule associated with orbits 1 and 6; (e) Schedule associated with orbits 2 and 4; (f)Schedule associated with orbits 2 and 5; (g) Schedule associated with orbits 2 and 6; (h)Schedule associated with orbits 3 and 4; (i)Schedule associated with orbits 3 and 5; (j) Schedule associated with orbits 3 and 6. 54 (a) (b) Figure3.13: Allclosedorbitfamiliesarecenteredondashedlinedefinedbyf 1 =f 3 . (a) As a family of closed orbits enclose smaller and smaller loops, they intersect at a point of tangency. The point of tangency must lie on the dashed line. (b) Examples of several closed loop orbits centered at different points down the line. 55 Chapter 4 The role of synergy and antagonism in designing multi-drug adaptive chemotherapy schedules Chemotherapeutic resistance via the mechanism of competitive release of resistant tumor cellsubpopulationsisamajorproblemassociatedwithcancertreatmentsandoneofthemain causesoftumorrecurrence. Often, chemoresistanceismitigatedbyusingmultidrugschedules (two or more combination therapies) that can act synergistically, additively, or antagonis- tically on the heterogeneous population of cells as they evolve. In this paper, we develop a three-component evolutionary game theory model to design two-drug adaptive schedules (timing and dose levels associated with C 1 (t) and C 2 (t)) that mitigate chemoresistance and delay tumor recurrence in an evolving collection of tumor cells with two resistant subpopu- lations: R 1 (sensitive to drug 1, resistant to drug 2), and R 2 (sensitive to drug 2, resistant to drug 1). A key parameter, e, takes us from synergistic (e > 0), to additive (e = 0), to antagonistic (e < 0) drug interactions. In addition to the two resistant populations, the model includes a population of chemosensitive cells, S that have higher baseline fitness but are not resistant to either drug. Using the nonlinear replicator dynamical system with a pay- off matrix of Prisoner’s Dilemma (PD) type (enforcing a cost to resistance), we investigate the nonlinear dynamics of the three-component system (S,R 1 ,R 2 ), along with an additional tumor growth model whose growth rate is a function of the fitness landscape of the tumor cell populations. We show that antagonistic drug interactions generally result in slower rates of adaptation of the resistant cells than synergistic ones, making them more effective in com- bating the evolution of resistance. We then design closed loops in the three-component phase 56 space by shaping the fitness landscape of the cell populations (i.e. altering the evolutionary stable states of the game) using appropriately designed time-dependent schedules (adaptive therapy), altering the dosages and timing of the two drugs using information gleaned from constant dosing schedules. We show that the bifurcations associated with the evolutionary stable states are transcritical, and we detail a typical antagonistic bifurcation that takes place between the sensitive cell populationS and theR 1 population, and a synergistic bifur- cation that takes place between the sensitive cell population S and the R 2 population for fixed values of C 1 and C 2 . These bifurcations help us further understand why antagonistic interactions are more effective at controlling competitive release of the resistant population than synergistic interactions in the context of an evolving tumor. 4.1 Introduction We study a mathematical model to explore the role of synergisitic vs. antagonistic multidrug interactions on an evolving population of cancer cells in a tumor. Our model builds on the (single drug) adaptive therapy model developed in [78] and is based on a replicator dynamical system of three (well-mixed) populations of cells: (i) sensitive cancer cells, S, that are sensitive to both drug 1 and drug 2, (ii) a resistant population, R 1 , that is sensitive to drug 1 but resistant to drug 2, and (iii) resistant population, R 2 , that is sensitive to drug 2 but resistant to drug 1. The replicator dynamical system governing the relative frequencies of (S,R 1 ,R 2 ) makes use of a 3×3 payoff matrixA of Prisoner’s Dilemma type [52, 85, 7, 119]. We introduce chemotherapeutic schedules using time-dependent dosing functions C 1 (t), C 2 (t) that alter the fitness levels of the sensitive cell population and the two resistant populations independently allowing us to shape the fitness landscape of the population of cells adaptively as the tumor evolves [40, 39, 107, 32, 37, 134, 34, 35]. Our general goal in this context is to design schedules that delay tumor recurrence (re-growth) due to competitive release of the resistant cell population [42, 17, 63, 58] and to quantify the role of synergistic and antagonistic drug interactions in this process. 57 To quantify the synergistic vs. antagonistic effects of the two toxins on the evolving populations, we use a parameter e which we alter. Whene> 0, the drugs interact synergis- tically, when e = 0, the drugs are additive, and when e < 0, they interact antagonistically. Level curves of the fitness profiles as a function of C 1 and C 2 are shown in figure 4.1 which can be compared with figure 1 of [50]. Details of how we define our fitness functions in the context of our replicator model will be described in §II. In the neutral case, additive drug interactions (figure 4.1(b)) accomplish exactly the kill rate that the sum of each of the two would accomplish acting independently. By contrast, synergistic drug interactions (figure 4.1(a)) bend the curves inward, indicating that the growth rates (fitness) are lowered more than the two dosages would accomplish independently, while antagonistic interactions (figure 4.1(c)) bend the curves outward, indicating that the growth rates are lowered less than the two dosages would accomplish independently. There is a large and dedicated literature on characterizing the interactions of many different combinations of toxins on static cell populations. As far as we are aware, the first comprehensive study of synergistic vs antagonistic effects was carried out by Bliss [14] in 1939, using joint probabilities, leading to a formula that is commonly called the Bliss index for drug interactions. A similar but slightly modified criterion was introduced by Loewe [71] and more recently developed further by Chou and collaborators [25]. These indices have all been used to help quantify the many different types of interactions that can occur with two or more toxins applied jointly in a static population of cells. In this context, it is common to assume that synergistic interactions are desirable in most circumstances, as a lower total dosage accomplishes the same kill rate as a higher dose would accomplish if the drugs acted independently. These kinds of studies have been used effectively to choose appropriate drug cocktails to individual patients by testing wide ranges of combinations on tissue samples obtained from patient tumors [110]. When the interacting population of cells are evolving, however, the relevant criteria become more complex. This is due to the fact that the subpopulations of cells respond differently to the different toxins applied, and as they respond, an ever-changing (adaptive) 58 Figure 4.1: (a) e > 0 showing synergistic profiles; (b) e = 0 showing additive profiles; (c) e< 0 showing antagonistic profiles. combination of toxins might be required to accomplish a given goal. Instead of necessarily killing the maximum number of cancer cells with the least amount of toxin, it is often the case that the goal becomes avoiding chemo-resistance and delaying unwelcome tumor recurrence to the maximum extent possible. A strategy called resistance management [94] is often advocated and occasionally implemented [10]. These kinds of strategies have been advocated and implemented in chemotherapy settings [42, 39, 134, 37, 122, 121, 9], but perhaps have been most elegantly and thoroughly carried out in a bacterial setting (since experimentsaremorepractical)byKishonyandcollaborators[50,76,3,100,11,57]whoeven discuss strategies that might reverse antibiotic resistance [11]. See also [128] for recent work discussing both microbial populations and cancer cells and [82, 131] for novel sequential therapy methods. In the context of evolving microbial populations [70], mutations occur frequently and it is important to consider not only pre-existing mutated subpopulations, but also mutations that occur as a result of the application of antibiotic agents. In the case of chemotherapeutic resistance in tumors, it is often assumed that resistant mutations occured before the application of treatment, hence it is common to separate the subpopulations into 59 sensitive and resistant subpopulations as we do in our deterministic model which does not include further mutations during treatment. See [16] for further discussions of these and related issues, and [53, 15] for discussions of the general approach of using evolutionary game theory in biology. Insection§IIweintroducethedetailsofthethree-componentreplicatordynamicalsystem that we use. Section §III describes the effects of constant chemotherapy schedules on the evolving populations, using the range of values 0≤ C 1 ≤ 1, 0≤ C 2 ≤ 1 (with a total dose upperthresholdC 1 +C 2 ≤ 1),alongwithourparametereoverarangeofpositive(synergistic) to negative (antagonistic) values. We describe in detail the transcritical bifurcations that occur and the tumor growth in response to the chemotoxins. In §IV we introduce adaptive time-dependent schedules (C 1 (t),C 2 (t)) along with the parametere with the goal of delaying tumorrecurrenceandwediscusstherateofadaptationofthecellpopulationsinthiscontext. Finally in §V we discuss the relevance of our model to the design of adaptive-therapy clinical trials. 4.2 A three-component replicator system 4.2.1 Why the prisoner’s dilemma game? To understand why the prisoner’s dilemma (PD) game is a useful paradigm for tumor growth, consider a two-component system, with a population of healthy cells H, and a population of sensitive cancer cells S, each competing for the best payoffs. The standard version of the prisoner’s dilemma payoff matrix is: A = a b c d = 3 0 5 1 , (4.1) where the first row corresponds to my payoffs if I choose to cooperate (associated with the healthy cell population (H)) and the second row corresponds to my payoffs if I choose to defect (associated with chemo-sensitive cells (S)). I am better off defecting than cooperating 60 if my opponent cooperates (first column), since my payoff is 5 instead of 3. I am also better off defecting than cooperating if my opponent defects (second column), since my payoff is 1 rather than 0. No matter what my opponent does, I am better off defecting. But my opponent uses the same logic and comes to the same conclusion - that she is also better off defecting, no matter what I do. The result is that we both defect and get a payoff of 1 (secondrow, secondcolumn)whichislessthanwewouldhavegottenhadwebothcooperated (payoff of 3). The fact that defect-defect is a sub-optimal Nash equilibrium (neither player is better off by making a different choice) is the essence of the prisoner’s dilemma game [52, 85] and the reason it is so widely used as the basis for modeling the evolution of cooperation in many different contexts. For cancer cell modeling, the healthy cells are the cooperators and the cancer cells are the defectors. Unlike other contexts in which game theory is used, cells are not making strategic decisions, instead, their strategy is encoded in their reproductive prowess, and selection is frequency dependent. In any mixed population ~ x = (x H ,x S ) T , 0≤x H ≤ 1; 0≤x S ≤ 1; x H +x S = 1, the fitness functions, ~ f = (f H ,f S ) T , associated with the two subpopulations are: ~ f = A~ x, (4.2) which in component form yields: f H = (A~ x) 1 = 3·x H + 0·x S , (4.3) f S = (A~ x) 2 = 5·x H + 1·x S , (4.4) while the average fitness of the total population is given by the quadratic form: hfi = ~ x T A~ x = 3x 2 H + 5x H x S +x 2 S ≥ 1. (4.5) The average fitness of the healthy state (x H ,x S ) = (1, 0) is given byhfi| (x H =1) = 3, while that of the cancerous state (x H ,x S ) = (0, 1) is given byhfi| (x S =1) = 1, which minimizes the 61 average fitness. For the static game, the cancerous state (x H ,x S ) = (0, 1)≡p ∗T is a strict Nash equilibrium sincep ∗T Ap ∗ >p T Ap ∗ , for allp [52]. We can then embed this static game into an evolutionary context using the replicator dynamical system [85]: ˙ x H = (f H −hfi)x H , (4.6) ˙ x S = (f S −hfi)x S , (4.7) from which (using values from eqn (4.1)) it is straightforward to show: ˙ x S =x S (1−x S )(2−x S ), (4.8) with fixed points at x S = 0, 1, 2. The cancerous state (x s = 1) then becomes an asymptoti- cally stable fixed point of the dynamical system (4.8) and an evolutionary stable state (ESS) of the system (4.6), (4.7) which serves to drive the system to the strict Nash equilibrium under the flow. The fact that this ESS also corresponds to the one with the lowest average fitness is an extra feature of the PD game. For any initial condition containing at least one cancer cell: 0<x S (0)≤ 1, we have: (i) x S → 1, x H → 0 as t→∞ (ii)hfi→ 1 as t→∞. The first condition (and the structure of the nonlinear equations) guarantees that the can- cer cell population will saturate at the carrying capacity of 1 in an S-shaped (logistic) growth curve, while the second guarantees that this asmptotically stable carrying capacity is sub-optimal, sincehfi| (x S =1) <hfi| (x H =1) . For these two reasons, the prisoner’s dilemma evolutionary game serves as a simple paradigm for tumor growth both in finite population models, as well as replicator system (infinite population) models [118, 119, 122, 114]. This two-component system alone, however, is not able to account for the evolution of resistance. 62 4.2.2 The three-component model The model we employ is a three-component replicator dynamical system for the three subpopulations of cells: (S,R 1 ,R 2 )≡ (x 1 ,x 2 ,x 3 ) ˙ x 1 = (f 1 −hfi)x 1 , (4.9) ˙ x 2 = (f 2 −hfi)x 2 , (4.10) ˙ x 3 = (f 3 −hfi)x 3 , (4.11) where each dependent variable represents the relative frequencies of cells, with x 1 +x 2 + x 3 = 1. In these equations, f i represents the fitness of subpopulation i = 1, 2, 3, while hfi represents the average fitness of all three subpopulations. These equations then give rise to the obvious interpretation that if a given subpopulation’s fitness is above(below) the average, it grows(decays) exponentially - reproductive prowess is directly associated with the deviation of the fitness of a subpopulation from the average fitness of the entire population. The fitness functions are frequency-dependent (i.e. non-constant), which couples eqns (4.9)-(4.11) nonlinearly. The subpopulation fitness f i (i = 1, 2, 3) function is given by: f i = 1−w i +w i (A~ x) i (4.12) A = a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 (4.13) whereA is a 3×3 payoff matrix which introduces the evolutionary game being played byx 1 , x 2 , andx 3 . In this paradigm, the sensitive population (x 1 ) are the defectors (higher fitness) and both groups of resistant cells (x 2 and x 3 ) are the cooperators (lower, but not equal fitness). We use 0≤ w i ≤ 1 as a parameter to determine the strength of selection in the 63 system. Whenw i ∼ 0, selection is relatively weak and the evolutionary game does not play a big role in the balance of the three subpopulations. Whenw i ∼ 1, selection is strong, and the game plays a bigger role. Both of those limiting cases have been discussed in the literature [129, 117], but the tumor response through the full range of values is poorly understood. We use the selection parametersw i as our mechanism to introduce chemotherapy, by way of altering the relative fitness of the three subpopulations: w 1 = w 0 (1−C 1 −C 2 −eC 1 C 2 ) (4.14) w 2 = w 0 (1−C 1 ) (4.15) w 3 = w 0 (1−C 2 ). (4.16) Here, C 1 is the chemotherapy parameter associated with drug 1, C 2 is that for drug 2, and e is our synergy (e> 0) vs. antagonism (e< 0) parameter, and we takew 0 = 0.1 which sets the timescale in our simulations. Level curves of the effective fitness functions defined by eqns (4.14)-(4.16) are shown in figure 4.1. The average fitness of the population is given by: hfi =f 1 x 1 +f 2 x 2 +f 3 x 3 . (4.17) The condition for the payoff matrix A to be of (PD) type is: a 21 <a 11 <a 22 <a 12 (4.18) a 31 <a 11 <a 33 <a 13 (4.19) a 32 <a 22 <a 33 <a 23 (4.20) 64 and for definiteness, we choose the specific values: A = 2 2.8 2.8 1.5 2.1 2.3 1.5 1.8 2.2 (4.21) For the static game, it is easy to show that the cancerous state (S,R 1 ,R 2 ) = (x 1 ,x 2 ,x 3 ) = (1, 0, 0)≡p ∗T is a strict Nash equilibrium in the absence of chemotherapy (C 1 = 0,C 2 = 0) and an ESS for the replicator system, as shown in the diagram of figure 4.2(a) where the entire triangular region is the basin of attraction for the S population. A seperate important ingredient in our model is our tumor-growth equation, which we take as: ˙ x tumor = (hfi−g)x tumor (4.22) where the growth(decay) of the tumor is a function of the average fitness associated with the tumor minus a constant background fitness level g, associated with the surrounding tissue (say healthy cells) and micro-environment [123]. When the average fitness level of the population of cancer cells is higher than g, the tumor grows, and when it is lower, it regresses. Our chemotherapy functions C 1 and C 2 largely control this complex dynamic by modifying the Nash equilibria and ESS’s of the system via the fitness function (4.12). 4.3 Constant chemotherapy First, we examine the dynamical system for different constant levels of the chemo- parameters C 1 , C 2 (0 ≤ C 1 ≤ 1; 0 ≤ C 2 ≤ 1; C 1 + C 2 ≤ 1) in the case of additive e = 0 interations, synergistic (e> 0) interactions, and antagonistic (e< 0) interactions. 4.3.1 Additive interactions e = 0 Figure 4.2 shows the panel of trajectories for three different scenarios. To start, figure 4.2(a) shows the trajectories with no chemotherapy – all trajectories lead to the S corner 65 which saturates the tumor. In these figures, solid circles denote the ESS states, while open circles denote the unstable states. Figure 4.2(b) shows trajectories withC 1 = 0,C 2 = 0.8. In this case, competitive release of the R 1 population allows it to take over the tumor, with all trajectories leading to the R 1 corner. For these parameter values, R 1 is the ESS (and strict Nash equilibrium) of the system. Figure 4.2(c) with C 1 = 0.8,C 2 = 0 depicts competitive release of the R 2 population, with all trajectories leading to the R 2 corner, showing that R 2 is the ESS of the system. In figure 4.2(d) we plot three different constant therapy schedules together, showing how they can intersect at different times. This gives the possibility of switching the therapies off and on at the intersection times in order to create a trajectory that stays in a closed loop and never reaches any of the corners – these closed loops represent scenarios in which the three subpopulations stay in perpetual competition, driven by a time-dependent schedule that shapes the fitness landscapes in such a way as to manage chemotherapeutic resistance. We will describe the systematic construction of these closed loops in §4.4. Figure 4.3 shows the tumor size (using eqn (4.22)) plotted logarithmically. The untreated tumor (C 1 = 0,C 2 = 0) grows exponentially, while each of the treated tumors initially show tumor regression up until t∼ 75, then tumor recurrence due to chemotherapeutic resistance (either of the R 1 or R 2 populations), an unwelcome common scenario. This scenario is clearly depicted in prostate cancer data sets shown in [120] and widely discussed in the clinical literature. Our goal is to show how specially designed adaptive multi-chemotherapies can push the tumor recurrence point further to the right on these plots. The full range of possible profiles are shown in figure 4.4. For certain ranges of chemo- dosing, there are mixed basins of attraction to each of the corners, hence multiple evolution- ary stable states of the system. In the top row, withC 1 = 0, the two critical values where the ESS’s bifurcate areC 2 = 1/3 (R 1 changes from unstable to stable via a transcritical bifurca- tion) andC 2 = 1/2 (S changes from stable to unstable via a transcritical bifurcation). In the left column, withC 2 = 0, the two critical values where the ESS’s bifurcate areC 1 = 7/18 (R 2 changes from unstable to stable via a transcritical bifurcation) andC 2 = 1/2 (S changes from 66 stable to unstable via a transcritical bifurcation). This panel shows the complexity associ- ated with the stability and basins of attraction of the three tumor subpopulations even in such a simple model, but also gives us the opportunity to exploit the inherent underlying dynamics associated with the different trajectories using piecewise constant chemo-dosing protocols. (a) (b) (c) (d) Figure 4.2: Depiction of evolutionary stable states (ESS) for different chemotherapy values. (a) With no chemotherapy, the tumor saturates to the S corner regardless of the initial make-up of the three subpopulations. (b) WithC 1 = 0, andC 2 = 0.8, competitive release of the resistant population R 1 drives all trajectories to the R 1 corner. (c) With C 1 = 0.8 and C 2 = 0, competitive release of the resistant population R 2 drives all trajectories to the R 2 corner. (d) Trajectories associated with three different constant combinations of C 1 andC 2 , depicting the overlap of the trajectories at different times. 67 Figure 4.3: Tumor growth curves (log-plots) for untreated and constant therapies. Tumor recurrence (dashed line) occurs at t∼ 75 dimensionless time units. 4.3.2 Non-additive interactions Using the panel in figure 4.2 as our guide, we focus in more detail on the valuesC 1 = 0.35, C 2 = 0.27, comparingtheeffectsofsynergisticinteractionsandantagonisticinteractionswith values e =−0.4,−0.3,−0.2,−01., 0, 0.1, 0.2, 0.3, 0.4 in figure 4.5. For strongly antagonistic interactions (e =−0.4,−0.3,−0.2), only the two subpopulationsS andR 2 compete for dom- inance (both ESS), with theR 1 population being an unstable state. For strongly synergistic interactions (e = 0.2, 0.3, 0.4), the two subpopulations R 1 and R 2 compete for dominance (both ESS), with the population of cells S sensitive to both drugs being an unstable state. It is the intermediate regime e =−0.1, 0, 0.1 that is the most interesting, and has all three evolutionary stable populations competing for dominance (all three ESS) with intertwined basins of attraction for each. We examine the details of the bifurcations that occur between the two-species co-existence and three-species co-existing states next. 68 Figure 4.4: Panel showing the full range of trajectories, ESS and their basins of attraction for constant additive therapies in the range 0≤C 1 ≤ 1, 0≤C 2 ≤ 1, C 1 +C 2 ≤ 1. Bifurcation values along top row are C 2 = 1/3, C 2 = 1/2. Bifurcation values along left column are C 1 = 7/18, C 1 = 1/2. Transcritical bifurcations Infigure4.6, wefocusontherangeofantagonisticvaluese =−0.21,−0.19,−0.175,−0.15 for fixed values of C 1 and C 2 . The relevant bifurcation that occurs at the critical value e (a) c =−310/1701 =−0.1822 takes place at the R 1 corner, when the fixed point R 1 = 1 goes from unstable (e<e (a) c ) to asymptotically stable (e>e (a) c ) in a transcritical bifurcation (exchange of stability). Through the bifurcation point, a stable fixed point (shown in figure 4.6(a) outside the triangle below R 1 ) moves up the R 1 −S side of the triangle (R 2 = 0), 69 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 4.5: The effect of synergy (e> 0) vs. antagonism (e< 0) on the ESS and the basins of attraction for a fixed choice of constant combination therapies C 1 = 0.35, C 2 = 0.27. (a) e =−0.4; (b) e =−0.3; (c) e =−0.2; (d) e =−0.1; (e) e = 0; (f) e = 0.1; (g) e = 0.2; (h) e = 0.3; (i) e = 0.4. and exchanges stability with the fixed point at the R 1 corner. In figure 4.7 we show the details of the collision of eigenvalues that takes place (figure 4.7(a)) and the process in the dS/dt vs. S plane (figure 4.7(b)-(e)). Figure 4.7(b) shows the classic transcritical bifurcation diagram (see [109]). When the level of antagonism is sufficiently large, there are only the two evolutionary stable states S and R 2 . In figure 4.8 we highlight the bifurcation that takes place in the synergistic regime around values e = 0, 0.1, 0.2, 0.3 for fixed values of C 1 and C 2 . The transcritical bifurcation occurs 70 at the critical valuee (s) c = 30/189 = 0.1587 takes place at theS corner, when the fixed point S = 1 goes from unstable (e>e (s) c ) to asymptotically stable (e<e (s) c ). In figure 4.9 we show the details of the bifurcation diagram that governs this synergistic transcritical bifurcation. When the level of synergism is sufficiently large, there are only the two evolutionary stable states R 1 and R 2 . In figure 4.10 we show the areas of three basins of attraction through the full range of values−0.3≤ e≤ 0.3. The basin areas begin to rapidly change in the antagonistic regime at e =−0.2 and, in general, show much more sensitivity to changes in e in the antogonistic regime than the synergistic regime. Figure 4.11 shows the average fitness curves (figure 4.11(a)) and tumor growth curves (figure 4.11(b)) through a range of values of e. Notice, in the beginning, the average fitness of the antagonistic case e =−0.3 is higher, but ends up lower over time. The tumor growth curve in this case follows the same trend - initially tumor growth is highest, then ends up lowest. Tumor recurrence time for the antagonistic case is also pushed later in time (t∼ 300) showing the effectiveness of antagonistic drug interactions in terms of managing resistance. 4.4 Designing adaptive schedules Many new possibilities can be created with time-dependent chemoschedulesC 1 (t),C 2 (t), if we monitor the balance of the subpopulations and adaptively make changes at judiciously chosen time points. The basic idea is shown in figure 4.2(d) where we see how the trajectories associated with different constant chemotherapy schedules cross. At any of the crossing times, it is possible to switch from one trajectory to another by switching the values ofC 1 or C 2 at those crossing times (typically termed bang-bang control [64]). This basic procedure allows us to design schedules that take us from any point A in the triangle to any other point B along the legs of a path that are separated by the switching times. Using multiple time-switching, we can also design trajectories that form closed loops and never converge to any of the corners. 71 (a) (b) (c) (d) Figure 4.6: The opening of the basin of attraction for R 1 via a transcritical bifurcation at e =−310/1701 =−0.1822, C 1 = 0.35, C 2 = 0.27. At the bifurcation point, the areas of the basins of attraction of the respective regions areS = 76.2%,R 2 = 13.8%. (a)e =−0.21; (b) e =−0.19; (c) e =−0.175; (d) e =−0.15. Figure 4.12 shows examples of how closed (piecewise differentiable) orbits are designed in practice. In each of the figures, point O is fixed, as is the untreated (blue) curve with C 1 = 0, C 2 = 0. Consider the loop OABO created figure 4.12(a), with additive interaction parameter e = 0. In traversing the OA leg, we use C 1 = 0.5, C 2 = 0.2. When the trajectory reaches point A, we switch to C 1 = 0, C 2 = 0, i.e. no therapy. When we reach point B, we switch toC 1 = 0.2,C 2 = 0.5 until we reach point O, and then we start the schedule again to traverse the same loop as many times as we desire. The dose density plot is shown in figure 72 (a) (b) (c) (d) (e) Figure 4.7: Four ways of depicting the antagonistic transcritical bifurcation. (a) Eigen- value collision that takes place as the two fixed points collide at the R 1 corner. The other eigenvalue remains negative for both fixed points; (b) Transcritical bifurcation diagram; (c) Pre-bifurcation phase plane; (d) Bifurcation phase plane; (e) Post-bifurcation phase plane. 4.12(b). Figures 4.12(c),(d) uses the same dosing values, but with e = 0.3 (synergistic). The loop in this case is larger (encloses more area), so there is a larger deviation in the subpopulations throughout the loop than there was for the additive case. Figure 4.12(e),(f) shows an example of an antagonistic e =−0.3 adaptive therapy loop. Table I shows the total doseD, time periodT, and average doseD/T for each case. The antagonistic adaptive loop delivers the smallest total dose over the shortest time-period to achieve one closed loop, whereas the synergistic adaptive loop delivers the smallest average dose over the loop, since the time-period is longest. Figure4.13showsthetumorgrowthcurvesforeachoftheadaptiveschedulesascompared with the untreated growth curve (exponential growth), and constant therapy curve (each 73 (a) (b) (c) (d) Figure 4.8: The closing of the basin of attraction for S via a transcritical bifurcation at e = 30/189 = 0.1587,C 1 = 0.35,C 2 = 0.27. At the bifurcation point, the areas of the basins of attraction of the respective regions are R 1 = 72.3%,R 2 = 17.7%. (a) e = 0; (b) e = 0.1; (c) e = 0.2; (d) e = 0.3. e -0.3 0 0.3 D 118.23 119.42 120.47 T 178.7 192.5 206.7 D/T 0.6616 0.6204 05828 Table 4.1: Total dose (D), Time period (T), and Average Dose (D/T) associated with adap- tive therapies with antagonistic, additive, and synergistic drug interactions. shows tumor regression followed by recurrence). In each case, the adaptive loop overcomes recurrence, with the antagonistic schedule minimizing the tumor re-growth leg (AB) of the schedule. 74 (a) (b) (c) (d) (e) Figure 4.9: Four ways of depicting the synergistic transcritical bifurcation. (a) Eigenvalue collision that takes place as the two fixed points collide at the S corner. The other eigen- value remains negative for both fixed points; (b) Transcritical bifurcation diagram; (c) Pre- bifurcation phase plane; (d) Bifurcation phase plane; (e) Post-bifurcation phase plane. A final metric that we use for comparisons is the rate of adaptation of the subpopulation i, defined by: α i (t) = 1 t Z t 0 x i (s)ds, (4.23) which is the time-average of subpopulation level x i (t). Figure4.14showstheratesofadaptationassociatedwiththetworesistantsubpopulations R 1 andR 2 duringthecourseoftheadaptiveschedules. Noticetherateofadaptationislowest for the antagonistic interaction, which is the main reason the tumor growth curve in figure 4.13(c) is most effective at controlling and delaying tumor recurrence as the cell population evolves. 75 Figure 4.10: Areas of the three basins of attraction as a function of the parameter e. Note the sensitivity of the S and R 1 regions in the antagonistic regime−0.2≤e≤−0.1. (a) (b) Figure 4.11: (a) Average fitness curves for e = 0.3, 0.15, 0,−0.15,−0.3; (b) Tumor growth curves (log-plot) for e = 0.3, 0.15, 0,−0.15,−0.3. 4.5 Discussion The role of synergistic vs. antagonistic combination drug interactions on the dynamical balance of evolving subpopulations of cells is not simple to characterize. In general terms, 76 (a) (b) (c) (d) (e) (f) Figure 4.12: Closed loop adaptive schedules OABO using three cycles for e = 0, 0.3,−0.3. (a)e = 0, OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (b) Corresponding adaptive schedule; (c) e = 0.3, OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (d) Corresponding adaptive schedule; (e) e =−0.3, OA: C 1 = 0.5,C 2 = 0.2, AB: C 1 = 0,C 2 = 0, BO: C 1 = 0.2,C 2 = 0.5; (f) Corresponding adaptive schedule. (a) (b) (c) Figure 4.13: Tumor growth curves for the adaptive schedules from figure 4.12 as compared with untreated growth (exponential) and constant schedule which eventually leads to tumor recurrence. (a) e = 0; (b) e = 0.3; (c) e =−0.3. antagonistic interactions are able to exert a more targeted and subtle range of influences on an evolving population than synergistic ones which, roughly speaking, cause the two drug combination to effectively act as one. This limits their flexibility for designing effective strategies to manage resistance, although might be useful in producing a bigger killing effect 77 (a) (b) Figure 4.14: Rates of adaptation for e = 0, 0.3,−0.3: (a) R 1 rate of adaptation; (b) R 2 rate of adaptation; (c) x 1 rate of adaptation; (d) x 2 rate of adaptation; with a smaller dose. As a by-product of the ability to deliver a more nuanced effect on the balance of cells, the sensitivity to small changes in the relative doses of the two drugs in an antagonistic setting seems to be higher than in a synergistic setting (as shown best in figure 4.10). This, perhaps, makes these interactions harder to control. Tumor recurrence times for the antagonistic drug interactions are delayed more effectively than for synergis- tic interactions, consistent with the fact that adaptation rates are slower for antagonistic interactions. Anatagonistic drug interactions provide more flexibiility in designing adaptive resistance management schedules. Additional features to our model could be the introduc- tion of mutations that might occur in response to the chemotherapy, which potentially could be handled using a mutational replicator system [85], or a finite-cell Moran process based model [118]. It might also be possible to analyze existing individual patient data and tumor response curves (as was done in [120]) to design and optimize better multi-drug strategies retrospectively, which is ongoing by our group. There are also ongoing adaptive multi-drug clinical trials at the Moffitt Cancer Center that show promise in prostate cancer patients. 78 Chapter 5 Optimal control of evolutionary games: The two-by-two case In this chapter, we focus on the optimal control problem in the replicator system with two species competing based on evolutionary games. The controlu(t) in this problem repre- sents an external/manual factor that alters the fitness landscapes, thus it alters the entries of the payoff matrix and changes the game type the species play. Additionally, we place an extra constraint on the integral of the control u(t), which indicates the external effect could only function for finite time. After the Pontryagin Maximum Principle is applied to the dynamical system according to an appropriate payoff function, we find that the stan- dard replicator equations and the adjusted replicator equation may have different behaviors. For example, if we want to minimize the final frequency of the defector in the Prisoner’s dilemma game, the standard replicator equation always has bang-bang control as optimal, while the adjusted replicator equation may have a very complex time-dependent function as the optimal solution. However, the adjusted replicator equation is a nonlinear re-scaling of the replicator equation. The comparison between two replicator systems highlights that adjusted replicator equations differs from the standard replicator equation when the timing is sensitive, especially in optimal control problems. The solution to this optimal control problem, therefore, helps to build the proper models for each individual problem in evolu- tionary games and it could contribute to the drug-scheduling, adaptive therapy, or resistance prevention when these models are applied to clinical data. 79 5.1 Introduction In an evolutionary process, the final outcome may not be ideal. It is desirable to develop a control theory framework for the evolutionary process in order to select for particular phenotypes or genotypes [28]. Then modeling and optimizing this whole process requires control and optimal control techniques, and there are a lot of applications reported, such as prey-predator system[101], drug dosage design and resistance prevention[55, 13], or coopera- tion networks[2]. Particularly, a lot of literature have discussed about optimal chemotherapy scheduling for cancer with differential equations[73, 111, 125, 79, 22], or stochastic birth death process[126, 26]. Thispaperisalsoinspiredbytheoptimizationforchemotherapy, butwekeeptheanalysis general, withoutincludingmuchdetailassociatedwithcancermodeling. Insection5.2weuse an evolutionary game model[108] to determine the fitness landscape of the two competing species. The control u(t) alters the fitness, thus it actually change the game type and the corresponding evolutionary stable strategy[108, 85, 52]. The system is then evolved with standard replicator equations or adjusted replicator equations. After that, we use the PontryaginMaximumprinciple[67,66,98,113]tofindouttheoptimalcontrolofthedynamic system in section 5.2.2. In section 5.3, we analyze the properties of the optimal control, and we find that standard replicator equations and adjusted replicator equations may have different behavior. With appropriate target function for optimization, the standard replicator will have bang-bang control to be optimal, while the adjusted replicator equations may have more complex time- dependent optimal control. Numerical simulations of a prisoner’s dilemma game[87, 85, 52] are provided in section 5.4. Because the Prisoner’s dilemma game has the defector as the Nash equilibrium[87] while the defector is not optimal for maximizing payoff, we use the Prisoner’s dilemma game as a cancer model[124] and controlling the defector is analog to controlling the tumor. Numerical results shows that standard replicator equation always has bang-bang optimal control to minimize the final defector. However, when the average fitness increaseswiththecontrolvalue, theadjustedreplicatorequationwillhaveaslowlyincreasing 80 function as the optimal control. The adjusted replicator equation has been proved to be the continuous limit of the Moran process[114, 115], and similar results has been reported[126] that high entropy performs the best with Moran process simulation. The reason of the difference is explained in section 5.45. High average fitness will weaken the effects of the control, and the bang-bang control in this case has the highest fitness when u(t) = 1, so it won’t be optimal any more. 5.2 model 5.2.1 evolutionary games and replicator dynamics Nowconsiderthereexisttwocompetingspecieswithx 1 andx 2 astheindividualfrequency, so the sum always equal to one x 1 +x 2 = 1. The selection is based on the fitness landscape of the system, and the individual fitness is determined by the 2-by-2 payoff matrix A: f i = (A~ x) i (5.1) A = a 11 a 12 a 21 a 22 (5.2) where ~ x = [x 1 ,x 2 ] T and i = 1, 2.Suppose u(t) is the control for the system, and it could be any external influence such as drugs for some disease, disinfection for the bacterial, foods for one single animal, etc. The value of the control u(t) is always positives but smaller than 1 at any time. 0≤u(t)≤ 1 (5.3) The fitness will be altered when the controlu(t) is none zero. Without losing generosity, we assume the altered selection tends to be in favor of x 1 , which means it will either reduce the fitness of x 2 or increase the fitness of x 1 . Suppose the change of fitness is proportional to the control and the coefficient is the linear combinations of the frequencies as shown in 81 eq.5.4.α 1 and α 2 are constant. There’s no constant in the coefficient because adding any constant is equivalent to adding it to both α 1 and α 2 . Δf = (α 1 x 1 +α 2 x 2 )u(t) (5.4) First if f 2 is reduced by Δf, then let’s denote the altered fitness as f 2 0 . Because of the change of the fitness landscape, the game type, indicated by the new payoff matrix A 0 , is varied as well. Reducing the fitness of the x 2 is equivalent to reducing the entries in the second raw. Similarly, if f 1 is increased, then entries in the first row of the payoff matrix will be increased as well. f 2 0 = f 2 − Δf⇔A 0 = a 11 a 12 a 21 −α 1 u a 22 −α 2 u (5.5) f 1 0 = f 1 + Δf⇔A 0 = a 11 +α 1 u a 12 +α 2 u a 21 a 22 (5.6) Once we construct the fitness landscape, the system is simulated with either the standard replicator equation or the adjusted replicator equation. It is worth pointing out that both of them can be transformed into an single variable ODE by eliminating x 1 . We will use this result to simplify our system later. ˙ x i = (f i −hfi)x i (5.7) First let’s consider the standard replicator equation in eq.5.7 wherehfi =f 1 x 1 +f 2 x 2 is the average fitness. It’s quite simple to combine the formula of the fitness with the replicator equations. In eq.5.8 the standard replicator equation is rewritten with single variable x 2 . 82 Here we only put one single equation because increasing f 2 and reducing f 1 by Δf lead to the same ordinary differential equations, thus the same trajectory. ˙ x 2 = (a 21 −a 11 −α 1 u(t))(1−x 2 ) 2 x 2 (5.8) +(a 22 −a 12 −α 2 u(t))(1−x 2 )x 2 2 (5.9) f 2 0 = f 2 − Δf⇔f 1 0 =f 1 + Δf (5.10) However, if we take the adjusted replicator equation in eq.5.11 for the simulation, we can find that increase f 2 or reducing f 1 are quite different because of the average fitness hfi in the denominator. In the next section we will discuss both for the adjusted replicator equation. ˙ x i = (f i −hfi) hfi x (5.11) f 2 0 = f 2 − Δf<f 1 0 =f 1 + Δf (5.12) Since the control u(t) goes against x 2 , if we aim at minimizing the population of x 2 and there’s no extra constraint on the control u(t), the optimal control will simply be u(t) = 1 for all time. In reality, for example about drug scheduling, the total dosage the patient can take with a period is limited, and the task for the doctor is to distribute drugs over time for the patient. Therefore the extra constraint in eq.5.13 is put on the integral of the control u(t). Z T 0 u(t)dt =D (5.13) Where D is a fixed constant and T is the terminal time. For simplicity, we denote x 2 as x as the single variable in the remaining sections. x ≡ x 2 (5.14) x 1 = 1−x (5.15) 83 5.2.2 Pontryagin Maximum Principle The Pontriagin Maximum(Minimum) Principle[67, 66, 98, 113] is a widely used method to solve optimal control problem. Assuming that our target is to minimize a general cost function in eq.5.16 minimize : Z T 0 L(x(t),u(t),t)dt +ϕ(x(T )) (5.16) The standard form for the Maximum Principle is applied for dynamic system with boundary value constraint, so we need to convert the integral constraint in eq.5.13 into the standard form. Eq.5.17 defines a new variabled(t) so that the constraint is transformed into an ODE with two boundary values. Together with the standard replicator equation the new dynamic system is constructed in eq.5.21. d (t) = Z t 0 u(s)ds (5.17) ˙ d (t) =u(t) (5.18) d (0) = 0 (5.19) d (T ) =D (5.20) ˙ X = F (X) (5.21) X = [x(t),d(t)] T (5.22) F (X) = [(f 2 −hfi)x,u(t)] T (5.23) If we use the adjusted replicator equation rather than the standard ones, we only need to update F (X): F (X) = [( (f 2 −hfi) hfi x,u(t)] T (5.24) 84 Then following the notations from [113], we construct the Control Theory Hamiltonian as eq.5.25. λ 1 andλ 2 are the co-state function, or in Physics they correspond to momentum. H(x,d,λ 1 ,λ 2 ,u) = λ T F (x) +L(x,u,t) (5.25) λ = [λ 1 ,λ 2 ] T (5.26) Assumingthatu ∗ (t)istheoptimalcontrolforthisproblem, withcorrespondingtrajectory x ∗ (t),d ∗ (t), thereexistco-statefunctionsλ 1 ∗ (t)andλ 2 ∗ (t)thatsatisfythecanonicalequation 5.27 ˙ x ∗ (t) = ∂H(x ∗ (t),d ∗ (t),λ 1 ∗ (t),λ 2 ∗ (t),u(t)) ∂λ 1 (5.27) ˙ d ∗ (t) = ∂H(x ∗ (t),d ∗ (t),λ 1 ∗ (t),λ 2 ∗ (t),u(t)) ∂λ 2 (5.28) ˙ λ 1 ∗ (t) =− ∂H(x ∗ (t),d ∗ (t),λ 1 ∗ (t),λ 2 ∗ (t),u(t)) ∂x (5.29) ˙ λ 2 ∗ (t) =− ∂H(x ∗ (t),d ∗ (t),λ 1 ∗ (t),λ 2 ∗ (t),u(t)) ∂d (5.30) with the corresponding boundary condition: x(0) = x 0 (5.31) d(0) = 0,d(T ) =D (5.32) λ 1 (T ) = ∂ϕ(x(T )) ∂x(T ) (5.33) At any time point, the optimal control u ∗ (t) will minimize the Control Theory Hamilto- nian as shown in eq.5.34 u ∗ (t) = arg min u(t) H(x ∗ (t),d ∗ (t),λ 1 ∗ (t),λ 2 ∗ (t),u(t)) (5.34) Now the optimization problem turns into the two-point boundary value problem, and the solution will be the optimal control and the corresponding trajectory. 85 5.3 Results NowconsiderthecasewherethefunctionLonlyhasexplicitdependenceonthetrajectory L = L(x(t)). We can use eq.5.34 to find some useful properties of the optimal control functionu ∗ (t). For standard replicator equation and adjusted replicator equation, we denote the Control Theory Hamiltonian as H std and H adj respectively in eq.5.35 and 5.36 H std = λ 1 (f 2 −hfi)x +λ 2 u +L(x(t)) (5.35) H adj = λ 1 (f 2 −hfi) hfi x +λ 2 u +L(x(t)) (5.36) In addition, thef 1 andf 2 are affine functions of the controlu(t) no matter how you alter the fitness, thus the average fitnesshfi is also an affine function ofu(t). As a result, so does the Control theory Hamiltonian H std for the standard replicator equation. To minimize H std , the optimal controlu ∗ (t) depends on the sign of ∂H std ∂u . If ∂H std ∂u 6= 0 holds almost everywhere, the then optimal control will be the bang-bang control. u ∗ (t) = 0 ∂H std ∂u > 0 1 ∂H std ∂u < 0 (5.37) For the case of the adjusted replicator equation, things are more complex and hard to see the minimal values ofH adj . Finding the minima is first step by taking the partial derivative of u(t) and setting it to 0 in eq.5.38: ∂H adj ∂u = λ 1 x(hfi ∂f 2 ∂u −f 2 ∂hfi ∂u )−λ 2 hfi 2 hfi 2 (5.38) = 0 (5.39) As mentioned in section 5.2 f 1 ,f 2 ,hfi are all affine functions of u(t), so the solution to the minima in eq.5.38 is actually solving a quadratic equation. Also the minima may not be the global minimal. The control is bounded between 0 and 1 so even if the minima falls within 86 the range, we still have to compare the minima with the boundary value at 0 and 1. In conclusion,it’s possible that the optimal controlu ∗ (t) is still bang-bang control, but for some cases it may be some value between 0 and 1, different from the standard replicator equation. We will give an example in section 5.4 5.4 discussion Now we consider the case of chemotherapy dosage scheduling for cancer. The game between the cancer cells and surroundings healthy cells is a Prisoner’s dilemma game[124] in eq.5.40. The cancer cells play as a defector and the healthy cells are the cooperator. In a Prisoner’s dilemma game the defector, thus the cancer cell, is an evolutionary stable strategy (ESS) and Nash Equilibrium (NE)[87, 85, 52], but the average fitness for all cancer cell will be lower than the one with all healthy cells. We want to find the optimal scheduling of the chemotherapy u ∗ (t) so that the frequency of cancer cells is minimized at the end of the therapy. The corresponding cost function is updated as eq.5.41 A = 3 0 4 1 (5.40) L = 0 (5.41) ϕ(x(T )) = x(T ) (5.42) 5.4.1 f 2 decrease Now we consider the altered game in eq.5.43, and the fitness of the cancer cells is reduced with α 1 = 5 and α 2 = 0.When u > 1/4, the game type changes from prisoner’s dilemma game to coordinate game, andx = 1 becomes an evolutionary stable strategy for the system. Before we find the optimal solution for this problem,I want to introduce two pre-determined therapy schedule. One is the Maximum tolerated Dose (MTD), which provides the patient 87 as high dose as possible in the beginning then stops the drug, and the Low dose Metronomic therapy, which gives the patient constant dosage throughout the treatment. A 0 = 3 0 4− 4u(t) 1 (5.43) Figure 5.1: Schema for the Maximum tolerated Dose (MTD) and Low Dose Metronomic (LMD) schedule. The integral of these two schedule are the same. Fig.5.1 is an example of the MTD versus LDM with T = 1,D = 0.3. The corresponding numerical simulations are in Fig.5.2. We can see that in this case the MTD works better than the LDM for both standard replicator equation and adjusted replicator equation. We know that unless singular control exists, bang-bang control will be optimal for standard replicator equation. After solving the two point boundary value problem from eq.5.27, the optimal control for both ODE system is the bang-bang control, and the schedule is exactly the same as the MTD therapy in Fig.5.1. The corresponding trajectory x ∗ (t), therefore, is just the same as the MTD curves in Fig.5.2(a) and Fig.5.2(b). 88 (a) (b) Figure 5.2: Simulation results with D = 0.3,T = 1,x(0) = 0.1 and payoff matrix in eq.5.43 (a) the corresponding trajectory for MTD and LDM simulated by standard replicator equa- tion.(b)the corresponding trajectory for MTD and LDM simulated by Adjusted replicator equation. Figure 5.3: Bang-bang optimal control for both standard replicator equation and adjusted replicator equation 89 5.4.2 f 1 increases However, as we discussed in section 5.3, the bang-bang control may not be optimal for adjusted replicator equation. Here we consider increasing the fitness of x 1 by the same amount as in section 5.4.1, and the new payoff matrix is shown in eq.5.44 with increased entries in the first row. A 0 = 3 + 4u(t) 0 4 1 (5.44) Because we have proved in section 5.2 that the increasing the first row or decreasing the second row are the same for standard replicator equation, the corresponding trajectories for MTD and LDM schedule should be identical to Fig.5.2(a). The optimal control for the standard replicator equation will just be bang-bang control, the same as Fig.5.3, so we will focus on the Adjusted replicator equation for this section. Fig.5.4 is the trajectories simulated with adjusted replicator equation with D = 0.3,T = 1,x(0) = 0.1. Contradictory to what we get in section 5.4.1, the LDM method works better than MTD schedule so the MTD is no longer the optimal control. Similar results simulated by Moran Process has been reported that high entropy schedule, correlated with LDM meth- ods, outperform other schedules. After we numerically solve the two point boundary value problem in eq.5.27 based on Pontryagin Maximum principle, we have the optimal control u ∗ (t) and trajectory x ∗ (t) in Fig.5.5 The optimal controlu ∗ (t) is non-trivial in this case as an slowly decreasing function, but the value is always close to 0.3 so it is more similar to the LDM schedule and that’s the reason why the LDM has better performance than MTD method, or the regular bang-bang control. 90 Figure 5.4: Simulation results with D = 0.3,T = 1,x(0) = 0.1 and payoff matrix in eq.5.44 with Adjusted replicator equation. (a) (b) Figure 5.5: (a)optimal control u ∗ (t) for adjusted replicator equation (b) the corresponding trajectory x ∗ (t) 91 5.4.3 Time re-scaling To convert the adjusted replicator equation to the standard replicator equation, we need the transformation[114] for the timing in eq.5.45 ds = dt hfi (5.45) dx i ds = (f i −hfi)x i (5.46) This transformation actually varies the time structure, so after that the terminal timeT and the integral value D are varied as well. Eq.5.47 provides the new terminal time as Ts and Eq.5.48givesthenewintegralunderthenewtimecoordinates. Wecanseethatiftheaverage fitnesshfi increase a lot when the control u(t) increases, then after the transformation the time range of highu(s) value will be shorten, the integral of dose valueu(s) is also reduced. For example, for bang-bang control, the new time range foru(s) = 1 will be relatively shorter but the new time range for u(s) = 0 will become longer. That’s why in section 5.4.2 the bang-bang control doesn’t work very well because increasing the control u(t) dramatically increases the average fitnesshfi. In the other hand, if the average fitnesshfi decreases with the control u(t), the bang-bang control could be enhanced. T s = Z T s 0 ds = Z T 0 dt hfi (5.47) D s = Z T s 0 u(s)ds = Z T 0 u(t)dt hfi (5.48) 5.5 Conclusion In this paper we solve the optimal control problem for 2-by-2 evolutionary game evolving with standard replicator equation or adjusted replicator equation. After we simplify the problem into the standard form, the optimal control is found by the Pontryagin Maximum principle. We take the prisoner’s dilemma game as the numerical example and the standard replicator equation always has bang-bang control as optimal to reduce the defector, while 92 the adjusted replicator equation, possibly Moran process as well, may have a no-trivial time dependent function as optimal control, especially when the average fitness of the system increases with the control value. This results take a deeper insight into the evolutionary game models, and it will potentially help to find the optimal drug scheduling as a medical problem. 93 Chapter 6 Optimal control of levitation by blowing Is it possible to levitate an object with constant blowing? We show that it is not possible to levitate the object from the origin to a pre-specified height x ∗ using constant blowing, because the point at which the downward gravitational force balances the upward blowing forceisanellipticfixedpointofthegoverningequations, sothereisnotransfertrajectorythat connectstheorigintox ∗ . Toovercomethis, wedesigntime-dependentblowingschedulesthat achieve the transfer borrowing on orbital transfer ideas developed in the orbital mechanics literature. Then we ask which of these schedules are optimal. We show that, generally, it is bang-bang (on-off) blowing schedules that achieve the transfer in minimal time, using minimal energy, and minimal air volume. For certain parameter values, however, there are more complicated blowing schedules that are optimal (with respect to energy) which can be designed using the Pontryagin Maximum Principle (PMP) and singular control. The model we describe is frictionless and treats the object (i.e. a ping-pong ball) as a point mass. 6.1 Introduction Is it possible to levitate a ping-pong ball from rest to a fixed height by constant blowing? At first glance, one might imagine the answer is as simple as supplying a blowing force from below that exactly balances the gravitational force pushing down at the height x ∗ where you want the ball to levitate as shown in figure 6.1. An experiment, such as the one shown here: https://www.youtube.com/watch?v=bCRjPFhlSYk might even confirm that it is achievable, until one examines the problem more closely. Let the ball start from rest, at 94 x = 0 where the blower is located. If one supplies a constant blowing force that is required to balance the ball at height x ∗ > 0, it turns out, surprisingly, that it is not possible to drive the ball to the position required for the forces to balance. We will show why, and describe how to overcome the problem by using non-constant blowing schedules designed appropriately. We formulate the problem as one in nonlinear dynamics and optimal control and show how some blowing schedules are better than others. The methods are borrowed from orbit transfer ideas developed initially in the orbital mechanics literature [60] in which intersecting phase space trajectories at different fixed parameter values are pieced together by switching from one parameter value to the other at precisely the time at which the trajectories intersect, creating a time-dependent control protocol that drives the system to a desirable state not achievable without the switching time. This idea is also currently being developed to design time-dependent chemotherapeutic schedules that achieve better results than standard protocols, making use of the replicator dynamical system as a tumor growth model in the context of evolution. In this paper, our goal is to use the simplest possible model to elucidate how the method can work in settings more commonly developed using more standard control theory ideas, and we show that the methods yield equivalent results. Using air to levitate objects (as opposed to acoustic or magnetic levitation) is a common method of conveyance and has been explored in micro-manufacturing settings as discussed in [59], but to our knowledge, the observation that non-constant blowing is required to levitate to a fixed height has not been discussed in the literature. Reference [24] discusses the design of a low-cost air levitation system for teaching control engineering. 6.2 The model The model we consider is the simple dimensionless nonlinear system: m¨ x =−mg +mgb(t) exp(−αx) (6.1) 95 Figure 6.1: Ping-pong ball levitated by blowing through a straw. The streaming air blows upwards along the positive x direction against the gravity. x ∗ marks the equilibrium posi- tion where the forces balance. For full video, see https://www.youtube.com/watch?v= bCRjPFhlSYk Here,b(t) is a dimensionless blowing parameter (proportion to amount of air per unit of time) that we will use to levitate the ball from rest (v(0) = 0), x(t) is the vertical height, m is mass, and α > 0 is a constant that controls the distance over which the blower’s influence is felt. The ball starts at position x(0) = 0 where the blower is located, as shown in figure 6.1 and the goal is to levitate it to height marked x ∗ . Our assumption in (6.1) is that the influence of the blowing decreases exponentially from the location of the blower (x = 0), but our results are not sensitive to that choice, as long as the influence decreases monotonically with height. We also assume there is no frictional loss of energy as the ball levitates and that the ball does not rotate, hence we treat it as a point mass.We emphasize that our model is 96 simple in that it does not take into account the detailed fluid mechanical forces and dynamics at the jet-ball contact point (that would create such effects as dissipation and rotation of the ball). Our model is just complex enough to show that (i) it is not possible to drive the ball from its resting state at the base to the levitated state at height x ∗ with constant blowing, and (ii) it is possible to drive the ball from its resting state to the levitated state by piecing together time-dependent blowing schedules with precise switching times determined from crossing phase-space trajectories. To levitate the ball, the upward blowing force must exceed the downward gravitational force at x = 0, which implies that the right hand side of eqn (6.1) be greater than zero (at x = 0). The condition for this is that b> 1. The governing Hamiltonian for the system is given by: H(x,v) = 1 2 mv 2 +mgx + mgb exp (−αx) α ≡E b (6.2) where we highlight the fact that b is a parameter we will vary in order to change energy values. Notice for fixed x and v, it is evident that E b 1 >E b 2 when b 1 >b 2 since eqn (6.2) is linear in b. The corresponding Hamiltonian equations of motion are: m ˙ x = ∂H ∂v =mv (6.3) m ˙ v = − ∂H ∂x =−mg +mgb exp(−αx). (6.4) For any given b, the forces balance (i.e. acceleration is zero) when: x≡x ∗ = ln (b) α (6.5) where the right hand side of eqn (6.1) is zero. With the condition b > 1, the equilibrium position is positive x ∗ > 0. 97 Thephasespacediagraminthe (x,v)plane, obtainedfromplottinglevelcurvesH(x,v) = E b =constant, is shown in figure 6.2(a) for different energy values. The elliptic fixed point where the forces balance correspond to: E ∗ b =H(x ∗ , 0) = mg α (ln(b) + 1). (6.6) At the origin where the blower is located, the energy value is: E 0 b =H(0, 0) = mgb α >E ∗ b . (6.7) We show in figure 6.2(b) curves E 0 b for increasing values of b. We also show in figure 6.2(c) the family of free-fall trajectories (b = 0) associated with energy values E 0 = 1 2 mv 2 +mgx which start at (x,v) = (0,v ∗ ). To compute the pointx max shown in figure 6.2(a), we follow the phase curve that passes through the origin to the pointx max , using the fact that the Hamiltonian is constant, which gives rise to the transcendental equation for x max : αx max =b(1− exp(−αx max )). (6.8) The left and right sides of eqn (6.8) are plotted in figure 6.2(d) where the intersection defines x max . Because x ∗ is an elliptic fixed point, it is not possible to traverse from the origin to this point on a phase curve, i.e. there is no constant energy trajectory that passes both through the origin andx ∗ . Instead, the levitating ball starting at the origin will reach height x max >x ∗ , then drop back down to the origin in a periodic orbit, always overshooting height x ∗ where we want to position the ball. For future use, notice the velocity formula from eqn (6.2) is: v =± v u u t 2 m " E b −mgx− mgb α exp(−αx) # . (6.9) 98 6.3 Time-dependent blowing To levitate the ball from the origin tox ∗ we need to drop the energy from its initial energy valueE 0 b to the final energyE ∗ b by altering our blowing parameter b. Consider the family of phase curves for two different choices of the blowing parameter b =b 1 >b 2 > 0, as shown in figure 6.3. The diagram depicts two different paths to traverse from pointO to pointx ∗ . The blue curve corresponds to energy value E b = E 0 b for value b = b 1 , giving rise to the elliptic balance point x ∗ | b 1 = ln(b 1 ) α . The red curve that intersects x ∗ | b 1 corresponds to energy value E b = E 0 b obtained by choosing blowing parameter b = b 2 < b 1 , where x max | b 2 = x ∗ | b 1 . To traverse the red elliptical path that connects O to the resting point x ∗ , we first use blowing parameter value b 2 > 1 until the ball arrives at x ∗ , say at time t =t 1 . Then we can balance it there if we instantaneously change to valueb 1 >b 2 att 1 , which traps the ball at the elliptic fixed point. We call this a two-step schedule since it requires choosing two different value of b with one switch-time t 1 . To calculate t 1 , we need to integrate eqn (6.9) from O to x ∗ : t 1 = Z ln(b 1 )/α 0 dx r 2 m h E 0 b 2 −mgx− mgb 2 α exp(−αx) i (6.10) A useful quantity to track is the energy: hEi = Z t f 0 E(t)dt (6.11) required to achieve a given transfer, where t f is the time the ball reaches x ∗ . In this case, we have: hEi = E 0 b 2 ·t 1 (6.12) = mgb 2 α ·t 1 . (6.13) 99 A second quantity we track is the area under the b(t) curve: hb(t)i = Z t f 0 b(t)dt (6.14) which we think of as the total volume of air used to transfer the ball through time t f . In this case, we have: hb(t)i =b 2 ·t 1 . (6.15) A second path that takes the ball from O to x ∗ traverses through point A. Consider the piecewise differentiable path O−→ A−→ x ∗ . On the first piece from O−→ A we launch the ball using blowing parameter valueb 1 until timet 1 when we reach pointA. At that time, we set b = 0 and let the ball free-fall along path A−→ x ∗ until time t 2 when we arrive at x ∗ . For this, we use the free-fall trajectory given by: x = − 1 2 gt 2 +v A t +x A (6.16) v = −gt +v A (6.17) where (x A ,v A ) mark the position and velocity coordinates at pointA in figure 6.3. We know from eqn (6.17) that v = 0 (maximium height) at time t = v A /g. Therefore for the total time to traverse the free-fall trajectory, we have: t 2 −t 1 =v A /g. (6.18) Here, we use: t 1 = Z x A 0 dx r 2 m h E 0 b 1 −mgx− mgb 1 α exp(−αx) i . (6.19) 100 Eliminating t in eqns (6.16),(6.17) gives the parabolic curve: x = − v 2 2g + v 2 A 2g +x A ! (6.20) beginning at (x A ,v A ) and ending at the point (x ∗ , 0), where: x ∗ = v 2 A 2g +x A ! = ln(b 1 ) α . (6.21) At time t 2 we then let: b =b 1 = exp " α v 2 A 2g +x A !# (6.22) to trap the ball at this resting point. We call this a three-step sequence since it requires three sequential values of b (b = b 1 ;b = 0;b = b 1 ) and two switching times t 1 and t 2 . The total energy along the trajectory is: hEi = E 0 b 1 ·t 1 +E 0 · (t 2 −t 1 ) (6.23) = mg α [b 1 ·t 1 + ln(b 1 )(t 2 −t 1 )], (6.24) where E 0 is the free-fall energy with b = 0, and (x,v) = (x ∗ , 0), i.e. E 0 = mgx ∗ . For this trajectory, we have: hb(t)i =b 1 ·t 1 + 0· (t 2 −t 1 ) =b 1 ·t 1 . (6.25) Path Time t f EnergyhEi Air volumehb(t)i O−→x∗ 1.036 21.83 2.227 O−→A−→x∗ 0.6662 16.84 0.7482 Table 6.1: Time, Energy and Air volume values for two trajectories shown in figure 6.3. The results are simulated using variables: b 1 = 6,α = 1m −1 ,g = 9.8m/s 2 ,m = 1kg. 101 It is clear that once we allow for a time-dependent blowing functionb(t), there exist many schedules that allow us to achieve the goal of moving the ball from the origin to the resting point x ∗ , each taking a different total time to achieve the transfer and each expending a different amount of energy and air. We show in Table 6.1 numerically computed values for the two different paths described in figure 6.3 which shows in all cases, the three-step path: O−→ A−→ x ∗ (which includes a free-fall segment) is more efficient than the more direct two-step path: O−→ x ∗ . A more comprehensive panel of schedules, using different fixed values of the parameter b, and the trajectories that achieve the transfer to the levitation point, is shown in figure 6.4. We are now is a position to formulate the problem as one in optimal control theory to ask, of all possible schedules, which one achieves the transfer fromO tox ∗ in minimal time? Which transfer uses minimal energy? Which transfer uses minimal amounts of air? 6.4 Optimal blowing schedules Consider the problem of finding time-dependent blowing schedules b(t) and the corre- sponding trajectories (x(t),v(t)) in the phase plane that accomplish the task of moving the ball from (0, 0) to (x ∗ , 0) in total time t f , subject to the equations: m¨ x = −mg +mgb(t) exp(−αx) (6.26) x(0) = 0; ˙ x(0) = 0 (6.27) x(t f ) = x ∗ ; ˙ x(t f ) = 0. (6.28) We use the fact thatx≥ 0, and we assume the blowing function is bounded from below and above: 0≤b(t)≤b max . 102 6.4.1 Minimum time strategy The transfer time t f is given by: t f = Z t f 0 dt = Z x ∗ 0 dx dx/dt = Z x ∗ 0 dx v . (6.29) To minimize t f , we need to minimize the area under 1/v(x) from 0 < x < x ∗ , as shown in eqn (6.29). We claim that a bang-bang strategy, in which we choose b = b max initially, until switching time t s at the point (x,v) = (x s ,v s ), when we choose a free-fall trajectory b = 0, accomplishes the minimization. To see this, consider the trajectories associated with the two extreme (constant) values b = 0 and b =b max as shown in figure 6.5(a). The lower (red) trajectory is given by the curve E 0 bmax . At the switching time t =t s , we switch values of the blowing parameter to b = 0, so that the ball travels on a free-fall trajectory E 0 to the final resting point (x ∗ , 0) at timet f . These curves define a sector in the first quadrant of the (x,v) plane as shown in the figure. We show in figure 6.5(b) the same curves in the 1/v vs. x plane and consider the areas under the curve using E 0 bmax from 0 < x < x s , and E 0 from x s < x < x ∗ . In the interval 0 < x < x s , for any (v(t),x(t)) associated with the schedule b(t), the following function will be used to derive the maximum velocity v max for each fixed x: d dt 1 2 v(t) 2 +gx (t) + b max g α e −αx(t) ! (6.30) = (b (t)−b max )vge −αx(t) (6.31) ≤ 0 (6.32) Thus the initial value gives the upper bound: 1 2 mv(t) 2 +mgx (t) + mb max g α e −αx(t) (6.33) ≤ 1 2 mv(0) 2 +mgx (0) + mb max g α e −αx(0) (6.34) = mb max g α =E 0 bmax (6.35) 103 From this, we know that v(t) is bounded by v≤ v u u t 2 m " E 0 bmax −gx− b max g α e −αx # =v max (6.36) Hence: Z xs 0 dx v max < Z xs 0 dx v . (6.37) Similarly, in the interval x s <x<x ∗ : d dt 1 2 v(t) 2 +gx (t) (6.38) = b (t)vg (6.39) ≥ 0 (6.40) Thus the final value gives the upper bound: 1 2 v(t) 2 +gx (t) (6.41) ≤ 1 2 v(t f ) 2 +gx (t f ) (6.42) = gx ∗ (6.43) From this, we know that v(t) is bounded by v≤ q 2 [gx ∗ −gx] =v max (6.44) Hence: Z x ∗ xs dx v max < Z x ∗ xs dx v . (6.45) 104 Therefore the minimum-time transfer is achieved using the bang-bang strategy. An example of an optimal time trajectory and schedule is shown in figure 6.6. 6.4.2 Minimum air volume strategy To minimize eqn (6.14), first consider that: ˙ v = dv(x) dx dx dt =v 0 (x)v (6.46) = −g +gb(t) exp(−αx), (6.47) which gives: b = v 0 (x)v +g g exp(−αx) . (6.48) From this, we have: Z t f 0 b(t)dt (6.49) = Z x f 0 b(t) dx v (6.50) = Z x f 0 v 0 (x)v +g g exp(−αx) dx v (6.51) = Z x f 0 v 0 (x) g exp(αx) + Z x f 0 exp(αx) v dx (6.52) The first term can be integrated by parts, and using the fact that v(0) =v(x f ) = 0 gives: Z t f 0 b(t)dt = Z x f 0 exp(αx) v dx (6.53) − Z x f 0 αv g exp(αx)dx. It is easy to see that maximizing v(x) will minimize both integrals, which gives the same result as the bang-bang strategy to minimize time (shown in figure 6.6). 105 6.4.3 Minimum energy strategy The problem of finding a minimum energy strategy is more complicated. We start with: Z t f 0 Hdt = Z t f 0 ( 1 2 v 2 +gx + bg exp(−αx) α )dt (6.54) = Z t f 0 ( 1 2 v 2 +gx + ˙ v +g α )dt (6.55) = Z x ∗ 0 ( 1 2 v 2 +gx + ˙ v +g α ) dx v (6.56) = Z x ∗ 0 ( 1 2 v 2 +gx + dv dx v +g α ) dx v (6.57) = Z x ∗ 0 ( 1 2 v 2 +gx + g α ) dx v + Z x ∗ 0 dv α (6.58) The last integal is simply v(x ∗ )−v(0) = 0, using the boundary conditions. So, Z t f 0 Hdt = Z x ∗ 0 f(v(x))dx, (6.59) f(v) = ( 1 2 v 2 +gx + g α )/v. (6.60) To minimize this integral, we consider each fixed, x, and minimize with respect to v(x), i.e. ∂f ∂v = 1 2 − gx + g α v 2 = 0 =⇒ v(x) = s 2gx + 2g α . (6.61) However, the value above may not be achievable for all x, e.g we have v(x = 0) = 0, so v(x) = q 2gx + 2g/α q 2gx + 2g/α≤v max v max q 2gx + 2g/α≥v max (6.62) According to figure 6.7, if the switching point (v s ,x s ) is on the right side of the curve v(x) = q 2gx + 2g α , then the extra transfer segment will exist connecting the b max trajectory with the b = 0 trajectory. This gives a sufficient condition: v s > s 2gx s + 2g α (6.63) 106 where x s = − ln(1− αx ∗ bmax ) α (6.64) v s = √ 2gx ∗ − 2gx s (6.65) After we simplify the equations: x s < x ∗ − 1 α 2 (6.66) Meanwhile, the inequality is the necessary condition as well because the equations: x s ≥ x ∗ − 1 α 2 (6.67) v max ≥ q 2gx + 2g/α (6.68) have no solution in the interval 0<x<x s . In summary eqn (6.66) is both a necessary and sufficient condition for the existence of the extra transfer segment shown in figure 6.7(a). 6.5 Pontryagin Maximum Principle The Pontryagin maximum (or minimum) principle (PMP) [91] is the most widely used method to find the optimal control scheme for certain dynamical systems. We will first give a statement of the PMP specifically for our problem and for proofs or more general versions, see reference [91, 66, 98, 113]. 107 Following the notation of [113],consider a general dynamical system with fixed intial condition X 0 , fixed final condition X f , and control b(t): ˙ X(t) = f(X(t),b(t)) (6.69) X(0) = X 0 (6.70) X(t f ) = X f (6.71) 0≤ b(t) ≤b max (6.72) And the target is to minimize some objective function: min b(t)∈[0,bmax] Z t f 0 L(X(t),b(t))dt. (6.73) First, we construct the control theory Hamiltonian function H by introducing the Lagrange multiplier vector: H (X(t),λ(t),b (t)) =L + T f (6.74) The dynamical system for this optimal control problem be same as eqn (6.3) and eqn (6.4) so we have: X = (x(t),v(t)) T (6.75) f = (v(t),−g +gb(t) exp(−αx(t))) T (6.76) Suppose = (λ x ,λ v ) T , the control theory Hamiltonian function for this problem will be: H =L +λ x v +λ v (−g +gb(t) exp(−αx)) (6.77) 108 Assuming that b ∗ (t) is our optimal control function, with the corresponding optimal trajectoryx ∗ (t),v ∗ (t), there exist functionsλ x ∗ ,λ v ∗ , which satisfy the canonical equations of Hamilton: ˙ x ∗ = ∂H (x ∗ ,v ∗ ,λ x ∗ ,λ v ∗ ,b ∗ ) ∂λ x ˙ v ∗ = ∂H (x ∗ ,v ∗ ,λ x ∗ ,λ v ∗ ,b ∗ ) ∂λ v ˙ λ ∗ x =− ∂H (x ∗ ,v ∗ ,λ x ∗ ,λ v ∗ ,b ∗ ) ∂x ˙ λ ∗ v =− ∂H (x ∗ ,v ∗ ,λ x ∗ ,λ v ∗ ,b ∗ ) ∂v . (6.78) The optimal controlb ∗ (t) will minimize the control theory Hamitonian H at any time point: b ∗ (t) = arg min b(t)∈[0,bmax] H (x ∗ ,v ∗ ,λ x ∗ ,λ v ∗ ,b(t)) (6.79) If we leave the final time t f as a free parameter, we need to impose one extra constraint on the final state: H (x ∗ (t f ),v ∗ (t f ),λ x ∗ (t f ),λ v ∗ (t f ),b ∗ (t f )) = 0. (6.80) Together with the initial condition in eqn (6.27) and final condition in eqn (6.28), the opti- mal control problem becomes two-point boundary value problem. For the three cases we considered in section 6.4, L = 1 minimize the time, L =b(t) minimizes the air volume and L =H minimize the energy respectively. It’s not hard to observe that for all the cases the control theory Hamitonian H is linear with the control b, so eqn (6.79) can be simplified based on the sign of ∂H ∂b : b ∗ (t) = 0 ∂H ∂b > 0 b max ∂H ∂b < 0 (6.81) This is the bang-bang schedule as calculated previously. The special case ∂H ∂b = 0 will be dealt with in the next section. If ∂H ∂b is non-zero almost everywhere, the optimal control b will be switching between 0 andb max , thus a bang-bang control. After solving the two-point 109 boundary value problem for L = 1 and L =b(t), the optimal bang-bang control is identical to section 6.4. 6.5.1 Singular control For the minimum energy case L =H, bang-bang control may not provide a solution to the two-point boundary value problem for the reason that ∂H ∂b = 0 will hold for at least some time interval in which case the Pontryagin minimum principle fails to yield the complete solution. For this, we need to implement what is called singular control [21]. The most simple method to solve this problem is to repeatly differentiate ∂H ∂b and set it to zero [21, 133]. First we know that ∂H ∂b = 0 (6.82) for some time interval, so set the first order and second order time derivatives of ∂H ∂b to 0 as well: d dt ∂H ∂b ! = 0 (6.83) d 2 dt 2 ∂H ∂b ! = 0. (6.84) Since the canonical equations (6.78) give the formulas of ˙ x ∗ , ˙ v ∗ , ˙ λ ∗ x , ˙ λ ∗ v , the time derivatives can be easily calculated with the chain rule. After solving the equations (6.82),(6.83) and (6.84), we have the following: λ ∗ v (t) = − 1 α (6.85) λ ∗ x (t) = −v ∗ (t) (6.86) b ∗ (t) = 2e αx ∗ (t) . (6.87) 110 Now we achieve the formulab(t). We will show that eqn (6.61) in section 6.4 (the transfer segment) is actually a singular arc. Take the time derivative of equation (6.61): d dt v(x) = d dt q 2gx + 2g/α = 2g 2 q 2gx + 2g/α dx dt = g q 2gx + 2g/α v =g (6.88) Together with eqn (6.1), the control b(t) for this segment will be: b(t) = ( ˙ v +g)e αx /g = 2e αx (6.89) which is identical to the singular arc in eqn (6.87) (shown as the transfer segment in figure 6.7(a)). 6.6 Discussion We have shown, using basic classical mechanics techniques along with orbit transfer ideas [60, 80], that in almost all parameter cases, it is a bang-bang control schedule that achieves the transfer of the ball to the levitation point x ∗ in minimum time, with minimum air volume, and expending minimum energy. There are, however, some parameter values for which a minimum energy transfer is only achieved with the use of an extra transfer segment that connects two bang-bang arcs. The transfer formula is obtained analytically and is an example of the use of singular control methods [21]. It might be of some interest to formulate similar optimal control problems associated with other forms of levitation, such as acoustic levitation or magnetic levitation. 111 It is worth noting that aside from choosing an appropriate time-dependent function b(t) to levitate the ball to the elliptic point, there are other mechanisms we could exploit. The simplest would be to add a small amount of dissipation to the model, proportional to the velocity: m¨ x =−mg +mgb(t) exp(−αx)−γ ˙ x. (6.90) Here γ is our dissipation parameter. We show in figure 6.8 a trajectory that ends up at x ∗ (as t→∞) after oscillating around it more and more tightly. 112 (a) (b) (c) (d) Figure 6.2: (a) Phase curves in the (x,v) plane for 5 different values of the energy E b . x max marks the maximum height that the ball can achieve on the energy curveE b =E 0 b , while the elliptic point x ∗ marks the levitation point where the forces balance E b = E ∗ b . (b) Energy curves E b = E 0 b going through the origin for four different values of b. (c) Phase curves of free fall (b = 0) with different values of energy. (d) x max marks the intersection of the two curves defined by the left and right hand sides of eqn (6.8). 113 Figure 6.3: Diagram showing two different paths from point O to point x ∗ . Orbits corre- sponding to energy E b = E 0 b for two different values, b = b 1 (blue) and b = b 2 < b 1 (red). Heightx ∗ represents the elliptic fixed point associated withb =b 1 . The value ofb 2 is chosen so that x max for the red orbit intersects x ∗ . The pink orbit corresponds to a free-fall tra- jectory b = 0 with an initial upward launch velocity v ∗ so that the trajectory intersects x ∗ . PointA, with coordinates (x A ,v A ), marks the intersection of the free-fall trajectory with the blue orbit. 114 (a) (b) (c) (d) (e) (f) (g) (h) Figure 6.4: Examples of piecewise constant blowing schedules that move the ball from the origin to the target point (x,v) = (x ∗ , 0). Parameter values are given by: b∈ [0, 6],α = 1m −1 ,g = 9.8m/s 2 ,x ∗ = log(2)meter,m = 1kg, and we choose T = 1.34s for time scale. (a) Two-step schedule: b(t) starts with b = αx ∗ e αx ∗ e αx ∗ −1 , then it switches to b = e αx ∗ . (b) The corresponding trajectories for the two-step schedule. Note that the ball will remain at equilibrium on second step, so the trajectory only has one segment. (c) Three-step schedule with blowing on the first segment, followed by free falling (b = 0) for the second segment. (d) The corresponding trajectory for (c). The trajectories will get closer and closer to vertical axis with the increase of initialb(t). (e) Examples of general three-step schedule, withb 1 for first segment and b 2 for the second. (f) The corresponding b(t) for the examples of general three step schedule. (g) A schematic of a four-step schedule, where t 1 t 2 and t 3 mark three switching points.(h) The corresponding trajectory for (g) 115 (a) (b) Figure 6.5: The schematic of the trajectory driven by a general blowing function b(t) as indicated in blue. For comparison, we show what we call the bang-bang strategy, starting with maximum blowing b max from t = 0 to t =t s , then free fall b = 0 from t =t s to t =t f . Any general trajectory associated withb(t) is trapped between the two segments of the other orbit. (a) (b) Figure 6.6: The bang-bang (on-off) strategy is optimal to attain the transfer in minimum time and using minimum air volume. (a) The actual trajectory for the transfer; (b) The time-dependent blowing schedule. 116 (a) (b) (c) (d) Figure 6.7: Minimum energy strategy.(a)phase space. Simulated with parameters b max = 20,α = 1m −1 ,x ∗ = ln(5)m,g = 9.8m/s 2 ,m = 1kg (b)corresponding schedule for (a) .(c)phase space simulated with parameters: b max = 5,α = 1m −1 ,x ∗ = ln(5)m,g = 9.8m/s 2 ,m = 1kg.(d) corresponding optimal control and it’s still bang-bang control. 117 Figure 6.8: Constant levitation with dissipation.The numbers for the simulation are b = 5,α = 1m −1 ,g = 9.8m/s 2 ,x ∗ =log(5)m,m = 1kg,γ = 1kg/s 118 Reference List [1] CA Aktipis, VSY Kwan, KA Johnson, SL Neuberg, and CC Maley. Overlooking evolution: a systematic analysis of cancer relapse and therapeutic resistance research. PloS One, 6(11):e26100, 2011. [2] Eitan Altman, Rachid ElAzouzi, Yezekeael Hayel, and Hamidou Tembine. An evolu- tionary game approach for the design of congestion control protocols in wireless net- works. In 2008 6th International Symposium on Modeling and Optimization in Mobile, Ad Hoc, and Wireless Networks and Workshops, pages 547–552. IEEE, 2008. [3] D Andersson, NQ Balaban, F Baquero, P Courvalin, P Glaser, U Gophna, R Kishony, S Molin, and T Tonjum. Antibiotic resistance: turning evolutionary principles into clinical reality. FEMS Microbiology Rev., 2020. [4] Dan I Andersson and Bruce R Levin. The biological cost of antibiotic resistance. Current opinion in microbiology, 2(5):489–493, 1999. [5] DI Andersson and BR Levin. The biological cost of antibiotic resistance. Current Opinion in Microbiology, 2:489–493, 1999. [6] Camille Stephan-Otto Attolini and Franziska Michor. Evolutionary theory of cancer. Annals of the New York Academy of Sciences, 1168(1):23–51, 2009. [7] Robert Axelrod, David E Axelrod, and Kenneth J Pienta. Evolution of cooperation among tumor cells. Proc. Natl. Acad. Sci., 103(36):13474–13479, 2006. [8] Katarina Bacevic, Robert Noble, Ahmed Soffar, Orchid Wael Ammar, Benjamin Boszonyik, Susana Prieto, Charles Vincent, Michael E Hochberg, Liliana Krasinska, and Daniel Fisher. Spatial competition constrains resistance to targeted cancer ther- apy. Nature communications, 8(1):1995, 2017. [9] D Basanta and ARA Anderson. Exploiting ecological principles to better understand cancer progression and treatment. Interface Focus, 3(4):20130020, 2013. [10] D Basanta, RA Gatenby, and ARA Anderson. Exploiting evolution to treat drug resis- tance: combination therapy and the double bind. Molecular pharmaceutics, 9(4):914– 921, 2012. [11] M Baym, LK Stone, and R Kishony. Multidrug evolutionary strategies to reverse antibiotic resistance. Science, 351(6268):DOI: 10.1126/science.aad3292, 2016. 119 [12] RA Beckman, GS Schemmarm, and CH Yeang. Impact of genetic dynamics and single- cell heterogeneity on the development of personalized non-standard medicine strategies for cancer. Proc. Natl. Acad. Sci., 109(36):14586–14591, 2012. [13] Sharon Bewick, Ruoting Yang, and Mingjun Zhang. Embedding evolutionary game theory into an optimal control framework for drug dosage design. In 2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pages 6026–6029. IEEE, 2009. [14] CI Bliss. The toxicity of poisons applied jointly. Annals of Applied Biology, 26:585–615, 1939. [15] K Bohl, S Hummert, S Werner, D Basanta, A Deutsch, S Schuster, G Theiben, and A Schroeter. Evolutionary game theory: molecules as players. Molecular BioSystems, 10(12):3044–3065, 2014. [16] I Bozic and MA Nowak. Timing and heterogeneity of mutations associated with drug resistance in metastatic cancers. Proc. Natl. Acad. Sci., 111(45):15964–15968, 2014. [17] I Bozic and MA Nowak. Resisting resistance. Ann. Rev. of Can. Bio., 1:203–221, 2017. [18] K Bratton, U Ledzewicz, and H Schattler. Modeling and control of heterogeneous tumors under chemotherapy. Biomath Comm., 1(1):1–2, 2014. [19] JS Brown and K Stankova. Game theory as a conceptual framework for managing insect pests. Current Opinion in Insect Sci., 21:26–32, 2017. [20] William L Brown and Edward O Wilson. Character displacement. Systematic Zoology, 5(2):49–64, 1956. [21] Arthur Earl Bryson. Applied Optimal Control: Optimization, Estimation and Control. Routledge, 2018. [22] Cécile Carrere. Optimization of an in vitro chemotherapy to avoid resistant tumours. Journal of Theoretical Biology, 413:24–33, 2017. [23] Gustavo F Carvalhal, Saima N Daudi, Donghui Kan, Dana Mondo, Kimberly A Roehl, Stacy Loeb, and William J Catalona. Correlation between serum prostate-specific antigen and cancer volume in prostate glands of different sizes. Urology, 76(5):1072– 1076, 2010. [24] J Chacon, J Saenz, L de la Torre, JM Diaz, and F Esquembre. Design of a low-cost air levitation system for teaching control engineering. Sensors, 17(2321), 2017. [25] TC Chou and DC Rideout. Synergism and Antagonism in Chemotherapy. Academic Press, NY, 1991. [26] Andrew J Coldman and JM Murray. Optimal control for a stochastic model of cancer chemotherapy. Mathematical biosciences, 168(2):187–200, 2000. 120 [27] Joseph H Connell. The influence of interspecific competition and other factors on the distribution of the barnacle chthamalus stellatus. Ecology, 42(4):710–723, 1961. [28] Charles Darwin and William F Bynum. The origin of species by means of natural selection: or, the preservation of favored races in the struggle for life. Penguin Har- mondsworth, 2009. [29] Luis A Diaz Jr, Richard T Williams, Jian Wu, Isaac Kinde, J Randolph Hecht, Jordan Berlin, Benjamin Allen, Ivana Bozic, Johannes G Reiter, Martin A Nowak, et al. The molecular evolution of acquired resistance to targeted egfr blockade in colorectal cancers. Nature, 486(7404):537–540, 2012. [30] Li Ding, Timothy J Ley, David E Larson, Christopher A Miller, Daniel C Koboldt, John S Welch, Julie K Ritchey, Margaret A Young, Tamara Lamprecht, Michael D McLellan, et al. Clonal evolution in relapsed acute myeloid leukaemia revealed by whole-genome sequencing. Nature, 481(7382):506–510, 2012. [31] Pedro M Enriquez-Navas, Yoonseok Kam, Tuhin Das, Sabrina Hassan, Ariosto Silva, Parastou Foroutan, Epifanio Ruiz, Gary Martinez, Susan Minton, Robert J Gillies, et al. Exploiting evolutionary principles to prolong tumor control in preclinical models of breast cancer. Science Translational Medicine, 8(327):327ra24–327ra24, 2016. [32] PM Enriquez-Navas, JW Wojtkowiak, and RA Gatenby. Application of evolutionary principles to cancer therapy. Cancer research, 75(22):4675–4680, 2015. [33] A Fisher, I Vazquez-Garcia, and V Mustonen. The value of monitoring to control evolving populations. Proc. Natl. Acad. Sci., 112(4):1007–1012, 2015. [34] J Foo and Michor F. Evolution of resistance to anti-cancer therapy during general dosing schedules. J. Theor. Bio., 263(2):179–188, 2010. [35] J Foo and Michor F. Evolution of acquired resistance to anti-cancer therapy. J. Theor. Bio., 355:10–20, 2014. [36] Jill Gallaher, Pedro M Enriquez-Navas, Kimberly A Luddy, Robert A Gatenby, and Alexander RA Anderson. Spatial heterogeneity and evolutionary dynamics modulate timetorecurrenceincontinuousandadaptivecancertherapies. Cancer Research, 2018. [37] Jill A Gallaher, Pedro M Enriquez-Navas, Kimberly A Luddy, Robert A Gatenby, and Alexander RA Anderson. Adaptive vs continuous cancer therapy: Exploiting space and trade-offs in drug scheduling. 2017. [38] RA Gatenby. A change of strategy in the war on cancer. Nature, 459(7246):508–509, 2009. [39] RA Gatenby and J Brown. The evolution and ecology of resistance in cancer therapy. Cold Spring Harbor Perspectives in Medicine, page doi: 10.1101/cshperspect.a033415, 2017. 121 [40] RAGatenby,ASSilva,RJGillies,andBRFrieden. Adaptivetherapy. Cancer Research, 69(11):4894–4903, 2009. [41] Robert A Gatenby, Robert J Gillies, and Joel S Brown. Of cancer and cave fish. Nature Reviews Cancer, 11(4):237, 2011. [42] RJ Gillies, D Verduzco, and RA Gatenby. Evolutionary dynamics of carcinogenesis and why targeted therapy does not work. Nature Reviews Cancer, 12(7):487–493, 2012. [43] GGomez, WSKoon, MWLo, JEMarsden, andSDRoss. Connectingorbitsandinvari- ant manifolds in the spatial restricted three-body problem. Nonlinearity, 17(5):1571, 2004. [44] F Gould. The evolutionary potential of crop pests. American Scientist, 79:496–507, 1991. [45] Peter R Grant. Convergent and divergent character displacement. Biological journal of the Linnean Society, 4(1):39–68, 1972. [46] Mel Greaves and Carlo C Maley. Clonal evolution in cancer. Nature, 481(7381):306– 313, 2012. [47] P Grover and SD Ross. Designing trajectories in a planet-moon environment using the controlled keplerian map. J. of Guidance, Control, Dyn., 32(2):436–443, 2009. [48] Erik Gullberg, Sha Cao, Otto G Berg, Carolina Ilbäck, Linus Sandegren, Diarmaid Hughes, and Dan I Andersson. Selection of resistant bacteria at very low antibiotic concentrations. PLoS pathogens, 7(7):e1002158, 2011. [49] Elsa Hansen, Robert J Woods, and Andrew F Read. How to use a chemotherapeutic agent when resistance to it threatens the patient. PLoS biology, 15(2):e2001110, 2017. [50] M Hegreness, N Shoresh, D Damiann, D Hartl, and R Kishony. Accelerated evolution of resistance in multidrug environments. Proc. Nat’l Acad. Sci., 105(37):13977–13981, 2008. [51] Keith S Hoek and Colin R Goding. Cancer stem cells versus phenotype-switching in melanoma. Pigment cell & melanoma research, 23(6):746–759, 2010. [52] Josef Hofbauer and Karl Sigmund. Evolutionary Games and Population Dynamics. Cambridge University Press, 1998. [53] S Hummert, K Bohl, D Basanta, A Deutsch, S Werner, G Theiben, A Schroeter, and S Schuster. Evolutionary game theory: cells as players. Molecular BioSystems, 10(12):3044–3065, 2014. [54] G Jamsen, R Gatenby, and CA Aktipis. Opinion: Control vs. eradication: Applying infectious disease strategies to cancer. Proc. Natl. Acad. Sci., 112(4):937–938, 2015. 122 [55] Hem Raj Joshi. Optimal control of an hiv immunology model. Optimal control appli- cations and methods, 23(4):199–213, 2002. [56] Kristel Kemper, Oscar Krijgsman, Paulien Cornelissen-Steijger, Aida Shahrabi, Fleur Weeber, Ji-Ying Song, Thomas Kuilman, Daniel J Vis, Lodewyk F Wessels, Emile E Voest, et al. Intra-and inter-tumor heterogeneity in a vemurafenib-resistant melanoma patient and derived xenografts. EMBO Molecular Medicine, 7(9):1104–1118, 2015. [57] S Kim, TD Lieberman, and R Kishony. Alternating antibiotic treatments constrain evolutionary paths to multidrug resistance. Proc. Nat’l Acad. Sci., 40:14494–14499, 2014. [58] NL Komarova and D Wodarz. Drug resistance in cancer: Principles of emergence and prevention. Proc. Natl. Acad. Sci., 102(27):9714–9719, 2005. [59] S Konishi, M Harada, Y Ogami, Y Daiho, Y Mita, and H Fujita. Measurements and analysis on characteristics of a conveyance system using air levitation. IEEE 6th International Conference on Emerging Technologies and Factory Automation, page DOI: 10.1109/ETFA.1997.616274, 1997. [60] WS Koon, MW Lo, JE Marsden, and SD Ross. Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics. Chaos, 10(427):doi:10.1063/1.166509, 2000. [61] K-A Kreuzer, P Le Coutre, O Landt, I-K Na, M Schwarz, K Schultheis, A Hochhaus, and B Dörken. Preexistence and evolution of imatinib mesylate-resistant clones in chronic myelogenous leukemia detected by a pna-based pcr clamping technique. Annals of Hematology, 82(5):284–289, 2003. [62] Pierre Laurent-Puig, Deniz Pekin, Corinne Normand, Steve K Kotsopoulos, Philippe Nizard, Karla Perez-Toralla, Rachel Rowell, Jeff Olson, Preethi Srinivasan, Delphine Le Corre, et al. Clinical relevance of kras-mutated subclones detected with picodroplet digitalpcrinadvancedcolorectalcancertreatedwithanti-egfrtherapy. Clinical Cancer Research, 2014. [63] O Lavi, MM Gottesman, and D Levy. The dynamics of drug resistance: A mathemat- ical perspective. Drug Res. Updates, 15(1-2):90–97, 2012. [64] U Ledzewicz and HM Schattler. Optimal bang-bang controls for a two-compartmental model in cancer therapy. J. Optimization Theory and Appl., 114(3):609–637, 2002. [65] U Ledzewicz, S Wang, H Schattler, N Andre, MA Heng, and E Pasquier. On drug resis- tance and metronomic chemotherapy: A mathematical modeling and optimal control approach. Math. Biosciences and Eng., 14(1):217–235, 2017. [66] Ernest Bruce Lee and Lawrence Markus. Foundations of optimal control theory. Tech- nical report, Minnesota Univ Minneapolis Center For Control Sciences, 1967. [67] S Lev. Pontryagin. mathematical theory of optimal processes, 1987. 123 [68] RC Lewontin. The units of selection. Annual Reviews of Ecol. Syst., 1:1–18, 1970. [69] Erez Lieberman, Christoph Hauert, and Martin A Nowak. Evolutionary dynamics on graphs. Nature, 433(7023):312–316, 2005. [70] M Lipsitch and BR Levin. The population dynamics of antimicrobial chemotherapy. Antimicrobial Agents and Chemotherapy, 41(2):363–373, 1997. [71] S Loewe. The problem of synergism and antagonism of combined drugs. Arzneim. Forsch., 3:285–290, 1953. [72] Alfred J Lotka. Contribution to the theory of periodic reactions. The Journal of Physical Chemistry, 14(3):271–274, 2002. [73] Rory B Martin. Optimal control drug scheduling of cancer chemotherapy. Automatica, 28(6):1113–1123, 1992. [74] Lauren MF Merlo, John W Pepper, Brian J Reid, and Carlo C Maley. Cancer as an evolutionary and ecological process. Nature Reviews Cancer, 6(12):924–935, 2006. [75] M Dror Michaelson, Stephane Oudard, Yen-Chuan Ou, Lisa Sengeløv, Fred Saad, Nadine Houede, Peter Ostler, Arnulf Stenzl, Gedske Daugaard, Robert Jones, et al. Randomized, placebo-controlled, phaseiiitrialofsunitinibplusprednisoneversuspred- nisone alone in progressive, metastatic, castration-resistant prostate cancer. Journal of Clinical Oncology, 32(2):76–82, 2013. [76] JBMichel, PYeh, RChait, RCMoellering, andRKishony. Druginteractionsmodulate the potential for evolution of resistance. Proc. Nat’l Acad. Sci., 105(39):14918–14923, 2008. [77] A Sorana Morrissy, Livia Garzia, David JH Shih, Scott Zuyderduyn, Xi Huang, Patryk Skowron, Marc Remke, Florence MG Cavalli, Vijay Ramaswamy, Patricia E Lindsay, et al. Divergent clonal selection dominates medulloblastoma at recurrence. Nature, 529(7586):351–357, 2016. [78] PK Newton and Y Ma. Nonlinear adaptive control of competitive release and chemotherapeutic resistance. Phys. Rev. E, 99(022404), 2019. [79] PK Newton and Y Ma. Nonlinear adaptive control of competitive release and chemotherapeutic resistance. Physical Review E, 99(2):022404, 2019. [80] PK Newton and Y Ma. Nonlinear adaptive control of competitive release and chemotherapeutic resistance. Phys. Rev. E, 99(022404), 2019. [81] D Nichol, P Jeavons, A G Fletcher, R A Bonomo, P K Maini, J L Paul, R A Gatenby, A R A Anderson, and J G Scott. Steering evolution with sequential therapy to prevent the emergence of bacterial antibiotic resistance. PLoS Comput Biol, 11(9):e1004493, 2015. 124 [82] D Nichol, P Jeavons, AG Fletcher, RA Bonomo, PK Maini, JL Paul, RA Gatenby, ARA Anderson, and JG Scott. Steering evolution with sequential therapy to prevent the emergence of bacterial antibiotic resistance. PLoS Comp. Bio., 11(9), 2015. [83] L Norton. A gompertzian model of human breast cancer growth. Can. Res., 48(24 Part 1):7067–7071, 1988. [84] L Norton and R Simon. Tumor size, sensitivity to therapy, and design of treatment schedules. Cancer Treatment Reports, 61(7):1307, 1977. [85] Martin A Nowak. Evolutionary Dynamics. Harvard University Press, 2006. [86] Peter C Nowell. The clonal evolution of tumor cell populations. Science, 194(4260):23– 28, 1976. [87] Martin J Osborne and Ariel Rubinstein. A course in game theory. MIT press, 1994. [88] MC Perry. The Chemotherapy Source Book. Lippincott Williams & Wilkins, 2008. [89] Daniel P Petrylak, Nicholas J Vogelzang, Nikolay Budnik, Pawel Jan Wiechno, Cora N Sternberg, Kevin Doner, Joaquim Bellmunt, John M Burke, Maria Ochoa de Olza, Ananya Choudhury, et al. Docetaxel and prednisone with or without lenalidomide in chemotherapy-naive patients with metastatic castration-resistant prostate cancer (mainsail): a randomised, double-blind, placebo-controlled phase 3 trial. The Lancet Oncology, 16(4):417–425, 2015. [90] Angela Oliveira Pisco, Amy Brock, Joseph Zhou, Andreas Moor, Mitra Mojtahedi, Dean Jackson, and Sui Huang. Non-darwinian dynamics in therapy-induced cancer drug resistance. Nature communications, 4:2467, 2013. [91] Lev Semenovich Pontryagin. Mathematical Theory of Optimal Processes. Routledge, 2018. [92] JE Prussing and BA Conway. Orbital Mechanics, 2nd Ed. Oxford University Press, 2013. [93] S Radu, WS Koon, MW Lo, JE Marsden, SD Ross, and RS Wilson. Halo orbit mission correction maneuvers using optimal control. Automatica, 38:571–583, 2002. [94] AF Read, T Day, and S Huijben. The evolution of drug resistance and the curious orthodoxy of aggressive chemotherapy. Proc. Natl. Acad. Sci., 108(2):10871–10877, 2011. [95] NicholasPRestifo,MarkEDudley,andStevenARosenberg. Adoptiveimmunotherapy for cancer: harnessing the t cell response. Nature Reviews Immunology, 12:269–281, 2012. [96] Catherine Roche-Lestienne and Claude Preudhomme. Mutations in the abl kinase domain pre-exist the onset of imatinib treatment. In Seminars in Hematology, vol- ume 40, pages 80–82. Elsevier, 2003. 125 [97] Alessandro Romanel, Delila Gasi Tandefelt, Vincenza Conteduca, Anuradha Jayaram, Nicola Casiraghi, Daniel Wetterskog, Samanta Salvi, Dino Amadori, Zafeiris Zafeiriou, Pasquale Rescigno, et al. Plasma ar and abiraterone-resistant prostate cancer. Science Translational Medicine, 7(312):312re10–312re10, 2015. [98] IM Ross. A Primer on Pontryagin’s Principle in Optimal Control, 2nd Ed. Collegiate Press, 2015. [99] SDRossandDJScheeres. Multiplegravityassists, capture, andescapeintherestricted three-body problem. SIAM J. Appl. Dyn. Sys., 6(3):576–596, 2007. [100] D Russ and R Kishony. Additivity of inhibitory effects in multidrug combinations. Nature Microbiology, 3:1339–1345, 2018. [101] Bean San Goh, George Leitmann, and Thomas L Vincent. Optimal control of a prey- predator system. Mathematical Biosciences, 19(3-4):263–286, 1974. [102] WH Sandholm, E Dokumacı, and F Franchetti. Dynamo: diagrams for evolutionary game dynamics. software, 2012. [103] RolandFSchwarz, CharlotteKYNg, SusannaLCooke, ScottNewman, JillianTemple, Anna M Piskorz, Davina Gale, Karen Sayal, Muhammed Murtaza, Peter J Baldwin, et al. Spatial and temporal heterogeneity in high-grade serous ovarian cancer: A phylogenetic analysis. PLoS Med, 12(2):e1001789, 2015. [104] H Segre, N DeMalach, Z Henkin, and R Kadmon. Quantifying competitive exclusion and competitive release in ecological communities: A conceptual framework and a case study. PloS One, 11(8):e0160798, 2016. [105] S Seton-Rogers. Chemotherapy: Preventing competitive release. Nature Rev. Can., 16(4):199–199, 2016. [106] Sreenath V Sharma, Diana Y Lee, Bihua Li, Margaret P Quinlan, Fumiyuki Taka- hashi, Shyamala Maheswaran, Ultan McDermott, Nancy Azizian, Lee Zou, Michael A Fischbach, et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell, 141(1):69–80, 2010. [107] AS Silva and RA Gatenby. A theoretical quantitative model for evolution of cancer chemotherapy resistance. Biology Direct, 5(25):1–25, 2010. [108] J Maynard Smith. Game theory and the evolution of fighting. On evolution, pages 8–28, 1972. [109] SH Strogatz. Nonlinear Dynamics and Chaos. Westview Press, 2015, 2nd Ed. [110] P Sudalagunta, MC Silva, RR Canevarolo, R Gatenby, Gillies R, R Baz, MB Meads, KH Shain, and AS Silva. A pharmacodynamic model of clinical synergy in multiple myeloma. EBioMedicine, 54:102716, 2020. 126 [111] George W Swan. Role of optimal control theory in cancer chemotherapy. Mathematical biosciences, 101(2):237–284, 1990. [112] Ian F Tannock, Ronald de Wit, William R Berry, Jozsef Horti, Anna Pluzanska, KimNChi, StephaneOudard, ChristineThéodore, NicholasDJames, IngelaTuresson, et al. Docetaxel plus prednisone or mitoxantrone plus prednisone for advanced prostate cancer. New England Journal of Medicine, 351(15):1502–1512, 2004. [113] Kok Lay Teo, C Goh, and K Wong. A unified computational approach to optimal control problems. 1991. [114] Arne Traulsen, Jens Christian Claussen, and Christoph Hauert. Coevolutionary dynamics: From finite to infinite populations. Physical Review Letters, 95(23):238701, 2005. [115] Arne Traulsen, Jens Christian Claussen, and Christoph Hauert. Coevolutionary dynamics in large, but finite populations. Physical Review E, 74(1):011901, 2006. [116] S Venkatesan and C Swanton. Tumor evolutionary principles: How intratumor het- erogeneity influences cancer treatment and outcome. In American Society of Clinical Oncology educational book. American Society of Clinical Oncology. Meeting, volume 35, page e141, 2016. [117] BWalshandMWBlows. Abundantgeneticvariation+strongselection=multivariate genetic constraints: A geometric view of adaptation. Annu. Rev. Ecol. Evol. Syst., 40:41–59, 2009. [118] J West, Z Hasnain, P Macklin, and PK Newton. An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity and gompertzian growth. SIAM Review, 58(4):716–736, 2016. [119] J West, Z Hasnain, J Mason, and PK Newton. The prisoner’s dilemma as a cancer model. Converg. Sci. Phys. Oncol., 2(3):035002, 2016. [120] J West, Y Ma, and PK Newton. Capitalizing on competition: An evolutionary model of competitive release in metastatic castration resistant prostate cancer treatment. J. Theor. Bio., 455:249–260, 2018. [121] J West, Y Ma, and PK Newton. Towards multi-drug adaptive therapy. Cancer Research, 80(7):1578–1589, 2020. [122] J West and PK Newton. Chemotherapeutic dose scheduling based on tumor growth ratesprovidesacaseforlow-dosemetronomichigh-entropytherapies. Cancer Research, pages doi: 10.1158/0008–5472.CAN–17–1120, 2017. [123] J West and PK Newton. Cellular interactions constrain tumor growth. Proc. Nat’l Acad. Sci., 116(6):1918–1923, 2019. 127 [124] Jeffrey West, Zaki Hasnain, Paul Macklin, and Paul K Newton. An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity driving gompertzian growth. SIAM Review, 58(4):716–736, 2016. [125] Jeffrey West, Yongqian Ma, and Paul K Newton. Capitalizing on competition: An evolutionary model of competitive release in metastatic castration resistant prostate cancer treatment. Journal of theoretical biology, 455:249–260, 2018. [126] Jeffrey West and Paul K Newton. Chemotherapeutic dose scheduling based on tumor growth rates provides a case for low-dose metronomic high-entropy therapies. Cancer research, 77(23):6717–6728, 2017. [127] Julia Wilkerson, Kald Abdallah, Charles Hugh-Jones, Greg Curt, Mace Rothenberg, Ronit Simantov, Martin Murphy, Joseph Morrell, Joel Beetsch, Daniel J Sargent, et al. Estimation of tumour regression and growth rates during treatment in patients with advanced prostate cancer: a retrospective analysis. The Lancet Oncology, 18(1):143– 154, 2017. [128] KB Wood, KC Wood, S Nishida, and P Cluzel. Uncovering scaling laws to infer multidrug response of resistant microbes and cancer cells. Cell Reports, 6:1073–1084, 2014. [129] B Wu, PM Altrock, L Wang, and A Traulsen. Universality of weak selection. Phys. Rev. E, 82:046106, 2010. [130] N Yoon, R Vander Velde, A Marusyk, and J G Scott. Optimal therapy scheduling based on a pair of collaterally sensitive drugs. Bull. Math. Bio., pages 1–34, 2018. [131] N Yoon, R Vander Velde, A Marusyk, and JG Scott. Optimal therapy scheduling based on a pair of collaterally sensitive drugs. Bull. Math. Bio., 80(7):1776–1809, 2018. [132] L You, JS Brown, F Thuijsman, JJ Cunningham, RA Gatenby, J Zhang, and K Stankova. Spatial vs. non-spatial eco-evolutionary dynamics in a tumor growth model. J. Theor. Bio., 435:78–97, 2017. [133] MI Zelikin and VF Borisov. Singular optimal regimes in problems of mathematical economics. Journal of Mathematical Sciences, 130(1):4409–4570, 2005. [134] J Zhang, JJ Cunningham, JS Brown, and RA Gatenby. Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer. Nature Comm., 8(1):1816, 2017. 128 Appendix A replicator dynamics of tumor growth Tounderstandwhytheprisoner’sdilemmaisausefulparadigmfortumorgrowthresulting from competing subpopulations, we focus on the 2× 2 case only. The standard version of the prisoner’s dilemma payoff matrix [85] in a 2× 2 setting in which healthy cells compete with cancer cells is: A = a b c d = 3 0 5 1 , (A.1) where the first row and column correspond to the payoffs associated with the cooperator (C) in the PD evolutionary game, and the second row and column correspond to the payoffs associated with the defector (D). In the simplest tumor growth paradigm in which a popu- lation of healthy cells competes with a population of cancer cells, the healthy cells are the cooperators, while the cancer cells are the defectors. In any mixed population~ x = (x C ,x D ) T , 0≤x C ≤ 1; 0≤x D ≤ 1; x C +x D = 1, the fitness functions, ~ f = (f C ,f D ) T , associated with the two subpopulations are: ~ f = A~ x, (A.2) 129 which in component form yields: f C = (A~ x) 1 = 3·x C + 0·x D , (A.3) f D = (A~ x) 2 = 5·x C + 1·x D , (A.4) while the average fitness of the total population is given by the quadratic form: hfi = ~ x T A~ x = 3x 2 C + 5x C x D +x 2 D ≥ 1. (A.5) The average fitness of the healthy state (x C ,x D ) = (1, 0) is given byhfi| (x C =1) = 3, while that of the cancerous state (x C ,x D ) = (0, 1) is given byhfi| (x D =1) = 1, which minimizes the average fitness. Tumor growth is then modeled as a 2× 2 evolutionary game governed by the replicator dynamical system: ˙ x C = (f C −hfi)x C , (A.6) ˙ x D = (f D −hfi)x D . (A.7) It is straightforward to show: ˙ x D = [(c−a)− (d−b)]x D (1−x D ) 1 1− d−b c−a −x D , (A.8) with fixed points at x D = 0, 1, (c−a) (c−a)−(d−b) . From this, we can conclude (using values from eqn (A.1)) that for any initial condition containing at least one cancer cell: 0<x D (0)≤ 1, we have: (i) x D → 1, x C → 0 as t→∞ (ii)hfi→ 1 as t→∞. Condition (i) guarantees that the cancer cell population will saturate, while condition (ii) guarantees that the saturated state is sub-optimal, sincehfi| (x D =1) <hfi| (x C =1) . For these 130 tworeasons, theprisoner’s dilemma evolutionary game serves as a simple paradigmfor tumor growth both in finite population models, as well as replicator system (infinite population) models [118, 119, 122, 114]. Note that the 2× 2 system alone it not able to account for the evolution of resistance. 131
Abstract (if available)
Abstract
In this thesis, we use a three by three replicator dynamical system with a payoff matrix of Prisoner’s Dilemma type to develop evolutionary game theory models of adaptive chemotherapy that avoid chemotherapeutic resistance in cancer. With the primary goal of controlling resistance and delaying tumor recurrence, we introduce and develop a trajectory design method from the orbital mechanics literature in order to construct families of closed loops, called ‘evolutionary cycles’, using adaptive time-dependent dose schedules. With time-dependent adaptive chemotherapy control, we attempt to alter the fitness landscape of the heterogeneous cell populations in order to delay tumor recurrence. We show that for a three-component population made up of healthy cells, sensitive cells, and resistant cells, this method outcompetes both maximum tolerated dose (MTD) schedules and low-dose metronomic (LDM) schedules, which are standard currently used protocols. We then consider a more complex system made up of two types of resistant cells and two drugs applied sequentially and in combination. In this case, the drugs can interact synergistically, additively, or antagonistically. We show how to develop multi-drug schedules in this context using a similar method of adaptive trajectory design. Generally speaking, we show that antagonistic drug interactions are more effective at delaying tumor recurrence than synergistic ones which effectively cause the two drugs to act as one, limiting the flexibility in schedule design. We then apply the Pontryagin Maximum Principle to two by two evolutionary games. We show examples of payoff matrices where bang-bang control is optimal, and others where singular control techniques need to be applied. As a final chapter, we use a Hamiltonian model along with the trajectory design method to show how to levitate a mass with specially designed time-dependent blowing schedules.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response
PDF
Stochastic multidrug adaptive chemotherapy control of competitive release in tumors
PDF
Integrative mathematical and computational approaches to investigate metastatic cancers
PDF
Reinforcement learning based design of chemotherapy schedules for avoiding chemo-resistance
PDF
Adaptive agents on evolving network: an evolutionary game theory approach
PDF
Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations
PDF
Neural network for molecular dynamics simulation and design of 2D materials
PDF
Stochastic and multiscale models for urban and natural ecology
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Multi-scale quantum dynamics and machine learning for next generation energy applications
PDF
Structure – dynamics – function analysis of class A GPCRs and exploration of chemical space using integrative computational approaches
PDF
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
Evolutionary genomic analysis in heterogeneous populations of non-model and model organisms
PDF
The role of translocator protein (TSPO) in mediating the effects of mitotane in human adrenocortical carcinoma cells
PDF
Novel techniques for analysis and control of traffic flow in urban traffic networks
PDF
Model based design of porous and patterned surfaces for passive turbulence control
PDF
Nonlinear dynamical modeling of single neurons and its application to analysis of long-term potentiation (LTP)
PDF
Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
Asset Metadata
Creator
Ma, Yongqian
(author)
Core Title
Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
07/31/2020
Defense Date
06/08/2020
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
cancer modeling,chemo-therapy scheduling,evolutionary games,nonlinear dynamics,OAI-PMH Harvest,optimal control,replicator dynamics
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Haas, Stephan (
committee chair
), Haselwandter, Christoph (
committee member
), Nakano, Aiichiro (
committee member
), Newton, Paul (
committee member
), Pahlevan, Niema (
committee member
)
Creator Email
myqian0202@gmail.com,yongqiam@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-353470
Unique identifier
UC11665935
Identifier
etd-MaYongqian-8836.pdf (filename),usctheses-c89-353470 (legacy record id)
Legacy Identifier
etd-MaYongqian-8836.pdf
Dmrecord
353470
Document Type
Dissertation
Rights
Ma, Yongqian
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
cancer modeling
chemo-therapy scheduling
evolutionary games
nonlinear dynamics
optimal control
replicator dynamics