University of Southern California Dissertations and Theses

Stochastic multidrug adaptive chemotherapy control of competitive release in tumors

(USC Thesis Other)

Stochastic multidrug adaptive chemotherapy control of competitive release in tumors

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content STOCHASTIC MULTI DRUG ADAPTIVE CHEMOTHERAPY CONTROL OF COMPETITIVE RELEASE IN TUMORS by Jiyeon Park A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Applied Mathematics) May 2022 Copyright 2022 Jiyeon Park TableofContents ListofTables iv ListofFigures v Abstract xvi Chapter1: Introduction 1 1.1 Evolutionary game theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 (Adjusted) Replicator dynamics in an innite population . . . . . . . . . . . . . . . 6 1.1.2 Stochastic models in a nite population . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2 Moran process in a nite population . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.1 Moran process: A birth-and-death process . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.2 Connection of Moran process to the adjusted replicator dynamics . . . . . . . . . . 15 1.3 Moran process and a chemotherapy concentration applied to the selection . . . . . . . . . 18 1.3.1 Prisoner’s dilemma game . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.3.2 Cost of resistance and competitive release . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.3 Fitness landscape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.4 Structure of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Chapter2: FixationprobabilityforMoranprocesseswithtwostrategies 29 2.1 Fixation probability and evolutionary stability . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.1.1 Computation of xation probability . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.1.2 Evolutionarily stable strategy in a nite population . . . . . . . . . . . . . . . . . . 32 2.2 Local approximation to xation probability under weak selection . . . . . . . . . . . . . . 34 2.2.1 Equivalence of stochastic processes under weak selection . . . . . . . . . . . . . . 35 2.2.2 Too limited validity of the local approximation . . . . . . . . . . . . . . . . . . . . 37 2.3 Bernstein approximation: global approximation to xation probability . . . . . . . . . . . . 42 Chapter3: Stochasticsingledrugchemotherapymodel 49 3.1 Fixation probability for the Moran process with three strategies . . . . . . . . . . . . . . . 53 3.2 Single drug adaptive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3 Single drug adaptive control for multiple evolutionary cycles . . . . . . . . . . . . . . . . . 69 3.4 Comparison of adaptive single drug chemotherapy schedule with standard clinical approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 ii Chapter4: Stochastictwodrugchemotherapymodel 79 4.1 S;R 1 ;R 2 multi drug model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2 Adaptive control of evolutionary cycles with additive multi drug schedule . . . . . . . . . 88 4.3 Adaptive control of evolutionary cycles with synergistic and antagonistic multi drug schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.3.1 Synergistic multi drug schedule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.3.2 Antagonistic multi drug schedule . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.4 Comparison of adaptive multi drug chemotherapy schedule with standard clinical approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Chapter5: Futuredirections 127 Bibliography 129 iii ListofTables 4.1 total dose (D), total time (T), and average dose (D/T) associated with adaptive chemother- apy schedules with antagonistic, additive and synergistic drug interactions during one evolutionary cycle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.2 the constant background tness (g) associated with adaptive chemotherapy schedules for antagonistic, additive and synergistic drug interactions. . . . . . . . . . . . . . . . . . . . . 104 iv ListofFigures 1.1 Defectors (D) playing a Prisoner’s dilemma game outcompete cooperators (C), grow in a S-shaped Gompertzian curve, and eventually saturate the entire population, eventually reaching a sub-optimal state according to the adjusted replicator system. For simulation, we setR = 3, S = 0, T = 5 andP = 1. (a) Initially nonzero defector subpopulations with the proportion, 0:01, keeps increasing over time in aS-shaped curve and eventually saturating the entire population; (b) The averaged tness function,<f >, is unfortunately minimized at that evolutionarily stable strategy,D, meaning thatx D = 1 is an asymptotically stable but sub-optimal state. . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1 The local approximation, N , in (2.12) to the rate of evolution,N 1 , of a singleTFT mutant under weak selection has limited validity with the upper threshold of selection strength, TRS, which is order ofO(N) for the repeated Prisoner’s dilemma game between TFT andALLD. (ParametersT = 5,R = 3,P = 1,S = 0 andn = 10 in all panels) (a)N = 10, TRS = 0:05567010; (b)N = 10 2 , TRS = 0:00437730; (c)N = 10 3 , TRS = 0:00042946; (d)N = 10 4 ,TRS = 0:00004287 . . . . . . . . . . . . . . . . . . . . 41 2.2 The Bernstein polynomials, B d (r 1 ), of degree, d, approximates the rate of evolution, r 1 :=N 1 , in (2.12) in the whole region of selection strength, being more accurate with the increase ofd. For simulation, we choose parameters,T = 5,R = 3,P = 1,S = 0 andn = 10, in all panels for the repeated Prisoner’s dilemma game betweenTFT and ALLD. (a)N = 10. Bothr 1 andB d (r 1 ) are increasing function on [0; 1] being greater than the rate of evolution,r w=0 1 :=N w=0 1 , under the neutral drift. Selection favorsTFT replacingALLD for all selection strength; (b)N = 10 2 . The same interpretation as for N = 10 is made but with a higher rate of evolution; (c)N = 10 3 . Selection favorsTFT replacingALLD for weak selection while it does not for strong selection according to either the exact or approximate functions withd > 1; (d)N = 10 4 . Selection opposes TFT replacingALLD for all selection strength according to the approximationB d (r 1 ) for alld, however, it does not if selection is extraordinarily weak according to the exact function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.3 The error, produced by the Bernstein polynomial,B d (r 1 ), in approximating the rate of evolution,r 1 =N 1 , is reduced as much as desired by increasing the degree,d, for each N although a relatively higher degree is required for an intermediate value ofN. Despite of its universality as a global approximation on [0; 1], it does not behave well particularly under weak selection, producing a relatively high error. (ParametersT = 5,R = 3, P = 1,S = 0 andn = 10 in all panels) (a)N = 10; (b)N = 10 2 ; (c)N = 10 3 ; (d)N = 10 4 47 v 3.1 The space set, N , of the Moran process withH,S andR is indicated as dots in the phase space,S 2 , whenN = 10. Each state, (H;S;R) = (i;Nij;j), is assigned to a lattice point inS 2 , wherei-numbering goes next to the side,HS, andj-numbering lies on the bottom of the side,SR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 The space set, N , of the Moran process withH,S andR gets denser in the phase space, S 2 , with the increase of the population size,N. For eachN,j N j is order ofO(N 2 ), being equal to (N+1)(N+2) 2 . (a)N = 10,j N j = 66; (b)N = 20,j N j = 231; (c)N = 30, j N j = 496; (d)N = 40,j N j = 861 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3 Local xation probability for the Moran process at the state, (H;S;R) = (0:8; 0:1; 0:1), approaches the proportion of area for the basin of attraction by the adjusted replicator dynamics as the population size, N, increases although unexpected features such as xation toH are observed for a smallN. (a)N = 100; (b)N = 500; (c)N = 1000 . . . . 55 3.4 The basins of attraction by the adjusted replicator dynamics with constant chemo concentration,C :=C(t), shows the smooth transition from the global attraction byS to the global absorption toR with the increase ofC, describing the competitive release of resistant subpopulations with high chemotherapy concentration. The system experiences bistability toS orR for an intermediate value ofC. (a)C 0:29; (b)C 0:37; (c) C 0:45; (d)C 0:53 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5 Fixation probability toR for the Moran process with constant chemotherapy concentra- tions,C :=C(t), shows the stochastic version of the competitive release of the resistant cells with high concentration for each population size,N. For a xedC, each region in a phase space with either a high or a low probability respectively approaches the basin of attraction to eitherR orS by the adjusted replicator system asN increases. (a)C 0:29, N = 100; (b)C 0:37,N = 100; (c)C 0:45,N = 100; (d)C 0:53,N = 100; (e) C 0:29,N = 500; (f)C 0:37,N = 500; (g)C 0:45,N = 500; (h)C 0:53, N = 500; (i)C 0:29,N = 1000; (j)C 0:37,N = 1000; (k)C 0:45,N = 1000; (l) C 0:53,N = 1000; (m)C 0:29,N = 2000; (n)C 0:37,N = 2000; (o)C 0:45, N = 2000; (p)C 0:53,N = 2000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.6 Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted replicator system for dierent constant chemotherapy values. (a) UnderC(t) 0, the competitive release of the sensitive subpopulations,S, to drug drives all trajectories to the S corner. (b) UnderC 0:7, the competitive release of the resistant subpopulationsR to drug drives all trajectories to theR corner. (c) Trajectories with two dierent constant chemotherapy combinations overlap at dierent times and generate a closed loop. . . . . . 60 3.7 Switching chemotherapy on and o at adequate times traps a trajectory associated with the adjusted replicator system within a closed loop. (a) The system, that starts atP , moves along a red line and arrives atQ when the high chemo dose,C(t) 0:7, is continuously administered duringT PQ = 19:36 unit time. When chemo is turned o during further T QP = 12:56 unit time, it returns to the initial point,P , eventually generating a closed loop,PQP ; (b) The tumor volume is controlled with the use of the adaptive schedule for an evolutionary cycle, experiencing tumor regression and recurrence. . . . . . . . . . . . . 61 vi 3.8 Realizations of multiple trajectories associated with the Moran process under adminis- tration of a constant chemotherapy show the ability of the stochastic system to behave similarly to what the adjusted replicator dynamics drive, getting closer as the population size increases. The Moran process, starting at a state near the cornerH withC(t) 0 (blue wiggled lines), possibly evolves and attains the homogeneous population of allS for eachN, having smoother trajectories asN increases and nally being similar to the deterministic trajectory (light blue line). Similarly, the Moran process starting at a state near the xed point (0; 0:6875; 0:3125) (red wiggled lines) is able to be driven to the state R withC(t) 0:7 for eachN. (a)C(t) 0,N = 1K; (b)C(t) 0,N = 5K; (a) C(t) 0,N = 10K; (a)C(t) 0,N = 50K; (e)C(t) 0:7,N = 1K; (f)C(t) 0:7, N = 5K; (g)C(t) 0:7,N = 10K; (h)C(t) 0:7,N = 50K . . . . . . . . . . . . . . . . 62 3.9 The spread of the distribution of points aroundQ (orP ) for 1; 000 realizations of the stochastic Moran process gets denser and demonstrates the shrunken randomness as the population size increases when each realization is under the administration of a constant chemo scheduleC(t) 0:7 (orC(t) 0) during half an evolutionary cycle,T PQ (or T QP ), since its exact start atP (orQ). (a)N = 1K; (b)N = 5K; (c)N = 10K; (d)N = 50K 63 3.10 The spread of the distribution of points in the principal axis coordinate system for 1; 000 realizations of the Moran process is, for a large population size, characterized as a multi Gaussian distribution aroundQ (orP ) when each realization is under the administration of a constant chemo scheduleC(t) 0:7 (orC(t) 0) during half an evolutionary cycle,T PQ (orT QP ), since its exact start atP (orQ). The mean frequency, S (or R ), of the sensitive (or resistant) subpopulations around the point,P (orQ), converges to the proportion ofS (orR) asN increases, with the decreasing semi-major axis, 1 , and semi-minor axis, 2 . (a)N = 1K, S = 0:1920, R = 0:2930; (b)N = 5K, S = 0:1223, R = 0:3452; (c)N = 10K, S = 0:1154, R = 0:3479; (d)N = 50K, S = 0:1099, R = 0:3590; (e)N = 1K, 1 = 0:2800, 2 = 0:1356; (f)N = 5K, 1 = 0:1401, 2 = 0:0386; (g)N = 10K, 1 = 0:1028, 2 = 0:0238; (h)N = 50K, 1 = 0:0460, 2 = 0:0093; (i)N = 1K, S = 0:8907, R = 0:0684; (j)N = 5K, S = 0:8987, R = 0:0606; (k)N = 10K, S = 0:8990, R = 0:0601; (l)N = 50K, S = 0:8992, R = 0:05989; (m)N = 1K, 1 = 0:0551, 2 = 0:0119; (n)N = 5K, 1 = 0:0224, 2 = 0:0050; (o)N = 10K, 1 = 0:0150, 2 = 0:0035; (p)N = 50K, 1 = 0:0067, 2 = 0:0016 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.11 The stochastic trajectory of one realization of the Moran process under the administration of the adaptive chemotherapy, associated with the adjusted replicator dynamics, develops a random walk along a lattice in a phase space,S 2 . The adaptive schedule is able to prevent the stochastic system from the saturation of cancer cells, even in a small population with N = 50. Starting atP = (0:04; 0:9; 0:06) (black dot), it moves along the lattice (red line) and reaches a neighborhood (green dot) ofQ = (0:53; 0:11; 0:36), duringT PQ evolution time. Turning o the chemo at that green dot, the stochastic system evolves (blue line) and eventually reaches a neighborhood (yellow dot) of the initial point,P , duringT QP . . . 66 vii 3.12 The averaged trajectory of 1; 000 realizations of the Moran process under the adaptive schedule, associated with the adjusted replicator system, during one evolutionary cycle (T PQ +T QP = 31:92 unit time) ts the corresponding deterministic trajectory for a large population size. The Moran process is likely to return nearly to the initial state with a high probability for a largeN even though the spread of the distribution of the points nearQ is still wide. (a) the distribution of the points associated with the adaptive chemo schedule forN = 10K; (b) the trajectory of one single realization of the Moran process withN = 10K; (c) the averaged trajectory of 1; 000 realizations of the Moran process withN = 10K; (d) the distribution of the points forN = 50K; (e) the trajectory of one single realization withN = 50K; (f) the averaged trajectory withN = 50K . . . . . . . . 67 3.13 The spread of the distribution of points for 1; 000 realizations of the stochastic Moran process withN = 50K is characterized as a multivariate Gaussian distribution, centered nearly at the initial point,P , when each realization is under the administration of the adaptive schedule during one evolutionary cycle (T PQ +T QP = 31:92 unit time). As shown in theSR and the principal axis coordinate system, the deviation is equally likely to each other in either directions. (a) the distribution of the points aroundP in theSR coordinate system; (b) the kernel density distribution in theSR coordinate system; (c) the distribution of the points in the principal axis coordinate system; (d) the kernel density distribution in the principal axis coordinate system . . . . . . . . . . . . . . . . . . . . . . 68 3.14 For eachN, the averaged trajectory of 1; 000 realizations of the stochastic Moran process shows that the saturation of cancer cells can be delayed until the 4th evolutionary cycle when each realization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start atP , although the tightness of the the averaged stochastic trajectory to the deterministic one lasts shorter with a smaller population size. (a) the trajectory of one single realization for the Moran process withN = 10K; (b) the averaged trajectory of 1; 000 realizations for the Moran process withN = 10K; (c) the trajectory of one single realization withN = 50K; (d) the averaged trajectory withN = 50K . . . . . . . . . . . 70 3.15 The spread of the distribution of terminal points (blue dots) aroundP in each round for 1; 000 realizations of the stochastic Moran process becomes wider as the number of rounds increases when each realization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start atP . (a)N = 10K, round 1; (b)N = 10K, round 2; (c)N = 10K, round 1; (d)N = 10K, round 4; (e)N = 50K, round 1; (f)N = 50K, round 2; (g) N = 50K, round 1; (h)N = 50K, round 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.16 As the number of rounds increases, the spread of the distribution of points aroundP for 1; 000 realizations of the stochastic Moran process forN = 50K becomes centered towards the homogeneous sensitive state and immediately looses the multivariate Gaussian-like distribution, where each realization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start atP . (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 1; (f) round 2; (g) round 3; (h) round 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 viii 3.17 The semi-major (and -minor) axis of the distribution of the points aroundP at the end of each evolutionary cycle for 1; 000 realizations of the Moran process overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the adaptive schedule associated with the adjusted replicator system during 8 evolutionary cycles since its exact start atP . (a)N = 10K; (b) N = 10K, the log-log t; (c)N = 50K; (d)N = 50K, the log-log t . . . . . . . . . . . . 73 3.18 Two standard clinical approaches, the maximum tolerated schedule (MTD) and the low-dose metronomic schedule (LDM), are designed in order to have the same total dose as the adaptive chemotherapy schedule associated with the deterministic replicator system has during 4 rounds. In each cycle of the length, 31:92, MTD delivers a drug at the highest concentration during the rst 13:552 unit time, followed by no chemo until the end of each cycle. On the other hand, a drug is continuously delivered during the whole cycles at as low concentration as 0:42456 for LDM. . . . . . . . . . . . . . . . . . . . . . . 74 3.19 The adaptive chemotherapy schedule associated with the replicator system is compared to the standard clinical approaches, MTD and LDM, in order to demonstrate its eciency in terms of delaying the time of saturation of tumor cells, which is attained before the rst round ends under either MTD or LDM schedules. The adaptive schedule beats the other two standard clinical chemo schedules since not only it prevents the system from converging to a cancerous state and but also the tumor size is thus controlled for 4 rounds on average. (Each of 1; 000 realizations of the Moran process evolves under the administration of the adaptive schedule, MTD or LDM independently during 4 evolutionary cycles since its exact start atP .) (a)N = 10K, one single realization; (b)N = 10K, the tumor size associated with the single trajectory; (c)N = 10K, the averaged trajectory; (d)N = 10K, the tumor size associated with the averaged trajectory; (e)N = 50K, one single realization; (f)N = 50K, the tumor size associated with the single trajectory; (g)N = 50K, the averaged trajectory; (h)N = 50K, the tumor size associated with the averaged trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.20 Among 1; 000 realizations of the Moran process, the rate that it is well controlled and returns back nearly to the initial state, neither gaining the full tumor volume nor having over 99% sensitive cells after 4 evolutionary cycles, is computed under the adaptive, MTD or LDM chemo schedules independently. The adaptive chemo schedule obviously is superior than other two schedules, having the nearly full rate whenN = 50K, while MTD has an extraordinarily small but positive success rate up to the 3 rounds. . . . . . . . 77 4.1 Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted replicator system for dierent constant chemotherapy values withe = 0. (a) UnderC 1 = 0 andC 2 = 0, the tumor saturates to theS corner regardless of the initial distribution of the three subpopulations. (b) UnderC 1 = 0:8 andC 2 = 0, the competitive release of the resistant subpopulationsR 2 to drug 1 drives all trajectories to theR 2 corner. (c) Under C 1 = 0 andC 2 = 0:8, the competitive release of the resistant subpopulationsR 1 to drug 2 drives all trajectories to theR 1 corner. (d) Trajectories with three dierent constant chemotherapy combinations of drug 1 and drug 2 overlap at dierent times and generate a closed loop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 ix 4.2 Realizations of multiple trajectories associated with the Moran process under administra- tion of dierent constant chemotherapy combinations of drug 1 and drug 2 (e = 0) show the ability of the stochastic system to behave similarly to what the adjusted replicator dynamics drive, getting closer as the population size increases. For example, the Moran process, starting at a state near the cornerR 2 withC 1 (t) 0,C 2 (t) 0 (blue wiggled lines), possibly evolves and attains the homogeneous population of allS for eachN, having smoother trajectories asN increases and nally being similar to the deterministic trajectory (light blue line). (a)C(t) 0,N = 1K; (b)C(t) 0,N = 5K; (a)C(t) 0, N = 10K; (a)C(t) 0,N = 50K; (e)C(t) 0:7,N = 1K; (f)C(t) 0:7,N = 5K; (g)C(t) 0:7,N = 10K; (h)C(t) 0:7,N = 50K . . . . . . . . . . . . . . . . . . . . . 85 4.3 Switching constant chemotherapy combinations of two drugs (e = 0) on and o at adequate times traps a trajectory, associated with the adjusted replicator system, within a closed loop, controlling a tumor size. (a)OP :C 1 = 0:8,C 2 = 0 duringT OP = 6:933 unit time,PQ:C 1 = 0,C 2 = 0 duringT PQ = 6:248 unit time,QO:C 1 = 0,C 2 = 0:8 during T QO = 6:324 unit time; (b) the trajectory treated according to the multi drug additive adaptive schedule and the untreated trajectory (pink) being driven to theS corner; (c) tumor size under untreated (pink) and the adaptive chemotherapy schedule. (For this numerical experiment, we takeg = 1:5519.) . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.4 The spread of the distribution of points aroundP (Q orO) for 1, 000 realizations of the stochastic Moran process gets denser and demonstrates the shrunken randomness as the population size increases when each realization is under the administration of a constant chemo combination C 1 (t) 0:8;C 2 (t) 0 (C 1 (t) 0;C 2 (t) 0 or C 1 (t) 0;C 2 (t) 0:8) withe = 0 during 1=3 evolutionary cycle,T OP (T PQ orT QO ), since its exact start atO (P orQ). (a)N = 1K; (b)N = 5K; (c)N = 10K; (d)N = 50K . 89 4.5 The spread of the distribution of points in the principal axis coordinate system for 1; 000 realizations of the Moran process is, for a large population size, characterized as a multivariate Gaussian distribution aroundP (Q orO) when each realization is under the administration of a constant chemotherapy combinationC 1 (t) 0:8;C 2 (t) 0 (C 1 (t) 0;C 2 (t) 0 orC 1 (t) 0;C 2 (t) 0:8) withe = 0 during 1=3 evolutionary cycle,T OP (T PQ orT QO ), since its exact start atO (P orQ). As the population size increases, both the semi-major axis, 1 and the semi-minor axis, 2 , decrease. (a)N = 1K, 1 = 0:0977, 2 = 0:0296; (b)N = 5K, 1 = 0:0460, 2 = 0:0127; (c)N = 10K, 1 = 0:0320, 2 = 0:0091; (d)N = 50K, 1 = 0:0138, 2 = 0:0040; (e)N = 1K, 1 = 0:0318, 2 = 0:0120; (f)N = 5K, 1 = 0:0140, 2 = 0:0085; (g)N = 10K, 1 = 0:0097, 2 = 0:0060; (h)N = 50K, 1 = 0:0044, 2 = 0:0027; (i)N = 1K, 1 = 0:0860, 2 = 0:0207; (j)N = 5K, 1 = 0:0402, 2 = 0:0093; (k)N = 10K, 1 = 0:0273, 2 = 0:0065; (l)N = 50K, 1 = 0:0125, 2 = 0:0030 . . . . . . . . . . . . . 90 4.6 The stochastic trajectory of one realization of the Moran process under the administration of the additive adaptive chemotherapy, associated with the adjusted replicator dynamics, as in Figure 4.3a develops a random walk along a lattice in a phase space,S 2 . The adaptive schedule is able to prevent the stochastic system from the saturation of cancer cells, even in a small population withN = 30. (a)N = 30; (b)N = 40; (c)N = 50 . . . . . . . . . . 91 x 4.7 The averaged trajectory of 1; 000 realizations of the Moran process under the additive adaptive schedule, associated with the adjusted replicator system, during one round (T OP +T PQ +T QO = 19:2653 unit time) ts the corresponding deterministic trajectory for a large population size withe = 0. The Moran process is likely to return nearly to the initial state with a high probability for a largeN even though the spread of the distribution of the points nearQ is still wide. (a) the distribution of the points associated with the adaptive chemo schedule forN = 10K; (b) the trajectory of one single realization of the Moran process withN = 10K; (c) the averaged trajectory of 1; 000 realizations of the Moran process withN = 10K; (d) the distribution of the points forN = 50K; (e) the trajectory of one single realization withN = 50K; (f) the averaged trajectory with N = 50K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.8 The spread of the distribution of points for 1; 000 realizations of the stochastic Moran process is characterized as a multivariate Gaussian distribution, centered nearly at the initial point, O, when each realization is under the administration of the multi drug adaptive chemo schedule associated with the adjusted replicator dynamics during one round (T OP +T PQ +T QO = 19:2653 unit time) withe = 0. The mean frequency, R 1 (or R 2 ), of the subpopulation,R 1 (orR 2 ), around the point,O, converges to the proportion ofR 1 (orR 2 ) asN increases, with the decreasing semi-major axis, 1 , and semi-minor axis, 2 . (a)N = 10K, R 1 = 0:3175 , R 2 = 0:0617; (b)N = 10K, 1 = 0:0453, 2 = 0:0099; (c)N = 50K, R 1 = 0:3177, R 2 = 0:0617; (d)N = 50K, 1 = 0:0211, 2 = 0:0045 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.9 Two realizations of the stochastic Moran process of sizeN = 10K, starting atO and evolving under the adaptive chemo schedule during 8 rounds withe = 0, show a great dierence in their tumor size as well as in their trajectories. (a) one realization; (b) corresponding tumor size to Figure 4.9a; (c) another realization; (d) corresponding tumor size to Figure 4.9c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.10 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the saturation of cancer cells can be delayed until the end of 8th round when each realization evolves under the administration of the multi drug adaptive chemo schedule, associated with the adjusted replicator system, withe = 0 during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.10a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.10c 96 4.11 The spread of the distribution of terminal points (green dots) aroundO in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug adaptive chemo schedule, associated with the adjusted replicator system, withe = 0 during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 . . . . . . . . . . . . . . . . . 96 xi 4.12 The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K, all starting atO and evolving under the multi drug adaptive chemo schedule with e = 0 during 8 rounds, shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , in the principal axis coordinate system as the number of rounds increases. (a) round 1, 1 = 0:0211, 2 = 0:0045; (b) round 2, 1 = 0:0319, 2 = 0:0068; (c) round 3, 1 = 0:0425, 2 = 0:0084; (d) round 4, 1 = 0:0561, 2 = 0:0096; (e) round 5, 1 = 0:0710, 2 = 0:0108; (f) round 6, 1 = 0:0871, 2 = 0:0116; (g) round 7, 1 = 0:1051, 2 = 0:0128; (h) round 8, 1 = 0:1243, 2 = 0:0136 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.13 The semi-major (and -minor) axis of the distribution of the points aroundO at the end of each round for 1; 000 realizations of the Moran process overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the multi drug adaptive chemo schedule associated with the adjusted replicator system withe = 0 during 8 rounds since its exact start atO. (a)N = 10K; (b) N = 10K, the log-log t; (c)N = 50K; (d)N = 50K, the log-log t . . . . . . . . . . . . 99 4.14 Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted replicator system for dierent constant chemotherapy values for each drug interaction. Under C 1 = 0 and C 2 = 0, the tumor saturates to the S corner regardless of the initial distribution of the three subpopulations. UnderC 1 = 0:5 andC 2 = 0:2, the competitive release of the resistant subpopulations,R 2 , to drug 1 drives all trajectories to theR 2 corner. UnderC 1 = 0:2 andC 2 = 0:5, the competitive release of the resistant subpopulations,R 1 , to drug 2 drives all trajectories to theR 1 corner. (a)e = 0; (b)e = 0:3; (c)e =0:3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.15 Deterministic evolution of subpopulations, S, R 1 andR 2 , by the adjusted replicator system generates a closed loop,OPQO, when it starts atO and evolves under the multi drug adaptive chemotherapy schedule during one round, resulting in the tumor size controlled. (a) Withe = 0, the total dose delivered to generate one evolutionary cycle is 17:7604 during 29:4440 unit time. OP : C 1 = 0:5,C 2 = 0:2 duringT OP = 14:226 unit time,PQ: C 1 = 0,C 2 = 0 duringT PQ = 4:072 unit time,QO: C 1 = 0:2,C 2 = 0:5 duringT QO = 11:146 unit time; (b) the corresponding deterministic trajectory; (c) the corresponding tumor size with the averaged background tness,g = 1:4527; (d) With e = 0:3, the total dose delivered to generate one evolutionary cycle is 17:7723 during 32:29 unit time. OP :C 1 = 0:5,C 2 = 0:2 duringT OP = 14:260 unit time,PQ:C 1 = 0, C 2 = 0 duringT PQ = 6:901 unit time,QO: C 1 = 0:2,C 2 = 0:5 duringT QO = 11:129 unit time; (e) the corresponding deterministic trajectory; (f) the corresponding tumor size withg = 1:4857; (g) Withe =0:3, the total dose delivered to generate one evolutionary cycle is 17:8150 during 26:8740 unit time.OP :C 1 = 0:5,C 2 = 0:2 duringT OP = 14:300 unit time,PQ:C 1 = 0,C 2 = 0 duringT PQ = 1:424 unit time,QO:C 1 = 0:2,C 2 = 0:5 duringT QO = 11:150 unit time; (h) the corresponding deterministic trajectory; (i) the corresponding tumor size withg = 1:4182 . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 xii 4.16 The spread of the distribution of terminal points (green dots) aroundO in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15a during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 . . . 105 4.17 The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug additive (e = 0) adaptive chemo schedule as in Figure 4.15a during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0156, 2 = 0:0076; (b) round 2, 1 = 0:0268, 2 = 0:0123; (c) round 3, 1 = 0:0432, 2 = 0:0172; (d) round 4, 1 = 0:0675, 2 = 0:0223; (e) round 5, 1 = 0:1069, 2 = 0:0288; (f) round 6, 1 = 0:1623, 2 = 0:0371; (g) round 7, 1 = 0:2233, 2 = 0:0443; (h) round 8, 1 = 0:2797, 2 = 0:0474 . . . . . . . . . . . . . . 107 4.18 1; 000 realizations of the stochastic Moran process with size N = 50K show that the saturation of cancer cells can be delayed when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15a during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.18a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.18c 108 4.19 The spread of the distribution of terminal points (green dots) aroundO in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15d during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 . . 109 4.20 The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug synergistic (e = 0:3) adaptive chemo schedule as in Figure 4.15d during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0166, 2 = 0:0082; (b) round 2, 1 = 0:0279, 2 = 0:0140; (c) round 3, 1 = 0:0445, 2 = 0:0200; (d) round 4, 1 = 0:0743, 2 = 0:0262; (e) round 5, 1 = 0:1231, 2 = 0:0344; (f) round 6, 1 = 0:1857, 2 = 0:0420; (g) round 7, 1 = 0:2496, 2 = 0:0467; (h) round 8, 1 = 0:3057, 2 = 0:0529 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 xiii 4.21 1; 000 realizations of the stochastic Moran process with size N = 50K show that the saturation of cancer cells can be delayed when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15d during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.21a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.21c 111 4.22 The spread of the distribution of terminal points (green dots) aroundO in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15g during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.23 The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug antagonistic (e =0:3) adaptive chemo schedule as in Figure 4.15g during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0146, 2 = 0:0073; (b) round 2, 1 = 0:0250, 2 = 0:0118; (c) round 3, 1 = 0:0374, 2 = 0:0158; (d) round 4, 1 = 0:0584, 2 = 0:0206; (e) round 5, 1 = 0:0906, 2 = 0:0267; (f) round 6, 1 = 0:1369, 2 = 0:0329; (g) round 7, 1 = 0:1942, 2 = 0:0369; (h) round 8, 1 = 0:2516, 2 = 0:0409 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.24 1; 000 realizations of the stochastic Moran process with size N = 50K show that the saturation of cancer cells can be delayed when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15g during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.24a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.24c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.25 The semi-major (and -minor) axis of the distribution of the points aroundP at the end of each evolutionary cycle for 1; 000 realizations of the Moran process withN = 10K or N = 50K overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the adaptive schedule associated with the adjusted replicator system during 8 evolutionary cycles since its exact start atO. (a)e = 0; (b)e = 0:3; (c)e =0:3; (d)e = 0, the log-log t; (e) e = 0:3, the log-log t; (f)e =0:3, the log-log t . . . . . . . . . . . . . . . . . . . . . . 116 xiv 4.26 The multiple additive (e = 0) standard clinical approaches are designed to have the same total dose, being equal to 17:7604, as the amount that is delivered during one round (29:4440 unit time) according to the adaptive chemo schedule in Figure 4.15a. (a) adaptive; (b)MTD 1 : the drug 1 is maximally administered during the beginningT MTD:C 1 = 9:3422 unit time while the drug 2 is delivered during the lastT MTD:C 2 = 8:4182 unit time. (c) MTD 2 : the drug 2 is maximally administered during the beginningT MTD:C 2 = 8:4182 unit time while the drug 1 is delivered at largest during the lastT MTD:C 1 = 9:3422 unit time. (d)LDM: both the drug 1 and drug 2 are constantly administered during the whole rounds at the level ofC 1 = 0:347287 andC 2 = 0:255905, respectively. . . . . . . . . . . . 118 4.27 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the saturation of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.27a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.27a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.27d; (f) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.27d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.28 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the saturation of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.28a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.28a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.28d; (f) the averaged tness ofS,R 1 and R 2 cells corresponding to Figure 4.28d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.29 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the saturation of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.29a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.29a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.29d; (f) the averaged tness ofS,R 1 and R 2 cells corresponding to Figure 4.29d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 xv Abstract This thesis is concerned with modeling tumor growth and chemotherapy response using stochastic evolu- tionary game theory models. In particular we develop stochastic N-cell models of a heterogeneous (mul- tiple cell types) tumor using a Moran process model with frequency dependent tness, which in the limit N!1 converges to the deterministic adjusted replicator dynamical system as its mean-eld limit. This limit and some of the details of our model and background literature are described in Chapter 1. Chapter 2 develops new results associated with the xation probability of cell types that do not necessarily depend on the assumption of weak selection, but are valid across the full range of selection strengths. Chap- ters 3 and 4 develop our main results on adaptive chemotherapy scheduling using our stochastic N-cell evolutionary game theory model, both single drug and multi-drug scheduling, with the goal of avoiding chemo-resistance via the mechanism of competitive release, a concept borrowed from the ecology litera- ture. The methods we develop are superior in measurable ways to more standard chemotherapy schedules currently in use at cancer centers all over the world, such as maximum tolerated dose (MTD) schedules and low-dose metronomic schedules (LDM). The nal chapter describes potential future research directions. xvi Chapter1 Introduction Cancer can be seen as a consequence of interactions of cells with their surrounding environment as well as each other all of which contribute to its development. One approach to modeling the co-evolution of cancer cell populations that make up a tumor is to apply evolutionary game theory, an approach ideally suited to model the interactions dictated by the various tness levels of the cells [37]. In the context of evolutionary game theory, tness means the success rate at reproduction of a cell whose genotype is passed from generation to generation. This information is quantied through a payo matrix that denes the game being played. Each group of cells is assumed to have a dierent success rate which depends not only on its own frequency in the population but also other types’ in a population. As the populations co-evolve, success breeds success, and failure spirals downward, which is the essence of the replicator dynamical system and reinforcement learning equations. From Darwinian evolutionary theory, we know that natural selection is the key mechanism of evolution of a co-evolving population. With no selection, the system proceeds under random drift. However, even under weak selection, some sub-populations can develop higher tness [12]. Eventually these frequency-dependent tness functions determine and guide the fate of groups of cells when it is compared to the averaged tness of the entire population. These ideas form the basis for the replicator dynamical system we use as our basic model, both in its mean-eld deterministic form and its stochastic nite population form. 1 Many people have attempted to mathematically model tumor evolution with the aim of controlling its growth at a theoretical/computational level in order to guide and motivate clinicians to try new strategies and test them in clinical trials. These eorts in turn gives rise to important insights that can aid clinicians in nding and exploring to methods and avenues of possible treatments for their patients. Chemother- apy is the most important clinical treatment that is executed for the purpose of alleviating/controlling cancer at a systemic level. When it is administered to a patient in standard treatment schedules, it often decreases the size of the tumor initially (tumor regression), but then can be followed by re-growth (tu- mor recurrence) due to chemo-resistance. A working hypothesis in this thesis is that designing adaptive chemotherapy schedules that change as the tumor evolves can sometimes be an eective tool in delaying tumor recurrence. We develop computational models in this thesis (based on evolutionary game theory) and test various hypothetical adaptive chemotherapy schedules for the purpose of designing schedules that combat chemo-resistance. Much like adaptive control theory, the chemo-dosing adapts to the current state of the tumor, which of course requires ecient monitoring of the sub-populations of cells making up the tumor.The idea in this approach is to focus more on controlling the tumor growth with an allowable increase in its size than on eradicating all present cancerous cells [18, 19]. To better understand the role of chemotherapy in the evolution of the tumor, we divide (coarse-grain) the tumor cells into three groups: healthy cells, chemo-sensitive cells, and chemo-resistant cells. Assuming that cancer is the outcome of interactions of only these three dierent types, modeling the interplay be- tween cancer and a drug is then equivalent to designing a tness landscape for each type by which selection dynamics result, having the concentration of a drug as a controller of the system. We interpret the tumor recurrence as a result of a continued drug use in the context ofthecompetitiverelease, which demonstrates that when two dierent types of groups are competing against each other for limited resources, currently subordinate subpopulations start proliferating and nearly take over the entire population when a catalyst that decreases the tness of presently dominant subpopulations is introduced into a population [11]. In 2 other words, the sensitive cell subpopulations are dominant when no drug is administered since they have higher tness, but when chemotherapy is introduced, it lowers the tness of the sensitive cell population allowing the resistant population to proliferate. That is, a high enough concentration of a drug acts as a catalyst that exchanges the tness of two dierent cancerous subpopulations competing each other and eventually the fate of the system. With this understanding that the drug concentration is a controller of the system, we let it enter the population equations through the selection of sensitive and/or resistant cells such that a high dose reduces the tness level of the sensitive cells. We adopt a Prisoner’s dilemma (PD) game to assign the payo that each cell receives from interaction with another. In its original form, the PD game is a two-player game between a cooperator and a defector from classical game theory, which is determined by what is called the Prisoner’s dilemma inequality [3]. It describes that the rationality of players who are selshly interested in maximizing their own rewards in playing a game eventually results in the sub-optimal scenario, both becoming defectors and receiving a lower payo than they would have gained if they cooperated. The PD game has received a great attention in evolutionary game theory since the emergence of cooperation, which is often seen in nature but cannot be understood if assuming rationality of participants as in classical game theory, can be explained when this game is played repeatedly by a number of individuals in a population along with the natural selection acting on [3, 21, 39, 42, 44, 55]. We model a payo matrix for healthy, sensitive and resistant cells such that healthy cells are cooperators whereas cancerous cells (resistant cells and sensitive cells) are defectors. For the payo from the interaction between resistant cells and sensitive cells, we take a cost of resistance into account [8, 13, 20]. In many studies of drug resistance, it has been shown that resistant cells are present in a population as a result of mutation long before any drug is administered, and there is a tness cost for cells becoming resistant. We adopt a PD game for the interaction between resistant and sensitive cells such that resistant cells are cooperators (lower tness) whereas sensitive cells are defectors (higher tness). The payo matrix we use assures growth of the sensitive cell population under low drug concentrations, 3 while allowing for the competitive release of resistant cells under drug concentrations at a high enough threshold level. In order to understand the co-evolutionary dynamics of these sub-groups in a nite population, we develop a Moran process model with three strategies, which is a discrete-time stochastic process [40, 51]. Knowing thatthe(adjusted)replicatordynamics, which assumes an innite size of population, is the asymp- totic dynamic (N!1) of the Moran process [53, 54] and that an adaptive chemotherapy schedule is able to control the growth of a tumor by allowing a slight increase in size at an acceptable level as much as desired by repeating several cycles when it is modeled with the replicator equations [36], we will apply the adaptive chemotherapy schedule associated with the adjusted replicator dynamics to the Moran process and investigate how long the application is successful in controlling tumor growth in a nite population or how badly this adaptation fails to be used as a proper treatment. Moreover, for comparative purposes, the same investigation is carried out for two typical clinical chemo-schedules that ensure the same amount of drug use as the adaptive schedule: the maximum tolerated dose (MTD) and the low-dose metronomic schedules (LDM). The natural extension of this study is obtained by generalizing the number of drugs. When more than one drug is used simultaneously, then some combinations act synergistically while others act antagonis- tically, compared to the case where each drug is independently delivered. The study of drug interaction has a long history and drug interactions can be mainly divided into three cases: additive, synergistic or antagonistic drug interactions [5, 9]. Following a drug interaction parameter introduced in [30] allows those three cases to be distinguishable from each other. We the model the eect of two drugs on cancerous populations by the Moran process with three sub-populations where there are two kinds of resistant cells to each of two drugs as well as sensitive cells to both drugs. A similar but more complicated investigation compared to a single drug model is achieved with this stochastic two drug chemotherapy model. This thesis is concerned with nding the optimal single drug and two-drug chemotherapy schedule that delays 4 the time of the saturation of cancerous cells in a population, and discuss its advantage compared to the standard clinical approaches. 1.1 Evolutionarygametheory Classical game theory, established by J. Von Neumann and O. Morgenstern in [35] in 1928, had been studied sporadically in the 19th century. It has received great attention and been applied to various areas from social science to ecology in the 20th century in an eort to understand how this decision-making model from interactions of rational players results in human behavior [33]. A classical game theory is originally evoked in a form of two-player games, where players participate in a game receiving a certain amount of payo according to the strategies of their choices, and are assumed to be rational and struggle for earning a maximal possible reward being selsh. In 1950, a strict Nash equilibrium, a solution to a game, was invented by J. Nash, which is a strategy with which a players cannot be beneted in payo by switching a strategy [34]. Game theory had not been signicantly considered valuable in biology until J. M. Smith founded the eld of evolutionary game theory by taking frequency dependence into account and applying it to evolution in order to modelDarwinism [12] in the early 1970’s in [46, 47, 48, 49] while similar attempts had been preceded in the 1930s by population geneticist such as R. A. Fisher, J. B. S. Haldane and S. Wright on the purpose of modelingrandomdrift with the neglect of frequency in a nite population. Evolutionary game theory is distinct from classical game theory although they share some of the same ingredients of a game such as the payo matrix that players or populations of players receive from interactions. However, being invented to describe the dynamics of subpopulations in a well-mixed reproductive population where each individual has the ability to reproduce its copy and a genotype is passed on to its ospring by heredity. In a population of cells where each cell has unique genotypic characteristics, cells are players of a game, genotypes are strategies, and tness or success rate of a cell at reproductions is payo. A group of cells with a higher expected tness produces more ospring whereas a group with a lower expected tness 5 replicates less. Unlike classical game theory, the focus of evolutionary game theory is not necessarily nding an optimal strategy with which a group of cells obtains a maximal payo, but at achieving a state in which a sub-population is stable against invasion by an innitesimally small increase of the frequency of a mutant population when a game is repeatedly played by a number of cells in a population. Thus, evolutionary game theory is described by a dynamical system, whereas classical game theory is described by a payo matrix. 1.1.1 (Adjusted)Replicatordynamicsinaninnitepopulation After the term "replicator" was rst invented by R. Dawkins in [14], the replicator equation has become the most popular model that describes the dynamics of coevolving populations in a well-mixed population assuming an innite size since its rst introduction by P. G. Taylor and L. B. Jonker in [52] in 1978 [23, 41, 45]. For a game with two strategies, it is dened by _ x A = (f A <f >)x A ; _ x B = (f B <f >)x B ; (1.1) where ~ x := (x A ;x B ) | is a frequency vector of A and Bsubpopulations such that x A +x B = 1, f A := (M~ x) 1 = ax A +bx B andf B := (M~ x) 2 = cx A +dx B are repectively the expected tness of A andBsubpopulations, and< f >:=~ x | M~ x = x A f A +x B f B is the averaged tness in the entire populations. Here,M is the payo matrix dening the interactions of the two (or more) sub-populations. It is a nonlinear deterministic ordinary dierential equation that describes the relative growth of competing subpopulations. It is clear from the frequency dependent equations in (1.1) not only that a certain type of subpopulations grows faster in density when it is abundant in a population but also that it grows at the rate of the relative dierence of its expected tness to the averaged tness,< f >, of the whole population. We only present the replicator equations for a game with two strategies but it is similarly dened for a 6 game with any nite number of strategies. Havingx A +x B = 1, the equations in (1.1) are equivalent to a single equation _ x A = (f A f B )(1x A )x A (1.2) which demonstrates that a certain type of subpopulations increases in frequency whenever it has a higher expected tness than other types of subpopulations. Deterministic evolutionary dynamic is often explained by the adjusted replicator equation instead of the regular replicator equation, which is dened by: _ x A = f A <f > <f > x A ; _ x B = f B <f > <f > x B ; (1.3) and withx A +x B = 1, is identical to _ x A = f A f B <f > (1x A )x A : (1.4) As long as the tness function,< f >, of the entire populations is positive, the equations (1.3) are not qualitatively dierent from equations (1.1), except for a dierence in time scale. The replicator equation is not the only one model that illustrates deterministic dynamical systems, and many studies in evolutionary game theory have been conducted with other models including ODEs, dierential inclusions, dierence equations and reaction-diusion system as summarized in [22]. However, it is worth deserving its fame when taking its equivalence to the Lotka-Volterra equation into account, which is a well-known model in ecology to describe the interactions of two populations over time, having a long history [6, 43]. 7 1.1.2 Stochasticmodelsinanitepopulation Although a deterministic model is often a convenient description of population dynamics, sometimes there is a need for a model describing nite population interactions which have inherent stochastic uctuations. D. Foster and P. Young suggested a continuous-time stochastic model that is designed by simply adding the Wiener process to the deterministic replicator equation, presuming that the white noise is congured by the Gaussian distribution in [17]. However, a discrete-time nite-space stochastic model is enough to characterize cell division and its inheritance from generation to generation as long as the model is equipped with tness. One of frequently used such models is a frequency-dependentWright-Fisherprocess although it was introduced a long ago by S. Wright and R. A. Fisher in 1930’s [16, 60] with the absence of frequency [57, 25]. Suppose that a population of sizeN consists of cells who has unique genotype betweenA and B. The number of cells in the subpopulationA then denes a state of the system, and it is assumed in this model to be determined byN independent Bernoulli trials where the success and failure probabilities are weighted by the expected tness ofA andB. In precise, the current state,i, changes to a state,j, in the next step with a probability,T i!j , T i!j = N j if A if A + (Ni)f B j (Ni)f B if A + (Ni)f B Nj (1.5) forj = 0; ;N. This sampling method, which is binomial for two strategies and multinomial for higher strategies, produces nonoverlapping generations and denes a discrete-time nite-space Markov chain with nontridiagonal transition matrix, meaning that every states is reachable from any states except for absorbing states,i = 0 ori = N. The power of this model is not only that it is built with a well-studied distribution, but also the population size is adjustable at all times without being xed so that the total population size,N, at the current step can jump to, say, 2N in the following step, which is benecial in 8 describing cell divisions. However, this process being a Markov chain with a nontridiagonal transition probability even in a constant population case challenges computational complexity for analytical study. Unlike the Wright-Fisher process not being a birth-and-death process, a birth-death process, so-called apairwisecomparisonprocess, has been studied in depth [2, 50, 54, 56, 58]. In this model, a randomly chosen individual updates its strategy to another randomly chosen individual’s with a bigger probability than 1=2 if the partner has a higher payo at each time step. Precisely, the probability,p Fermi , that anAindividual replaces aBindividual is given by p Fermi = 1 1 +e (f A f B ) ; (1.6) where it is represented as a function of tness dierence, f := f A f B , and 0 is the strength of natural selection corresponding to the random drift when = 0. Thus, an update is more likely to occur as the tness dierence, f, gets bigger. The fact that the number of individuals of a certain type can change at most by1 in one time step and the total population size remains xed in evolution denes a Markov chain with a tridiagonal transition probability. For this pairwise comparison model, the transition probability,T + i , that the system increases the number ofAindividuals fromi toi + 1 and the transition probability,T i , that it decreases the size ofAindividuals fromi toi 1 are given by T i = i N Ni N 1 1 +e (f A f B ) : (1.7) The advantages of this model is that the selection intensity parameter,, is unbounded and hence it can be set as arbitrary strong as needed while either weak or strong selections are of major interests to be analyzed. Also, the exponential form of probabilities in (1.6) and its expression only through the tness dierence, f, helps to obtain fruitful analytical observations such as xation probability and xation time, which are of key interests for stochastic models and we will be discussed in Chapter 2. This pairwise 9 comparison process is much more known asFermiprocess with the recognition of use of the Fermi function, stemming from physics, for probability of switching strategies in (1.6). Having a specic function for switching probability (1.6), the Fermi process is one special case of a local update model although it generalizes selection intensity to an unbounded interval. A local update model instead has a switching probability,p local update , given by p local update = 1 2 + w 2 f A f B f max ; (1.8) where f max is the maximum possible tness dierence and 0 w indicates the selection intensity. With this local update model, the transition probabilities similar to (1.7) are given by T i = i N Ni N 1 2 w 2 f A f B f max : (1.9) The local update model has been mainly studied in 2 2 game especially in the comparison to other stochastic processes [54, 57] and it experienced a generalization in the number of strategies of a game in [53]. At the expense of complexity of functions of interest such as transition or xation probability, it can also be generalized by modifying the switching or imitating probability, shown in (1.8), such that it possesses the same quality that the tter grows more rapidly with a higher probability of switching and that a bigger tness dierence leads to a higher updating probability. Dening the updating probability by a function,g(f), of both selection and tness dierence such thatg 0 (f) is nondecreasing in f, a generalized local update model was discussed in [61]. These local update models were discussed at a group level as well as an individual level in the context of spatial dynamics in [44] with the consideration that cells are more interacting with their neighbors than ones in a far distance. 10 1.2 Moranprocessinanitepopulation Among many stochastic processes describing co-evolving nite populations, the one that has attracted many researchers in evolutionary game theory is the frequency-dependentMoranprocess since it was de- veloped by M. A. Nowaketal. in [40] in 2004 [1, 2, 10, 29, 31, 32, 40, 51, 53, 54]. A local update model that only requires local information, which is the tness dierence of two randomly chosen individuals, for an update in each evolution time step is advocated by some researchers in the sense that cell interactions also depend on the physical distance and it can be then easily extended to the study of spatial dynamics. How- ever, the computational simplicity, leading to a wide range of analytical observations besides numerical results while capturing the same amount of essence of description for selection dynamics acting on asex- ual cell divisions as other processes do, is enough to have fascinated and convinced scientists the use of the frequency-dependent Moran process. The Moran process was originally introduced by an Australian statistician, P. A. P. Moran, in [31, 32] in 1958 to express population dynamics, especially random drift, in a nite population in the absence of frequency in describing tness. That is, the model in his introduction assumed that the each individual has a constant tness that is never aected by the number of its own groups’ or other groups’ tness during evolution. However, M. A. Nowak et al. insisted that how fast a mutant reproduces its ospring takes the frequencies of all dierent types into account and they modied it to be a frequency-dependent model. 1.2.1 Moranprocess: Abirth-and-deathprocess The Moran process is a birth-and-death model that describes population dynamics in a well-mixed popula- tion, in which all individuals compete and interact with each other for limited resources, without imposing a population structure. For a basic Moran process, one individual is chosen and reproduces its ospring of the same type with a probability that is proportional to its expected tness, and one individual, possi- bly the same one, is randomly chosen and eliminated in one evolution time step. This reproduction and 11 elimination in one time step thus results in a constant total population size during the whole evolution though the number of subpopulations are ever changing until it get absorbed to a homogeneous state. With this model, the number of at most two dierent types of subpopulations are aected and each group either increases or decreases at most by 1 in one time step, without both increasing or both decreasing, while the number of other types keeps unchanged. For example, consider a population of sizeN with two genotypes,A andB, in which every pair of two cells in one generation interact each other equally likely and their tness is evaluated according to the payo matrix,M, in (1.1). A state of the Moran process is determined by the number of subpopulations of one type, say typeA, since a constant total population size is maintained during the whole evolution process, and then the number of cells in the subpopulation, B, is automatically given by the dierence of the number of cells in the subpopulation,A, from the total population size,N. Knowing that the state of the Moran process can change at most by 1 in one time step, this denes a discrete-time 1dimensional Markov chain on a nite space,f0; 1; ;Ng, with a (N + 1) (N + 1) tridiagonal transition probability matrix. To be precise, leti be the number of individuals in the subpopu- lation,A, at the current step. LetT + i be the transition probability that the state,i, increases by 1 toi + 1, and T i the transition probability that the state decreases to i 1 from i. These are the only possible nonzero transition probabilities with the probability, T 0 i = 1T + i T i , that the system stays at the current state under the assumption of no mutation. Bothi = 0 andi = N are two absorbing states that correspond to a homogeneousBpopulation and a homogeneousApopulation, respectively, and these obviously lead toT + 0 = 0 andT N = 0, or in other words,T 0 0 = 1 andT 0 N = 1. At all other internal states, the transition probability,T i , at the state,i, is computed by considering the birth rate and the random 12 death. The birth rate is assumed to be proportional to the expected tness,f A i andf B i , ofA andB which interpret the expected payos, A i and B i , that are given by excluding the self-interaction: A i = a(i 1) +b(Ni) N 1 ; B i = ci +d(Ni 1) N 1 : (1.10) As introduced as a key mechanism by Darwin, it is common for the Moran process that natural selec- tion whose intensity parameter,w, ranges between 0 and 1 goes through tness function in a way that whenw = 0, it reduces to the background tness which is usually normalized to 1 for all types while whenw = 1, it is exactly the same as the expected payos in (1.10). Thus,w = 0 corresponds to neutral drift where all subpopulations of dierent types are equally likely at reproduction whereasw = 1 is where the payo, on which the researchers’ preliminary knowledge and background about the components of a game is reected, has the greatest inuence on the selection dynamic. Fitness functions are then dened by: f A i := 1w A +w A A i = 1w A +w A a(i 1) +b(Ni) N 1 ; f B i := 1w A +w A B i = 1w B +w B ci +d(Ni 1) N 1 ; (1.11) wherew A andw B are selection intensities ofA andB, respectively. Finally, the birth rates ofA and Bsubpopulations at the state,i, are given by if A i if A i +(Ni)f B i and (Ni)f B i if A i +(Ni)f B i , and these help to compute the transition probabilities for the Moran process: T + i = if A i if A i + (Ni)f B i Ni N ; T i = (Ni)f B i if A i + (Ni)f B i i N (1.12) 13 fori2f1; ;N 1g. It allows the Moran process to have an explicit (N + 1) (N + 1) tridiagonal transition probability,T , such as: T = 0 B B B B B B B B B B B B B B @ 1 0 0 0 0 0 T 1 T 0 1 T + 1 0 0 0 . . . . . . . . . . . . . . . . . . . . . 0 0 0 T N1 T 0 N1 T + N1 0 0 0 0 0 1 1 C C C C C C C C C C C C C C A : (1.13) Due to the facts that the tness functions in (1.11) become the expected payos in (1.10) whenw = 1 and that these tness functions are factors of transition probabilities in (1.12), it is necessary to assume that bothf A i andf B i are nonnegative for alli’s. These are satised if either selection intensities are small or all entries of a payo matrixP in (1.1) are nonnegative. It is common to implicitly assume thata;b;c;d 0 for the Moran process when the process with the whole range of selection is aimed to be analyzed. Whenever no one is tter than another at all states, that is, f A i = f B i for alli’s , we haveT + i = T i and neither genotype is favored for reproduction at all states. This is the case to which the neutral drift case belongs. Also, it is a random walk inZ + 0 though not simple since it still does not have a stationary property, in fact,T i 6=T j fori6=j. Thus, many known theories in a random walk inZ + 0 are ready-to-use, allowing fruitful analytical results available. Along the studies of evolutionary game theory by modeling a population with the Moran process since its reinterpretation as a frequency-dependent model, it has experienced generalizations in many aspects. For example, imposing a structure into a population by placing all individuals at vertices, and giving a weight to a directed edge, which represents a probability that an adjacent vertex is replaced by the ospring of the vertex from which the edge ies on a graph, a stochastic model can be designed to a various extents by varying either structures and weights as discussed in [29] in 2005. From their design, a 14 homogeneous population illustrated by the Moran process is considered as a particular case of a complete graph on which all edges have the same weight. They deduced that evolution is greatly dependent on the underlying structure and characterized graphs for which the quality of xation coincides with that of a homogeneous population when it is modeled by the Moran process. On the other hand, P. M. Altrock and A. Traulsen proposed another generalization of the Moran process by introducing tness-dependency in death rate as well as in birth rate for 2 2 game in [1] in 2009. According to their design, a cell is chosen for elimination not at random but at a probability proportional to inverse tness, and the transition probabilities,T i , in (1.12) are then modied as: T + i = ie +w b A i ie +w b A i + (Ni)e +w b B i (Ni)e w d B i ie w d A i + (Ni)e w d B i ; T i = (Ni)e +w b B i ie +w b A i + (Ni)e +w b B i ie w d A i ie w d A i + (Ni)e w d B i ; (1.14) wherew b 0 andw d 0 are the selection intensities at birth and death, respectively, and the tness is dened in the form of an exponential function. Leaving these for further studies, we model a nite well-mixed population using the tness-dependent Moran process in this thesis. 1.2.2 ConnectionofMoranprocesstotheadjustedreplicatordynamics Evolutionary game theory, which classically began with the replicator dynamics equations for an innite population, was later developed for a nite (stochastic) population. How the stochastic process for a nite population is related to the replicator system for an innite population as the population size increases was a natural question in this eld. A. Traulsen et al. rst proved in [54] in 2005 for 2 2 game that the replicator dynamic is the expected system of a local update stochastic model dened with a switching probability in (1.8) in the population size limit whereas the adjusted replicator dynamic is the mean system for the Moran process. They generalized this result for a general stochastic process with an arbitrary nite number of strategies, where mutation is also allowed in [53] in 2006. We briey discuss the derivation of 15 the adjusted replicator system from the Moran process in the population size limit for a 2 2 game and the maximal selection intensity in this section since it captures the essence of the argument in its simplest form, and the generalization to higher dimension is straightforward. The Moran process with two strategies is described at a microscopic level in terms of the probability, P i (), that the system is in the state,i, at time, with the transition probabilities: P i ( + 1)P i () =P i1 ()T + i1 +P i+1 ()T i+1 P i ()T i P i ()T + i ; (1.15) wherei is the number of agents in subpopulationA in a population of sizeN as before. This equation (1.15) is called the master equation and it describes the net change in probability from state,i, in one time step from time,. Letting x =i=N; (1.16) t ==N; (1.17) T (x) =T i ; (1.18) the equation in (1.15) is rewritten in terms of the probability density function, (x;t) := NP i (), as follows: (x;t +N 1 )(x;t) =(xN 1 ;t)T + (xN 1 ) +(x +N 1 ;t)T (x +N 1 ) (x;t)T (x)(x;t)T + (x): (1.19) 16 For a largeN, the Taylor expansion atx andt for both the probability density function,, and the transition probabilities,T , up to the rst order term inN 1 gives rise to the Kolmogorov forward equation: @ @t (x;t) = @ @x [a(x)(x;t)] + 1 2 @ 2 @x 2 [b 2 (x)(x;t)] (1.20) with the drift term,a(x) =T + (x)T (x), and the diusion term,b(x) = p (T + +T )=N. Applying the It^ o calculus on(x;t), the equation (1.20) is equivalent to a stochastic dierential equation, called the Langevin equation: dX t =a(X t )dt +b(X t )dB t ; (1.21) whereB t is the one dimensional standardWienerprocess andX t is the state of the system at timet. Noting that the diusion term,b(X t ), vanishes with the rate of 1= p N asN!1, the limiting system is solely determined by the drift term,a(X t ), as follow: _ x = lim N!1 (T + i T i ): (1.22) By equations in (1.12), it becomes _ x = lim N!1 f A i f B i (if A i + (Ni)f B i )=N Ni N i N = f A f B <f > (1x)x (1.23) and nally the adjusted replicator equation in (1.4) is recovered in the limit ofN from the master equation for the Moran process withx =x A . This derivation is available in [53] even when the Moran process is applied to a population having more than two strategies, sayd strategies. However, this stochastic process withd > 2 strategies is no longer a random walk on (Z + 0 ) d1 . In the same manner, the Fermi process dened by (1.6) and the local 17 update process dened by (1.8) are also shown to converge to the regular replicator system with a constant factor that only re-scales the time. 1.3 Moran process and a chemotherapy concentration applied to the selection Chemotherapy is a standard systemic clinical approach and is ubiquitous in various cancers as a treatment to alleviate disease by generally focusing on killing as many malignant cells as possible in order to reduce the size of tumor. The goal of this thesis is to develop adaptive chemotherapy protocols that prevent the deterministic adjusted replicator system from saturating the population with cancer cells, in the spirit of [36] and [59], then to test these strategies with the stochastic Moran process model, both using the single drug and multi-drug protocols. We will then highlight the advantages of using the adaptive chemotherapy schedules with respect to avoiding chemo-resistance of the tumor. Our approach is to view tumor growth as the outcome of the interactions of three subpopulations of cells: healthy cells (H), sensitive cancer cells (S) and resistant cancer cells (R) to a drug. The model we employ is the 2dimensional Moran process for a population of sizeN with three cell types,H,S and R. Since the Moran process keeps a constant population size during the whole evolution, the frequency of the dierent cell types determines the state of the system. We specify the number of healthy cells and the number of resistant cells as independent variables, so 2dimensional process, and leti be the size of healthy subpopulations andj the size of resistant subpopulations. Then for eachN, the set N :=f(i;j)2Z + 0 Z + 0 j0iN; 0jN;i +jNg (1.24) 18 denes the space of states for the Moran process with three strategies, and the number of elements, j N j = (N+1)(N+2) 2 , of N increases in the order ofO(N 2 ). When we refer to a state in N in (1.24), we interchangeably use (H;R) = (i;j) or (H;S;R) = (i;Nij;j). We assume that a nite population is well-mixed and every pair of cells in this population interacts and each cell receives a payo from that interaction according to a 3 3 payo matrix,A: A = 0 B B B B B B @ H S R H a 11 a 12 a 13 S a 21 a 22 a 23 R a 31 a 32 a 33 1 C C C C C C A : (1.25) From those equally likely interactions, a cell of typeX,X2fH;S;Rg, receives the expected payo, X i;j , at the state, (i;j), excluding self-interaction as it does for the 1dimensional Moran process in (1.10): H i;j = a 11 (i 1) +a 12 (Nij) +a 13 j N 1 ; S i;j = a 21 i +a 22 (Nij 1) +a 23 j N 1 ; R i;j = a 31 i +a 32 (Nij) +a 33 (j 1) N 1 : (1.26) The expected tness function,f X i;j , is dened, similarly to (1.11), by f X i;j = 1w X +w X X i;j ; (1.27) wherew X 2 [0; 1] is the selection intensity of the type,X. Modeling the role of a drug on a population committed to tumor regression/recurrence is then equivalent to shaping tness functions based on the scientic understandings about the relation between dierent types of cells when a drug is on or o. In 19 order to complete the description of our model, we introduce some key concepts in the following two subsections. 1.3.1 Prisoner’sdilemmagame A Prisoner’s dilemma (PD) game is a two-player-game, originally framed in classical game theory by M. Flood and M. Dresher while working at the RAND Corporation in 1950. What makes this game adilemma game, in the context of classical game theory, is that rational decisions to maximize their own payo leads them to choose a strategy that instead provides them with a lower payo than they would have gained if they chose the alternative. Later, the game was used in an iterated form to capture the essence of the emergence of cooperation [21, 38, 39, 42, 44, 55]. Being hard to be understood in the framework of classical game theory, this altruism is easily seen in nature, especially in animal worlds. For example, in some species, it is often seen that one gives an alarm for their peers when it encounters a predator in order to let them hide while exposing itself to danger. It is also seen in global climate-change that all countries take advantage from maintaining a stable climate, but any single country often hesitates to regulate CO 2 emission, thinking of maintaining this behavior as being more benecial to itself than it would be if all countries decided to reduce CO 2 emission. This scenario has been dubbed ’the tragedy of the commons’. Achievable altruism has been deeply studied in evolutionary theory with the variety of generalization from iteration to the number of strategies. In a classical Prisoner’s dilemma (PD) two-player-game, each player can choose to cooperate (C) or defect (D) in playing a game each round. Given a 2 2 payo matrix M = 0 B B @ C D C R S D T P 1 C C A = 0 B B @ C D C 3 0 D 5 1 1 C C A ; (1.28) 20 each receivesR when both cooperate whereas each receivesP when both defect. When two players with dierent strategies encounter, a cooperator receivesS and a defector receivesT . What makes this game in (1.28) a Prisoner’s dilemma game is the Prisoner’s dilemma inequality: T >R>P >S: (1.29) Under this assumption, all rational players will end up playingD since it gives a higher compensation in either cases when a partner cooperates or defects. This is called a Nash equilibrium. However, they would have received a higher payo, R, each if they cooperated. Any single combination of R, S, T , and P that satises the inequality in (1.29) denes a Prisoner’s dilemma game, but the canonical choice is when R = 3,S = 0,T = 5 andP = 1 as given in (1.28). In evolutionary game theory, strategies are interpreted as genotypes which are selected and passed to the next generation by natural selection depending on their frequencies in a population. Consider a innite population of cooperators (C) and defectors (D). Letx C be the frequency of cooperators andx D that of defectors, then we havex C +x D = 1. Letting the selection intensities,w C andw D , ofC andD be equal to 1 so that the payo matrix,A, has the most inuence on the tness functions dened similarly to (1.11) but in continuous way, the tness functions,f C andf D , are given by f C = (M~ x) 1 = 3x C + 0x D ; f D = (M~ x) 2 = 5x C + 1x D ; (1.30) where~ x = (x C ;x D ) | . Then the averaged tness function,<f >, in the entire population becomes <f >:=~ x | M~ x = 3x 2 C + 5x C x D +x 2 D (1.31) and it is straightforward to show that<f > has the minimum value, 1. 21 (a) (b) Figure 1.1: Defectors (D) playing a Prisoner’s dilemma game outcompete cooperators (C), grow in a S-shaped Gompertzian curve, and eventually saturate the entire population, eventually reaching a sub- optimal state according to the adjusted replicator system. For simulation, we setR = 3,S = 0,T = 5 and P = 1. (a) Initially nonzero defector subpopulations with the proportion, 0:01, keeps increasing over time in aS-shaped curve and eventually saturating the entire population; (b) The averaged tness function, < f >, is unfortunately minimized at that evolutionarily stable strategy,D, meaning thatx D = 1 is an asymptotically stable but sub-optimal state. Employing the adjusted replicator equation in (1.4) gives rise to _ x C = (1 +x C )(1x C )x C <f > (1.32) which implies that the deterministic system has xed points atx C = 0; 1. It is easy to check thatx C = 0 is asymptotically stable. Figure 1.1 illustrates that the proportion of defectors in an innite population, that initially takes only 1% of the total population, increases along a S-shaped Gompertzian curve and reaches plateaus at 1 [27], but the averaged tness function,< f > is minimized atx C = 0 being equal to 1. It is this Gompertzian growth law that makes the PD game useful as a tumor growth model. In fact, the system is driven to the stable homogeneous state of purely defectors whenever there initially exists the nonzero number of defectors, but that is a sub-optimal state. These two characteristics (Gompertzian growth and sub-optimal tness) allow us to adopt a PD game for describing interactions between healthy and cancerous cells in our model by letting healthy cells be cooperators and cancer cells defectors. 22 1.3.2 Costofresistanceandcompetitiverelease As mentioned earlier, one main factor that makes chemotherapy challenging is tumor recurrence followed by the initial tumor regression due to adaptive ability of a population, leading it to form resistance to a drug. This chemo-resistance is often explained by the ecological mechanism of the competitive release that was originally termed by J. H. Connell in [11] in 1961. The phenomenon unfolds when two species compete for a limited resource, with one species out competing the other, thereby suppressing the growth of the weaker species (i.e. the one with lower tness). When a change in the environment occurs that kills the dominating species, it releases the weaker one from the competition, thereby allowing it to grow. For competing cancer cells, in the absence of chemotherapy, the dominant population are the chemo-sensitive cells (higher tness), and the sub-dominant ones are the chemo-resistant ones (lower tness due to the cost of resistance). With chemotherapy, the sensitive cells are mostly killed o, which releases the resistant cells from the competition and allows them to grow, rendering the tumor resistant. Though this term had been supported later by a lot of in vivo experiments mainly in the eld of chemical ecology, cancer can also be modeled based on this notion in order to describe the chemotherapeutic response onto a population of healthy, sensitive and resistant cells. It is known that resistant cells to a drug as well as sensitive cells typically exist even before chemother- apy is carried out as a result of pre-existing mutations. Typically, however, resistant cells are present in a population only in small numbers due to the cost of resistance: it is inherently expensive for cancer cells to maintain resistance. In other words, resistant cancer cells pay a tness cost in the form of lowering their reproductive success. It has been demonstrated in laboratory experiments that resistant cells are less t and less ourishing than sensitive cells in the absence of chemotherapy [8, 13, 20]. The application of chemotherapy, however, then lowers the tness of the sensitive cell population, and the resistant cells begin to ourish. We will adopt these two notions to model the chemotherapeutic resistance of a tumor population. 23 1.3.3 Fitnesslandscape We quantify the previous notions in our choice for the payo matrix,A, in (1.25) to shape tness landscape of all types of cells together with drug concentration. We let the concentration, 0C(t) 1, of a drug, which is a function of time indicating the drug intensity rather than the absolute amount, enter the cell population through selection functions,w H ,w S andw R , in a following way: w H (t) =w 0 ; w S (t) =w 0 (1C(t)); w R (t) =w 0 ; (1.33) wherew 0 is a constant which we will set equal to 1 for computational convenience since it can be scaled out without changing qualitative features. This denition in (1.33) allows the selection functions of all types to be time-dependent functions and quantitatively establishes that the concentration function,C(t), becomes a controller of the system. When it is turned on, it lowers the tness of the sensitive cell population, leaving the others unchanged. A Prisoner’s dilemma game is an adequate model to describe the relations between healthy cells and cancerous cells in a way that under no chemotherapy, healthy cells are cooperators whereas cancerous cells are defectors since it possesses two important features when it is interpreted by the replicator system in the evolutionary game theory perspective: (i) initially nonzero portion of defectors eventually saturates the entire population, and (ii) that asymptotically stable homogeneous state of defectors is sub-optimal. Thus we impose the payo matrix,A, to satisfy the Prisoner’s dilemma inequalities in (1.29) such that a 21 >a 11 >a 22 >a 12 ; (1.34) a 31 >a 11 >a 33 >a 13 ; (1.35) 24 where (1.34) denes a PD game between healthy cells and sensitive cells whereas (1.35) denes a PD game between healthy cells and resistant cells. In addition, having a higher tness, sensitive cells are abundant under no treatment though a small number of resistant cells are present. However, when a drug is delivered, it lowers the tness of the sensitive cell population which, in turn, allows the resistant population to ourish (competitive release). In other words, the delivery of a drug, which is an articial interference to a population, plays the role of catalyst for the proliferation of resistant cells by lowering the tness of sensitive cells, releasing their competing agents (resistant cells), and allowing those resistant cells to use more resources, proliferate and eventually survive. These features, that are termed as the cost of resistance and the competitive release, convince us to adopt a PD game for interactions between sensitive cells and resistant cells where in the absence of treatment, resistant cells are cooperators whereas sensitive cells are defectors. We design the payo matrix,A, to satisfy the following prisoner’s dilemma inequality a 23 >a 33 >a 22 >a 32 (1.36) as well as a 21 >a 31 (1.37) to guarantee that sensitive cells have a higher tness than resistant cells at all time, that is, regardless of the distribution of a population when chemotherapy is absent. Also, the fact that a drug administration at a high dose causes sensitive cells to be less t than resistant cells is reected by requiring that the tness function, f R , of resistant cells is greater than the tness function, f S j C(t)1 , of sensitive cells under a constant high concentration. Knowing from (1.33) and (1.27) that the tness of sensitive cells is almost 25 equal to the background tness, which is set to be 1, when the constant high chemo-dose is nearly equal to 1, this is obtained if we further set a 31 ;a 32 ;a 33 1 (1.38) with all three not being equal to 1 simultaneously. For two drug model, we will design the tness of all types in an analogous manner based on PD games, the cost of resistance and competitive release by letting drug concentration functions enter through selection functions. This will be discussed further in Chapter 4. 1.4 Structureofthethesis In a nite population evolutionary game theory, what has received great attention by researchers are the state that a single mutant successfully enters a homogeneous population and the state that that single mutant survives and eventually takes over the whole population. The former is termed as invasion and the later is termed as xation. Between two, what distinguishes a stochastic process from a deterministic system in modeling a co-evolving population is the xation since it quanties a chance of observing a new scenario driven by natural selection which would not be seen if it is modeled by the replicator system assuming an innite population. Many researchers have devoted eort to nding either xation probability or xation time with many dierent stochastic models, resulting in fruitful analytical results. In Chapter 2, we briey go through the derivation of the exact formula of xation probability for stochastic processes with 2 strategies. In classical evolutionary game theory assuming an innite popu- lation, the evolutionary stable strategy (ESS) is well established as a state that survives and persists, hence is of great interest. However, this concept underwent modication in a nite population and we mention this evolutionary stable strategy (ESS N ) in a nite population, which is dened in terms of both xation probability and invasion. We also refer to both one well-known approximation to the rate of evolution, 26 which is the multiplication of xation probability and population size, under weak selection, and how the denition ofESS N is sharpened with the use of this approximation. We point out that this approximation, in the form of the truncated Taylor expansion with the degree of 1, is only valid under an extraordinarily narrow weak selection window with the example of repeated Prisoner’s dilemma game, and precisely give the upper threshold of those weak selection. We also show that the regular linear Taylor approximation becomes a better approximation over those limited region of weak selection. Moreover, we suggest a global approximation to the rate of evolution usingBernsteinpolynomial, which is dened on the whole selection interval, and explains its eciency on a strong selection limit. In Chapter 3, we adopt the frequency-dependent Moran process model with three strategies in order to model the role of single drug chemotherapy acting on a nite population of healthy (H), sensitive (S) and resistant (R) cells, all participating in tumor development. We shape the tness functions of cells as discussed in section 1.3.3. Knowing that the xation probability for the Moran process with three strategies does exist but is, so far, only expressed as a solution to a linear system, proved by E. M. Ferreira and A. G. M. Neves in [15], where the coecient matrix is order ofO(N 4 ) for the population sizeN, we numerically compute xation probability atj N j many points distributed in phase space. Visualizing these xation probabilities asN increases, we show how the stochastic version of competitive release converges to the basins of attraction for the adjusted replicator equation. Moreover, as a second part, the adaptive chemotherapy schedule associated with the adjusted replicator equations, developed in [36] but with the regular replicator equations, is applied to the stochastic Moran process to examine the chemotherapeutic response. It is obvious that applying this chemotherapy schedule to the stochastic process will fail in trapping the wiggly trajectory in a closed loop though it is designed to do for the expected deterministic system. It is for the randomness, driven by the nite population size, leads not only one single realization but also the averaged system of many realizations to settle on a neighborhood instead of the exact initial state after each evolutionary cycle. We quantify how badly this schedule works or how fast its use on the 27 Moran process is not justied as increasing the number of cycle. We later evaluate this adaptive schedule in terms of tumor volume, that is quantied as a sum of cancerous cells, by comparing it to the standard clinical approaches that have the same amount of the total dose of toxin: the maximum tolerated dose schedule (MTD) and the low-dose metronomic schedule (LDM). In Chapter 4, we generalize Chapter 3 by introducing an additional drug, where the two drugs act asymmetrically, meaning that they have dierent response rates. Following [30], we equip our model with a drug interaction parameter,1 e 1, in order to distinguish antagonistic (e < 0), additive (e = 0) and synergistic (e> 0) drug interactions from each other. We consider a nite population that is composed of sensitive (S) cells to both drugs, a group of cells (R 1 ) that are sensitive to drug 1 but resistant to drug 2, and a group of cells (R 2 ) that are sensitive to drug 2 but resistant to drug 1. Shaping the tness functions of these three types of cancer cells analogously to but relatively more complex than that of a single drug model, the evolution of a population of cancer cells is modeled by the Moran process with three strategies. After we develop an adaptive chemotherapy schedule associated with the adaptive replicator equations for each drug interaction case that use up the similar total dose of toxin from both drugs, we apply each schedule to the Moran process. For the evaluation of these schedules, we separately adopt the tumor-growth model, which is the nonlinear deterministic dierential equation, to describe tumor volume. With these two models, we study how successfully the adaptive chemotherapy schedule associated with the adjusted replicator dynamic works and how long the application of this chemotherapy can be justied, just as done in Chapter 3. Evaluation of this adaptive two drug schedule is then made either across or within drug interactions by comparing it to two standard clinical approaches such as MTD and LDM having the same total dose. Finally in Chapter 5, we will discuss future potential research directions which would generalize our models in several important ways. 28 Chapter2 FixationprobabilityforMoranprocesseswithtwostrategies Begun by M. A. Nowak et al. in modeling a nite population that play 2 2 game with the frequency- dependent Moran process, xation has been intensively studied from nding its probability to nding the time that the underlying stochastic process takes until its occurrence, across many stochastic processes [2, 29, 40, 51, 56]. Fixation refers to a state that a single (or, in general, a nite number of) mutant is favored by selection replacing resident subpopulations and eventually takes over the entire population. The reason why this concept attracted many researchers is that this is what makes stochastic processes distinguished from deterministic models and what quanties the randomness driven by niteness in population size. For example, as shown in [51] for a 2 2 game, the Moran process for a nite population drives eight selection scenario whereas the deterministic replicator system for an innite population leads to four selection ones, and those extended categorization was made based on xation probability. Also, the well- established concept,evolutionarystability, in an innite population was modied by M. A. Nowaketal. in [40] by introducing xation probability as an extra condition besides the original invasion condition when adopting this concept into a nite population. 29 2.1 Fixationprobabilityandevolutionarystability Fixation probability has an exact formula for a 2 2 game though complex, and it was obtained by ex- ploiting that 1 dimensional Moran process with two strategies can be considered as a random walk on a subset ofR and using well-known ready-to-use theorems. In fact, the closed form is available for any birth-and-death processes as well as the Moran process. However, in general, the Moran process withd many strategies is not a random walk inR d1 whend 3, and this makes the computation challenging for a higher dimension. As far as we know, E. M. Ferreira and A. G. M. Neves were the rst who recently proved the existence of the xation probability ford = 3, and they represented it as a solution to a sys- tem of linear equations by formalizing the xation probability as a solution of discretized 2dimensional Dirichlet problem for the Laplace equation in [15] in 2020. We will discuss more about it in the next chap- ter, but for a while we mainly focus on the xation probability for the Moran process with two strategies in this chapter. 2.1.1 Computationofxationprobability Consider the Moran process with two strategies,A andB, for a nite population of sizeN, all interacting equally likely each other and receiving a payo according to the payo matrixM in (1.1). We dene a state of the Moran process by the number,i, of cells in the subpopulation,A,i = 0; 1; ;N. For eachi, let A i be the probability thati mutants of typeA replace allBmutants and take over the population, and similarly B i forB. Since all states of the Moran process are transient, the stochastic system eventually achieves the homogeneity of a population in either allA’s or allB’s, leading to B i = 1 A i . Thus it is sucient to compute A i for a 2 2 game. We will use A i and i interchangeably for convenience. 30 If mutation is not allowed, the Moran process with two strategies has two absorbing states,i = 0 and i =N, at which xation probability is complete as follows: 0 = 0; N = 1: (2.1) For an intermediate states,i = 1; ;N 1, the state can have a jump only to eitheri + 1 ori 1 besides staying unchanged. Each jump is made at a probability,T + i andT i , respectively. Then it allows to have a following recurrence relation in terms of xation probability, i , and transition probabilities,T i : i =T i i1 +T 0 i i +T + i i+1 ; (2.2) where the probability,T 0 i , that the stochastic system remains unchanged is equal toT 0 i = 1T + i T i . Dening variables, i := i i1 and i :=T i =T + i , the equation (2.2) is re-written: 1 = 1 ; i+1 = i i (2.3) fori = 1; ;N 1, and this is equivalent to i+1 = 1 i Y k=1 i : (2.4) Summing (2.4) overi gives rise to P N i=1 i = ( 1 0 ) + ( 2 1 ) + + ( N N1 ) = N 0 = 1 and it leads to 1 = 1 1 + P N1 j=1 Q j k=1 k : (2.5) 31 Using the same summation of (2.4) but up to as less number of terms as desired, it implies i = 1 (1 + P i1 j=1 Q j k=1 k ), and the general xation probability, i , ofA starting from the state,i, is driven as follows: i = 1 + P i1 j=1 Q j k=1 k 1 + P N1 j=1 Q j k=1 k (2.6) with equating P i1 j=1 Q j k=1 k to 0 wheni = 1. The derivation of xation probability in (2.6) has not made use of the underlying stochastic process except for that it is a birth-and-death process so that a state can move at most by1 in one time step. For this reason, the expression in (2.6) is valid for any 1dimensional birth-and-death processes with no mutation assumption. We end this subsection by specifying this xation probability for the Moran process for the later use. From (1.12), the variable, i := T i =T + i , can be explicitly written in terms of tness functions ofA andB for the Moran process, and the xation probability in (2.6) is nothing but i = 1 + P i1 j=1 Q j k=1 f B k =f A k 1 + P N1 j=1 Q j k=1 f B k =f A k : (2.7) 2.1.2 Evolutionarilystablestrategyinanitepopulation Corresponding to the Nash equilibrium, which is a strategy with which no plyer improves its reward by switching to another, in classical evolutionary game theory, one of the main interests in evolutionary game theory is to nd a strategy, pure or mixed, persistent to natural selection through time in a coevolutionary population. Such strategy is what is called an evolutionarily stable strategy (ESS), originally invented by M. Smith in [46, 47, 48, 49]. Assuming an innite population,ESS is a strategy that for an innitesimally large population playing this strategy, an innitesimally small population playing a dierent strategy is 32 not able to invade into the nearly homogeneous population. Precisely in a 22 game with a payo matrix M in (1.1), a strategy,B, is said to be evolutionarily stable if either (i)d>b or (2.8) (ii)d =b and a<c: (2.9) This condition is made by letting the tness of the subpopulation, A, smaller than that of B when the population is nearly composed of individuals of typeB, so that natural selection opposes the spread of innitesimally small fractions ofA in innitely large populations ofB. As evolutionary game theory developed and the needs for more realistic models describing a nite coevolving population arose, this concept of evolutionary stability had been challenged to be adopted to a nite population. D. Foster and P. Young suggested a new concept of stability, which they called stochas- tically stable equilibrium (SSE), for a stochastic system that is designed using the Wiener process added to the deterministic replicator equations in [17] in 1990. By deningSSE as a long-run viability, they convinced that attractors in dynamical system are not sucient to be an alternative due to the dependence of the limiting behavior of the system on the initial distribution. In fact, they denedSSE to be a state at which the stochastic system asymptotically stays within its small neighborhood with a positive probability as the stochastic eects decreases, where the decreased eects is measured as the shrinking variance of the Wiener process. Recognizing that a single mutant, successfully invaded a resident population, is possibly able to take over the entire nite population, M. A. Nowak et al. modied the concept ofEES, that M. Smith had suggested, in [40] in a way that selection opposes both invasion and xation, and began to be widely used 33 in this eld. Denoting it byESS N , a strategy,B, is said to be evolutionarily stable in a nite population of sizeN for a 2 2 game in (1.1) if followings are satised: (i)b(N 1)<c +d(N 2); (2.10) (ii) A 1 < 1 N : (2.11) The condition in (2.10) is equivalent to f A 1 < f B 1 and guarantees that selection opposes A invading a resident population ofB whereas the condition in (2.11) describes thatA is not favored taking over the entire population by natural selection. The quantity, 1 N , in (2.11) is the xation probability ofA when the selection is absent for both types such that the coevolutionary population evolves exclusively under the neutral drift. Knowing k = 1 in (2.5) underw A = w B = 0, this quantity is easily computed, and we denote it by w=0 1 to specify no selection intervention. 2.2 Localapproximationtoxationprobabilityunderweakselection In Darwinian evolution theory, all species of organisms competing each other for the same limited re- sources survive and reproduce being favored by the natural selection if tter, where how tter is slightly dierentiated from each other compared to the neutral drift. For this reason, the weak selection limit got a great attention, and this is where abundant analytical results were obtained. Fixation probability for a birth-and-death process has an exact formula for a 2 2 game, but its complexity drove many researchers to consider an approximation instead under weak selection. For a two-player-game in (1.1), one well- known local approximation to the rate of evolution,N 1 , which is multiplying the xation probability by 34 the population size, was suggested in [40] in the form of the rst order truncated Taylor expansion with the shared selectionw :=w A =w B as follows: N A 1 1 1 (N)w=6 ; (2.12) where =a + 2bc 2d and = 2a +b +c 4d. We denote the rational function ofw in the right hand side of (2.12) by N or by N (w) to specically the independent variable,w. 2.2.1 Equivalenceofstochasticprocessesunderweakselection Using the local approximation in (2.12) or equivalently A 1 1 [1 (N)w=6]N (2.13) under weak selection, the second condition ofESS N in (2.11) for a nite population is sharpened, and the denition ofESS N under weak selection in Section 2.1.2 is replaced by: (i)b(N 1)<c +d(N 2); (2.14) (ii)N <: (2.15) M. A. Nowaketal. in [40] linked the second condition to an unstable internal xed point of a coordina- tion game, which is dened with two constraints,a>c andb<d. According to the replicator equations, this game has an internal equilibrium,x := db abc+d such that the deterministic system is driven to the homogeneous state of allB’s whenever the initial frequency ofA is less thanx whereas it is attracted by 35 x = 1 (allA’s) if the initial proportion ofA exceedsx . It is straightforward that for a large but niteN, the inequality (2.15) is equivalent toa + 2b<c + 2d, and this can be re-written in terms ofx as x > 1=3: (2.16) Surprisingly, this says that theA mutant is opposed replacing the resident populationB by the natural selection if the frequency ofA is higher than 1=3 whereas favored otherwise. This is what is termed as a ’1=3law’, and it is the key condition that is used to show the equivalence of any two stochastic processes. The ’1=3 law’ for the Moran process was derived by M. A. Nowak et al. in [40] in 2004, and to this model the equivalence of the Wright-Fisher model in [25], of the local update process in [57], of the Fermi process in [58, 61], and of the generalized Moran process in [1] were shown in the sense that all of these processes hold the ’1=3 law’ under weak selection, where its derivation for the Moran process exploits the local approximation to xation probability in (2.13). In fact, ’1=3 law’ was proven in [28] to be valid for any processes under weak selection as long as it is in the domain of Kingman’s coalescence, and all the stochastic processes mentioned earlier belong to this set, thus being equivalent to each other. However, the violation of ’1=3 law’ was also discussed when the extra condition,Nw >> 1, is introduced in [57], and when the xation probability is approximated by a Taylor expansion with a higher degree than 1 in [61]. We end this subsection by discussing the relationship between the traditionESS for an innite popula- tion andESS N for a nite population under the weak selection limit with two extreme cases of population size. ForN = 2, the conditions, (2.14) and (2.15), forESS N reduce simply tob < c, and it implies that 36 ESS N is a neither necessary nor sucient condition forESS, referring to (2.8) and (2.9). On the other hand, whenN is nite but large, the conditions, (2.14) and (2.15), are re-written by (i)b<d; (2.17) (ii)x < 1=3; (2.18) meaning thatESS N is a sucient but not necessary condition forESS. In summary,ESS N recovers the traditionalESS asN increases, and this is compatible with the fact that a stochastic process converges to the deterministic (adjusted) replicator equations in the increase ofN. 2.2.2 Toolimitedvalidityofthelocalapproximation Despite of the height of its fame as a local approximation under weak selection, the approximation, N , in (2.12) to the rate of evolution,N 1 , in the form of the linear truncated Taylor expansion has its validity on an extraordinarily small subinterval of [0; 1]. In fact, it is limited from above and the upper threshold varies depending on the population size, N. This limited validity is due to the fact that for eachN the local approximation, N (w), is a rational function in selection,w, and this takes negative values on the majority part of the whole interval [0; 1] of selection, being inappropriate to be used in approximating a positive function. Taking that the rate of evolution, N 1 , is a positive function and ranges 0 toN into account, an approximation is acceptable at least if it shares the same range, [0;N]. Forcing 0 N N gives rise to an upper threshold ofw wheneverN > 0: w 6(N 1) (N)N : (2.19) Denoting the right hand side of (2.19) by TRS N , we note that the upper threshold, TRS N , is order of O(1=N 2 ) and it shrinks to 0 asN increase if> 0. 37 In order to see how ineective the local approximation, N , is, we come back to a Prisoner’s dilemma game dened for two players in (1.28) with the Prisoner’s dilemma inequality in (1.29). First, we note that the strategy,D, is an evolutionarily stable strategy for both a nite and innite population. Recalling (2.8) and (2.9), it is obvious thatD is an evolutionary stable strategy for an innite population sinceP >S. For a nite population modeled by the Moran process, it is straightforward to obtain (2.10), and we consider the factor, k , of products in (2.6) for the second condition ofESS N in (2.11). For the Moran process, we have k =f D k =f C k , and in the strongest selection temperature for bothA andB it is equivalent, for a PD game, to k = Tk +P (Nk 1) R(k 1) +S(Nk) : (2.20) Using the PD inequality in (1.29), it is easy to check k > 1 and it implies that C 1 < 1 N , meaning that aC mutant is opposed by selection replacingD and taking over the residentD populations. As a result,D is also an evolutionarily stable strategy for a nite population. Thus, it is hard to understand the emergence of cooperation if a PD game is played only one time. For this reason a repeated PD game for a nite/innite number of times was studied extensively in the eld of both classical and evolutionary game theory with the desire for the descriptive paradigm of the emergence of cooperation. In a repeated PD game, players choose a strategy betweenC andD each round and a series of strategies during the set number of times establishes a strategy for this game. For example,"Alwayscooperate(ALLC)" is a strategy that a player cooperates all the time whereas"Alwaysdefect(ALLD)" is a strategy with which a player joins a game by defecting all the time. Obviously, a player withALLD outcompetes another with ALLC, receiving a higher payo each round. In fact,ALLD is a strict Nash equilibrium. Besides these, one strategy of direct reciprocity for a repeated PD game is what is known to be Tit-For-Tat (TFT), with which a player cooperates at the rst round and does whatever its opponent did in the previous round from the second round. In other words, a player with TFT is just as generous as the opponent is but never forgives once the opponent betrays. This strategy, TFT , was rst invented by A. Rapoport and 38 won the iterated PD game tournament held by R. Axelrod in 1984 [3] though it was only the best among those of the participants. Later, the winner was replaced by a strategy taking a probability into account, calledGenerous Tit-For-Tat (GTFT), and a player withGTFT forgives with a positive probability and keeps cooperating even if an opponent defected in the previous round. We consider an iterated PD game with the strategies,TFT andALLD, duringn rounds. Then the payo matrix betweenTFT andALLD is as follows: 0 B B @ TFT ALLD TFT nR S + (n 1)P ALLD T + (n 1)P nP 1 C C A : (2.21) It is straightforward in an innite population thatALLD is anESS in an innite population sincenP is always greater thanS + (n 1)P for alln sinceS <P . However, ifnR>T + (n 1)P is additionally assumed, or equivalently n> TP RP ; (2.22) thenTFT also becomes anESS. The condition, (2.22), is established by forcing eachTFT -individual to receive a higher expected payo,nR, than the expected payo,T +(n1)P , that eachALLD-individual receives whenTFT is abundant in a population. With this extra condition that requires a large enough number of repetition of a PD game, the invasion of innitesimally smallALLD-subpopulation is opposed by selection to an innitesimally largeTFT -subpopulation. In this homogeneous state of allTFT ’s, all individuals cooperate and the ultimate piece is achieved since no one defects in a population, and this supports the emergence of cooperation in an innite population. We now turn our gear to a nite population under weak selection, particularly on [0; 0:01]. Unlike being anESS in an innite population, the strategy,ALLD, is no longer anESS N for almost everyN under 39 weak selection. In order to specify the critical value of the population size over which the evolutionarily stability ofB is not guaranteed under weak selection, we consider two conditions in (2.14) and (2.15). The condition in (2.14) is satised ifALLD has a higher expected tness thanTFT when there is only one TFT mutant in a population, which can be written asf TFT 1 < f ALLD 1 . Assuming the shared selection strength,w :=w TFT =w ALLD , forTFT andALLD, it is then equivalent to show: [S + (n 1)P ](N 1)<T + (n 1)P +nP (N 2); (2.23) and this is rearranged to (PS)(N 2) + (TS)> 0: (2.24) The Prisoner’s dilemma inequality in (1.29) allows the inequality in (2.24) to hold for allN 2. Hence, selection opposes a singleTFT mutant invading the residentALLD population for allN. In fact, this invasion inability of a singleA-individual to theB-populations holds regardless of selection strength as long asTFT andALLD share the same selection strength. On the other hand, the xation condition in (2.14) is equivalent to [nR + 2ST (n 1)P ]N < 2nR +S +T 2(n + 1)P (2.25) under weak selection. With the choice of parameters,T = 5,R = 3,P = 1,S = 0 andn = 10, (2.25) leads toN < 43=14. Selection favors a singleTFT mutant taking over the entire population for allN but N = 2 andN = 3. Thus the strategy,ALLD, is not anESS N ifN 4 under weak selection, and this is the new scenario that the deterministic system does not result. Recall that the xation condition in (2.14) ofESS N is obtained from adopting the local approximation, N , in (2.12) to the rate of evolution,N A 1 , under weak selection. For the numerical purpose, we plot the 40 (a) (b) (c) (d) Figure 2.1: The local approximation, N , in (2.12) to the rate of evolution,N 1 , of a singleTFT mutant under weak selection has limited validity with the upper threshold of selection strength, TRS, which is order ofO(N) for the repeated Prisoner’s dilemma game betweenTFT andALLD. (ParametersT = 5, R = 3, P = 1, S = 0 and n = 10 in all panels) (a) N = 10, TRS = 0:05567010; (b) N = 10 2 , TRS = 0:00437730; (c)N = 10 3 ,TRS = 0:00042946; (d)N = 10 4 ,TRS = 0:00004287 exact value,N TFT 1 (=N 1 ), in red line together with its approximation, N , in blue line for the repeated Prisoner’s dilemma game under weak selection, precisely on [0; 0:01], forN = 10; 10 2 ; 10 3 ; 10 4 with the choice of T = 5, R = 3, P = 1, S = 0 and n = 10 (Fig.2.1). Comparing to the rate of evolution, N w=0 1 (= 1) (black solid line), in the random drift case,N 1 in Fig.2.1 shows that the xation ofTFT is favored by selection under all weak selection on [0; 0:01] for all N = 10; 10 2 ; 10 3 ; 10 4 though how fast aTFT mutant takes the entire population depends on the total population size,N. However, more importantly, the well-known approximation, N , that performs well on [0; 0:01] for an extraordinarily small population size such asN = 10 becomes an invalid estimate on an increasing subinterval of [0; 0:01] asN increases since N approximatesN 1 by a negative number oncew passes the vertical asymptote. The exact upper threshold, TRS N (blue dashed line), of weak selection below which the use of N is 41 acceptable was computed for each N in (2.19) and is also added to Figure 2.1. It is clearly shown that TRS N shifts towards 0 as N increases, making N work on a too tight window on which the regular linear Taylor expansion performs even better. 2.3 Bernsteinapproximation:globalapproximationtoxationprobability As seen in the previous section, the local approximation, N , to the rate of evolution,N 1 , is often used under weak selection though it has limited validity, particularly for a largeN. However, unlike the one in [56] for the Fermi process using the fact that its tness function is designed to cover any nonnegative selection strength, none of global approximations on [0; 1] for the Moran process has been discussed so far. We suggest one global approximation by using the Bernstein polynomial that was originally constructed by S. Bernstein in [4] in 1912 and used to prove Weierstrass approximation theorem, which states that a continuous function on a closed interval can be uniformly approximated by a series of polynomials. Given a function,f(x), dened on the closed interval, [0; 1],theBernsteinpolynomialofdegreedofthe functionf(x) is dened by: d X k=0 f( k d ) d k x k (1x) dk : (2.26) Denoting (2.26) byB d (f), it is proven that iff is continuous on [0; 1], thenB d (f) converges uniformly to f on [0; 1] asd!1 with an error bound: jf(x)B f d (x)j (1 + 1 4 d 2 )!(f;d 1=2 ); (2.27) where the modulus,!(f;d 1=2 ), of continuity off is dened by !(f;d 1=2 ) = sup x;y2[0;1];jxyj< fjf(x)f(y)jg: (2.28) 42 Expanding the formula in (2.26), it is easy to check thatB d (f), which is a linear combination of polynomials of degree at mostd, is the weighted average of functionsx;x 2 ; ;x d . From its construction, obtaining the global approximation,B d (f), precisely requires to knowd + 1 many function values at selection values, w = k d , k = 0; 1; ;d, that equally partition the interval, [0; 1], into d subintervals whereas getting the local approximation, N , is based on the stronger dierentiability condition ofN 1 but at only one locationw = 0. For a nite population of individuals playing a two-player game dened in (1.1), we denote byr i the rate of evolution, N i , for i = 1; ;N 1. Then this approximation object, N i , can be globally approximated on [0; 1] by the Bernstein polynomials in (2.26) as follow: B d (r i ;w) := d X k=0 r i ( k d ) d k w k (1w) dk : (2.29) Furthermore, if the nite population is modeled by the Moran process, then it presumes that all entries of a payo matrix are positive since it otherwise brings in a negative transition probability. This positiveness guarantees that the tness functions in (1.11) for the Moran process are continuous on [0; 1], and so is i in (2.7). This completes the continuity ofr i , and we have that for eachi,B d (r i ) uniformly converges tor i asd!1. We visit again the repeated Prisoner’s dilemma game between TFT and ALLD, dened in (2.21), where this game is repeated more than TP RP times to guarantee that selection opposes invasion ofALLD into a homogeneousTFT -population, having the emergence of cooperation in mind. In order to model a nite population playing this game by the Moran process, it is necessary to assume that all entries of 43 (a) (b) (c) (d) Figure 2.2: The Bernstein polynomials,B d (r 1 ), of degree,d, approximates the rate of evolution,r 1 :=N 1 , in (2.12) in the whole region of selection strength, being more accurate with the increase ofd. For simu- lation, we choose parameters,T = 5,R = 3,P = 1,S = 0 andn = 10, in all panels for the repeated Prisoner’s dilemma game betweenTFT andALLD. (a)N = 10. Bothr 1 andB d (r 1 ) are increasing func- tion on [0; 1] being greater than the rate of evolution,r w=0 1 :=N w=0 1 , under the neutral drift. Selection favors TFT replacing ALLD for all selection strength; (b) N = 10 2 . The same interpretation as for N = 10 is made but with a higher rate of evolution; (c)N = 10 3 . Selection favorsTFT replacingALLD for weak selection while it does not for strong selection according to either the exact or approximate functions withd > 1; (d)N = 10 4 . Selection opposesTFT replacingALLD for all selection strength according to the approximationB d (r 1 ) for alld, however, it does not if selection is extraordinarily weak according to the exact function. (2.21) are positive to allow the whole range of selection for this selection dynamics. Using the Prisoner’s dilemma inequality in (1.29), it is equivalent to P > 0; S + (n 1)P > 0: (2.30) 44 With the choice,T = 5 ,R = 3,P = 1,S = 0 andn = 10, the inequalities in (2.30) are satised and the uniform convergence of the Bernstein polynomial,B d (r i ), to the rate of evolution of a singleTFT - mutant is implied for this repeated Prisoner’s dilemma game. With this selection of parameters, we plot the Bernstein approximation, B d (r 1 ), with a set of degrees, d = 1; 2; 4; 5; 10; 20; 25; 50, for each N = 10; 10 2 ; 10 3 ; 10 4 along with both the exact function ofN 1 in black solid line and the rate of evolution, N w=0 1 (= 1), under the neutral drift in black dashed line (Fig.2.2). For eachN, it is clearly seen thatB d (r 1 ) generally better approximatesr 1 asd increases and it is also supported by the reduced error,jB d (r i )r i j (Fig.2.3). WhenN = 10 4 ,N 1 is nearly equal to 0 on the majority of selection except for weak selection on which it slightly increases nearw = 0 and is followed by a sudden drop to almost 0. On the other hand, all the Bernstein polynomials,fB d (r 1 )g d , take the value, 1, atw = 0 and they monotonically decrease to 0 asw increases though the higherd, the faster decrease. These approximations lead us to interpret that a singleTFT -mutant is opposed by selection replacing theALLD-mutants regardless of selection intensity whereas it is implied by the exact function,N 1 , that selection opposes the xation ofTFT on [0; 1] but on an extraordinarily small weak selection window. This interpretation is compatible with the fact that the deterministic adjusted replicator system has a bistability for this game so that the system is driven to allALLD’s if the initial frequency ofTFT does not exceeds a critical proportion. Precisely, assuming the same nonzero selection intensity forTFT and ALLD, the deterministic system is driven toTFT as long asf TFT >f ALLD , wheref TFT f ALLD for this game is given by f TFT f ALLD = [nRx + (S + (n 1)P )(1x)] [(T + (n 1)P )x +nP (1x)] = [(n(RP ) (TP )) + (PS)]x (PS) (2.31) 45 withx :=x TFT denoting the frequency ofTFT in an innite population. Using the Prisoner’s dilemma inequality and the assumption on the number of iteration in (2.22), both (n(RP ) (TP )) + (PS) andPS are positive, and we havef TFT >f ALLD if and only if x> PS n(RP ) (TP )) + (PS) =: ^ x (2.32) and we denote the internal equilibrium by ^ x. Our choice of parameters gives rise to ^ x = 1 17 . The initial state,i = 1, in a nite populations withN = 10 4 corresponds to the state,x 0 = 10 4 , in an innite popu- lation. The initial state,x 0 , is denitely less than ^ x, and it leads the deterministic system to allALLD’s by (2.32). Noting that the adjusted replicator equation is the asymptotic system of the Moran process in the large population size limit, it implies that for a nite but large population selection is not likely to favor a singleTFT -mutant taking over the entire population if the initial proportion ofTFT is less than the critical value, ^ x = 1=17. This is what is shown in Figure.2.2d. Every Bernstein approximation,B d (r 1 ), of any degree is strictly less thanN w=0 1 (= 1) for allw6= 0, meaning that selection opposes the xation of TFT . However, things interesting take place with the exact function,N 1 , on an absolutely tiny interval includingw = 0. On that extremely small window, selection slightly favorsTFT replacingALLD, and this is a new feature that is driven by stochasticity from the nite population size and was not captured by the deterministic system. This new feature continues on a wider weak selection interval as N decreases due to the stronger stochasticity. WhenN = 10 3 , r 1 is greater than 1 on weak selection, [0; 0:43), whereas it is less than 1 on strong selection, (0:43; 1]. Thus, the xation ofTFT is reached with the extended weak selection intensity with this smaller population size,N = 10 3 , compared to that withN = 10 4 (Fig.2.2c). This can be similarly translated byB d (r 1 ) as long asd> 2 so thatB d (r 1 ) crosses the horizontal line ofN 1 (= 1), with the more accurate interpretation for a higherd. Eventually, whenN becomes as small asN = 10 or 46 N = 10 2 , we observe that selection favors the saturation ofTFT , though the rate is higher withN = 10 2 , according to eitherr 1 andB d (r 1 ). (a) (b) (c) (d) Figure 2.3: The error, produced by the Bernstein polynomial,B d (r 1 ), in approximating the rate of evolu- tion,r 1 =N 1 , is reduced as much as desired by increasing the degree,d, for eachN although a relatively higher degree is required for an intermediate value ofN. Despite of its universality as a global approxima- tion on [0; 1], it does not behave well particularly under weak selection, producing a relatively high error. (ParametersT = 5,R = 3,P = 1,S = 0 andn = 10 in all panels) (a)N = 10; (b)N = 10 2 ; (c)N = 10 3 ; (d)N = 10 4 In summary, the Bernstein polynomial,B d (r 1 ), plays the role of a global approximation to the rate of evolution,r 1 = N 1 , on [0; 1]. Its performance improves asd increase for eachN so that the error can be reduced as much as desired, and the major renement is achieved for strong selection. However, the weakness of this approximation is on weak selection, especially nearw = 0, where it produces a relatively bigger error compared to that on strong selection intensity. Despite of this bigger error nearw = 0, this error is overall smaller for either a small or a large population size, as shown in Figure 2.2a and Figure 2.2d. On the other hand, the error remains big nearw = 0 even with a higher degree for an intermediate value ofN, particularly forN at which selection dynamic for xation changes at some mild selection strength in 47 (0; 1) as shown in Figure 2.2c. In fact, the numerical simulation for the repeated Prisoner’s dilemma game with the same choice of parameters results that the regular linear Taylor expansion outperforms both the Bernstein approximation,B d (r 1 ), and the local approximation, N , on the weak selection window over which N is valid with the constraint in (2.19). 48 Chapter3 Stochasticsingledrugchemotherapymodel For our stochastic single drug chemotherapy model for a nite population with health, sensitive and re- sistant cells by the Moran process, we specify the payo matrix, A, in (1.25) satisfying (1.34)-(1.38) as follows: A = 0 B B B B B B @ H S R H 3 1:5 1:5 S 4 2 2:8 R 3:9 1 2:2 1 C C C C C C A (3.1) Each state, (i;j), of the Moran process represents a population withi healthy cells,Nij sensitive cells andj resistant cells, and it is assigned to a lattice point in a phase space that is 2simplex,S 2 . Each vertex of the phase space represents a homogeneous population and is denoted by its type, and at a point in each side ofS 2 the population consists of only two types of subpopulations.inumbering goes next to the side,HS, andjnumbering lies on the bottom of the side,SR (Fig.3.1). For example, WhenN = 10, there are 66 distinct states, and 30 of them consist of less than three dis- tinct subpopulations and they are assigned to boundary points (black dots in Fig.3.1) inS 2 . The remaining states correspond to the population with no vanishing subpopulations and they are mapped to interior points (red dots in Fig.3.1). 49 Figure 3.1: The space set, N , of the Moran process withH, S andR is indicated as dots in the phase space,S 2 , whenN = 10. Each state, (H;S;R) = (i;Nij;j), is assigned to a lattice point inS 2 , wherei-numbering goes next to the side,HS, andj-numbering lies on the bottom of the side,SR. The space set, N , dened in (1.24), has exactly (N+1)(N+2) 2 elements, and more precisely this set is a disjoint union of a set of (N2)(N1) 2 interior points and a set of 3N boundary points. The number of interior points is order ofO(N 2 ) whereas the number of boundary points is order ofO(N). The points gradually expand asN increases and the state space of an innite population, which we denote by 1 , eventually lls in the whole 2simplex,S 2 (Fig.3.2). Lettingx H andx R the frequencies of healthy cells and resistant cells, a state, (H;S;R) = (x H ; 1x H x R ;x R )2 1 of the deterministic system is assigned to the point inS 2 that is distant away fromSR byx H and fromSH byx R . (a) (b) (c) (d) Figure 3.2: The space set, N , of the Moran process withH,S andR gets denser in the phase space,S 2 , with the increase of the population size,N. For eachN,j N j is order ofO(N 2 ), being equal to (N+1)(N+2) 2 . (a)N = 10,j N j = 66; (b)N = 20,j N j = 231; (c)N = 30,j N j = 496; (d)N = 40,j N j = 861 Since the Moran process is a discrete-time and nite-space Markov chain that each type of subpopu- lations changes its size by at most 1 in one time step, a state, (i;j), only moves to one of its neighboring 6 lattice points in a phase space or stays remained, each with a certain probability. This microscopic proba- bilistic movement denes the Moran process with three strategies to be a random walk in N , not inR 2 . 50 DeneT XY i;j to be the transition probability at (i;j) that the number ofX subpopulations increases by 1 whereas the number ofY subpopulations decreases by 1, whereX;Y 2fH;S;Rg. For example, recalling that a state is determined by the number ofH andR subpopulations,T HS i;j is the transition probability that the state, (i;j), moves to (i + 1;j) in one time step. LetT cons i;j be the probability that the system at the state, (i;j), stays remained in one time step. Then all, not necessarily zero, 7 transition probabilities are given by: T HR i;j = if H i;j N w i;j j N ; T HS i;j = if H i;j N w i;j Nij N ; T RH i;j = jf R i;j N w i;j i N ; T RS i;j = jf R i;j N w i;j (Nij) N ; T SH i;j = (Nij)f S i;j N w i;j i N ; T SR i;j = (Nij)f S i;j N w i;j j N ; T cons i;j = 1 (T HR i;j + +T SR i;j ); (3.2) where the weighted population size, N w i;j , is dened byN w i;j = if H i;j + (Nij)f S i;j +jf R i;j and the expected tness functions, f X i;j , are dened by the equation in (1.27) with the time-dependent selection functions in (1.33). Then the Moran process withH,S andR eventually evolves according to an initial state and these transition probabilities in (3.2). Together with the stochastic Moran process, we also consider the deterministic adjusted replicator equations for an innite population playing the same game in (3.1) since this is the asymptotic system of the Moran process in the large population size limit and hence we may have an intuition for the dynamics 51 driven by the Moran process for a nite populations with a large population size. The adjusted replicator dynamic is governed by the following equations: _ x H = f H <f > <f > x H ; _ x S = f S <f > <f > x S ; _ x R = f R <f > <f > x R ; (3.3) where the expected tness functions at~ x := (x H ;x S ;x R ) | are given, with the time-dependent selections functions in (1.33), by: f H = 1w H +w H (A~ x) 1 ; f S = 1w S +w S (A~ x) 2 ; f R = 1w R +w R (A~ x) 3 : (3.4) It is analytically proven that when a constant chemo dose is administered over time, the system has asymptotically stable states: (x H ;x S ;x R ) = (0; 1; 0) wheneverC(t) < 1 3 and (x H ;x S ;x R ) = (0; 0; 1) whenever C(t) > 1 2 . The sensitive populations dominate the total population for constant low dose chemotherapy whereas the system is led to the all resistant populations for constant high dose chemother- apy regardless of an initial distribution. A low chemo dose acts as a main tool to remove the competing agents, resistant cells, helping sensitive cells proliferate in a population by making the environment ben- ecial for them. Similarly, a high chemo dose allows resistant cells to win in a competition with sensitive cells and eventually take over the whole population. That is, the chemo dose plays the role of controller of the deterministic system, causing the competitive release between sensitive cells and resistant cells. For an intermediate value ofC(t), the system has bistability to cancerous states, and the initial state determines its fate to either S or R. In other words, the proportion of area of states in a phase space 52 eventually converging to the allS state is equal to 1 forC < 1 3 , and it starts gradually decreasing in the increase ofC(t) and nally reaches 0 forC > 1 2 . The proportion of area of states in a phase space being attracted to the allR state is exactly the reverse increasing on [0; 1]. Moreover, the adjusted replicator system has an internal xed point 3(48C + 47) 2(32C + 13) 7 2 ; 48C + 47 4(32 + 13) ; 9 2 7(48C + 47) 4(32C + 13) (3.5) for 19 48 <C < 0:625, and the system spirals out from this internal xed point. 3.1 FixationprobabilityfortheMoranprocesswiththreestrategies In this section, we want to examine how the stochasticity driven by the niteness of population size distorts the adjusted replicator dynamics in terms of xation probability. In other words, the question is to study at what probability the Moran process with three strategies leads the nite population to be xated at the all sensitive populations when it stats at (i;j), and similarly to the all resistant populations. To this end, we denote by X i;j the xation probability toX 2fH;S;Rg at the state (i;j). We will dropX in X i;j when the specication of type is not needed. At a boundary state, a population has at most two dierent types of cells and the revival of the eliminated type never happens since no mutation is assumed in our model. Thus, the further evolution since the arrival at a boundary state is then described by the Moran process with two strategies, for which the xation probability is ready to be used in (2.6). Thus the question reduces to nding the xation probability only at interior states, for instance, at all colored dots in Figure 3.2. 53 At an interior point, (i;j)2int( N ) :=f(i;j)j1iN 2; 1jN 2;i +jN 1g, we have a following recursion in terms of the xation probability, i;j , and transition probabilities in (3.2): (1T cons i;j ) i;j =T HR i;j i+1;j1 +T RH i;j i1;j1 +T HS i;j i+1;j +T SH i;j i1;j +T SR i;j i;j1 +T RS i;j i;j+1 : (3.6) Using the lexicographic order,r(i;j), (1; 1); (1; 2); ; (1;N 2); (2; 1); (2; 2); ; (2;N 3); ; (N 3;1); (N 3; 2); (N 2; 1); (3.7) we assign to each state, (i;j)2int( N ), the order: r(i;j) = (i 1)(2N 2i) 2 +j: (3.8) Equating i;j to(r(i;j)), the recursion in (3.6) is re-written in a matrix form: EX =F; (3.9) whereE is a (N1)(N2) 2 (N1)(N2) 2 square matrix and X = [(r(1; 1));(r(1; 2)); ;(r(N 2; 1)] T (3.10) 54 is a (N1)(N2) 2 1 column vector. Note that all elements ofX are the xation probabilities, i;j ’s, at (i;j)2 int( N ), listed using the order r(i;j) in (3.8), and the elements in the r(i;j) th row of E are precisely determined by the transition probabilities to either its surrounding six states such as (i 1;j), (i;j1), (i1;j1) or the same state, (i;j). The entries ofF are completely determined by the xation probabilities at the boundary states. By solving the linear systemEX =F in (3.9), we may compute the xation probability onint( N ). The proof and the detailed explanation for the idea of getting the xation probability for the Moran process with three strategies can be found in [15]. We do not go over it but we will numerically study the stochastic eect of a nite population on asymptotic dynamics using this ndings. In order to see a local eect, we x a state with the proportion, (H;S;R) = (0:8; 0:1; 0:1), where healthy cells are abundant in a total population. At this state, we numerically compute the xation prob- abilities to H, S or R when chemotherapy is administered at a constant dose in [0; 1] (Fig.3.3). Unlike what is analytically shown that the adjusted replicator dynamics has two attracting state,S orR, depend- ing on the constant drug concentration, the Moran process evolves and may eventually reach any one of homogeneous states,H,S,R, with a nonzero probability for any value ofC(t), particularly for a smallN. (a) (b) (c) Figure 3.3: Local xation probability for the Moran process at the state, (H;S;R) = (0:8; 0:1; 0:1), ap- proaches the proportion of area for the basin of attraction by the adjusted replicator dynamics as the population size,N, increases although unexpected features such as xation toH are observed for a small N. (a)N = 100; (b)N = 500; (c)N = 1000 55 The stochasticity introduced to a nite population by its size allows the stochastic system to have nonzero xation probability to H as well as to eitherS or R. For example, in Figure 3.3a, the xation probability toH (green line) at the state, (0:8; 0:1; 0:1), forN = 100 is overall 10% beating the xation probability to S (blue line) for some 1=3 < C < 1=2, with which the adjusted replicator dynamic is asymptotically driven to either all sensitive state or all resistant state depending on the drug intensity. However, for a xed suchC, this xation probability toH at this state tends to decrease asN increase and the xation toH seems almost impossible whenN reaches 500. That is, this new scenario, that is not expected by the adjusted replicator dynamic, such as the absorption of the system toH is captured with a smallN, but it disappears shortly with the increase ofN. Also, on [1; 1=3], the xation probability toR (red line) is nonzero and increases, eventually reaching nearly 0:5 aroundC = 1 3 forN = 100. However, the adjusted replicator system converges to the all sensitive states for C < 1 3 regardless of the initial distribution, and the area of the basin of attraction toR is in fact equal to zero (orange). Similarly, with the increase ofN, the xation probability toR decreases to zero for eachC on [0; 1 3 ], and it numerically supports the convergence of the Moran process to the adjusted replicator system again. For comparison, the proportion of areas on which attraction to eitherS orR is reached by the adjusted replicator system is presented along with the xation probabilities in each panel of Figure 3.3. All of these is literally the universal property throughout all states in a population though the xation probability was examined at one state with the proportion, (0:8; 0:1; 0:1). We will now focus on the global xation probability by xing some values ofC’s instead: C(t) 0:29; 0:37; 0:45; 0:53. For a while, we point out numerical constraint on the size of a population for the computation of global xation proba- bilities. In order to numerically compute the xation probability at all states requires to solve the linear system in (3.9), where the coecient matrix,E, has the size, (N1)(N2) 2 (N1)(N2) 2 . However, the matrix,E, is sparse of which each row has at most 7 nonzero elements since the Moran process with three strategies is able to change its current state, (i;j), to one of 6 nearby states, (i;j), 56 (i;j), (i;j), or to stay unchanged at (i;j) in one time step. In fact, ther(i;j)th row ofE, where (i;j) corresponds to a red dot in Figure 3.2, has exactly 7 nonzero elements, the row corresponding to a black dot has 5 nonzero elements, and the row corresponding to a green dot has 3 nonzero elements. This sparsity helps to reduce the question to solving the linear equation in terms of the sparse matrix ofE, where the number of elements of the sparse matrix,E, is order ofO(N 2 ) while it is order ofO(N 4 ) for the full matrix,E. For a numerical solution, we use an iterative method, the least squared method (lsqr), with the error tolerance, 10 6 . This setting allows us to obtain the global xation probability forN up to 2; 000 in a reasonable computation time. The computation time dramatically increases inN and it successfully computed the numerical global xation probability forN = 2; 000 within a day, but it fails during a week forN = 2; 500. (a) (b) (c) (d) Figure 3.4: The basins of attraction by the adjusted replicator dynamics with constant chemo concentration, C :=C(t), shows the smooth transition from the global attraction byS to the global absorption toR with the increase ofC, describing the competitive release of resistant subpopulations with high chemotherapy concentration. The system experiences bistability toS orR for an intermediate value ofC. (a)C 0:29; (b)C 0:37; (c)C 0:45; (d)C 0:53 For the purpose of comparison with the Moran process, we investigate the basin of attraction by the adjusted replicator dynamics withC = 0:29; 0:37; 0:45; 0:53 with which the dierent key dynamics are shown (Fig.3.4). The basin of attraction toS andR are respectively depicted in blue and red regions. The asymptotically dominant subpopulation of typeS withC 0:29 starts reducing its attraction region with the increase ofC. The basin of attraction toR withC 0:37 gradually expands its region from states with fewer sensitive cells in a population towards the all sensitive state withC 0:53. We indicate the separatrix in a black solid line for eachC. Meanwhile, the deterministic system has an internal xed point 57 with the proportion, (0:2555; 0:6259; 0:1186), withC 0:45 from which the system spirals out. In fact, this internal xed points exists for allC such that 19 28 <C < 0:625. (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) Figure 3.5: Fixation probability toR for the Moran process with constant chemotherapy concentrations, C := C(t), shows the stochastic version of the competitive release of the resistant cells with high con- centration for each population size,N. For a xedC, each region in a phase space with either a high or a low probability respectively approaches the basin of attraction to eitherR orS by the adjusted replicator system asN increases. (a)C 0:29,N = 100; (b)C 0:37,N = 100; (c)C 0:45,N = 100; (d) C 0:53, N = 100; (e)C 0:29, N = 500; (f)C 0:37, N = 500; (g)C 0:45, N = 500; (h) C 0:53,N = 500; (i)C 0:29,N = 1000; (j)C 0:37,N = 1000; (k)C 0:45,N = 1000; (l) C 0:53,N = 1000; (m)C 0:29,N = 2000; (n)C 0:37,N = 2000; (o)C 0:45,N = 2000; (p) C 0:53,N = 2000 In Figure 3.5, the xation probability toR at every state is described in a way that blue indicates the xation toR never happens whereas red indicates that the xation surely occurs. WithC 0:29, the 58 Moran process has no chance to reach the xation toR at the majority of states, but the system at states where sensitive cells are rare in a population evolves and the resistant cells nally dominate the entire population with a high probability for a small population sizeN = 100 (Fig.3.5a). As expected, the region of positive xation probability toR gradually shrinks asN increases toN = 2; 000 (Fig.3.5a,3.5e, 3.5i, 3.5m). WhenC 0:37, the area of states in which the xation probability toR is reached with a high probability is near the side,HR, where there are relatively a small number of sensitive cells in a population. This region is included in the basin of attraction toR of the adjusted replicator dynamics, colored in red in Figure 3.4b, for eachN. It expands towards the separatrix, and the area with non zero xation probability toR gets sharper and sharper around the separatrix in the increase ofN (Fig.3.5b, 3.5f, 3.5j, 3.5n). When C 0:45, as the adjusted replicator system has an internal xed point from which the system spirals out, the Moran process also shows that the high and low probability of the xation toR is obtained in each side of the basin boundary, and the region of the neutral xation probability gets narrower around the basin boundary asN increase (Fig.3.5c, 3.5g, 3.5k, 3.5o). WithC 0:53, the almost sure xation toR in the total phase space is hardly met even when a population size is as big asN = 2000, however, the region where the xation is certainly attained gets enlarged towards the state,S, in the increase ofN (Fig.3.5d, 3.5h, 3.5l, 3.5p). The xation probability to R for the Moran process overall shows the asymptotically similar pat- terns to the basins of attraction toR by the adjusted replicator dynamics for a largeN for eachC(t) 0:29; 0:37; 0:45; 0:53. Thus it is rarely reached at a state outside the basin boundary though it is still pos- sible at a low probability depending on the mixture of the subpopulations. However, when the population size is small so that the stochasticity is strong, somewhat dierent dynamic from the deterministic system is seen mainly around the basin boundary withC 0:37; 0:45. 59 3.2 Singledrugadaptivecontrol Setting the chemotherapy concentration, C := C(t), constant in time, it was found that the adjusted replicator system has two attractor states depending on the constant value of the chemo concentrationC. With the choice ofC(t) 0, the trajectories of the adjusted replicator dynamics starting at states near the all healthy state show their eventual convergence to the all sensitive state in a phase space. In fact, the basin of attraction toS is the whole region of a phase space except for two states: the all healthy state and the all resistant state. WhenC(t) 0:7, the trajectories of the deterministic system starting at the states near the equilibrium, 0; 9C3 4C+2 ; 55C 4C+2 ) = (0; 0:6875; 0:3125 , show that the system is driven to the all resistant state, and this absorption is made at all state but bothS andR (Fig.3.6). More interestingly, overlapping these two families of dynamics with two dierent constant chemo concentrations produce several closed loops (Fig.3.6c). By scheduling a chemotherapy as a step function betweenC(t) 0 and C(t) 0:7, it is possible to permanently trap the deterministic system in a closed loop and prevent it from the absorption to a cancerous state as under a constant chemo dose. We pick, for numerical experiment, one closed loop and denote the intersection near theS corner byP and the intersection near the line,HR, byQ, whereP = (0:04; 0:9; 0:06) andQ = (0:53; 0:11; 0:36) (Fig.3.7a). (a) (b) (c) Figure 3.6: Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted replica- tor system for dierent constant chemotherapy values. (a) UnderC(t) 0, the competitive release of the sensitive subpopulations,S, to drug drives all trajectories to theS corner. (b) UnderC 0:7, the competi- tive release of the resistant subpopulationsR to drug drives all trajectories to theR corner. (c) Trajectories with two dierent constant chemotherapy combinations overlap at dierent times and generate a closed loop. 60 In order to design a chemotherapy schedule precisely, letT PQ be the time it takes that the point,P , evolves by the adjusted replicator equations and reaches the point,Q, and letT QP be the time that the system takes to get from Q to P . Administering a high chemo C(t) 0:7 during T PQ unit time, the adjusted replicator system starting at pointP moves along the red line arriving atQ. Then, turning the chemo o for T QP unit time, the point Q is driven by the system back to the point P along the blue line (Fig. 3.7a). We call the time,T PQ +T QP , that it takes to return to the initial state one evolutionary cycle. The evolutionary cycle varies with the choice of two states and thus the chemo schedule needs to be carefully designed with the deep investigation of the patient’s medical status in the clinical purpose. Also for the later use, we call the time,T PQ orT QP , a half evolutionary cycle. In our deterministic model, each half cycle isT PQ = 19:36 andT QP = 12:56 and hence one evolutionary cycle is equal to 31:92 unit time. (a) (b) Figure 3.7: Switching chemotherapy on and o at adequate times traps a trajectory associated with the adjusted replicator system within a closed loop. (a) The system, that starts atP , moves along a red line and arrives atQ when the high chemo dose,C(t) 0:7, is continuously administered duringT PQ = 19:36 unit time. When chemo is turned o during furtherT QP = 12:56 unit time, it returns to the initial point, P , eventually generating a closed loop, PQP ; (b) The tumor volume is controlled with the use of the adaptive schedule for an evolutionary cycle, experiencing tumor regression and recurrence. Continuing this chemotherapy schedule, what is called an adaptive schedule, for innitely many evo- lutionary cycles, the adjusted replicator system starting atP returns back toP everyT PQ +T QP unit time. In addition, the initially high volume of tumor decreases forT PQ since the high chemo helps to re- move the sensitive cells in a population at the beginning, however, the drug administration during longer than necessary period results in the regrowth of tumor slightly before it reachesT PQ . After turning the 61 chemo o, the tumor volume keeps increasing and nally recovering the initial tumor volume when one evolutionary cycle is over (Fig.3.7b). Using the adaptive schedule for multiple cycles causes the tumor size to oscillate and be controlled by preventing the system from the absorption to eitherS orR corners. (a) (b) (c) (d) (e) (f) (g) (h) Figure 3.8: Realizations of multiple trajectories associated with the Moran process under administration of a constant chemotherapy show the ability of the stochastic system to behave similarly to what the adjusted replicator dynamics drive, getting closer as the population size increases. The Moran process, starting at a state near the cornerH withC(t) 0 (blue wiggled lines), possibly evolves and attains the homogeneous population of allS for eachN, having smoother trajectories asN increases and nally being similar to the deterministic trajectory (light blue line). Similarly, the Moran process starting at a state near the xed point (0; 0:6875; 0:3125) (red wiggled lines) is able to be driven to the stateR withC(t) 0:7 for each N. (a)C(t) 0,N = 1K; (b)C(t) 0,N = 5K; (a)C(t) 0,N = 10K; (a)C(t) 0,N = 50K; (e) C(t) 0:7,N = 1K; (f)C(t) 0:7,N = 5K; (g)C(t) 0:7,N = 10K; (h)C(t) 0:7,N = 50K We will apply this adaptive schedule to the Moran process withH,S,R and see if how badly our model behaves for one evolutionary cycle. To this end, we rst investigate the stochastic trajectory against the deterministic trajectory in a phase space with the use of a constant chemo schedule forN = 1K, 5K, 10K, 50K, where the stochastic movement is determined by the transition probabilities in (3.2). We are able to nd one realization of the Moran process starting at a state near theH corner and achieving the xation toS with no chemo administered for eachN. One realization of the Moran process starting at a state near the xed point, (0; 0:6875; 0:3125), and eventually converging to R with the constant use of a drugC(t) 0:7 is also obtained. For both constant drug schedules, the stochastic trajectories get 62 wiggled around the deterministic trajectories and become smoother and smoother in the increase ofN nally being similar to the deterministic ones for a largeN (Fig.3.8). However, it is just one realization of the Moran process as an example and a stochastic system can only be understood in terms of its distribution in general. Recalling the relation in (1.17) between the evolution step,, for the Moran process and the evolution time,t, for the deterministic adjusted replicator system through which the convergence of the Moran process to the deterministic system is driven, we simulate 1; 000 individual Moran processes all starting atP during the half evolution steps,T PQ N, with C(t) 0:7 and indicate the terminal points of all realizations as red dots for eachN = 1K, 5K, 10K and 50K (Fig.3.9). (a) (b) (c) (d) Figure 3.9: The spread of the distribution of points aroundQ (orP ) for 1; 000 realizations of the stochastic Moran process gets denser and demonstrates the shrunken randomness as the population size increases when each realization is under the administration of a constant chemo scheduleC(t) 0:7 (orC(t) 0) during half an evolutionary cycle,T PQ (orT QP ), since its exact start atP (orQ). (a)N = 1K; (b)N = 5K; (c)N = 10K; (d)N = 50K It is obvious that they form a distribution near the point,Q, at which the deterministic ends its half cycle starting fromP , and the deviation becomes smaller in the increase ofN getting denser aroundQ. When they are plotted in theSR coordinate system, the mean frequency, S , of the sensitive subpopulations at the terminal points around Q for 1; 000 realizations approaches the proportion of S at Q, which is precisely equal to 0:11 (Fig.3.10a-3.10d). It means that asN increase, the mean frequency S , which is the sample mean, becomes a less biased estimator for the proportion ofS atQ. The mean frequency, R , of the resistant subpopulations aroundQ also becomes a more accurate estimator for the frequency ofR at the state,Q, being equal to 0:36. One thing notable from the distribution inSR coordinate system is that 63 the deviation in the direction ofR is in general smaller than in the direction ofS for eachN though both are getting smaller in the increase ofN, and it is because the expected deterministic system diverges near P moving towardsH and then eventually towardsR withC 0:7 while it converges toS withC 0 as shown in Figure 3.6. When the terminal points aroundQ are plotted in the principal axis coordinate system, we see that not only both the semi-major and the semi-minor axis decrease inN but also it forms a multi Gaussian distribution for a largeN (Fig.3.10e-3.10h). We separately simulate 1; 000 individual Moran process all starting atQ during the evolution steps, T QP N, withC(t) 0 and indicate the terminal points of each realization in blue dots. It shows the similar pattern that the rst half cycle withC(t) 0:7 results in terms of the shrinking distribution. The mean frequency, S , of the sensitive subpopulations at the terminal points of the1; 000 realizations of the Moran process aroundP converges to the frequency ofS atP being equal to 0:9 while the mean frequency, R , of the resistant subpopulations aroundP approaches the proportion ofR atP being equal to 0:06 asN increases. However, unlike that the terminal points from the rst half evolutionary cycle are more widely distributed along theS axis in theSR coordinate system, the end pints from the second half cycle are equally likely distributed along both axes, and it is because the expected deterministic system converges toS underC(t) 0 (Fig.3.10i-3.10p). In addition, when they are plotted in the principal axis coordinate system, we see that the multi Gaussian-like distribution is achieved with the smallerN compared to the rst half evolution cycle, and both the semi-major and the semi-minor axis are apparently smaller for each N. Now, we extend the evolution steps from a half evolutionary cycle up to the full one cycle, (T PQ + T QP )N. It means that we apply the adaptive chemotherapy schedule associated with the deterministic system to the Moran process in the following way: each realization of the Moran process starts atP under a high constant chemo dose,C(t) 0:7, duringT PQ N evolution steps, and once it reaches a neighborhood ofQ in the last step, we turn the high chemo o and the Moran process evolves from that neighborhood 64 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) Figure 3.10: The spread of the distribution of points in the principal axis coordinate system for 1; 000 realizations of the Moran process is, for a large population size, characterized as a multi Gaussian distri- bution aroundQ (orP ) when each realization is under the administration of a constant chemo schedule C(t) 0:7 (orC(t) 0) during half an evolutionary cycle,T PQ (orT QP ), since its exact start atP (or Q). The mean frequency, S (or R ), of the sensitive (or resistant) subpopulations around the point,P (or Q), converges to the proportion ofS (orR) asN increases, with the decreasing semi-major axis, 1 , and semi-minor axis, 2 . (a)N = 1K, S = 0:1920, R = 0:2930; (b)N = 5K, S = 0:1223, R = 0:3452; (c)N = 10K, S = 0:1154, R = 0:3479; (d)N = 50K, S = 0:1099, R = 0:3590; (e)N = 1K, 1 = 0:2800, 2 = 0:1356; (f) N = 5K, 1 = 0:1401, 2 = 0:0386; (g) N = 10K, 1 = 0:1028, 2 = 0:0238; (h) N = 50K, 1 = 0:0460, 2 = 0:0093; (i) N = 1K, S = 0:8907, R = 0:0684; (j)N = 5K, S = 0:8987, R = 0:0606; (k)N = 10K, S = 0:8990, R = 0:0601; (l)N = 50K, S = 0:8992, R = 0:05989; (m)N = 1K, 1 = 0:0551, 2 = 0:0119; (n)N = 5K, 1 = 0:0224, 2 = 0:0050; (o)N = 10K, 1 = 0:0150, 2 = 0:0035; (p)N = 50K, 1 = 0:0067, 2 = 0:0016 duringT QP N evolution steps, which expectedly arrive nearP . According to the convergence of the system, it is expected that for a largeN, this return is achievable with a high probability, however, it is 65 hard to occur in a small population due to the strong stochasticity though not impossible. We show one example of the stochastic trajectory with a small population of size,N = 50, along the lattice together with the deterministic trajectory in a phase space when a drug is administered according to the adaptive chemotherapy schedule. Figure 3.11: The stochastic trajectory of one realization of the Moran process under the administration of the adaptive chemotherapy, associated with the adjusted replicator dynamics, develops a random walk along a lattice in a phase space,S 2 . The adaptive schedule is able to prevent the stochastic system from the saturation of cancer cells, even in a small population withN = 50. Starting atP = (0:04; 0:9; 0:06) (black dot), it moves along the lattice (red line) and reaches a neighborhood (green dot) ofQ = (0:53; 0:11; 0:36), duringT PQ evolution time. Turning o the chemo at that green dot, the stochastic system evolves (blue line) and eventually reaches a neighborhood (yellow dot) of the initial point,P , duringT QP . The strong stochasticity allows the trajectory of the Moran process starting atP underC(t) 0:6 to jump back and forth around the deterministic trajectory along the lattice in a phase space and terminates its movement at a point far away fromQ. However, it is also the strong stochasticity that leads the terminal point withC(t) 0:7 to a point nearP afterT QP N evolution steps underC(t) 0 preventing the system from the absorption to a cancerous state. Again, it is one realization and we simulate 1; 000 realizations of the Moran process under the adaptive chemo schedule to understand it in distribution withN = 10K and 50K, and see how largeN enables the stochastic system to resume its initial state after one evolutionary cycle with a high probability. The distribution of the points nearQ shows the similar pattern that 1; 000 realization of the Moran process during the rst half cycle results since it is just one another sample of the same size. However, the distribution of the points nearP after the full one cycle is more broad compared to the case during the 66 (a) (b) (c) (d) (e) (f) Figure 3.12: The averaged trajectory of 1; 000 realizations of the Moran process under the adaptive sched- ule, associated with the adjusted replicator system, during one evolutionary cycle (T PQ +T QP = 31:92 unit time) ts the corresponding deterministic trajectory for a large population size. The Moran process is likely to return nearly to the initial state with a high probability for a largeN even though the spread of the distribution of the points nearQ is still wide. (a) the distribution of the points associated with the adaptive chemo schedule forN = 10K; (b) the trajectory of one single realization of the Moran process withN = 10K; (c) the averaged trajectory of 1; 000 realizations of the Moran process withN = 10K; (d) the distribution of the points forN = 50K; (e) the trajectory of one single realization withN = 50K; (f) the averaged trajectory withN = 50K second half evolutionary cycle shown in Figure 3.9 for eachN (Fig. 3.12a, 3.12d). It is obviously because each realization starts at one of red dots nearQ for the second half cycle under the adaptive chemo schedule unlike it starts at the exact point,Q, for Figure 3.9, and this drives the system to settle down at a point further away fromP resulting in a broad distribution around. Clearly, the distribution around the initial stateP is of the greater interest due to the fact that the randomness works in a cooperative way so that the system fallen far away fromQ after the rst half evolution cycle might be able to end its evolution near the stateP in the end of the full one round. Thus, we will investigate the distribution only around the initial point when it comes to the full evolu- tion cycle or longer. The distribution of points aroundP in theSR coordinate system particularly when 67 N = 50K shows the deviation in either direction is equally likely forming a multivariate Gaussian distri- bution with the mean frequency, S = 0:8973, of the sensitive subpopulations and the mean frequency, R = 0:06125, of the resistant subpopulations (Fig.3.13). These play the role of the estimators for the proportion, 0:9, ofS and the proportion, 0:06, ofR at the state,P , being less accurate compared to the case with the half evolutionary cycle discussed in Figure 3.9. (a) (b) (c) (d) Figure 3.13: The spread of the distribution of points for 1; 000 realizations of the stochastic Moran process withN = 50K is characterized as a multivariate Gaussian distribution, centered nearly at the initial point, P , when each realization is under the administration of the adaptive schedule during one evolutionary cycle (T PQ +T QP = 31:92 unit time). As shown in theSR and the principal axis coordinate system, the deviation is equally likely to each other in either directions. (a) the distribution of the points around P in theSR coordinate system; (b) the kernel density distribution in theSR coordinate system; (c) the distribution of the points in the principal axis coordinate system; (d) the kernel density distribution in the principal axis coordinate system We also show the whole trajectory of the Moran process under the adaptive schedule for one evolution- ary cycle in Figure 3.12. One realization of the process moves closely to the closed loop around the points, P andQ, for eachN while the trajectory passes the point,Q, during the rst half cycle withN = 10K but it does not reach close enoughQ withM = 50K (Fig.3.12b, 3.12e). However, the averaged trajectory of the Moran process withN = 50K denitely shows the closer t to the deterministic trajectory (Fig.3.12c, 3.12f). This implies that the adaptive chemo schedule helps the stochastic system recover the initial state after one evolutionary cycle at a high probability for a largeN. 68 3.3 Singledrugadaptivecontrolformultipleevolutionarycycles Now, we want to continue applying this adaptive chemotherapy schedule to the Moran process but during multiple evolutionary cycles, and see how early the recovery of the approximate initial state is broken for a large population. To be precise, the Moran process initially starting atP evolves under the adaptive chemo schedule for one evolutionary cycle and arrives at a point in the neighborhood ofP . Then at that terminal state of the rst cycle, we turn the chemo on again to start the second round of the adaptive chemo schedule for the next evolutionary cycle, and continue doing these steps for a few more rounds. That is, the process in the current round starts at the terminal state in the previous round for one evolutionary cycle. This is a natural extension of the application of the adaptive schedule to the Moran process, but it is not straightforward or easy to answer the question due to the facts obtained from its expected system. The adjusted replicator system is sensitive in the initial state as well as in the chemotherapy concen- tration: initially close two points can have completely dierent fate at the end under the same constant chemo schedule, particularly either if they are located across the basin boundary for a xedC such that 1 3 < C < 1 2 or if they are in the neighborhood of an unstable xed point of the system. Also the basin boundary of the deterministic system experiences the smooth transition in a phase space as the constant drug concentration,C, which acts as a parameter, increases from 1=3 to 1=2. Hence the system starting at a xed initial state can be lead to a dierent homogeneous state with a slight change inC if the change causes the initial point to be located on the other side of the basin boundary. For these reasons and with the randomness of the stochastic system added, the fate of the Moran process in the following round is indeterminate even though the Moran process under the adaptive schedule returns close enough to the initial state in the current round. Here, we consider the adaptive schedule during 4 cycles for the numerical experiment withN = 10K and 50K. The stochastic trajectory during the rst, second, third and fourth round is plotted in a red, green, blue and pink line, respectively while the deterministic trajectory is shown in light blue in Figure 69 (a) (b) (c) (d) Figure 3.14: For each N, the averaged trajectory of 1; 000 realizations of the stochastic Moran process shows that the saturation of cancer cells can be delayed until the 4th evolutionary cycle when each real- ization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start atP , although the tightness of the the averaged stochastic trajectory to the deterministic one lasts shorter with a smaller population size. (a) the trajectory of one single realization for the Moran process withN = 10K; (b) the averaged trajectory of 1; 000 real- izations for the Moran process withN = 10K; (c) the trajectory of one single realization withN = 50K; (d) the averaged trajectory withN = 50K 3.14. Also, we mark the terminal point in each round in the starred shape, colored accordingly. One realization of the Moran process starting atP under the adaptive schedule withN = 10K shows that it moves along the red line in the rst round as expected by the deterministic system (light blue line) and arrives at the red point nearP before it terminates its life at the all sensitive state (Fig.3.14a). The second round starts at that red point nearP , and it moves along the green line being o the expected trajectory and returns to the green point nearP in the end of the second round. The third round starts at the red point nearP , and it arrives at the blue point nearP after moving along the blue line. Unfortunately, the blue point at which the third round ends represents the state where one subpopulation is rare in the total population, and it possesses a high potential that the system is driven to theSH line where the population is composed of only two subpopulations of the type,S andH. Once it enters a state of two subpopulations, the continuation of the adaptive chemo schedule eventually causes the system to converge to a cancerous state. A similar scenario is shown in the single realization of the Moran process withN = 50K (Fig.3.14c). It moves closely around the closed loop associated with the deterministic system up to the second round, however, it starts being far o the expected orbit and nally ends at the pink point nearP at the end of the 70 fourth round. Then it is from now on laid in the similar situation to the end of the third round in Figure 3.14a withN = 10K. Nevertheless, it is not that disappointing to apply the adaptive chemo schedule when the averaged trajectory for 1; 000 realizations of the Moran process is taken into account. On average, the adaptive schedule keeps cancerous cells from xating the total population for 4 evolutionary cycles. The Moran process starting at P with N = 10K under the adaptive schedule approximately re-attains the initial stateP in the end of each round though it becomes further and further fromP as the number of rounds increases and it does not behave similarly to the deterministic system from the third round (Fig.3.14b). However, the distance fromP to the averaged point in the end of each round gets smaller asN increases, and it results that the Moran process behaves similarly to the deterministic system in its trajectory much longer (Fig.3.14d). (a) (b) (c) (d) (e) (f) (g) (h) Figure 3.15: The spread of the distribution of terminal points (blue dots) around P in each round for 1; 000 realizations of the stochastic Moran process becomes wider as the number of rounds increases when each realization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start at P . (a) N = 10K, round 1; (b) N = 10K, round 2; (c)N = 10K, round 1; (d)N = 10K, round 4; (e)N = 50K, round 1; (f)N = 50K, round 2; (g)N = 50K, round 1; (h)N = 50K, round 4 In fact, the points nearP for 1; 000 realizations of the Moran process in the end of each round are distributed wider and wider as the number of rounds increases for eachN = 10K and 50K (Fig.3.15). 71 Obviously, the relatively thick cluster aroundP is maintained in each round withN = 50K compared to N = 10K. Finally, the cluster around P starts blurry from the third round with N = 10K and in the fourth round withN = 50K (Fig.3.15). As they lose clustering aroundP , many of them end their 4 rounds of evolution near the all sensitive state meaning that the xation ofS eventually reached with a high probability. This becomes more apparent to see when the points are plotted in the SR coordinate system. For N = 50K, the range of the distribution of S and R are equally likely up to the last round (Fig.3.16). However, they get more clustered towards the point (1; 0) in theSR coordinate system as the number of rounds increases, where (1; 0) represents the all sensitive state. Finally they are no longer clustered around P in the last round and it means that the sensitive mutants nally xates the total population with a high probability after 4 evolutionary cycle. Also, the multivariate Gaussian distribution is already broken in the third round, and both the semi-major and the semi-minor axis keeps increasing in the number of rounds. (a) (b) (c) (d) (e) (f) (g) (h) Figure 3.16: As the number of rounds increases, the spread of the distribution of points aroundP for 1; 000 realizations of the stochastic Moran process forN = 50K becomes centered towards the homogeneous sensitive state and immediately looses the multivariate Gaussian-like distribution, where each realization evolves under the administration of the adaptive schedule, associated with the adjusted replicator system, during 4 evolutionary cycles since its exact start atP . (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 1; (f) round 2; (g) round 3; (h) round 4 72 Extending the number of the evolutionary cycles up to 8 rounds, we clearly see the overall increase of these values in the number of rounds for eachN = 10K and 50K while they are overall smaller with N = 50K due to the thicker cluster aroundP (Fig.3.17). Fitting the discrete data in the log-log plot shows the power-law dependence of the semi-major (minor) axis on the number of rounds. Precisely, the tted curve for the semi-major (minor) axis of the distribution of the points aroundP for 1; 000 Moran processes starting atP as a function of the number of rounds in the log-log scale is following: forN = 10K 1 0:0505n 1:3269 ; 2 0:0071n 1:8753 ; (3.11) and forN = 50K 1 0:0149n 1:8519 ; 2 0:0025n 2:0424 ; (3.12) where 1 is the semi-major axis, 2 is the semi-minor axis, andn is the number of rounds. (a) (b) (c) (d) Figure 3.17: The semi-major (and -minor) axis of the distribution of the points around P at the end of each evolutionary cycle for 1; 000 realizations of the Moran process overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the adaptive schedule associated with the adjusted replicator system during 8 evolutionary cycles since its exact start atP . (a)N = 10K; (b)N = 10K, the log-log t; (c)N = 50K; (d)N = 50K, the log-log t 73 3.4 Comparison of adaptive single drug chemotherapy schedule with standardclinicalapproaches Among several clinical chemotherapy schedules, we consider two frequently used standard schedules: the maximum tolerated dose schedule (MTD) and the low-dose metronomic schedule (LDM). The former is the highest dose of a drug by which any signicantly unacceptable side eects are not produced while the later is the low dose of a drug that is continuously given for a nite time. We compare the eectiveness of the adaptive chemotherapy schedule to those standard clinical chemo schedules in terms of either the prevention ability of the system from converging to a cancerous state or the control ability of its induced tumor size. Figure 3.18: Two standard clinical approaches, the maximum tolerated schedule (MTD) and the low-dose metronomic schedule (LDM), are designed in order to have the same total dose as the adaptive chemother- apy schedule associated with the deterministic replicator system has during 4 rounds. In each cycle of the length, 31:92, MTD delivers a drug at the highest concentration during the rst 13:552 unit time, followed by no chemo until the end of each cycle. On the other hand, a drug is continuously delivered during the whole cycles at as low concentration as 0:42456 for LDM. For our comparison, each chemo schedule is subject to having the same total dose over a nite time. For a drug concentration function,C(t), dened on [0;T ] forT > 0, the total dose,D(T ), delivered during T unit time is dened by D(T ) := Z T 0 C(t)dt: (3.13) 74 LetC Adaptive (t), C MTD (t) andC LDM (t) be the drug concentration functions of the adaptive, MTD and LDM schedules, respectively, at time,t. We assume that Z T 0 C Adaptive (t)dt = Z T 0 C MTD (t)dt = Z T 0 C LDM (t)dt; (3.14) whereT = T PQ +T QP = 31:92 for the numerical purpose. Then it is straight forward that the highest dose is delivered until 13:552 unit time followed by no drug administration until a cycle ends for MTD, and the constant dose as low asC LDM (t) 0:42456 is continuously administered during 31:92 unit time, the whole cycle, for LDM to have same total dose each other in one evolutionary cycle. These schedules are to be repeated as many times as the number of the evolutionary cycles increases (Fig.3.18). As shown in Figure 3.14, the adaptive schedule works eectively for 4 rounds by turning chemo on and o at adequate times, delaying the competitive release and inhibiting the xation of cancerous cells in a population. For a denser population, the trajectory was better controlled on average, tightly moving around the closed loop associated with the deterministic system during longer evolutionary cycles. Dene a tumor size as the sum of frequencies of cancerous cells, that is, the sum of sensitive cells and resistant cells. Note that a tumor size is a function in time ranging on [0; 1] where 0 indicates healthy state with no tumor cells in a population while 1 indicates that no hope expected since a population is merely composed of sensitive and resistant cells. Then, the tumor size under the adaptive schedule is well controlled on average during 4 cycles never reaching 1, the full tumor state. In our model, since the Moran process starts atP where the dominant cells are of the sensitive type, the tumor size at the state, P , is as big as 0:96. As drug is administered at a high dose forT PQ N evolution steps on the Moran process withN = 50K, starting atP , in the beginning of each cycles, tumor size reduces nearly about 0:5 with the slight increase at the end. This slight increase in tumor size origins from the resistance to drug as a result of a long drug administration 75 (a) (b) (c) (d) (e) (f) (g) (h) Figure 3.19: The adaptive chemotherapy schedule associated with the replicator system is compared to the standard clinical approaches, MTD and LDM, in order to demonstrate its eciency in terms of delaying the time of saturation of tumor cells, which is attained before the rst round ends under either MTD or LDM schedules. The adaptive schedule beats the other two standard clinical chemo schedules since not only it prevents the system from converging to a cancerous state and but also the tumor size is thus controlled for 4 rounds on average. (Each of 1; 000 realizations of the Moran process evolves under the administration of the adaptive schedule, MTD or LDM independently during 4 evolutionary cycles since its exact start atP .) (a)N = 10K, one single realization; (b)N = 10K, the tumor size associated with the single trajectory; (c) N = 10K, the averaged trajectory; (d)N = 10K, the tumor size associated with the averaged trajectory; (e)N = 50K, one single realization; (f)N = 50K, the tumor size associated with the single trajectory; (g) N = 50K, the averaged trajectory; (h)N = 50K, the tumor size associated with the averaged trajectory causing the re-growth of tumor. When chemo is shortly turned o, tumor size keeps increasing and almost recovers the full tumor state inT QP N evolution steps. It is true that the tumor size consequently increases on average compared to when the chemo just starts being turned o. However, the composition of the population totally changes to the mostly sensitive cells so that it is ready to sensitively react to a drug. As the second cycle begins, the tumor size undergoes down and up in the similar pattern on average, and continues until the last round (Fig.3.19h). On the other hand, the two standard clinical approaches show hopeless results for a large N. The maximum tolerated dose schedule shortly drives the system to states where the resistant cells are abundant by immediately causing the resistance to a drug and bringing the re-growth of tumor. After the chemo is turned o, the tumor size nally reaches its full volume and the system permanently enters the state 76 where a healthy cell no longer exists even in the end of the rst evolutionary cycle. Later, the control of the chemo dose just changes its state between the more sensitive and the more resistant, but the tumor size never reduces eternally on average. The low-dose metronomic schedule does not seem to be a promising method to attempt to cure a caner once a patient’s state is already full of cancer cells, that the initial point P represents in our model. According to the adjusted replicator dynamics, the population starting atP is likely to end up being the all sensitive with the constant concentration,C(t) 0:42456, becauseP is located in the basin of attraction toS. Thus, the Moran process starting atP is more likely to converge toS with a higher probability asN increases. It implies that not only the adaptive schedule but also the clinical approach must be applied after carefully investigating the patient’s state. We summarizing this section by presenting the eectiveness of the adaptive schedule for the Moran process starring atP during 4 evolutionary cycles. In order to quantify the ecacy, we dene that one realization of the Moran process starting atP successfully ends its evolution during 4 evolutionary cycles if it neither gains the full tumor volume on the way nor has over 99% sensitive cells in the end. Figure 3.20: Among 1; 000 realizations of the Moran process, the rate that it is well controlled and returns back nearly to the initial state, neither gaining the full tumor volume nor having over 99% sensitive cells after 4 evolutionary cycles, is computed under the adaptive, MTD or LDM chemo schedules independently. The adaptive chemo schedule obviously is superior than other two schedules, having the nearly full rate whenN = 50K, while MTD has an extraordinarily small but positive success rate up to the 3 rounds. 77 Then we see that the adaptive schedule works better for eachN since it has the highest rate compared to the clinical standard schedules (Fig.3.20). Also, asN increases, its success approximately increases in the form of an exponential function inN: 0:9823e 0:0001N + 1; (3.15) whereN is the population size. However, LDM unfortunately brings no positive fate to the population even when the population is small, suggesting the dierent approach if the patient is at the state,P . The maximum tolerated dose is also not promising at this state with nearly zero chance. When the population size is bigger than 20K, it shows no chance with 1; 000 realizations, however, that results in a small being smaller than 10% but positive probability that the system is successfully not driven to a cancerous state for a smallN. 78 Chapter4 Stochastictwodrugchemotherapymodel In the previous chapter, we studied the single drug chemotherapeutic response to a nite population of competitively interacting healthy, sensitive and resistant cell when modeled by the Moran process. In gen- eral, many clinical approaches use more than one toxin simultaneously. Drugs may interact in a way that they result in a greater cumulative eect than the combined individual eects from the multiple indepen- dent drug use. On the other hand, drug mixture may cause a consequently smaller eect than the summed individual eects. The former is termed as synergism whereas the later as antagonism, where additivity refers to the individual summation of the eects from the independent drug uses. Both synergism and antagonism have been studied and quantied for about a century since the ignition by C. I. Bliss in [5] in 1939, and more recently in [9]. In this thesis, we adopt the drug interaction parameter, e2 [1; 1], quantied by Y. Ma and P. K. Newton in [30] in a way thate < 0 indicates antagonistic drug interactions,e = 0 refers to the additive interactions, ande> 0 stands for synergistic interactions. We model two drug chemotherapeutic response by the Moran process in a similar spirit to the single drug model, however, what makes these two models distinguishable is the drug interaction parameter,e. In this framework, synergism is considered the most promising drug interactions to cure cancer in the sense that it kills as many cancer cells using a smaller amount of toxins as drugs in isolation would do as a sum with a higher amount of dose. However, drug 79 administration continuously changes tness of cells in all type in a coevolving population, and scheduling chemotherapy based on a track of distribution of all cell types seems to be more realistic in controlling tumor volume. 4.1 S;R 1 ;R 2 multidrugmodel The model that we employ is the Moran process with all cancerous cells of three dierent types: (i) sensitive cancer cells (S) that are sensitive to both drug 1 and drug 2, (ii) resistant cancer cells (R 1 ) that are sensitive to drug 1 but resistant to drug 2, and (iii) resistant cancer cells (R 2 ) that are sensitive to drug 2 but resistant to drug 1. We letC i (t) be the concentration of drugi,i = 1; 2, in time,t, ranging from 0 and 1 with an additional condition,C 1 (t) +C 2 (t) 1, for allt. These concentration functions together with the drug interaction parameter,e, enters though selection functions,w X ;X2fS;R 1 ;R 2 g, and distorts the tness functions as follows: w S (t) =w 0 (1C 1 (t)C 2 (t)eC 1 (t)C 2 (t)); w R 1 (t) =w 0 (1C 1 (t)); w R 2 (t) =w 0 (1C 2 (t)); (4.1) wherew 0 is constant that scales out and is set to be equal to 1 for the computational purpose. Similarly to the equation in (1.33), the selection function of sensitive cells is dened to be inversely proportional to the concentration,C i (t), of whatever drugs to which they act sensitively, in the additive drug interaction case. The value ofe in (4.1) accounts for a greater or smaller cumulative eects for synergism and antagonism, compared to the sum of single eects from independent drug use. Thus, for a xede, the concentration functions,C i (t), act as controllers of the stochastic system and allow us to schedule chemotherapy on the purpose of controlling tumor growth with the allowance of acceptable tumor increase. 80 Then these selection functions shapes the tness landscape, dened similarly as in (1.27) but withX2 fS;R 1 ;R 2 g, with the payo matrix wherei counts the number ofS subpopulations whereasj counts the number ofR 2 subpopulations. We assume a pairwise Prisoner’s dilemma game for three subpopulations with a 3 3 payo matrix,A: A = 0 B B B B B B @ S R 1 R 2 S a 11 a 12 a 13 R 1 a 21 a 22 a 23 R 2 a 31 a 32 a 33 1 C C C C C C A : (4.2) Considering the cost of resistance and the competitive release,S cells are set to be defectors whileR i ’s are cooperators, just as in the single drug model withH,S andR. Also, it is obviously more interesting to let each pair ofR 1 andR 2 individuals play an asymmetric game, and we design the payos between R 1 andR 2 in a way thatR 1 cells are dominant subpopulations in a long term according to the adjusted replicator dynamics. Thus the payo matrix,A, forS,R 1 andR 2 are assumed to satisfy the Prisoner’s dilemma inequalities: a 21 <a 11 <a 22 <a 12 ; a 31 <a 11 <a 33 <a 13 ; a 32 <a 22 <a 33 <a 23 (4.3) with the extra conditions that enables sensitive cells to have a higher expected tness than that of each resistant cells: a 13 >a 23 ; a 13 >a 33 : (4.4) 81 In addition, noting that the selection-dependent expected tness function, f X i;j , is equal to 1 when w X = 0, we force all entries ofA, which correspond to payos that resistant cells receives, to be bigger than the background tness, 1, so that it allows the competitive release of resistant cells when the sum of two drugs is intense and the tness of sensitive cells are as little as 1: a 21 ;a 22 ;a 23 1; a 31 ;a 32 ;a 33 1: (4.5) For the numerical purpose, we specify the payo matrix,A, that is assumed to satisfy the inequalities in (4.3) - (4.5) as follows: A = 0 B B B B B B @ S R 1 R 2 S 2 2:8 2:8 R 1 1:5 2:1 2:3 R 2 1:5 1:8 2:2 1 C C C C C C A : (4.6) For each (i;j), the Moran process can change its state into 6 neighborhood in N , dened in (1.24), in one time step, besides the inactivity at the current state. Each transition probabilities are as follows: T SR 2 i;j = if S i;j N w i;j j N ; T SR 1 i;j = if S i;j N w i;j Nij N ; T R 2 S i;j = jf R 2 i;j N w i;j i N ; T R 2 R 1 i;j = jf R 2 i;j N w i;j (Nij) N ; T R 1 S i;j = (Nij)f R 1 i;j N w i;j i N ; T R 1 R 2 i;j = (Nij)f R 1 i;j N w i;j j N ; (4.7) 82 andT cons i;j = 1 (T SR 2 i;j + +T R 1 R 2 i;j ), where the weighted population,N w , is dened byN w =if S i;j + (Nij)f R 1 i;j +jf R 2 i;j . These microscopic probabilistic movement eventually determine the evolution of cancer cells. In measuring tumor size,x tumor , following the notation in [30], an independent important tool aside from the stochastic Moran process is the tumor-growth equation: _ x tumor = (<f >g)x tumor ; (4.8) whereg is the constant background tness, and< f > is the averaged tness in an innite population. Precisely, letting~ x := (x S ;x R 1 ;x R 2 ) | be the frequency vector ofS;R 1 ;R 2 subpopulations, the averaged tnessf X ofX,X2fS;R 1 ;R 2 g is given by: f S = 1w S +w S (A~ x) 1 ; f R 1 = 1w R 1 +w R 1 (A~ x) 2 ; f R 2 = 1w R 2 +w R 2 (A~ x) 3 (4.9) withx S +x R 1 +x R 2 = 1. The nonlinear averaged tness function,< f >, in the entire population is dened to be <f >=f S x S +f R 1 x R 1 +f R 2 x R 2 : (4.10) First, we consider situation in which both drugs are constantly administered, and those two drugs are additive, that is,e = 0. We assume an innite size of a population for a while in order to understand the expected dynamic, the adjusted replicator dynamic: 83 _ x S = f S <f > <f > x S ; _ x R 1 = f R 1 <f > <f > x R 1 ; _ x R 2 = f R 2 <f > <f > x R 2 : (4.11) The system has a dierent set of the evolutionary stable states depending on the concentration com- bination of two drugs. It is numerically depicted in Figure 4.1 for three dierent constant chemotherapy combinations of drug 1 and drug 2. When no drug is delivered, that is,C 1 = 0 andC 2 = 0, the tumor saturates to the S corner regardless of the initial distribution of three subpopulations (Fig.4.1a). When C 1 = 0:8 and C 2 = 0, the high chemo dose of drug 1 causes the competitive release of the resistant subpopulationsR 2 to drug 1 and all trajectories are driven to theR 2 corner (Fig.4.1b). UnderC 1 = 0 and C 2 = 0:8, the competitive release and of the resistant subpopulationsR 1 to drug 2 is caused as a result of administering high chemo dose of drug 2 and all trajectories are instead driven to theR 1 corner (Fig.4.1c). (a) (b) (c) (d) Figure 4.1: Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted repli- cator system for dierent constant chemotherapy values with e = 0. (a) Under C 1 = 0 and C 2 = 0, the tumor saturates to theS corner regardless of the initial distribution of the three subpopulations. (b) UnderC 1 = 0:8 andC 2 = 0, the competitive release of the resistant subpopulationsR 2 to drug 1 drives all trajectories to theR 2 corner. (c) UnderC 1 = 0 andC 2 = 0:8, the competitive release of the resistant subpopulationsR 1 to drug 2 drives all trajectories to theR 1 corner. (d) Trajectories with three dierent constant chemotherapy combinations of drug 1 and drug 2 overlap at dierent times and generate a closed loop. Knowing that the adjusted replicator dynamic is the limiting system of the Moran process, it is expected that the stochastic Moran process whose fate is determined by the transition probabilities in (4.7) shows 84 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 4.2: Realizations of multiple trajectories associated with the Moran process under administration of dierent constant chemotherapy combinations of drug 1 and drug 2 (e = 0) show the ability of the stochastic system to behave similarly to what the adjusted replicator dynamics drive, getting closer as the population size increases. For example, the Moran process, starting at a state near the corner R 2 with C 1 (t) 0,C 2 (t) 0 (blue wiggled lines), possibly evolves and attains the homogeneous population of all S for eachN, having smoother trajectories asN increases and nally being similar to the deterministic trajectory (light blue line). (a)C(t) 0,N = 1K; (b)C(t) 0,N = 5K; (a)C(t) 0,N = 10K; (a) C(t) 0,N = 50K; (e)C(t) 0:7,N = 1K; (f)C(t) 0:7,N = 5K; (g)C(t) 0:7,N = 10K; (h) C(t) 0:7,N = 50K more similar trajectory to the deterministic one as the population size increases for the same constant chemotherapy schedule of two drugs. According to the relation between the evolutionary step,, for the Moran process of size,N, and the evolution time,t, for the deterministic system as in (1.17), the number of steps that is under a chemotherapy schedule is determined. For example, a constant combination of drug 1 and drug 2 are administered duringtN steps for the stochastic Moran process if the same concentration of drugs are delivered fort unit time for the deterministic replicator system. One realization of the Moran process is presented in Figure 4.2 under three drug combinations of two drugs: (i)C 1 = 0;C 2 = 0, (ii) 85 C 1 = 0:8;C 2 = 0, or (iii)C 1 = 0;C 2 = 0:8 for eachN = 1K, 5K, 10K and 50K. For a largeN = 50K, it is relatively easy to nd a stochastic trajectory that moves tightly to the deterministic trajectory for all three drug combinations. Each realization under C 1 = 0;C 2 = 0 often ends its journey at the S corner which is the evolutionary stable strategy of the deterministic system. It is true for the other two combinations. For a smallN, it is not impossible to obtain a trajectory that moves as the deterministic trajectory does though it is wiggled and o the solid deterministic line and that arrives at theESS of the deterministic system. However, it hardly takes place in general for as smallN as 1K that the xation to thatESS occurs with a positive probability for all drug combinations. Moreover, overlapping gures, 4.1a, 4.1b and 4.1c, we see that trajectories with three dierent constant chemotherapy combinations of drug 1 and drug 2 intersect at dierent times and eventually generate a closed loop in a phase space (Fig.4.1d). On the numerical purpose, we xO(0:617, 0:321, 0:062),P (0:44, 0:13, 0:43) andQ(0:807, 0:051, 0:142), and consider a deterministic trajectory starting atO, governed by the replicator dynamic under an adaptive schedule. This chemotherapy schedule delivers drug 1 at the constant concentration,C 1 = 0:8, with no drug 2 forT OP := 6:933 unit time, turns both drugs o for T PQ := 6:248 unit time, and administers drug 2 at the constant concentration, C 2 = 0:8, with drug 1 missing forT QO := 6:324 unit time successively (Fig.4.3a). We call the sum-up time,T OP +T PQ +T QO = 19:2653, one evolutionary cycle, and use the term 1=3 cycle to refer to each piece of time in a cycle for convenience. As can be seen in the gure 4.3b, under this adaptive chemotherapy schedule the adjusted replicator system starting at O moves along a red line under C 1 = 0:8 and C 2 = 0 and arrives at P after T OP unit time. Turning chemotherapy o atP , the system begins to have more sensitive cells in a population and reachesO, moving along a blue line, afterT PQ unit time sinceT OP . Again turning chemotherapy on being underC 1 = 0 andC 2 = 0:8 atO, it nally evolves along a green line and reaches backO in T QO unit time afterT OP +T PQ , generating a closed loop,OPQO, for one evolutionary cycle. When no 86 (a) (b) (c) Figure 4.3: Switching constant chemotherapy combinations of two drugs (e = 0) on and o at adequate times traps a trajectory, associated with the adjusted replicator system, within a closed loop, controlling a tumor size. (a)OP : C 1 = 0:8,C 2 = 0 duringT OP = 6:933 unit time,PQ: C 1 = 0,C 2 = 0 during T PQ = 6:248 unit time,QO:C 1 = 0,C 2 = 0:8 duringT QO = 6:324 unit time; (b) the trajectory treated according to the multi drug additive adaptive schedule and the untreated trajectory (pink) being driven to theS corner; (c) tumor size under untreated (pink) and the adaptive chemotherapy schedule. (For this numerical experiment, we takeg = 1:5519.) drug is administered, the system in the stateO starts rapidly having more sensitive cells and sensitive cells nearly saturates in a population taking 99:56% within one evolutionary cycle (Fig.4.3b). It is then clear that if both drugs are administered according to the adaptive chemotherapy schedule for longer than 1 cycle, even when 2 cycles, the saturation of cancerous cells is preventable while it would if no drugs were delivered for the same time period. The ecacy of the application of the adaptive chemotherapy schedule can be seen in terms of tu- mor volume that is given by the equation in (4.8). It is straightforward that the following equation of an exponential form solves that tumor-growth equation: x tumor (t) =x tumor (0) exp( Z t 0 (<f >g) dt): (4.12) Also, it is clear that the solution in (4.12) isT -periodic ifg is set to be as follows: g = Z T 0 (<f >g) dt: (4.13) 87 For the simulation of our specic model with the payo matrix,A, in (4.6), we setT = T OP +T PQ + T QO to be one cycle andg = R T 0 (<f >g) dt as in (4.13) so that the tumor volume returns back to the initial size after one evolutionary cycle when the adaptive chemotherapy schedule is applied. Using the Euler method to quantify g followed by some manual adjustment, such g turns out to be equal to 1:5519 for the additive drug interaction case. Setting the initial tumor size to be 1 for simplicity, we see that the rst 1=3 cycle with a high dose of drug 1 reduces the tumor size, the second 1=3 cycle with no drug used allows the regrowth of tumor, and the last 1=3 cycle with a high dose of drug 2 mitigates the cancerous condition and resumes the initial tumor size. Continuing the drug use repeatedly along with this adaptive chemotherapy schedule theoretically controls both the state of patients’ condition by letting their corresponding deterministic trajectories trapped in aS,R 1 ,R 2 phase space and the tumor volume going through a series of decreases and an increase. This is a great advantage of designing a chemotherapy schedule compared to when it is untreated. In fact, the tumor size explosively increases if untreated, and the initial tumor size begin 1 reaches over 10 even in the rst 1=3 evolutionary cycle while it decreases reaching nearly 0 when the adaptive chemotherapy schedule is applied (Fig.4.3c). 4.2 Adaptive control of evolutionary cycles with additive multi drug schedule Now we apply the adaptive chemotherapy schedule (Fig.4.3a) associated with the deterministic adjusted replicator dynamics to the Moran process for a few cycles, and observe for a largeN how fast the stochas- ticity does not justify its application in an additive (e = 0) multi drug case by increasing the number of rounds. As seen in Figure 4.2 for constant drug use, the stochastic trajectory evolve around the deter- ministic trajectory and it is xated to the correspondingESS of the deterministic system with nonzero probability. However, it is just one realization and does not represent the overall behavior of the stochastic 88 process. Like we did in the previous chapter, we simulate 1; 000 realizations of the Moran process in order to understand its behavior through their distribution and their sample mean. (a) (b) (c) (d) Figure 4.4: The spread of the distribution of points aroundP (Q orO) for 1, 000 realizations of the stochastic Moran process gets denser and demonstrates the shrunken randomness as the population size increases when each realization is under the administration of a constant chemo combinationC 1 (t) 0:8;C 2 (t) 0 (C 1 (t) 0;C 2 (t) 0 orC 1 (t) 0;C 2 (t) 0:8) withe = 0 during 1=3 evolutionary cycle,T OP (T PQ or T QO ), since its exact start atO (P orQ). (a)N = 1K; (b)N = 5K; (c)N = 10K; (d)N = 50K First, for eachN = 1K, 5K, 10K and 50K, we apply the adaptive chemotherapy schedule to 1; 000 realizations of the Moran process of size N during 1=3 evolutionary cycle when each starts at a given point, O, P or Q. The Moran process starting at O under a high constant drug 1 schedule, C 1 = 0:8 andC 2 = 0 duringT OP N evolutionary steps ends its journey around the pointP where the adjusted replicator system settles down after the rst 1=3 evolutionary cycle. Each of terminal points of 1; 000 realizations is indicated as a red dot and they form a cluster nearP in Figure 4.4. Likewise, for the second 1=3 evolutionary cycle, each starts at a xed pointP and evolves as untreated duringT PQ N evolutionary steps. All of them arrives at a point nearQ, forming a blue cluster aroundQ. For the third 1=3 cycle, all the realizations start atQ and their destiny under a high constant drug 2 schedule,C 1 = 0 andC 2 = 0:8, is shaped aroundO as a green cloud. It infers that the tumor size at the last evolutionary step of each 1=3 cycle driven by the Moran process would be more or less than one associated with the deterministic system (Fig.4.3c). Obviously, all distributions around O, P and Q are more spread for a small N and they get more centered at the corresponding point as N increases. We zoom in all these clouds and plot them in the principal axis coordinate system (Fig.4.5). We observe that for each xed constant chemotherapy schedule, 89 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 4.5: The spread of the distribution of points in the principal axis coordinate system for 1; 000 real- izations of the Moran process is, for a large population size, characterized as a multivariate Gaussian distri- bution aroundP (Q orO) when each realization is under the administration of a constant chemotherapy combinationC 1 (t) 0:8;C 2 (t) 0 (C 1 (t) 0;C 2 (t) 0 orC 1 (t) 0;C 2 (t) 0:8) withe = 0 during 1=3 evolutionary cycle, T OP (T PQ or T QO ), since its exact start at O (P or Q). As the population size increases, both the semi-major axis, 1 and the semi-minor axis, 2 , decrease. (a)N = 1K, 1 = 0:0977, 2 = 0:0296; (b)N = 5K, 1 = 0:0460, 2 = 0:0127; (c)N = 10K, 1 = 0:0320, 2 = 0:0091; (d) N = 50K, 1 = 0:0138, 2 = 0:0040; (e)N = 1K, 1 = 0:0318, 2 = 0:0120; (f)N = 5K, 1 = 0:0140, 2 = 0:0085; (g)N = 10K, 1 = 0:0097, 2 = 0:0060; (h)N = 50K, 1 = 0:0044, 2 = 0:0027; (i) N = 1K, 1 = 0:0860, 2 = 0:0207; (j)N = 5K, 1 = 0:0402, 2 = 0:0093; (k)N = 10K, 1 = 0:0273, 2 = 0:0065; (l)N = 50K, 1 = 0:0125, 2 = 0:0030 they are distributed following a multi Gaussian with the decreasing semi-major and semi-minor axes in the increase ofN. More interestingly, the semi-major and semi-minor axis nearQ are signicantly small compared to those around eitherP orO, being approximately three times less. It is true that the expected deterministic system converges, for all three constant drug combinations, toS underC 1 = 0;C 2 = 0,R 2 underC 1 = 0:8;C 2 = 0, andR 1 underC 1 = 0;C 2 = 0:8. However, the expected terminal point around 90 Q after the second 1=3 cycle is much closer to itsESS, theS corner, and it leads to the smaller semi-major (-minor) axis. (a) (b) (c) Figure 4.6: The stochastic trajectory of one realization of the Moran process under the administration of the additive adaptive chemotherapy, associated with the adjusted replicator dynamics, as in Figure 4.3a develops a random walk along a lattice in a phase space,S 2 . The adaptive schedule is able to prevent the stochastic system from the saturation of cancer cells, even in a small population withN = 30. (a)N = 30; (b)N = 40; (c)N = 50 We now extend the 1=3 cycle to one full cycle during which the adaptive chemotherapy schedule is applied to the Moran process who starts atO. All gures in the gure 4.6 visualize the random walk of the stochastic process along grids in a phase space for a smallN asN = 30, 40 and 50. Iterating states of the process that starts atO underC 1 = 0:8 andC 2 = 0, it travels along red segments while either arriving one of the nearest 6 grid points or staying at the same point in one step, and it nally reaches a pink point afterT OP N evolutionary steps. Turning both drugs o at that pink points drives the process to move along blue segments and arrive at a light blue point inT PQ N steps. When a high chemo schedule of drug 2 asC 1 = 0 andC 2 = 0:8 is constantly delivered, the light blue point initiate its move along a green line, terminating its journey at a light green point inT QO N steps. For these smallN, the Moran process rarely return to the neighborhood of the initial point,O, with a high probability. However, with several trials, at least one realization that shows the approximate return to the original point under the adaptive chemotherapy schedule was made. Increasing the population size, it is more likely that the trajectory of the stochastic process under the adaptive schedule is shaped tightly to the closed loop, OPQO, which the associated deterministic 91 (a) (b) (c) (d) (e) (f) Figure 4.7: The averaged trajectory of 1; 000 realizations of the Moran process under the additive adaptive schedule, associated with the adjusted replicator system, during one round (T OP +T PQ +T QO = 19:2653 unit time) ts the corresponding deterministic trajectory for a large population size withe = 0. The Moran process is likely to return nearly to the initial state with a high probability for a largeN even though the spread of the distribution of the points nearQ is still wide. (a) the distribution of the points associated with the adaptive chemo schedule forN = 10K; (b) the trajectory of one single realization of the Moran process withN = 10K; (c) the averaged trajectory of 1; 000 realizations of the Moran process withN = 10K; (d) the distribution of the points forN = 50K; (e) the trajectory of one single realization withN = 50K; (f) the averaged trajectory withN = 50K system generates for one full evolutionary cycle, still being o from it. We examine how well or badly the adaptive schedule allows the stochastic system to resume its initial frequency after one evolutionary cycle for 1; 000 realizations of the Moran process. Each realizations starts atO and both the second and the third 1=3 evolutionary cycles start at the point where their previous 1=3 cycle ends. Our simulation of 1; 000 realizations of the Moran process is given with the population sizeN = 10K and 50K in Figure 4.7. For eachN, each simulation starts atO and it proceeds underC 1 = 0:8 andC 2 = 0 for the rst 1=3 evolutionary cycle, arriving at one of red points nearP (Fig.4.7a,4.7d). At the red point, both drugs starts being taken o and it is led to one of blue points nearO by no drug strategy at the end of the second 1=3 cycle. The blue point is nally transferred to one of green point 92 nearO underC 1 = 0 andC 2 = 0:8 after the last 1=3 cycle. Indicating trajectories along grids under the dierent constant drug combinations in red, blue and green in order, the trajectory of one such realization is given in Figure 4.7b and the gure 4.7e for N = 10K and N = 50K, respectively. Both stochastic trajectories are drawn around the deterministic trajectory,OPQO, though it ts better the closed loop, OPQO, with the bigger population size. In order to understand the general behavior of the stochastic process with the xed size, simulation of 1; 000 such independent Moran process is carried out and it produces the distribution of points aroundO, P andQ. It is obviously expected that the spread of the distribution of points is denser with the bigger population as can be seen in Figure 4.7a and Figure 4.7d. However, even with this bigger deviation in distribution whenN = 10K, the averaged trajectories in Figure 4.7c of those 1; 000 realizations of the Moran process already well represents its limiting system. One other thing to point out in Figure 4.7a and Figure 4.7d is when each of these is compared to Figure 4.4c and Figure 4.4d with the same population size. Points in Figure 4.7a and Figure 4.7d are more widely spread aroundQ andO. It is because, for example, the distribution of points nearQ in Figure 4.4 is obtained when the adaptive chemotherapy schedule is applied to the Moran process that starts at the xed point, P , in order to observe the short term inuence duringT PQ N evolutionary steps. However, Figure 4.7 is given with the more extended perspective that it is for understanding the longer term inuence of the chemotherapy schedule, thus, the stochastic process that ends at one of red dots nearP starts at that point for the next 1=3 cycle instead of beginning at the exact point,P . This perturbation in the initial state for the second or the third 1=3 cycle results in the greater deviation of the distribution of points aroundQ or O. This is numerically better supported in terms of the semi-major (-minor) axis of the distribution of points aroundO when they are plotted in the principal axis coordinate system as shown in Figure 4.8. With N = 10K, the semi-major axis is equal to 1 = 0:0273 and the semi-minor axis is equal to 2 = 0:0065 93 (a) (b) (c) (d) Figure 4.8: The spread of the distribution of points for 1; 000 realizations of the stochastic Moran process is characterized as a multivariate Gaussian distribution, centered nearly at the initial point,O, when each realization is under the administration of the multi drug adaptive chemo schedule associated with the adjusted replicator dynamics during one round (T OP +T PQ +T QO = 19:2653 unit time) withe = 0. The mean frequency, R 1 (or R 2 ), of the subpopulation, R 1 (orR 2 ), around the point, O, converges to the proportion ofR 1 (orR 2 ) asN increases, with the decreasing semi-major axis, 1 , and semi-minor axis, 2 . (a)N = 10K, R 1 = 0:3175 , R 2 = 0:0617; (b)N = 10K, 1 = 0:0453, 2 = 0:0099; (c)N = 50K, R 1 = 0:3177, R 2 = 0:0617; (d)N = 50K, 1 = 0:0211, 2 = 0:0045 for the distribution of points nearO when the third piece of the adaptive chemotherapy schedule is applied to the Moran process starting atQ as shown in Figure 4.5k while 1 = 0:0453 and 2 = 0:0099 at the end of the full one cycle of the adaptive schedule as in Figure 4.8b. It is true forN = 50K that the deviation of the distribution of points at the end of one full cycle of the adaptive chemotherapy schedule is bigger than one in the last 1=3 cycle due to the accumulated randomness in the initial points in the last two 1=3 cycles. Moreover, the distribution of points aroundO driven by the application of the adaptive chemother- apy schedule during one full evolutionary cycle closely follows the multi Gaussian centered closely toO (Fig.4.8b, 4.8d). In fact, the mean frequency, R 1 , of theR 1 subpopulations is equal to R 1 = 0:31749 and the mean frequency, R 2 , of theR 2 subpopulations is equal to R 2 = 0:061724 forN = 10K while R 1 = 0:31766 and R 2 = 0:061714 forN = 50K (Fig.4.8a,4.8c). Knowing the proportion of subpopu- lations, (S;R 1 ;R 2 ) = (0:617; 0:321; 0:062), atO, it is shown that the center of the distribution is closely located toO with the reduced norm for the bigger population size. Furthermore, plotting the points in the R 1 R 2 coordinate system shows that theR 1 axis almost coincides with the major principal axis and theR 2 axis with the minor principal axis, explaining that the continuous delivery of drug 2 in high concentration 94 in the last 1=3 evolutionary cycle mostly forces the trade-o between the sensitive cellsS and the resistant cellsR 1 to drug 2. (a) (b) (c) (d) Figure 4.9: Two realizations of the stochastic Moran process of sizeN = 10K, starting atO and evolving under the adaptive chemo schedule during 8 rounds withe = 0, show a great dierence in their tumor size as well as in their trajectories. (a) one realization; (b) corresponding tumor size to Figure 4.9a; (c) another realization; (d) corresponding tumor size to Figure 4.9c We have just seen that the stochastic Moran process with as a big population size,N, asN = 10K or 50K returns closely to the initial state on average in one evolutionary cycle of the adaptive chemother- apy schedule. It implies that the deduced tumor size using the tumor-growth equation in (4.8) is overall controlled undertaking a series of increases and decreases as the tumor size associated with the adjusted replicator system drops during the rst 1=3 cycle, rises up when untreated for the second 1=3 cycle, and reduces during the third 1=3 cycle, nally recovering the initial tumor size at the end of the one full evolu- tionary cycle. In order to answer the justication of the use of drugs according to the adaptive chemother- apy schedule, we continued this adaptation for 8 evolutionary cycles to the Moran process of the size N = 50K. It starts atO for the rst round, and each round starts at the state where the previous round ends. For example, one realization of the Moran process starting atO terminates at one of green points at the end of one evolutionary cycle, which we call a round interchangeably. Then, the second round starts at that point instead of the exact point,O, and the following round begins in the same manner. Two realizations of such application to the Moran process of sizeN = 10K whose trajectories do not converge to any of cancerous states within 8 cycles are presented in Figure 4.9. Figure 4.9a and Figure 4.9c show that their trajectories maintain a triangle-like shape in each cycle without being xated to a 95 (a) (b) (c) (d) Figure 4.10: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed until the end of 8th round when each realization evolves under the administration of the multi drug adaptive chemo schedule, associated with the adjusted replicator system, withe = 0 during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.10a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.10c cancerous state though they closely evolve around the deterministic trajectory,OPQO, only up to rst a few rounds. Nonetheless, the tumor size associated with Figure 4.9a does not have a monotonically increasing maximum value of each cycle and gets doubled in the 7th round while one associated with Figure 4.9d incredibly soaks being over 500 in the same round. Likewise, each realization of the stochastic process results in apparently dierent dynamic and understanding stochastic process in distribution is needed. (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.11: The spread of the distribution of terminal points (green dots) around O in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug adaptive chemo schedule, associated with the adjusted replicator system, withe = 0 during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 96 One another realization of the Moran process with a bigger population sizeN = 50K whose trajectory moves more closely around the deterministic trajectory,OPQO, for the entire cycles is given in Figure 4.10a along with the associated tumor size in Figure 4.10b. For this bigger population size, the maximum of the tumor size in each cycle gradually increases as the number of rounds increases for this single simu- lation, and this pattern turns out to be true for the averaged trajectory of 1; 000 realizations of the Moran process and its related tumor size whenN = 50K (Fig.4.10c,4.10d). The trajectory under the adaptive chemotherapy schedule gradually shifts towards theS corner on average and it leads to the rise in the maximum of the tumor size attained at the end of the second 1=3 cycle as the number of cycle increases. Finally, the tumor size that is initially equal to 1 reaches over 15 within 8 cycles though it is reduced to under 5 at the end of the 8th cycle. Focusing on the spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of size N = 50K, we indicate them as green points at the end of each cycle (Fig.4.11). They form a cloud each cycle, and it gets less thick as the number of rounds increases and some of them mostly consisted of sensitive cells. It is clear that the spread is relatively small inR 2 subpopulations but the randomness is more reected in the proportion ofS orR 1 subpopulations. This spread is measured through the semi-major and the semi-minor axis by plotting the points in each cycle in the principal axis coordinate system (Fig.4.12). The spread overall follows the multi Gaussian distribution at the end of the rst round, but it gradually starts being distorted from the following round, having the center of the cloud o the most dense part. The semi-major axis which is equal to 1 = 0:0211 in the rst round increases roughly by 6 times in 8 cycles while the semi-minor axis being equal to 2 = 0:0045 in the rst round increases by 3 times. The more rapid rate of increase of the semi-major axis in the number of rounds is also described in Figure 4.13a and Figure 4.13c for the population sizesN = 10K and N = 50K, respectively, where darker dots represent the semi-major axis and lighter green dots stand for the semi-minor axis. As seen in these raw data, both metrics increase nonlinearly in the number of rounds, 97 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.12: The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K, all starting atO and evolving under the multi drug adaptive chemo schedule withe = 0 during 8 rounds, shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , in the principal axis coordinate system as the number of rounds increases. (a) round 1, 1 = 0:0211, 2 = 0:0045; (b) round 2, 1 = 0:0319, 2 = 0:0068; (c) round 3, 1 = 0:0425, 2 = 0:0084; (d) round 4, 1 = 0:0561, 2 = 0:0096; (e) round 5, 1 = 0:0710, 2 = 0:0108; (f) round 6, 1 = 0:0871, 2 = 0:0116; (g) round 7, 1 = 0:1051, 2 = 0:0128; (h) round 8, 1 = 0:1243, 2 = 0:0136 denitely with the fact that both of them are bigger with a smaller population size. For each population size, we could see that the semi-major (-minor) axis has the power-law dependency on the number of rounds, and the tted curves to the semi-major axis, 1 , and the semi-minor axis, 2 , withN = 10K are given below in the log-log scale: 1 0:0412n 0:7782 ; 2 0:0105n 0:4967 ; (4.14) wheren is the number of rounds. Similarly, withN = 50K, 1 0:0186n 0:8602 ; 2 0:0046n 0:5235 : (4.15) 98 (a) (b) (c) (d) Figure 4.13: The semi-major (and -minor) axis of the distribution of the points aroundO at the end of each round for 1; 000 realizations of the Moran process overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the multi drug adaptive chemo schedule associated with the adjusted replicator system withe = 0 during 8 rounds since its exact start atO. (a)N = 10K; (b)N = 10K, the log-log t; (c)N = 50K; (d)N = 50K, the log-log t From these equations in (4.14) and (4.15), we check that the semi-major axis has a bigger increase rate than the semi-minor axis does for eachN, meaning that the spread of the distribution of points nearO at the end of each cycle is made more rapidly in the direction of the major principal axis than towards the both direction of the minor principal axis. It is also observed that for eachi the increase rate of i with N = 50K is bigger in the log-log scale, however, i is a dominating function withN = 10K during 8 evolutionary cycles in the raw scale. In other words, despite of the fact that the distribution becomes more rapidly disseminated with a bigger population sizeN = 50K, the absolute spread withN = 50K is not as larger as one withN = 10K due to the stronger randomness from a smaller population size. 4.3 Adaptivecontrolofevolutionarycycleswithsynergisticandantagonistic multidrugschedules The multi drug adaptive chemotherapy schedule helps delay the saturation of cancerous cells on average for the 1; 000 realizations of the Moran process and it keeps tumor size well controlled compared to the case under no treatment when two drugs are additive for a few evolutionary cycles. In fact, the averaged trajectory circles around the closed loop,OPQO, that the adjusted replicator system associated with the adaptive schedule produces during one evolutionary cycle when the adaptive schedule is designed as if 99 each drug concentration function,C i (t), is a step function between two dierent levels, 0:8 and 0, across over time (Fig.4.3a). In this section, we have a similar question: how long and how well/badly the adaptive schedule overall helps postpone the direct increase of tumor volume, which would be attained if untreated, when a nite population is modeled through the Moran process with a large enough size and, more im- portantly, when two drugs act either synergistically or antagonistically. The drug interaction parameter, e, ranging from1 to 1 is what it determines if the action of two drugs are antagonistic, additive or syn- ergistic. Recall that two drugs are additive whene = 0, two drugs interact synergistically whene > 0 and they interact antagonistically whene < 0. For simulation, we x the value ofe as suche = 0:3 for synergistic ande =0:3 for antagonistic drug interactions. (a) (b) (c) Figure 4.14: Deterministic trajectories describe the evolutionary stable states (ESS) of the adjusted repli- cator system for dierent constant chemotherapy values for each drug interaction. UnderC 1 = 0 and C 2 = 0, the tumor saturates to theS corner regardless of the initial distribution of the three subpopula- tions. UnderC 1 = 0:5 andC 2 = 0:2, the competitive release of the resistant subpopulations,R 2 , to drug 1 drives all trajectories to theR 2 corner. UnderC 1 = 0:2 andC 2 = 0:5, the competitive release of the resistant subpopulations,R 1 , to drug 2 drives all trajectories to theR 1 corner. (a)e = 0; (b)e = 0:3; (c) e =0:3 We also consider a new combination of two drugs of constant concentration such asC i (t) 0:2 or C i (t) 0:5 instead of the combination of C i (t) 0 and C i (t) 0:8 since this combination leads to trajectories of the adjusted replicator dynamics traversing more widely in a phase space over time when they start at the same set of initial points for each value ofe that we set, and hence it allows us to more easily come up with an shared initial distribution at which the adaptive chemotherapy schedule starts its play. 100 We present the numerical result for additive drug interaction case as well as synergistic and antagonistic cases for comparison under this new combination of drug concentration. The numerical simulation describes evolutionary stable states of the adjusted replicator system for dierent constant chemotherapy combinations of drug 1 and drug 2 in Figure 4.14. For all the values of drug interaction parameter,e, such ase = 0;0:3, the deterministic system is driven to theS corner regardless of the initial distribution of the three subpopulations when no treatment is carried out. The tumor saturates to theR 2 corner when a higher concentration of drug 1 is constantly administered, precisely whenC 1 = 0:5 andC 2 = 0:2, and this results from the competitive release of the resistant subpopulationsR 2 to drug 1. In the same manner, the competitive release of the resistant subpopulations R 1 to drug 2 drives all trajectories to theR 1 corner when a higher concentration of drug 2 is continuously delivered. With the xed drug use, trajectories are distinguished by having a dierent unstable xed point for eache though the deterministic system have the unique evolutionary stable state regardless of the value ofe considered. For example, the unstable xed point composed of onlyS andR 1 subpopulations, which is indicated as a red open dot in Figure 4.14, underC 1 = 0:5 andC 2 = 0:2 represents a state that consists of moreS subpopulations whene = 0:3. This keeps consisting of moreR 1 subpopulations ase decreases, and it eventually merges into the homogenous state,R 1 , before it reachese =0:3. Similarly, the unstable xed point with noR 1 population, which is marked as a green open dot, underC 1 = 0:5 andC 2 = 0:2 refers to a state having a higher proportion ofS subpopulations whene = 0:3 while a smaller frequency ofS whene =0:3. Basically, the drug interaction parameter,e, plays the role of relocation of an unstable xed point, not changing the nature of the system in terms of the evolutionary stability with this new drug combination. It is also seen in Figure 4.14 that for eache, trajectories with three dierent constant chemotherapy combinations of drug 1 and drug 2 overlap at dierent times and generate a clear-cut closed loop. Fix- ing a pointO(0:812; 0:123; 0:065), we consider two drug adaptive chemotherapy schedules by which the 101 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 4.15: Deterministic evolution of subpopulations,S,R 1 andR 2 , by the adjusted replicator system generates a closed loop,OPQO, when it starts atO and evolves under the multi drug adaptive chemother- apy schedule during one round, resulting in the tumor size controlled. (a) Withe = 0, the total dose de- livered to generate one evolutionary cycle is 17:7604 during 29:4440 unit time.OP :C 1 = 0:5,C 2 = 0:2 duringT OP = 14:226 unit time,PQ: C 1 = 0,C 2 = 0 duringT PQ = 4:072 unit time,QO: C 1 = 0:2, C 2 = 0:5 duringT QO = 11:146 unit time; (b) the corresponding deterministic trajectory; (c) the corre- sponding tumor size with the averaged background tness,g = 1:4527; (d) Withe = 0:3, the total dose delivered to generate one evolutionary cycle is 17:7723 during 32:29 unit time.OP :C 1 = 0:5,C 2 = 0:2 duringT OP = 14:260 unit time,PQ: C 1 = 0,C 2 = 0 duringT PQ = 6:901 unit time,QO: C 1 = 0:2, C 2 = 0:5 during T QO = 11:129 unit time; (e) the corresponding deterministic trajectory; (f) the cor- responding tumor size with g = 1:4857; (g) With e =0:3, the total dose delivered to generate one evolutionary cycle is 17:8150 during 26:8740 unit time. OP : C 1 = 0:5,C 2 = 0:2 duringT OP = 14:300 unit time, PQ: C 1 = 0, C 2 = 0 during T PQ = 1:424 unit time, QO: C 1 = 0:2, C 2 = 0:5 during T QO = 11:150 unit time; (h) the corresponding deterministic trajectory; (i) the corresponding tumor size withg = 1:4182 102 adjusted replicator dynamic that starts at the exactO generates a closed loop,OPQO, in its own evolu- tionary cycle for eache as shown in Figures 4.15b, 4.15e, 4.15h. For all these values ofe, drug delivery underC 1 = 0:5 andC 2 = 0:2 causes the deterministic system to traverse a red line and reachesP inT OP unit times. When it is switched toC 1 = 0 andC 2 = 0 atP , the system reachesQ moving along a blue line inT PQ unit times. Then drug starts being delivered in the concentrationC 1 = 0:2 andC 2 = 0:5, and it results in the arrival of the system atO inT QO unit times. Note that the initial stateO is directed to the S corner for all drug interaction cases if it is chemotherapeutically untreated as seen along with a closed loop,OPQO. However, the two drug adaptive chemotherapy schedule creates a closed orbit in the trajectory of the deterministic system during one evolutionary cycle avoiding the saturation of tumor though the time periods,T OP ,T PQ orT QO , under a piece of xed constant drug combination all dier from each other for alle. These schedules with the exact time periods for each legs ofOP ,PQ andQO are precisely given for each case in Figures 4.15a, 4.15d, 4.15g. We agree that eache shares the same starting pointO but neither P nor Q as well as neither T OP , T PQ nor T QO . However, since no confusion arises, we use the same notations for those points at which the chemotherapy schedule changes drug concentration, and specify if needed. e 0:3 0 0:3 D 17:8150 17:7604 17:7723 T 26:8740 29:4440 32:2900 D=T 0:6629 0:6032 0:5504 Table 4.1: total dose (D), total time (T), and average dose (D/T) associated with adaptive chemotherapy schedules with antagonistic, additive and synergistic drug interactions during one evolutionary cycle. Our attempt to come up with closed loopsOPQO that are created using, ideally, the same total dose during one evolutionary cycle for alle’s numerically ended up being with closed orbits created using not a xed but similar amount of total dose in dierent total times instead. All these quantities such as total 103 time,T (one evolutionary cycle), total dose,D, and the averaged dose,D=T , used for our simulation for one evolutionary cycle is given in Table 4.1 following the same notation in [30]. When two drugs are additive, the adaptive schedule (Fig.4.15a) uses 17:7604 dose in total from both drugs causing the system to be trapped in the closed loopOPQO during one evolutionary cycle, 29:4440. When two drugs interact synergistically, the schedule (Fig.4.15d) uses 17:7723 in total and the associated trajectory creates a bigger closed loop during the longer time period, 32:2900 unit times. In contrast, using 17:8150 in total when two drugs interact antagonistically, the schedule (Fig.4.15g) generates the smallest loop during the shortest time period, 26:8740. What is understood in Table 4.1 is that using similar total dose being equal to about 17:78, the closed loop is generated distinguishably faster when two drugs interact antagonistically and this leads to the highest averaged total dose during one evolutionary cycle. e 0:3 0 0:3 g 1:4182 1:4527 1:4857 Table 4.2: the constant background tness (g) associated with adaptive chemotherapy schedules for an- tagonistic, additive and synergistic drug interactions. Another observation is made in terms of tumor size for all drug interactions. Tumor size that is assumed to be equal to 1 initially grows or decays according to the tumor-growth equation in (4.8), where the constant background tness,g, is dened by the equation in (4.13). These background tnessg’s associated with the adaptive chemotherapy schedules in Figure 4.15a, 4.15d, 4.15g are computed in Table 4.2 for each drug interaction, and each of these values enable tumor size to recover the initial volume at the end of one evolutionary cycle going through a series of rises and declines. What is observed in all cases in common is that tumor size initially reduces when a higher chemo dose of drug 1 is administered, rises when no chemotherapeutic treatment is done, and decreases again when a higher chemo dose of drug 2 is delivered, reaching back the initial volume in one evolutionary cycle. That is, tumor size is controlled when the two drug adaptive chemotherapy schedule is applied during one evolutionary cycle and hence as many times as 104 desired. This is a great advantage of the use of the adaptive schedule since tumor size dramatically soaks if no drugs are administered. In fact, tumor size under no drug surpasses its maximum in the adaptive therapy scenario even earlier than the rst 1=3 cycle ends for all drug interaction cases. Though tumor is controlled in size as a result of the adoption of the adaptive chemo schedule, it is signicantly dierent in deviation. Compared to the additive drug interaction case, the variation in tumor size is greater when two drugs interact synergistically while it is smaller when two drugs interact antagonistically, and this is compatible with the size of the closed loop,OPQO. The bigger the closed loop is, the greater tumor size deviates during one cycle. (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.16: The spread of the distribution of terminal points (green dots) around O in each round for 1; 000 realizations of the stochastic Moran process with size N = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15a during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 Now, we apply the adaptive chemotherapy schedule associated with the deterministic adjusted repli- cator system to the Moran process of a large but nite population size N, N = 50K, starting at O, and gure out up to how many of 8 cycles this adoption is valid by comparing averaged trajectories and tumor sizes of 1; 000 realizations of the Moran process for antagonistic, additive and synergistic drug 105 interactions represented with parameters, e = 0:3, 0, 0:3, respectively. For the purpose of compar- ison, we will rst look into the additive drug interaction case. Note that the one evolutionary cycle, T := T OP +T PQ +T QO = 29:4440 unit times, associated with the adaptive chemo schedule in Fig- ure 4.15a whene = 0 corresponds to the evolutionary steps, := TN = 1; 472; 200, for the Moran process withN = 50K. In 1; 472; 200 evolutionary steps, 1; 000 realizations of the Moran process that start atO form a distribution of points nearO, at which each of them starts its following run during the same number of evolutionary steps. Applying the adaptive chemo schedule during 8 cycles in this manner, we get 8 individual spreads of distribution of points aroundO from each round, and those points are depicted as green dots in Figure 4.16. Since the adaptive schedule that we are applying is associated with the deterministic schedule that both starts and ends atO, the gradual larger deviation in distribution in the increase of the number of runs from the continued adoption of this schedule is denitely expected not only by the stochastic structure itself but also by the fact that the initial point of each round is not the exactO but a neighborhood from the second round. Plotting the points nearO in the principal axis coordinate system and computing the semi-major axis, 1 , and semi-minor axis, 2 , for each run gives more comprehensive understanding (Fig.4.17). Applying the adaptive schedule one round guarantees its use for the later round since the spread of the distribution of the terminal points aroundO is closely centered at the initial point,O, and, in addition it follows the multi Gaussian distribution. Until the 3rd round, the application of this chemo schedule seems to suitably control the stochastic dynamic allowing the system to return back to the starting point on average although both the semi-major axis and the semi-minor axis increase as the number of rounds goes. However, the distribution starts being o the initial point,O, from the 4th round and the spreads are not anymore densely disseminated around its mean with the increasing distance fromO, losing the Gaussian distribution form. 106 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.17: The spread of the distribution of points aroundO for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug additive (e = 0) adaptive chemo schedule as in Figure 4.15a during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0156, 2 = 0:0076; (b) round 2, 1 = 0:0268, 2 = 0:0123; (c) round 3, 1 = 0:0432, 2 = 0:0172; (d) round 4, 1 = 0:0675, 2 = 0:0223; (e) round 5, 1 = 0:1069, 2 = 0:0288; (f) round 6, 1 = 0:1623, 2 = 0:0371; (g) round 7, 1 = 0:2233, 2 = 0:0443; (h) round 8, 1 = 0:2797, 2 = 0:0474 One of those 1; 000 realizations of the Moran process to which the adaptive chemo schedule is applied during 8 cycles is randomly selected and its full trajectory is shown in Figure 4.18a. This realization chosen shows a similar trajectory to the associated deterministic trajectory with the adjusted replicator dynamic, depicted as a light blue line, up to the second round, but it keeps moving towards theR 1 corner as rounds go while maintaining a triangle-shaped trajectory in each cycle. We adopted the tumor-growth equation in (4.8) to describe the growth and the decay of tumor volume, and the associated tumor size with that single trajectory is given in Figure 4.18b along with the tumor size associated with the adjusted replicator dynamic in light blue. With this large population sizeN = 50K, the tumor size driven by the stochastic process almost coincides with the deterministic one during the rst evolutionary cycle. However, the maximum tumor size achieved at the end of the second piece, that is no drug part, keeps increasing and it nally reaches over 30 in the 5th cycle though the tumor size reduces below 15 in the following cycle. 107 (a) (b) (c) (d) Figure 4.18: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15a during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.18a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.18c This pattern that the single realization has is similarly reected in the averaged ones of 1; 000 real- izations while both scale and size is reduced. Precisely, the averaged trajectory also traverse tightly near the closed loop,OPQO, for the rst few rounds but each leg of the trajectory in each round starts being apart from the deterministic one, moving towards to the R 1 corner (Fig.4.18c). This is similar to what the single realization shows in Figure 4.18a, however, deviation is much smaller on average. The smaller deviation is also captured in the averaged tumor size that it maximally reaches only up to 10 until the 6th evolutionary cycle while it increases to over 25 in the last round. On average, tumor size increases but the rate of increase is a lot more mild compared to Figure 4.18b. 4.3.1 Synergisticmultidrugschedule Under the same constant combination of two drugs, i.e.C 1 = 0:5;C 2 = 0:2 duringT OP ,C 1 = 0;C 2 = 0 during T PQ and C 1 = 0:2;C 2 = 0:5 during T QO unit times, we evaluate the adaptive chemotherapy applied to the Moran process of sizeN = 50K when two drug interact synergistically. Recall from Figure 4.15d that for synergistic drug interactions (e = 0:3),T OP = 14:260,T PQ = 6:901,T QO = 11:129, and thus the associated one evolutionary cycle,T , in this case is equal toT = 32:29 and the corresponding 108 evolutionary steps,, for the Moran process withN = 50K is then equal to = 1; 614; 500. Compared to the additive drug interaction case, this schedule has a longer cycle and hence more evolutionary steps. (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.19: The spread of the distribution of terminal points (green dots) around O in each round for 1; 000 realizations of the stochastic Moran process with size N = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15d during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 Again we consider 1; 000 realizations of the Moran process that start atO. Applying the two drug adaptive chemo schedule for 8 evolutionary cycles to these individual stochastic systems provides the spread of the distribution of points aroundO at the end of each cycle and those points are depicted as green dots in Figure 4.19. Similar patterns that the additive drug interaction cases possess are shown for synergistic drug interactions. Expected by both the stochasticity driving the terminal point of the system to settle down on a neighborhood and the continued application of the adaptive chemotherapy associated with the closed loop,OPQO, at that inexact point each cycle, the points aroundO get more widely spread as the number of rounds increases. Moreover, for a xed cycle, the distribution of the points vary less in the frequency of theR 2 subpopulations but more in eitherS orR 1 subpopulations. This is where the stochastic competitive release ofR 1 as a result of the continuous delivery of a higher dose of drug 1 is reected. 109 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.20: The spread of the distribution of points around O for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug synergistic (e = 0:3) adaptive chemo schedule as in Figure 4.15d during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0166, 2 = 0:0082; (b) round 2, 1 = 0:0279, 2 = 0:0140; (c) round 3, 1 = 0:0445, 2 = 0:0200; (d) round 4, 1 = 0:0743, 2 = 0:0262; (e) round 5, 1 = 0:1231, 2 = 0:0344; (f) round 6, 1 = 0:1857, 2 = 0:0420; (g) round 7, 1 = 0:2496, 2 = 0:0467; (h) round 8, 1 = 0:3057, 2 = 0:0529 Unlike the additive drug interactions, it is observed for a xed round that the synergistic drug interac- tions bring in the wider distribution than the additive drug interactions entail. This deviation is quantied in terms of the semi-major and semi-minor axis by plotting the points aroundO in the principal axis co- ordinate system. As it can be seen by comparing Figure 4.17 and Figure 4.20, both the semi-major axis, 1 , and the semi-minor axis, 2 , are dominating in all cycles when two drugs interact synergistically. How- ever, any signicantly new feature in the form of distribution is not shown. It starts loosing the multi Gaussian structure from the 4th round and the area where points are heavily distributed gets further away from the point,O, as the additive drug interactions result. Trajectories of one randomly chosen realization and the averaged of 1; 000 samples during 8 cycles are given in Figure 4.21. One sample in Figure 4.21a gets o much the expected trajectory quite early from the second round, and the associated tumor size explosively increases, reaching over 300 even in the 4th 110 (a) (b) (c) (d) Figure 4.21: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15d during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.21a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.21c round (Fig.4.21b). When the observation is made on the sample average, then the mean trajectory becomes mild in the scope of movement from the closed loop,OPQO, maintaining its move around the loop closely until the 4th round and expanding its journey towards theR 1 corner from the next round, which is the observed feature for the additive drug interactions. Although the additive and synergistic drug interactions do not dier much in their trajectories and the distribution of the terminal points of each round except for the deviation, things change when it comes to tumor size. Recall that tumor size is modeled by the tumor-growth equation in (4.8) where the constant background tness,g, enters. For the numerical purpose, we set dierent values ofg for each drug inter- action as in Table 4.2 and this allows the tumor size associated with the adjusted replicator equation to be periodic with the period,T , where the one evolutionary cycle,T , also varies depending on the type of drug interactions. For both the additive and synergistic drug interactions, tumor size is well controlled with a slight increase in the maximum size for the rst 3 to 4 rounds and the maximal tumor size keep increasing as the number of rounds rises. Compared to its own maximal tumor size associated with the deterministic system, the additive drug interaction causes approximately 8:2 times increase while the maximal tumor size is increased by approximately 37:23 times for the synergistic drug interaction. What it infers is that 111 synergistic drug interactions allow the use of two drug chemotherapy schedule for a shorter time than additive drug interactions permit. 4.3.2 Antagonisticmultidrugschedule Evaluation of the two drug adaptive chemotherapy schedule composed of the same constant combination of two drugs, i.e. C 1 = 0:5;C 2 = 0:2 duringT OP ,C 1 = 0;C 2 = 0 duringT PQ andC 1 = 0:2;C 2 = 0:5 duringT QO unit times, is to be made with 1; 000 realizations of the Moran process of sizeN = 50K for antagonistic drug interactions. With the choice of the drug interaction parameter,e =0:3, and those constant combination of drugs, the time period for each leg is given in Figure 4.15g asT OP = 14:300, T PQ = 1:424 andT QO = 11:150, which in turn results in the total time period,T , being equal toT = 26:8740. With the selected sizeN = 50K of the Moran process, it gives rise to the evolutionary steps, , being equal to = 1; 343; 700, and it is the shortest number of steps among all drug interactions considered. (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.22: The spread of the distribution of terminal points (green dots) around O in each round for 1; 000 realizations of the stochastic Moran process with sizeN = 50K becomes wider as the number of rounds increases when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15g during 8 rounds since its exact start atO. (a) round 1; (b) round 2; (c) round 3; (d) round 4; (e) round 5; (f) round 6; (g) round 7; (h) round 8 112 Indicating the last points of each run for 1; 000 individual simulations of the Moran process as green dots visualizes the spread of the distribution of points aroundO in Figure 4.22. The distribution obviously gets more widely disseminated in the rise of the number of rounds as similarly shown for other drug interaction cases. It is also apparently seen that the antagonistic drug interaction preserves thick mass aroundO for a longer cycle in comparison with the additive drug interaction. (a) (b) (c) (d) (e) (f) (g) (h) Figure 4.23: The spread of the distribution of points around O for 1; 000 realizations of the stochastic Moran process of sizeN = 50K shows the increase of both the semi-major axis, 1 , and the semi-minor axis, 2 , as the number of rounds increases from 1 to 8 in the principal axis coordinate system when each realization starts atO and evolves under the multi drug antagonistic (e =0:3) adaptive chemo schedule as in Figure 4.15g during 8 evolutionary cycles. Though forming a multivariate Gaussian distribution nearly aroundO at the beginning few rounds, the spread gets further away from the initial point,O, as the adaptive schedule is repeated. (a) round 1, 1 = 0:0146, 2 = 0:0073; (b) round 2, 1 = 0:0250, 2 = 0:0118; (c) round 3, 1 = 0:0374, 2 = 0:0158; (d) round 4, 1 = 0:0584, 2 = 0:0206; (e) round 5, 1 = 0:0906, 2 = 0:0267; (f) round 6, 1 = 0:1369, 2 = 0:0329; (g) round 7, 1 = 0:1942, 2 = 0:0369; (h) round 8, 1 = 0:2516, 2 = 0:0409 When the points are plotted in the principal axis coordinate system in Figure 4.23, we observe that the distribution in the form of multi Gaussian structure is maintained up to the third cycle similarly to two other drug interaction cases, but the distance between the area where points are heavily crowded and the exact point,O, slowly increases in this case. In fact, the point,O, in the principal axis coordinate system is located at the border line in the last round for both the additive and synergistic cases when the plots are smoothened using the kernel density estimation method while Figure 4.23h shows the heavy mass is 113 still close to the point,O, even in the 8th evolutionary cycle. In addition, the deviation of the distribution around O for the antagonistic drug interaction turns out to be smallest in all cycles. It must be due to that administering a similar amount of the total dose, the adaptive chemotherapy schedule results in the smallest closed loop for the expected system, the adjusted replicator dynamic, in its one evolutionary cycle when drugs interact antagonistically as shown in Figure 4.15. Thus the stochastic Moran process has less chance to reect its randomness on its path as it evolves. The exact values of the semi-major axis, 1 , and the semi-minor axis, 2 , are given on top of each plot in Figure 4.23. (a) (b) (c) (d) Figure 4.24: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, as in Figure 4.15g during 8 rounds since its exact start atO. (a) one realization; (b) tumor size corresponding to Figure 4.24a; (c) the averaged trajectory of 1; 000 realizations; (d) the averaged tumor size corresponding to Figure 4.24c The trajectories of one realization chosen at random and the mean for the 1; 000 realizations of the Moran process are provided in Figure 4.24. One realization of the stochastic system in Figure 4.24a evolves as much as there are moreR 1 subpopulations at the end of 8 cycles, showing the early development of the group of the R 1 subpopulations as a result of the competitive release caused by the higher dose of drug 2. However, the stochastic Moran process can only be understood in terms of a realistic estimate, and we take it to be the sample mean since it is of a nite but large population size. The trajectory in Figure 4.24c corresponding to the sample mean moves tightly to the deterministic path until the 5th round and gradually starts being o towards to theR 1 corner, but not as much as the one realization shows. Except that the averaged trajectory under the antagonistic drug interaction generates a smaller open-triangle-like 114 shape in each cycle than under two other drug interaction cases, it shares the same pattern with other two cases. However, the considerable advantage of the antagonistic drug interactions is proven when we focus on tumor size. With the choice of the constant background tness,g, as in Table 4.2, tumor size is not only controlled relatively longer with little increase in its maximum until the 4th round but also the increase that starts from the 5th round is much slower compared to either the additive or synergistic drug interactions. Precisely, the maximal tumor size obtained at the end of the second 1=3 cycle in each round is increased by approximately 3:56 times in 8 evolutionary cycles, and this is a considerably smaller rate if we note that the rates were 8:2 and 37:23 for the additive and synergistic drug interactions, respectively. The fact that the maximum tumor size even in the last round is less than only 6 is also notable. Set the threshold of tumor size to be equal to, say, 10, then the two drug adaptive chemotherapy schedule can be applied during 5 cycles for the additive drug interaction as seen in Figure 4.18d but it can be only applied during one cycle for the synergistic drug interaction as veried in Figure 4.21d. Figure 4.24d nally justies that the most powerful drug interaction in controlling tumor size when applying the adaptive chemo schedule is the antagonistic one since tumor size most slowly grows and it grows much less than it is allowed. As mentioned before, one of reasons that the antagonistic drug interaction when applying the adaptive chemotherapy schedule is that the spread of the distribution of points aroundO keeps centered nearO with the smallest deviation as the number of rounds increases, making the adoption of the adaptive chemo schedule for the next round more feasible. The measures, the semi-major and semi-minor axis, to explain the deviation in the distribution in each round are depicted in Figure 4.25 to summarize that the largest deviation is obtained with the synergistic interaction while the smallest one with the antagonistic drug interaction. We also present the tted curves of both 1 and 2 in log-log scale in Figure 4.25d, 4.25e, 4.25f 115 (a) (b) (c) (d) (e) (f) Figure 4.25: The semi-major (and -minor) axis of the distribution of the points around P at the end of each evolutionary cycle for 1; 000 realizations of the Moran process withN = 10K orN = 50K overall increases in the number of rounds, showing the power-law dependency, where each realization evolves under the administration of the adaptive schedule associated with the adjusted replicator system during 8 evolutionary cycles since its exact start atO. (a)e = 0; (b)e = 0:3; (c)e =0:3; (d)e = 0, the log-log t; (e)e = 0:3, the log-log t; (f)e =0:3, the log-log t for the additive, synergistic and antagonistic interaction, respectively, and all of them show the power-law dependency in the number,n, of rounds as follows withN = 10K: e=0 1 0:0295n 1:2504 ; e=0 2 0:0165n 0:9114 ; e=0:3 1 0:0337n 1:2108 ; e=0:3 2 0:0175n 0:9904 ; e=0:3 1 0:0261n 1:2214 ; e=0:3 2 0:0139n 0:9734 : (4.16) 116 Similarly, withN = 50K: e=0 1 0:0115n 1:4387 ; e=0 2 0:0069n 0:9152 ; e=0:3 1 0:0120n 1:4752 ; e=0:3 2 0:0077n 0:9216 ; e=0:3 1 0:0107n 1:3981 ; e=0:3 2 0:0068n 0:8525 : (4.17) 4.4 Comparison of adaptive multi drug chemotherapy schedule with standardclinicalapproaches The evaluation of the adaptive two drug chemotherapy schedule for a few cycles with the xed amount of the total dose was made in the previous section across the dierent types of drug interactions when it is applied to the Moran process that starts at a xed point. It was concluded that the therapy associated with the closed loop,OPQO, turns out to be the most ecient when two drugs interact antagonistically in the sense that it gives rise to the smallest trajectory and the lowest increase in tumor size during the whole cycles. On the other hand, when two drug interact synergistically, it involves the largest trajectory with the greatest deviation in the spread of the distribution of point around the initial point,O. On top of that, it provides the signicantly rapid increase in tumor size compared to two other drug interactions. In this section, xing a drug interaction, we evaluate the adaptive two drug chemotherapy schedule by comparing it to the most frequently used standard clinical schedules: the maximum tolerated dose schedule (MTD) and the low-dose metronomic schedule (LDM). We already considered these two clinical schedules in Section 3.4 with a single drug chemotherapy model, where we had only one controller of the system which is the concentration function of the single drug with a constraint 0C(t) 1 for allt. In 117 two drug chemotherapy model, we have two controllers,C 1 (t) andC 2 (t), that are concentration functions of drug 1 and drug 2, and these functions have come with one more constraint C 1 (t) +C 2 (t) 1 (4.18) for allt as well as 0C i (t) 1. The comparison of two drug chemotherapy schedules is made by xing both a type of drug interaction, more precisely a value of the drug interaction parameter,e, and the same total dose delivered over a nite time, where the total dose, D(T ), of two drugs delivered duringT unit time is similarly dened to the equation in (3.13) by: D(T ) := Z T 0 C 1 (t) +C 2 (t)dt: (4.19) (a) (b) (c) (d) Figure 4.26: The multiple additive (e = 0) standard clinical approaches are designed to have the same total dose, being equal to 17:7604, as the amount that is delivered during one round (29:4440 unit time) according to the adaptive chemo schedule in Figure 4.15a. (a) adaptive; (b)MTD 1 : the drug 1 is maximally administered during the beginningT MTD:C 1 = 9:3422 unit time while the drug 2 is delivered during the lastT MTD:C 2 = 8:4182 unit time. (c)MTD 2 : the drug 2 is maximally administered during the beginning T MTD:C 2 = 8:4182 unit time while the drug 1 is delivered at largest during the lastT MTD:C 1 = 9:3422 unit time. (d)LDM: both the drug 1 and drug 2 are constantly administered during the whole rounds at the level ofC 1 = 0:347287 andC 2 = 0:255905, respectively. Our two drug adaptive chemotherapy schedule in Figure 4.15a associated with the closed loop,OPQO, in Figure 4.5b for the additive drug interaction (e = 0) has the total doseD = 17:7604 for one evolutionary cycle which we setT to be. Having in mind that MTD delivers the highest dose of a drug by which any signicantly unacceptable side eects are not caused, we may naturally come up with two MTD schedules 118 with the same total dose, D, duringT = 29:4440 unit time forcing them to satisfy the equation (4.18). One MTD, which we denote byMTD 1 , is to deliver the highest dose of drug 1 rst until the total dose of drug 1 underMTD 1 meets half the total dose of drug 1 under the adaptive schedule, that isD=2. For the restD=2, the highest dose of drug 2 is delivered at the end of one cycle during the time over which the total dose of drug 2 underMTD 1 is equal toD=2. LetT MTD:C1 andT MTD:C2 be the time period over which drug 1 and drug 2, respectively, is maximally delivered so that the total dose of both drugs during one evolutionary cycle is equal toD. ThenT MTD:C1 = 9:3422,T MTD:C2 = 8:4182 for the additive drug interaction, andMTD 1 is well described in Figure 4.26b. Another MTD schedule is exactly the same as MTD 1 , except that the delivery order of drugs is switched so that drug 2 is delivered rst duringT MTD:C2 unit times and drug 1 is administered at the end of the cycle duringT MTD:C1 unit times, having a pause between. We denote this schedule byMTD 2 and it is depicted in Figure 4.26c. For both MTD schedules, we simply put the delivery times of two drugs apart, at the beginning and at the end, so that the area of the maximal concentration is not overlapped dur- ing one cycle in order to make sure the constraint in (4.18). However, it is true that we may design innitely many MTD schedules by adjusting time window as long as two maximal windows are not overlapped. LDM is one other typical clinical approach where the low dose is continuously delivered during a nite time period. It is straight forward to design LDM by simply averaging the total dose of each drug under the adaptive schedule for the entire time period. From the adaptive schedule in Figure 4.15a for the additive drug interaction, the constant low dose of drug 1 turns out to be 0:347287 and the one of drug 2 is lower and equal to 0:255905. This LDM schedule for the additive interaction is given in Figure 4.26d. Thus we so far have 4 two drug chemotherapy schedules that have the same total dose for a xed one evolutionary 119 cycle in hand to compare when two drugs are additive. Denoting the type of schedule as superscript for the concentration functions,C i (t), of drugi, all the schedules in Figure 4.26 satises D = Z T 0 C Adaptive 1 (t) +C Adaptive 2 (t)dt = Z T 0 C MTD 1 1 (t) +C MTD 1 2 (t)dt = Z T 0 C MTD 2 1 (t) +C MTD 2 2 (t)dt = Z T 0 C LDM 1 (t) +C LDM 2 (t)dt (4.20) , whereT =T OP +T PQ +T QO = 29:4440 andD = 17:7604. Note that we omit indicating here the type of drug interaction since the comparison to be made is within a xed type and the construction of MTD’s and LDM is designed in the same manner to satisfy the equation in (4.20) with simply dierent values of D andT for each drug interaction as determined in Table 4.2. Recall from Figure 4.18c,4.18d that applying the adaptive chemotherapy schedule to the 1; 000 realiza- tions of the Moran process of the population size, N = 50K, successfully controls both trajectory and tumor size for the rst few rounds on average when two drugs are additive as the well designed schedule eventually plays the role of delaying the competitive release ofR 1 orR 2 . Although the stochasticity drives the system gradually to be o the expected route in the direction towards theR 1 corner resulting in the growth of tumor in size, its increase is below 30 during 8 evolutionary cycles. These trajectory and tumor size for the adaptive schedule are depicted in pink in Figure 4.27. On the other hand, when MTD schedules are applied, the rst maximal dose of a drug helps the system shortly build up the resistant subpopulations even though the duration of the maximal dose is relatively shorter than the rst 1=3 piece of the adaptive schedule. ForMTD 1 which administer the drug 1 rst, the resistantR 2 subpopulations to drug 1 starts quickly proliferating in a population and it nearly overtakes the entire population at the end of the delivery of the maximal dose even in the rst round. However, a 120 (a) (b) (c) (d) (e) (f) Figure 4.27: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug additive (e = 0) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.27a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.27a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.27d; (f) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.27d pause, that is no drug, between two maximal dose administration of dierent drugs allows the system to step back towards theS corner for a while. When the administration of the maximal dose of drug 2 begins, the resistantR 1 subpopulations to drug 2 starts growing and the population nally has a higher frequency ofR 1 than others at the end of the rst round. This pattern continues as the number of rounds increases, but the absent mutation kills the whole sensitive subpopulations early and the population becomes full of resistant populations. In fact, the maximal dose administration of a drug let the selection functions of sensitive cells to the drug used vanish while it lets resistant cells to the drug used have the strongest selection regardless of the value of e according to the equation (4.1) and the payo matrix, A, in (4.6). Then the expected tness of those sensitive cells become less tter, being equal to 1. This is reected to 121 our stochastic simulation under MTD.MTD 1 administers drug 2 in the last and it leads to the proliferation of the resistant subpopulationsR 1 to that drug 2 as can be seen from the terminal point of the last cycle marked as a blue dot in Figure 4.27a, 4.27d. The whole trajectory corresponding toMTD 1 is described in a blue line in Figure 4.27. MTD 2 acts exactly in an opposite way since it interchanges the order of administration of two drugs. The maximal dose of drug 2 at the beginning leads the dynamic of the stochastic system to nearly theR 1 corner since the highest dose of drug 2 rapidly help the resistant subpopulationsR 1 to drug 2 develop. We noticed that this nearly xation toR 2 is reached even in a shorter timeT MTD:C2 = 8:4182 thanMTD 1 results in the almost full growth ofR 2 subpopulations forT MTD:C1 = 9:3422 unit times. It is hard to tell that one drug is stronger or less eective than another, but it is related to the distance from the starting point,O, in fact the initial distribution,O, is closer to theR 1 corner in a phase space. ThenMTD 2 takes a pause for a while and it causes the sensitive subpopulations to recover some, and when the highest dose of drug 1 begin its delivery, the development of the resistant subpopulations ofR 2 to drug 1 is achieved though it is smaller in extent as there are more sensitive subpopulations. Again, this pattern continues for the following rounds, leading the stochastic system to eventually have more and moreR 1 subpopulations on its way. The trajectories or tumor size, of both one single simulation and the averaged one, corre- sponding toMTD 2 are shown in green in Figure 4.27. The dynamics determined by both MTD schedules are dierent in directions especially in the growth and decay of resistant subpopulations, however, they are almost the same in their induced tumor size as can be seen in Figure 4.27e. The rapid development of the resistant subpopulations and the dramatic decrease in the number of sensitive cells brings out the explosive increase in tumor size even before the rst pause ends though tumor size is initially decreased sensitively reacting to the drug that is rst delivered due to the high dose. The low-dose metronomic schedule (LDM) shows somewhat dierent dynamic from both the adaptive and the maximum tolerated dose schedules. This constant schedule guides the stochastic system to the 122 (a) (b) (c) (d) (e) (f) Figure 4.28: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug synergistic (e = 0:3) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.28a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.28a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.28d; (f) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.28d xation toS on average and the associated trajectory is given in Figure4.27d as a red line. For one real- ization chosen at random also shows the direct approach to theS corner. In fact,S is an asymptotically stable xed point and trajectories converge to theS corner regardless of initial distributions as long as both drugs are delivered at a low dose. The exact description of theESS is explained in [30] for additive drug interactions. In detail,C 1 = 7=18 is a bifurcation value where the xed pointR 2 changes its stabil- ity from unstable to stable, andS changes its stability from stable to unstable atC 1 = 1=2. For drug 2, C 2 = 1=3 is a bifurcation value andR 1 changes from an unstable to a stable xed point. AtC 2 = 1=2,S changes its stability from stable to unstable. InLDM schedule,C 1 = 0:347261 andC 2 = 0:255905 , and the initial point,O, is laid in the basin of attraction forS for this combination. As a result, the averaged 123 trajectory of 1; 000 realizations of the Moran process of a large population size shows the convergence to the evolutionarily stable strategyS under LDM. Though the populations is xated to one of cancerous state,S, the associated tumor size surprisingly reduces during 8 cycles as shown in red in Figure 4.27e. It seems to conict with the exponentially in- creasing tumor size associated with the adjusted replicator system when no drug is used shown in Figure 4.15c. However, recall that tumor size is modeled using the tumor-growth equation in (4.8) by which it grows whenever the averaged tness, < f >, ofS, R 1 andR 2 in the entire population is greater than the constant background tness,g, that is set to be equal tog = 1:4527, and that< f > is a function of selectionsw S ,w R 1 andw R 2 as in (4.1) whereC i ’s enters a population through. Thus, a dierent set ofC i ’s gives rise to a dierent value of<f > even at a xed state and this possibly brings out a totally opposite destiny in tumor size. We computed the averaged tness,<f >, for the averaged Moran process during 8 evolutionary cycles and it shows that<f > is less tter than the constant background tness,g, during the entire time period as seen in Figure 4.27f and this allows the tumor size under LDM to decrease. We have investigated how successfully 4 dierent two drug chemotherapy schedules work when they are applied to the Moran process that starts a xed point,O, for the additive drug interaction. We conclude that bothMTD’s work the worst showing the immediate nearly full development of resistant subpopula- tions and the extraordinarily rapid growth rate of tumor size too early in the whole 8 cycles. The adaptive schedule associated with the deterministic adjusted replicator system is successful in controlling both the trajectory and tumor size by preventing the competitive release of a resistant subpopulation for the rst few cycles. Depending on how much big threshold of tumor size we set, this schedule can be applied longer or shorter. However, it does not work as eectively as LDM does which completely kills out tumor. What LDM does is that the continuous low dose administration helps lower the averaged tness of can- cerous cells,S,R 1 andR 2 , making the tness other factors such as healthy cells or any environment that contribute to tumor development tter. 124 (a) (b) (c) (d) (e) (f) Figure 4.29: 1; 000 realizations of the stochastic Moran process with sizeN = 50K show that the satu- ration of cancer cells can be delayed longer on average compared to being under either LDM and MTD’s when each realization evolves under the administration of the multi drug antagonistic (e =0:3) adaptive chemo schedule, associated with the adjusted replicator system, during 8 rounds since its exact start atO. (a) one single realization; (b) tumor size corresponding to Figure 4.29a; (c) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.29a; (d) the averaged trajectory of all realizations; (e) the averaged tumor size corresponding to Figure 4.29d; (f) the averaged tness ofS,R 1 andR 2 cells corresponding to Figure 4.29d Two MTD’s and LDM schedules can be easily constructed for synergistic and antagonistic drug inter- actions referring to its own adaptive schedule associated with eachOPQO and total dose,D, during one evolutionary cycle,T . Then applying those 4 schedules to the same number of realizations of the Moran process with size,N = 50K is processed in the same manner for eache = 0:3;0:3. We present a set of trajectory, tumor size, and the averaged tness,<f >, for both one randomly chosen realization and the averaged of 1; 000 individuals in Figure 4.28 for the synergistic interaction and in Figure 4.29 for the an- tagonistic drug interaction. Regardless of the drug interaction parameter,e, the maximum tolerated dose schedules turn out to be the worst in ecacy of controlling tumor size and the adaptive chemo schedule is successful for a few cycles while it may be applied longer if two drug interact antagonistically but shorter 125 if synergistically. LDM wins all other three chemo schedules delivering the same total dose for 8 cycles since the low dose of both drugs at that starting point,O, helps either noncancerous cells or environment tter and hence it gradually reduces tumor size. However, remember that one point in a phase space corresponds to a unique status of a patient. Thus the change in the initial status of a patient leads to a dierent starting point instead ofO, then to a dierent closed loop, a new associated adaptive chemo schedule, and a dierent MTD or LDM. The result in this thesis proves LDM as the most powerful chemotherapeutic method at one point, and the conclusion may vary depending on the given initial stage. 126 Chapter5 Futuredirections Originating from this thesis, generalization of our model can be made in many dierent ways. First, for both single and two drug chemotherapy models, we assumed for simplicity that no further mutation is allowed as the system evolves since the population initially attains the minimally required ingredients (e.g. resistant cells) by preexisting mutation. However, mutation is one of key factors in the Darwinian evolutionary theory along with heredity and natural selection, and there are a lot of literatures that describe an evolving population allowing mutation during the whole evolution process. Allowing mutation during evolution process also gives rise to qualitatively dierent but desired results in asymptotic behavior such as the emergence of cooperation as it was shown for both a nite and innite population by in [7, 24, 26]. Thus, it can be our one possible generalization of our model to introduce a mutation rate to each group of cells and let them evolve according to either the replicator-mutator equations for an innite populations or the Moran process with mutation for a nite population, and examine if the introduction of mutation helps us have more promising results in terms of the maximum number of cycles during which the application of the adaptive schedule to the Moran process is successful in regulating tumor volume. Second, for two drug chemotherapy model, tumor volume is designed using the tumor-growth equa- tion in (4.8) where the background tness, that captures the information of healthy cells or surrounding environment, is set to be constant. However, when drug concentration enters the total population that 127 contributes tumor development, not only it changes the tness of cancer cells but also whatever ingredi- ents that inuence on forming the background tness also adapt, making the background tness uctuate rather than stay constant. Representing the background tness as a function ofC 1 (t) andC 2 (t) is more re- alistic and generalizes our model. This should be designed along with clinical and laboratory observations about reproduction rate of healthy cells and surrounding environment. The last but not the least important interest for the next study stems from the evaluation of the adap- tive two drug chemotherapy model where the low-dose metronomic schedules turns out to be the best chemotherapeutic strategy for the Moran process that starts at one exact point,O, regardless of drug in- teractions in the sense that it was the most ecient in reducing the tumor volume. However, in reality, the dierent patients’ initial status are assigned to all dierent points fromO and from each other in a phase space, and the conclusion of the evaluation of 4 discussed chemotherapy schedules might vary depending on initial states. For that reason, we would like to visualize the area of initial points in a phase space, at which each of those schedules work the best. 128 Bibliography [1] P. M. Altrock and A. Traulsen. “Deterministic evolutionary game dynamics in nite populations”. In: Physical Review E 80.1 (2009), p. 011909. [2] P. M. Altrock and A. Traulsen. “Fixation times in evolutionary games under weak selection”. In: New Journal of Physics 11 (2009), p. 013012. [3] R. Axelrod. The evolution of cooperation. New York, USA: Basic Books, 1984. [4] S. Bernstein. “Démonstration du théorème de Weierstrass fondée sur le calcul des probabilités (Proof of the theorem of Weierstrass based on the calculus of probabilities)”. In: Communications of the Kharkov Mathematical Society 13.1 (1912), pp. 1–2. [5] C. I. Bliss. “The toxicity of poisons applied jointly”. In: Annals of Applied Biology 26.3 (1939), pp. 585–615. [6] I. M. Bomze. “Lotka-Volterra equation and replicator dynamics: A two-dimensional classication”. In: Biological Cybernetics 48 (1983), pp. 201–211. [7] I. M. Bomze and R. Burger. “Stability by mutation in evolutionary games”. In: Games and Economic Behavior 11.2 (1995), pp. 146–172. [8] “Changes associated with the development of resistance to imatinib (STI571) in two leukemia cell lines expressing p210 Bcr/Abl protein”. In: Cancer 100.7 (2004), pp. 1459–1471. [9] T. C. Chou and D. C. Rideout. Synergism and antagonism in chemotherapy. San Diego, CA, USA: Academic Press, 1991. [10] J. C. Claussen and A. Traulsen. “Non-Gaussian uctuations arising from nite populations: Exact results for the evolutionary Moran process”. In: Physical Review E 71.2 (2005), p. 025101. [11] J. H. Connell. “The inuence of interspecic competition and other factors on the distribution of the barnacle Chthamalus stellatus”. In: Ecology 42.4 (1961), pp. 710–723. [12] C. R. Darwin. On the origin of species. London, UK: John Murray, 1859. 129 [13] J. D. Davidson and et al. “An increase in the expression of ribonucleotide reductase large subunit 1 is associated with gemcitabine resistance in non-small cell lung cancer cell lines”. In: Cancer Research 64.11 (2004), pp. 3761–3766. [14] R. Dawkins. The selsh gene. Oxford, UK: Oxford University Press, 1976. [15] E. M. Ferreira and A. G. M. Neves. “Fixation probabilities for the Moran process with three or more strategies: general and coupling results”. In:JournalofMathematicalBiology 81 (2020), pp. 277–314. [16] R. A. Fisher. The genetical theory of natural selection. Oxford, UK: Clarendon Press; Oxford University Press, 1930. [17] D. Foster and P. Young. “Stochastic evolutionary game dynamics”. In: Theoretical Population Biology 38.2 (1990), pp. 219–232. [18] R. A. Gatenby. “A change of strategy in the war on cancer”. In: Nature 459 (2009), pp. 508–509. [19] R. A. Gatenby, A. S. Silva, R. J. Gillies, and B. R. Frieden. “Adaptive therapy”. In: Cancer Research 69.11 (2009), pp. 4894–4903. [20] A. Gautam and G. Bepler. “Suppression of lung tumor formation by the regulatory subunit of ribonucleotide reductase”. In: Cancer Research 66.13 (2006), pp. 6497–6502. [21] C. Hauert, S. D. Monte, J. Hofbauer, and K. Sigmund. “Replicator dynamics for optional public good games”. In: Journal of Theoretical Biology 218.2 (2002), pp. 187–194. [22] J. Hofbauer and K. Sigmund. “Evolutionary game dynamics”. In: Bulletin of the American Mathematical Society 40 (2003), pp. 479–519. [23] J. Hofbauer and K. Sigmund. Evolutionary games and population dynamics. Cambridge, UK: Cambridge University Press, 1998. [24] L. A. Imhof, D. Fudenberg, and M. A. Nowak. “Evolutionary cycles of cooperation and defection”. In: Proceedings of the National Academy of Sciences of the United States of America 102.31 (2005), pp. 10797–10800. [25] L. A. Imhof and M. A. Nowak. “Evolutionary game dynamics in a Wright-Fisher process”. In: Journal of mathematical biology 52 (2006), pp. 667–681. [26] L. A. Imhof and M. A. Nowak. “Stochastic evolutionary dynamics of direct reciprocity”. In: Proceedings of the Royal Society B: Biological Sciences 277.1680 (2010), pp. 463–468. [27] W. S. Kendal. “Gompertzian Growth as a consequence of tumor heterogeneity”. In: Mathematical Biosciences 73.1 (1985), pp. 103–107. [28] S. Lessard and V. Ladret. “The probability of xation of a single mutant in an exchangeable selection model”. In: Journal of Mathematical Biology 54 (2007), pp. 721–744. 130 [29] E. Lieberman, C. Hauert, and M. A. Nowak. “Evolutionary dynamics on graphs”. In: Nature 35 (2005), pp. 312–316. [30] Y. Ma and P. K. Newton. “Role of synergy and antagonism in designing multidrug adaptive chemotherapy schedules”. In: Physical Review E 103.3 (2021), p. 032408. [31] P. A. P. Moran. “Random process in genetics”. In: Mathematical Proceedings of the Cambridge Philosophical Society 54.1 (1958), pp. 60–71. [32] P. A. P. Moran. The Statistical Processes of Evolutionary Theory. Oxford, UK: Clarendon Press; Oxford University Press, 1962. [33] O. Morgenstern and J. V. Neumann. Theory of games and economic behavior. Princeton, NJ, USA: Princeton University Press, 1953. [34] J. F. Nash. “Equilibrium points in n-person games”. In: Proceedings of the National Academy of Sciences of the United States of America 36.1 (1950), pp. 48–49. [35] J. V. Neumann. “Zur Theorie der Gesellschaftsspiele [On the Theory of Games of Strategy]”. In: Mathematische Annalen [Mathematical Annals] (in German) 100 (1928), pp. 295–320. [36] P. K. Newton and Y. Ma. “Nonlinear adaptive control of competitive release and chemotherapeutic resistance”. In: Physical Review E 99.2 (2019), p. 022404. [37] M. A. Nowak. Evolutionary dynamics: Exploring the equations of life. Cambridge, MA, USA: Harvard University Press, 2006. [38] M. A. Nowak and R. M. May. “Evolutionary games and spatial chaos”. In: Nature 359 (1992), pp. 826–829. [39] M. A. Nowak and R. M. May. “The spatial dilemmas of evolution”. In: International Journal of Bifurcation and Chaos 3.1 (1993), pp. 35–78. [40] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg. “Emergence of cooperation and evolutionary stability in nite populations”. In: Nature 428 (2004), pp. 646–650. [41] M. A. Nowak and K. Sigmund. “Evolutionary dynamics of biological games”. In: Science 303.5659 (2004), pp. 793–799. [42] M. A. Nowak and K. Sigmund. “The evolution of stochastic strategies in the Prisoner’s dilemma”. In: Acta Applicandae Mathematicae (in Netherlands) 20 (1990), pp. 247–265. [43] K. M. Page and M. A. Nowak. “Unifying evolutionary dynamics”. In: Journal of Theoretical Biology 219.1 (2002), pp. 93–98. [44] F. C. Santos, M. D. Santos, and J. M. Pacheco. “Social diversity promotes the emergence of cooperation in public goods games”. In: Nature 454 (2008), pp. 213–216. 131 [45] P. Schuster and K. Sigmund. “Replicator dynamics”. In: Journal of Theoretical Biology 100.3 (1983), pp. 533–538. [46] J. M. Smith. Evolution and the theory of games. Cambridge, UK: Cambridge University Press, 1982. [47] J. M. Smith. On evolution. Edinburgh, Scotland: Edinburgh University Press, 1972. [48] J. M. Smith. “The theory of games and the evolution of animal conicts”. In: Journal of Theoretical Biology 47.1 (1974), pp. 209–221. [49] J. M. Smith and G. R. Price. “The logic of animal conicts”. In: Nature 246 (1973), pp. 15–18. [50] G. Szabo and C. Toke. “Evolutionary prisoner’s dilemma game on a square lattice”. In: Physical Review E 58.1 (1998), p. 69. [51] C. Taylor, D. Fudenberg, A. Sasaki, and M. A. Nowak. “Evolutionary game dynamics in nite populations”. In: Bulletin of Mathematical Biology 66 (2004), pp. 1621–1644. [52] P. D. Taylor and L. B. Jonker. “Evolutionary stable strategies and game dynamics”. In: Mathematical Biosciences 40.1–2 (1978), pp. 145–156. [53] A. Traulsen, J. C. Claussen, and C. Hauert. “Coevolutionary dynamics in large, but nite populations”. In: Physical Review E 74.1 (2006), p. 011901. [54] A. Traulsen, J. C. Claussen, and C. Hauert. “Coevolutionary dynamics: from nite to innite populations”. In: Physical Review Letters 95.23 (2005), p. 238701. [55] A. Traulsen and M. A. Nowak. “Evolution of cooperation by multilevel selection”. In:Proceedingsof the National Academy of Sciences of the United States of America 103.29 (2006), pp. 10952–10955. [56] A. Traulsen, M. A. Nowak, and J. M. Pacheco. “Stochastic dynamics of invasion and xation”. In: Physical Review E 74.1 (2006), p. 011909. [57] A. Traulsen, J. M. Pacheco, and L. A. Imhof. “Stochasticity and evolutionary stability”. In: Physical Review E 74.2 (2006), p. 021905. [58] A. Traulsen, J. M. Pacheco, and M. A. Nowak. “Pairwise comparison and selection temperature in evolutionary game dynamics”. In: Journal of Theoretical Biology 246.3 (2007), pp. 522–529. [59] J. West, Y. Ma, and P. K. Newton. “Capitalizing on competition: An evolutionary model of competitive release in metastatic castration resistant prostate cancer treatment”. In: Journal of Theoretical Biology 455 (2018), pp. 249–260. [60] S. Wright. “Evolution in Mendelian populations”. In: Genetics 16.2 (1931), pp. 97–159. [61] B. Wu, P. M. Altrock, L. Wang, and A. Traulsen. “Universality of weak selection”. In: Physical Review E 82.4 (2010), p. 046106. 132

Abstract (if available)

Abstract This thesis is concerned with modeling tumor growth and chemotherapy response using stochastic evolutionary game theory models. In particular we develop stochastic N-cell models of a heterogeneous (multiple cell types) tumor using a Moran process model with frequency dependent fitness, which in the population size limit converges to the deterministic adjusted replicator dynamical system as its mean-field limit. Briefly describing this limit and some of the details of our model and background literature, we develop new results associated with the fixation probability of cell types that do not necessarily depend on the assumption of weak selection, but are valid across the full range of selection strengths. We also develop our main results on adaptive chemotherapy scheduling using our stochastic N-cell evolutionary game theory model, both single drug and multi-drug scheduling, with the goal of avoiding chemo-resistance via the mechanism of competitive release, a concept borrowed from the ecology literature. The methods we develop are superior in measurable ways to more standard chemotherapy schedules currently in use at cancer centers all over the world, such as maximum tolerated dose (MTD) schedules and low-dose metronomic schedules (LDM).

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Reinforcement learning based design of chemotherapy schedules for avoiding chemo-resistance

PDF

Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance

PDF

Integrative mathematical and computational approaches to investigate metastatic cancers

PDF

Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response

PDF

Parameter estimation problems for stochastic partial differential equations from fluid dynamics

PDF

Adaptive agents on evolving network: an evolutionary game theory approach

PDF

Well posedness and asymptotic analysis for the stochastic equations of geophysical fluid dynamics

PDF

Unique continuation for parabolic and elliptic equations and analyticity results for Euler and Navier Stokes equations

PDF

Certain regularity problems in fluid dynamics

PDF

Stability analysis of nonlinear fluid models around affine motions

PDF

Development of a colorectal cancer-on-chip to investigate the tumor microenvironment's role in cancer progression

PDF

Asymptotic expansion for solutions of the Navier-Stokes equations with potential forces

PDF

Mach limits and free boundary problems in fluid dynamics

PDF

Regularity problems for the Boussinesq equations

PDF

Provable reinforcement learning for constrained and multi-agent control systems

PDF

Global existence, regularity, and asymptotic behavior for nonlinear partial differential equations

PDF

Learning and control in decentralized stochastic systems

PDF

Stochastic Variational Inference as a solution to the intractability problem of Bayesian inference

PDF

Model-based approaches to objective inference during steady-state and adaptive locomotor control

PDF

Modeling and simulation of circulating tumor cells in flow. Part I: Low-dimensional deformation models for circulating tumor cells in flow; Part II: Procoagulant circulating tumor cells in flow

Asset Metadata

Creator Park, Jiyeon (author)

Core Title Stochastic multidrug adaptive chemotherapy control of competitive release in tumors

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Applied Mathematics

Degree Conferral Date 2022-05

Publication Date 01/14/2022

Defense Date 12/08/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cancer evolution,evolutionary game theory,Moran process,OAI-PMH Harvest,optimal control theory,Prisoner's dilemma game,replicator dynamics

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kukavica, Igor (committee chair), Newton, Paul K. (committee member), Ziane, Mohammed (committee member)

Creator Email jiyeon.park66@gmail.com,parkjiye@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC110520281

Unique identifier UC110520281

Legacy Identifier etd-ParkJiyeon-10346

Document Type Dissertation

Format application/pdf (imt)

Rights Park, Jiyeon

Type texts

Source 20220124-usctheses-batch-908 (batch), 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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

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

Repository Email cisadmin@lib.usc.edu