Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response
(USC Thesis Other)
Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
UNIVERSITY OF SOUTHERN CALIFORNIA Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response by Jerey West A dissertation presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY (Mechanical Engineering) December 2017 Dedication Soli Deo gloria i Acknowledgements Nobel laureate Hermann Muller was right in saying, \If you prefer an `academic life' as a retreat from reality, do not go into biology. This eld is for a man or woman who wishes to get even closer to life." With this in mind, I rst thank my advisor Paul Newton, who beck- oned me to join the crusade against cancer using math. It has been immensely rewarding to work with you. I thank my committee members Aiichiro Nakano and Niema Pahlevan. I also wish to thank helpful collaborators and mentors such as (but not limited to), Paul Macklin, the Kuhn lab (Peter Kuhn, Jorge Nieva), Stacey Finley, the CAMM group (David Agus, Shannon Mumenthaler, Dan Ruderman). I thank my fellow comrades Jeremy Ma- son, Zaki Hasnain, and Yongian Ma for their patient feedback, helpful suggestions, and kind friendship. ii Contents Dedication i Acknowledgements ii List of Figures vi Abstract viii 1 Introduction 1 1.1 Cancer as an evolutionary and ecological process . . . . . . . . . . . . . . . 3 1.2 Chapter Summaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Emergent Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 The prisoner's dilemma as a cancer model 9 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 The prisoner's dilemma evolutionary game . . . . . . . . . . . . . . . . . . . 13 2.4 A tumor growth model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4.1 Mutations and heritability . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4.2 The tness landscape . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.3 Heterogeneity drives growth . . . . . . . . . . . . . . . . . . . . . . . 25 2.5 Simulated drug dosing strategies and therapeutic response . . . . . . . . . . 28 2.6 Mathematical modeling and tumor analytics . . . . . . . . . . . . . . . . . . 32 3 An evolutionary model of tumor cell kinetics and the emergence of molec- ular heterogeneity driving Gompertzian growth 36 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.3 Description of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.1 The Moran birth-death process . . . . . . . . . . . . . . . . . . . . . 42 3.3.2 The prisoner's dilemma payo matrix . . . . . . . . . . . . . . . . . 43 3.3.3 The tness landscape . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.4 Passenger and driver mutations . . . . . . . . . . . . . . . . . . . . . 49 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 iii Contents iv 3.4.1 Gompertzian tumor growth and three growth regimes . . . . . . . . 52 3.4.2 Heterogeneity and growth via statistical mechanics . . . . . . . . . . 56 3.4.3 Quantitative measures of tumor heterogeneity and growth . . . . . . 58 3.4.4 Dynamic phylogenetic trees and evolution of tness . . . . . . . . . 62 3.4.5 A comparison of early vs. late therapy . . . . . . . . . . . . . . . . . 65 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4 Chemotherapeutic dose scheduling based on tumor growth rates: the case for low dose metronomic high entropy therapies 70 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.1 Administration of metronomic chemotherapy . . . . . . . . . . . . . 76 4.2.2 The classic tumor regression laws . . . . . . . . . . . . . . . . . . . . 78 4.2.2.1 Skipper's Laws . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.2.2.2 Norton-Simon Hypothesis . . . . . . . . . . . . . . . . . . . 79 4.2.3 The Implications of the Norton-Simon Hypothesis . . . . . . . . . . 81 4.3 Materials and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.1 Chemotherapeutic agents alter the tness landscape . . . . . . . . . 83 4.3.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3.3 Dose concentration versus dose density . . . . . . . . . . . . . . . . . 89 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.4.1 Quantifying chemotherapeutic strategies via entropy metric . . . . . 91 4.4.2 LDM versus MTD chemotherapies . . . . . . . . . . . . . . . . . . . 94 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5 Harnessing the evolutionary cost of chemotherapeutic resistance by shap- ing the tness landscape of a tumor 100 5.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.2.1 Pre-existing resistance . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.2.2 Using evolutionary principles to model chemotherapy . . . . . . . . 106 5.3 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.3.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.3.2 Replicator equation dynamics . . . . . . . . . . . . . . . . . . . . . . 114 5.3.3 The pairwise prisoner's dilemmas . . . . . . . . . . . . . . . . . . . . 116 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.4.1 Treating the strategy: competitive release . . . . . . . . . . . . . . . 119 5.4.2 Treating the game: adaptive therapy . . . . . . . . . . . . . . . . . . 122 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6 High Performance Computing Techniques 128 6.1 Hybrid MPI + OpenMP Parallelization . . . . . . . . . . . . . . . . . . . . 129 6.2 Error in Hybrid MPI + OpenMP Parallelization . . . . . . . . . . . . . . . 130 Contents v 6.3 Speedup and Eciency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4 Performance Tunability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 7 Future Work 136 7.1 Chemotherapy scheduling and metastasis . . . . . . . . . . . . . . . . . . . 137 7.2 Control theory methods to minimize stable tumor burden . . . . . . . . . . 137 7.3 Heterogeneity drives tumor growth . . . . . . . . . . . . . . . . . . . . . . . 138 7.4 Integrating preclinical and clinical trial data . . . . . . . . . . . . . . . . . . 139 Bibliography 140 List of Figures 2.1 Schematic of the Moran Process . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Emergence of Gompertzian growth via selection . . . . . . . . . . . . . . . . 19 2.3 Markov Point Mutation Diagram . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Tumor tness drives tumor growth . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Moran Process t to Gompertzian Growth Data . . . . . . . . . . . . . . . 25 2.6 Sample Dendritic Phylogenetic Tree . . . . . . . . . . . . . . . . . . . . . . 27 2.7 Eects of varied dose density for early-stage, mid-stage, and late-stage ther- apies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.8 Growth-dependent tumor regression . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Stochastic Moran birth-death process . . . . . . . . . . . . . . . . . . . . . . 44 3.2 Fitness as a function of the selection parameter ww H w C . . . . . . . 50 3.3 Markov Point Mutation Diagram . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4 Gompertzian equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5 Moran birth-death process with selection . . . . . . . . . . . . . . . . . . . . 61 3.6 Tumor initiation prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.7 Comparison of stochastic Moran birth-death process, Gompertzian, and Shannon entropy growth curves . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.8 Emergence of genetic heterogeneity . . . . . . . . . . . . . . . . . . . . . . . 64 3.9 Simulated therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.10 A ow chart of the Moran process with selection and mutation algorithm . 69 4.1 Chemotherapy is a selective agent that alters the tness landscape of cells . 85 4.2 Classical Tumor Regression Laws . . . . . . . . . . . . . . . . . . . . . . . . 87 4.3 Response of murine tumors to 5-Fluorouracil (5-FU) treatment with model best-t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.4 Diminishing returns of dose escalation compared to linear relationship of dose density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Shannon entropy as an index to compare treatment strategies . . . . . . . . 92 4.6 High entropy, LDM-like chemotherapies outperform low entropy MTD-like chemotherapies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.7 High entropy strategies lead to an increase in tumor regression . . . . . . . 97 5.1 Clonal evolution of competitive release . . . . . . . . . . . . . . . . . . . . . 106 5.2 Schematic of competitive release in a tumor . . . . . . . . . . . . . . . . . . 109 vi List of Figures vii 5.3 Dynamics of competitive release under continuous therapy . . . . . . . . . . 112 5.4 Fitness landscape before and during therapy . . . . . . . . . . . . . . . . . . 114 5.5 The eect of growth rate, cost, and fractional resistance under continuous therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.6 The eect of dose on tumor relapse and progression free survival under continuous therapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.7 Dynamic phase portraits before and during chemotherapy . . . . . . . . . . 123 5.8 Adaptive therapy strategy to control resistant population . . . . . . . . . . 125 6.1 Schematic diagram of the Hybrid MPI + OpenMP Parallelization technique. 130 6.2 Strong scalability of Hybrid Parallelization . . . . . . . . . . . . . . . . . . . 131 6.3 Error in the Hybrid MPI + OpenMP primary tumory growth simulation. . 133 Abstract Tumor development is an evolutionary process in which a heterogeneous population of cells with dierent growth capabilities compete for resources in order to gain a proliferative advantage. We describe a cell-molecular based evolutionary mathematical model of tumor development driven by a stochastic Moran birth-death process, in which cancer cells and healthy cells compete for dominance in the population. Each are assigned payos according to a prisoner's dilemma evolutionary game where the healthy cells are the cooperators and the cancer cells are the defectors. The cells in the tumor carry molecular information in the form of a numerical genome which we represent as a four-digit binary string used to dierentiate cells into 16 molecular types. The binary string is able to undergo stochastic point mutations that are passed to a daughter cell after each birth event. The value of the binary string determines the cell tness, with lower t cells (e.g. 0000) dened as healthy phenotypes, and higher t cells (e.g. 1111) dened as malignant phenotypes. The model is used to explore key emergent features associated with tumor development, including tumor growth rates (e.g. Gompertzian growth) as it relates to intratumor molecular heterogeneity (e.g. heterogeneity is a driver of tumor growth). Next, the model is extended to include tumor regression due to chemotherapy, allowing for the design of chemotherapeutic \strategies," or schedules, tailored to dierent tumor growth characteristics. Using the Shannon entropy as a novel tool to quantify dosing strategies, we contrast maximum tolerated dose (MTD) strategies as compared with low dose, high density metronomic strategies (LDM) for tumors with dierent growth rates, concluding that LDM strategies can outperform MTD strategies, especially for fast growing tumors Abstract ix that thrive on long periods of unhindered growth without chemotherapy. The model shows good agreement with tumor growth data for unperturbed tumor growth and regression under drug therapy. The nal section of this thesis is focused on the pre-existing resistant clones present in a tumor at the start of treatment, which remains a major problem in cancer therapeutics today. Tumor relapse and the development of chemotherapeutic resistance is now thought largely to be a consequence of the mechanism of \competitive release" of pre-existing resistant cells in the tumor which are selected for growth after chemotherapeutic agents attack the sub-population of chemo-sensitive cells which had previously dominated the collection of competing sub-clones in the tumor. We explain the important parameters (cost of resistance and initial fraction of resistance) in anticipating the evolutionary adaptations of the tumor in order to design therapies that exploit or mitigate the harmful eects of potential future adaptions. We introduce an \adaptive therapy" concept of utilizing information in the dynamical phase plane to target the growth of the resistant population indirectly, by introducing periods of drug holidays. Chapter 1 Introduction In this thesis, a stochastic birth-death mathematical model is developed that describes the complex tumor kinetics arising from a single malignant cell leading to the evolution of the tness of a tumor and its subsequent growth. The model incorporates cell reproduction and death, mutations, natural selection, and genetic drift. The mathematical framework developed will be used as a platform for testing chemotherapeutic strategies to mitigate the development and growth of the tumor. These models are developed in order to better understand the complex interplay of several biological principals associated with tumor growth that drive the non-equilibrium dynamics of clonal subpopulations in a tumor and how to exploit these dynamics for the optimization of drug schedules that mitigate chemo- therapeutic resistance. This thesis is careful to give the historical context of mathematical tumor modeling and a defense of the need for the model to include evolutionary principals when used in the design and optimization of chemotherapeutic drug schedules. 1 Introduction 2 A fundamental observation of cancer is the upregulation of cell division among malignant cells which leads to the increase in the size of a tumor [1]. This observation led to an explosion of mathematical models with the goal of addressing questions of how tumor size increases over time (see [2] for one summary of historical development of mathematical models of tumor growth). The simplest model of cell division (whereby each cell repeatedly passes through the cell cycle and divides in an unconstrained fashion) is the exponential growth model, implying a constant doubling time for tumors [3]. While shown to be a good model during early tumor growth, long term slowed growth was rst addressed in a seminal paper by Laird [4]. A tumor increasing in size contains a proliferating region of a roughly constant width, which means a diminishing fraction of proliferating cells. This fraction doubles at a constant rate, but the overall tumor growth rate is modeled as decreasing at a constant growth rate, known as Gompertzian growth [5]. Later, these growth models were extended to give insight into the eects of chemotherapy [3, 6, 7]. These predictive drug models were pioneered in a clinical setting by Skipper et. al. [8] who stated that the relationship between dose and tumor cytotoxicity is linear-log (i.e. exponential decay) [6]. Skipper et al. were the rst to develop a set of theoretical laws governing the behavior (and imply the design) of chemotherapy schedules in cancer in the late 1970's [9]. This understanding was later rened by Norton and Simon [7] in the context of Gompertzian growth where the targeted proliferating fraction changes (i.e. decreasing) over time. Introduction 3 1.1 Cancer as an evolutionary and ecological process Yet cancer is an evolutionary and ecological process [10{13]. Ecology is the study of the dynamics of competition within a population of interacting individuals, each with various levels of tness [10]. In ecological interactions, the competing individuals can be classied into many dierent types of relationships including competition, parasitism, predation, mutualism, commensalism and amensalism [10]. One of these relationships observed in nature, cooperation (or mutualism) had puzzled evolutionary biologists for many years [14]. As Nowak eloquently explains in his evolutionary dynamics book [15], \In fact, for biology the problem is as old as evolution itself. Evolutionary progress, the construction of new features, often requires the cooperation of simpler parts that are already available. . . . Animals cooperate to form social structures, groups and societies. Worker bees risk their lives to defend the beehive. Moreover they forgo their own reproductive potential and instead raise the ospring of another individual, the queen. In some bird species, helpers assist the parents in feeding their young. Humans cooperate on a large scale, giving rise to cities, states and countries. Cooperation allows specialization. Nobody needs to know everything. But cooperation is always vulnerable to exploitation." Nowak extensively studied the interplay between cooperation and defection as seen repeatedly in nature and published a book describing methods of applying evolutionary game theory to various biological interactions [15]. Tumor progression can be described by the same phenomena of cooperators (i.e. healthy cells) overtaken by the invasive and exploiting malignant cancer cells. In tumor progression, there is evidence of competition, predation, parasitism, and mutualism between co-evolving Introduction 4 clones [10]. In order to model these evolutionary traits in cancer, many theoretical biol- ogists (including work done in this thesis) have followed Nowak's approach and turned to a coevolutionary game theory framework to study cancer progression (see [14, 16{20]). Evolutionary game theory provides a framework for analyzing contests (competition) be- tween various species in a population as `strategies' and provides mathematical tools to predict the prevalence of each species over time based on the strategies. Game theory was originally developed by John von Neumann and Oscar Morgenstern in order to predict the optimal strategies in competition [21]. The rules of the game govern the payo of each player, where the strategy chosen by each player is the approach on how they play the game. Game theory can be extended to model the evolution of interactions between multiple players in a population competing according to predened strategies in a frame- work called coevolutionary game theory [22]. In this research, we begin with a well-studied game known as the Prisoner's Dilemma [23] to model cell-scale competition and kinetics of a malignant tumor. Studying cancer as a disease of clonal evolution has major implications on tumor progres- sion, prevention and therapy [24, 25]. From an ecological perspective, tumourigenesis is the \process by which the homoeostasis that characterizes a healthy tissue is disrupted either via changes in the tissue microenvironment, or by an invading species (some bacteria and viruses are known to be able to lead to tumourigenesis), or by a local invasion (a resident species producing a brand new one as a result of mutations)" [26]. As such, the need for tumor growth models to include concepts of evolutionary forces such as genetic drift with heritable mutations and natural selection operating on a tness landscape that is in uenced Introduction 5 by tumor microenvironment and the interactions between competing cell types is well doc- umented [27, 28]. Selection will in uence the rates of proliferation and survival, which cause the population of cells within a tumor to progress toward more invasive, metastatic, or therapeutic resistant cell types [29, 30]. As mentioned previously, cancer is an evolutionary process which is driven by mutational events and genetic diversication typically arising via waves of clonal and sub-clonal ex- pansions [10, 25]. As a result, the chance of success is low for all but the most well-designed and tailored therapeutic strategies combating tumor growth. In fact, some tumors develop an enhanced ability to resist future therapeutic attacks [31]. In order to avoid drug resis- tance, chemotherapy has the prospect for greatest impact during the initial stages of tumor progression, before the tumor cell population has been selected for growth and survival, the tumor size is small, and genetic diversication is low [32]. Therapies can work to alter the competition between cancerous and healthy cells by altering the tness landscape, select for genetic stability, optimize therapy timing, and introduce multi-drug therapies [10]. Math modeling has already been used successfully to better understand the evolutionary nature of therapeutic resistance and relapse in cancer (e.g. [33{35]). This thesis will seek to develop a model of drug resistance and rened therapeutic techniques in the context of primary tumor growth with the future goal of extending some of the same principles to the deadly metastatic phase of cancer. Introduction 6 1.2 Chapter Summaries In chapter 2 of this thesis, the mathematical model is introduced and general features such as neutral drift, the tness landscape [10, 11], and genetic heterogeneity [36] are described. The Prisoner's dilemma is introduced as a model of cell-cell interactions (the clonal and subclonal expansion of cancer subpopulations that are competing with the surrounding healthy tissue) and evolution (heritable genetic passenger or driver mutations). A dis- cussion of the model within the context of previously developed mathematical hypotheses such as Gompertzian tumor growth (unperturbed growth) and the Norton-Simon cell kill hypothesis (kinetic resistance to chemotherapy) is presented. In chapter 3, the model is developed in more detail with a special focus on the role of molecular heterogeneity in the growth of the tumor. We present a statistical mechanics model which focuses on a measure of tumor heterogeneity driving tumor growth. We introduce a metric based on Shannon entropy [37] for quantifying molecular heterogeneity in the Prisoner's dilemma evolutionary model. Finally, we highlight the need for early detection of cancer in the high growth fraction (highly susceptible to chemotherapeutic drugs) and low diversity growth regimes. In chapter 4, the model is used to design and test various chemotherapeutic dosing `strate- gies' including current paradigms used in clinical practice such maximum tolerated dose (MTD) and low dose metronomic chemotherapies (LDM) [38, 39]. All strategies are sim- ulated using the stochastic model, and the Shannon entropy is used as a novel tool to quantify the eectiveness of each strategy in reducing the tumor size. Key results show that LDM chemotherapy strategies can outperform MTD strategies [40]. The advantage Introduction 7 is magnied for fast growing tumors that thrive on long periods of unhindered growth between doses. The eectiveness of LDM-like chemotherapy strategies is negligible after a single cycle but grows after each subsequent cycle of repeated chemotherapy. In chapter 5, the model is extended to include resistant cell types and the evolution- ary phenomenom of competitive release is discussed [41]. Important parameters such as initial fraction of resistance must be noted prior to therapy and traditional metrics like progression-free survival (PFS) are often misleading when describing the eectiveness of a therapy which ultimately fails due to resistance upon recurrence [42]. Finally, a proof of concept of adaptive therapy using the phase portrait to inform bang-bang control is introduced. In chapter 6, computational methods used to solve stochastic Moran process models for large cell numbers are explained, and chapter 7 outlines future work beyond the scope of this thesis and exciting directions for future research. 1.3 Emergent Features The famous mathematician and historian of science Jacob Bronowski once dened science as \the organization of our knowledge in such a way that it commands more of the hidden potential in nature." A model is nothing more than a clever organization of knowledge that is used for testing scientic hypotheses. The concept of emergence, or emergent features, describes the properties that result from developing a model formed of its many distinct parts interacting as whole, complex system [43]. While the emergent features are not explicitly modeled or observed within the subsystems or constituent parts of the overall Introduction 8 system, they are viewed as novel properties, functions or behaviors that result by the com- plex interactions between subsystems. It is useful to study emergent features in modeling biological systems to understand the in uence of individual, subsystem parameters on the system-level behavior. Emergent features can also be used as a form of model validation; a model that realizes some biological phenomena on a macro-scale can help uncover the basic mechanisms driving that phenomena [44]. In this thesis, several emergent features resulting from the evolutionary underpinning of tumor kinetics are highlighted. For example, the evolutionary game of the Prisoner's dilemma reveals the benets of early detection and early therapy (see chapter 2). The emergence of molecular heterogeneity results in Gompertzian-like growth (see chapter 2). The novel development of Shannon entropy as a tool to quantify drug dosing strategies leads to an emergence of superior performing, high-entropy strategies for fast growing tumors (see chapter 4). The development of models that describe the cost of resistant mutations leads to the counter-intuitive result that drug therapy holidays outperform continuous drug therapy schedules (see chapter 5). The last section will be devoted to a discussion of future goals of continuing work. We believe mathematical models of the type developed in this thesis will play an increasingly important role in identifying and implementing novel chemotherapeutic strategies and will join the arsenal of tools that clinicians will turn to as they continue to develop strategies to ght cancer. Chapter 2 The prisoner's dilemma as a cancer model Jerey West, Zaki Hasnain, Jeremy Mason, Paul K. Newton The Prisoner's Dilemma as a Cancer Model [45] Convergent Science Physical Oncology 2:3, 035002 (2016) 9 The prisoner's dilemma as a cancer model 10 2.1 Abstract Tumor development is an evolutionary process in which a heterogeneous popu- lation of cells with dierent growth capabilities compete for resources in order to gain a proliferative advantage. What are the minimal ingredients needed to recreate some of the emergent features of such a developing complex ecosys- tem? What is a tumor doing before we can detect it? We outline a math- ematical model, driven by a stochastic Moran process, in which cancer cells and healthy cells compete for dominance in the population. Each are assigned payos according to a prisoner's dilemma evolutionary game where the healthy cells are the cooperators and the cancer cells are the defectors. With point mu- tational dynamics, heredity, and a tness landscape controlling birth and death rates, natural selection acts on the cell population and simulated `cancer-like' features emerge, such as Gompertzian tumor growth driven by heterogeneity, the log-kill law which (linearly) relates therapeutic dose density to the (log) probability of cancer cell survival, and the Norton{Simon hypothesis which (linearly) relates tumor regression rates to tumor growth rates. We highlight the utility, clarity, and power that such models provide, despite (and because of ) their simplicity and built-in assumptions. The prisoner's dilemma as a cancer model 11 2.2 Introduction Cancer is an evolutionary process taking place within a genetically and functionally hetero- geneous population of cells that trac from one anatomical site to another via hematoge- nous and lymphatic routes [11, 46{49]. The population of cells associated with the primary and metastatic tumors evolve, adapt, proliferate, and disseminate in an environment in which a tness landscape controls survival and replication [10]. Tumorigenesis occurs as the result of inherited and acquired genetic, epigenetic and other abnormalities accumu- lated over a long period of time in otherwise normal cells [50, 51]. Before we can typically detect the presence of a tumor, the cells are already competing for resources in a Darwinian struggle for existence in tissues that progressively age and evolve. It is well established that the regenerative capacity of individual cells within a tumor, and their ability to trac multi-directionally from the primary tumor to metastatic tumors all represent signicant challenges associated with the ecacy of dierent cancer treatments and our resulting abil- ity to control systemic spread of many soft-tissue cancers [52, 53]. Details of the metastatic and evolutionary process are poorly understood, particularly in the subclinical stages when tumors are actively developing but not yet clinically visible [54]. It could be argued that in order to truly understand cancer progression at the level in which quantitative predictions become feasible, it is necessary to understand how genetically and epigenetically heteroge- neous populations of cells compete and evolve within the tumor environment well before the tumor is clinically detectable. Additionally, a better understanding of how these pop- ulations develop resistance to specic therapies [55, 56] might help in developing optimal strategies to attack the tumor, slow disease progression, or maintain it at a stable level. The prisoner's dilemma as a cancer model 12 Evolutionary game theory is perhaps the best quantitative framework for modeling evolu- tion and natural selection. It is a dynamic version of classical game theory in which a game between two (or more) competitors is played repeatedly, giving each participant the ability to adjust their strategy based on the outcome of the previous string of games. While this may seem like a minor variant of classical (static) game theory, as developed by the math- ematicians von Neumann and Morgenstern in the 1940's [57], it is not. Developed mostly by the mathematical biologists John Maynard Smith and George Price in the 1970s [22, 58] and Martin Nowak and Karl Sigmund [59, 60] more recently, this dynamic generalization of classical game theory has proven to be one of the main quantitative tools available to evolutionary biologists (if coupled with a tness landscape) whose goal is to understand natural selection in evolving populations. In this biological context, a strategy is not nec- essarily a deliberate course of action, but an inheritable trait [61]. Instead of identifying Nash equilibria, as in the static setting [62, 63], one looks for evolutionary stable strategies (ESS) and xation probabilities [59, 64] of a subpopulation. This subpopulation might be traced to a specic cell with enhanced replicative capacity (high tness), for example, that has undergone a sequence of mutations and is in the process of clonally expanding [24]. A relevant question in that case is what is the probability of xation of that subpopulation? More explicitly, how does one subpopulation invade another in a developing colony of cells? One game in particular, the Prisoner's Dilemma game, has played a central role in cancer modeling (as well as other contexts such as political science and economics) [14, 16{18, 60, 64{78]. It was originally developed by Flood, Dresher and Tucker in the 1950s as an example of a game which shows how rational players might not cooperate, even if it seems to be in their best interest to do so. The evolutionary version of the Prisoner's The prisoner's dilemma as a cancer model 13 Dilemma game has thus become a paradigm for the evolution of cooperation among a group of selsh individuals and thus plays a key role in understanding and modeling the evolution of altruistic behavior [14, 65]. Perhaps the best introductory discussion of these ideas is found in Dawkins' celebrated book, The Selsh Gene [79]. The framework of evolutionary game theory allows the modeler to track the relative frequencies of competing subpopulations with dierent traits within a bigger population by dening mutual payos among pairs within the group. From this, one can then dene a tness landscape over which the subpopulations evolve. The tness of dierent phenotypes is frequency dependent and is associated with reproductive prowess, while the `players' in the evolutionary game compete selshly for the largest share of descendants [64, 78]. Our goal in this article is provide a brief introduction to how the Prisoner's Dilemma game can be used to model the interaction of competing subpopulations of cells, say healthy, and cancerous, in a developing tumor and beyond. 2.3 The prisoner's dilemma evolutionary game An evolutionary game between two players is dened by a 2 x 2 payo matrix which assigns a reward to each player (monetary reward, vacation time, reduced time in jail, etc.) on a given interaction. Let us call the two players A and B. In the case of a prisoner's dilemma game between cell types in an evolving population of cells, let there be two subpopulations of cell types which we will call `healthy', and `cancerous'. We can think of the healthy cells as the subpopulation that is cooperating, and the cancer cells as formerly cooperating cells that have defected via a sequence of somatic driver The prisoner's dilemma as a cancer model 14 mutations. Imagine a sequence of `games' played between two cells (A and B) selected at random from the population, but chosen in proportion to their prevalence in the population pool. Think of a cancer-free organ or tissue as one in which a population of healthy cells are all cooperating, and the normal organ functions are able to proceed, with birth and death rates that statistically balance, so an equilibrium healthy population is maintained (on average). Now imagine a mutated cell introduced into the population with enhanced proliferative capability as encoded by its genome as represented as a binary sequence of 0's and 1's carrying forward its genetic information (which is passed on to daughter cells). A schematic diagram associated with this process is shown in Figure 1. We can think of this cancer cell as a formerly cooperating cell that has defected and begins to compete against the surrounding population of healthy cells for resources and reproductive prowess. From that point forward, one can imagine tumor development to be a competition between two distinct competing subpopulations of cells, healthy (cooperators) and cancerous (defectors). We are interested in the growth rates of a `tumor' made up of a collection of cancer cells within the entire population, or equivalently, we are interested in tracking the proportion of cancer cells,i(t), vs. the proportion of healthy cells, Ni(t), in a population ofN cells comprising the simulated tissue region. To quantify how the interactions proceed, and how birth/death rates are ultimately as- signed, we introduce the 2 x 2 prisoner's dilemma payo matrix: A = 0 B B @ a b c d 1 C C A = 0 B B @ 3 0 5 1 1 C C A : (2.1) The prisoner's dilemma as a cancer model 15 What denes a prisoner's dilemma matrix are the inequalities c>a>d>b. The chosen values in (2.1) are relatively standard, but not unique 1 . The essence of the prisoner's dilemma game is the two players compete against each other, and each has to decide what best strategy to adopt in order to maximize their payo. This 2 x 2 matrix assigns the payo (e.g. reward) to each player on each interaction. My options, as a strategy or, equivalently, as a cell type, are listed along the rows, with row 1 associated with my possible choice to cooperate, or equivalently my cell type being healthy, and row 2 associated with my possible choice to defect, or equivalently my cell type being cancerous. Your options are listed down the columns, with column 1 associated with your choice to cooperate (or you being a healthy cell), and column 2 associated with your choice to defect (or you being a cancer cell). The analysis of a rational player in a prisoner's dilemma game runs as follows. I do not know what strategy you will choose, but suppose you choose to cooperate (column 1). In that case, I am better o defecting (row 2) since I receive a payo of 5 instead of 3 (if I also cooperate). Suppose instead you choose to defect (column 2). In that case, I am also better o defecting (row 2) since I receive a payo of 1 instead of 0 (if I were to have cooperated). Therefore, no matter what you choose, I am better o (from a pure payo point of view) if I defect. What makes this game such a useful paradigm for strategic interactions ranging from economics, political science, biology, and even psychology [58, 65, 78] is the following additional observation. You will analyze the game in exactly the same way I did (just switch the roles of me and you in the previous rational analysis), so you will also decide to defect no matter what I do. The upshot if we both defect is that we will each receive a payo of 1, instead of each receiving a 1 A general investigation of how the values in the PD payo matrix aects evolutionary dynamics of the subpopulations is addressed in [80]. The prisoner's dilemma as a cancer model 16 payo of 3 if we had both chosen to cooperate. The defect-defect combination is a Nash equilibrium [62, 63], and yet it is sub-optimal for both players and for the system as a whole. Rational thought rules out the cooperate-cooperate combination which would be better for each player (3 points each) and for both players combined (6 points). In fact, the Nash equilibrium strategy of defect-defect is the worst possible system wide choice, yielding a total payo of 2 points, compared to the cooperate-defect or defect-cooperate combination, which yields a total payo of 5 points, or the best system-wide strategy of cooperate-cooperate yielding a total payo of 6 points. The game becomes even more interesting if it is played repeatedly [78], stochastically [77], and with spatial structure [81] with each player allowed to decide what strategy to use on each interaction so as to accumulate a higher payo than the competition over a sequence of N games. In order to analyze this kind of an evolving set-up, a tness function must be introduced based on the payo matrix A. Let us now switch our terminology so that the relevance to tumor cell kinetics becomes clear. When modeling cell competition, one has to be careful about the meaning of the term `choosing a strategy'. Cells do not choose a strategy, but they do behave in dierent ways depending on whether they are normal healthy cells cooperating as a cohesive group, with birth and death rates that statistically balance, or whether they are cancer cells with an overactive cell division mechanism (as triggered by the presence of oncogenes) and an underactive `break' mechanism (as triggered by the absence of tumor suppressor genes) [49]. In our context, it is not the strategies that evolve, as cells cannot change type based on strategy (only based on mutations), but the prevalence of each cell type in the population is evolving, with the winner identied as the sub-type that rst saturates in the population. The prisoner's dilemma as a cancer model 17 2.4 A tumor growth model Consider a population of N cells driven by a stochastic birth-death process as depicted in Figure 2.1, with red cells depicting cancer cells (higher tness) and blue cells depicting healthy cells (lower tness, but cooperative). We model the cell population as a stochastic Moran process [82] of N cells, `i' of which are cancerous, `Ni' of which are healthy. If each cell had equal tness, the birth-death rates would all be equal and a statistical balance would ensue. At each step, a cell is chosen (randomly but based on the prevalence in the population pool) and eliminated (death), while another is chosen to divide (birth). If all cells had equal tness, the birth/death rates of the cancer cells would be i=N, while those of the healthy cells would be (Ni)=N. With no mechanism for introducing a cancer cells in the population, the birth/death rates of the healthy cells would be equal, and no tumor would form. Now, introduce one cancer cell into the population of healthy cells, as shown in Figure 2.1a. At each step, there would be a certain probability of this cell dividing (P i;i+1 ), being eliminated (P i;i1 ), or simply not being chosen for either division or death (P i;i ). Based on this random process, it might be possible for the cancer cells to saturate the population, as shown by one simulation in Figure 2.2 depictingN = 1000 cells, with initiallyi = 1 cancer cell, andNi = 999 healthy cells. However, the growth curve would not show any distinct shape (Figure 2.2 (black)), and might well become extinct after any number of cell divisions, as opposed to reaching saturation. But we emphasize that without mutational dynamics, heritability, and natural selection operating on the cell population, the shape of the growth curve would look random, and we know this is not how tumors tend to grow [4, 83]. By The prisoner's dilemma as a cancer model 18 (a) (b) Figure 2.1: Schematic of the Moran Process | (a) The number of cancer cells, i, is dened on the state space i = 0; 1;:::;N where N is the total number of cells. The cancer population can change at most by one each time step, so a transition exists only between state i and i 1, i, and i + 1. (b) During each time step, a single cell is chosen for reproduction, where an exact replica is produced. With probability m (0m 1), a mutation may occur. contrast, Figure 2.2 (red) shows a Gompertzian growth curve starting with exponential growth of the cancer cell subpopulation, followed by linear growth, ending with saturation. The growth rate is not constant throughout the full history of tumor development, but after an initial period of exponential growth, the rate decelerates until the region saturates with cancer cells. The basic ingredients necessary to sustain Gompertzian growth seem to be: an underlying stochastic engine of developing cells, mutational dynamics, heritability, and a tness landscape that governs birth and death rates giving rise to some sort of natural selection. The prisoner's dilemma as a cancer model 19 Figure 2.2: Emergence of Gompertzian growth via selection | Random drift (black) plotted for a single simulation of 10 3 cells for 4 10 4 generations shows no par- ticular shape. A single simulation of the Moran process (red) with selection (w = 0:5) and mutations (m = 0:1) gives rise to the characteristic S-shaped curve associated with Gompertzian growth. 2.4.1 Mutations and heritability Each of the N cells in our simulated population carries with it a discrete packet of infor- mation that represents some form of molecular dierences among the cells. In our model, we code this information in the form of a 4-digit binary string from 0000 up to 1111, giving rise to a population made up of 16 distinct cell types. At each discrete step in the birth- death process, one of the digits in the binary string is able to undergo a point mutation [50, 84], where a digit spontaneously ips from 0 to 1, or 1 to 0, with probability p m . The mutation process is shown in Figure 2.1, while a mutation diagram is shown in Figure 2.3 in the form of a directed graph. This gure shows the possible mutational transitions that can occur in each cell, from step to step in a simulation. A typical simulation begins with a population of N healthy cells, all with identical binary strings 0000. The edges on the directed graph represent possible mutations that could occur on a given step. The rst 11 The prisoner's dilemma as a cancer model 20 binary string values (0-10) represent healthy cells in our model that are at dierent stages in their evolutionary progression towards becoming a cancer cell (the exact details of this genotype to phenotype map do not matter much). Mutations strictly within this subpop- ulation are called passenger mutations as the cells all have the same tness characteristics. The rst driver mutation occurs when a binary string reaches value 11-15. The rst cell that transitions from the healthy state to the cancerous state is the renegade cell in the population that then has the potential to clonally expand and take over the population. How does this process occur? Figure 2.3: Markov Point Mutation Diagram | Diagram shows 16 genetic cell types based on 4-digit binary string and the eect of a point mutation on each cell type. Blue indicates healthy cell type (0000 | 1010), red indicates cancerous cell type (1011 | 1111). Black arrows indicate passenger mutations (healthy to healthy or cancer to cancer), red arrows indicate driver mutations (healthy to cancer). 2.4.2 The tness landscape At the heart of how the Prisoner's Dilemma evolutionary game dictates birth and death rates which in turn control tumor growth, is the denition of cell tness. Let us start by The prisoner's dilemma as a cancer model 21 laying out the various probabilities of pairs of cells interacting and clearly dening payos when there are i cancer cells, and Ni healthy cells in the population. The probability that a healthy cell interacts with another healthy cell is given by (Ni 1)=(N 1), whereas the probability that a healthy cell interacts with a cancer cell is i=(N 1). The probability that a cancer cell interacts with a healthy cell is (Ni)=(N 1), whereas the probability that a cancer cell interacts with another cancer cell is (i 1)=(N 1). In a xed population of N cells, with i cancer cells, the number of healthy cells is given by Ni. The average payo of a single cell ( H ; C ), is dependent on the payo matrix value weighted by the relative frequency of types in the current population: H i = a(Ni 1) +bi N 1 (2.2) C i = c(Ni) +d(i 1) N 1 (2.3) Here, a = 3, b = 0, c = 5, d = 1 are the entries in the Prisoner's Dilemma payo matrix (2.1). For the Prisoner's dilemma game, the average payo of a single cancer cell is always greater than the average payo for a healthy cell (Figure 2.4c). With the invasion of the rst cancer cell, the higher payo gives a higher probability of survival when in competition with a single healthy cell. Selection acts on the entire population of cells as it depends not on the payo, but on the eective tness of the subtype population. The eective tness of each cell type (f H , The prisoner's dilemma as a cancer model 22 f C ) is given by the relative contribution of the payo of that cell type, weighted by the selection pressure: f H i = 1w +w H i (2.4) and the tness of the cancer cells as: f C i = 1w +w C i (2.5) The probability of birthing a new cancer cell depends on the relative frequency (random drift) weighted by the eective tness, and the death rate is proportional to the relative frequency. The transition probabilities can be written: P i;i+1 = if C i if C i + (Ni)f H i Ni N (2.6) P i;i1 = (Ni)f H i if C i + (Ni)f H i i N (2.7) P i;i = 1P i;i+1 P i;i1 ; P 0;0 = 1; P N;N = 1: (2.8) In the event of the rst driver mutation, the rst cancer cell is birthed. At the beginning of the simulation, the eective tness of the healthy population is much greater than the The prisoner's dilemma as a cancer model 23 tness of the cancer population (Figure 2.4b). This is because although the single cancer has a higher payo than any of the healthy cells, the number of healthy cells far outnumber the single cancer cells. That single cancer cell initiates a regime of explosive high growth and the tness of the cancer population steadily increases. Cancer cells are continually competing with healthy cells and receiving a higher payo in this regime (compare the payo entries of a cancer cell receiving c = 5 vs a healthy cell receiving b = 0). At later times, growth slows because cancer cells are competing in a population consisting mostly of other cancer cells. The payo for a cancer cell is dramatically lower when interacting with a cancer cell (observe the payo entry of both cancer cells receiving d = 1 when interacting). As the cancer population grows, the payo attainable decreases and growth slows. In addition, the average tness of the total population steadily declines because each interaction derives less total payo, from c +b = 5 to d +d = 1. It is precisely the payo structure of the Prisoner's Dilemma matrix that produces this declining average population tness as the cancer cells saturate the population. Although they receive higher payos than healthy cells on pairwise interactions, these cancer-healthy interactions mostly take place early on in the evolution of the tumor. As the cancer cells take over the population, most of the interactions take place between pairs of cancer cells (i.e. they eventually start competing with each other) causing the population tness to decline. This complex process of competition among cell types and survival of subpopulations, where defection is selected over cooperation, produces a Gompertzian growth curve shown in Figure 2.5, and compared with a compilation based on a wide range of data rst shown in [4, 83]. It is now well established that tumor cell populations (and other competing populations, such as bacteria and viral populations) generally follow this growth pattern, The prisoner's dilemma as a cancer model 24 (a) (b) (c) Figure 2.4: Tumor tness drives tumor growth | (a) The average of 25 stochastic simulations (N = 1000 cells,w = 0:5,m = 0:1) is plotted for 20,000 cell divisions to show the cancer cell population (defectors) saturating. The pink lines delineate the regions of tumor growth (dened by the maximum and minimum points of the second-derivative of i(t)). (b) Fitness of the healthy population, cancer population, and total population plotted for the range cancer cell proportion. (c) Average payo of a single healthy cell, cancer cell, and all cells plotted for the range cancer cell proportion. although the literature is complicated by the fact that dierent parts of the growth curve have vastly dierent growth rates [4, 83], and it is nearly impossible to follow the growth The prisoner's dilemma as a cancer model 25 of a population of cancer cells in vivo from the rst cancer cell through to an entire tumor made up of O(10 9 10 12 ) cells. Growth rates are typically measured for a short clinical time period [4, 83], and then extrapolated back to the rst renegade cell, and forward to the fully developed tumor population. Figure 2.5: Moran Process t to Gompertzian Growth Data | The mean and deviation of 25 stochastic simulations (N = 10 3 cells,w = 0:7,m = 0:3) is overlaid on data from a \normalized" Gompertzian [4, 83]. Values for m and w were chosen by implementing a least-squares t to the data over a range of m (0m 1), and w (0w 1). Pink lines delineate regions of growth (dened by the maximum and minimum points of the second-derivative of i(t)). 2.4.3 Heterogeneity drives growth Insights into the process by which growth rates vary and conspire to produce a Gom- pertzian shape can be achieved by positing that growth is related to molecular and cellular heterogeneity of the developing population [36, 48, 85]. Indeed, an outcome of the model is that molecular heterogeneity (i.e. the dynamical distribution of the 4-digit binary string The prisoner's dilemma as a cancer model 26 0000|1111 making up the population of cells) drives growth. Consider entropy [86, 87] of the cell population as a measure of heterogeneity: E(t) = N X i=1 p i log 2 p i (2.9) (here, log is dened as base 2). The probability p i measures the proportion of cells of typei, withi = 1;:::; 16 representing the distribution of binary strings ranging from 0000 to 1111. We typically course-grain this distribution further so that cells having strings ranging from 0000 up to 1010 are called `healthy', while those ranging from 1011 to 1111 are `cancerous'. Then growth is determined by: dn E dt =E(t) (2.10) It follows from (3.9) that the cancer cell proportionn E (t) can be written in terms of entropy as: n E (t) = Z t 0 E(t)dt (2.11) This relationship between growth of the cancer cell population and entropy is pinned down and detailed in [82]. We consider it to be one of the key emergent features of our simple model. A typical example of the emergence of genetic heterogeneity in our model system is shown in the form of a phylogenetic tree in Figure 2.6. This particular tree is obtained via The prisoner's dilemma as a cancer model 27 a simulation of only 30 healthy phenotypic cells (0000), which during the course of a simulation expand out (radially in time) to form a much more heterogeneous population of cells at the end of the simulation. In our model, the genetic time-history of each cell is tracked and the population can be statistically analyzed after the simulation nishes. Figure 2.6: Sample Dendritic Phylogenetic Tree | Sample dendritic phylogenetic tree tracking point mutations as time extends radially, depicting the emergence of molec- ular heterogeneity. The tree shows a simulation of 30 cells (all with genetic string 0000 at the beginning of the simulation) with strong selection (w = 1, m = 0:2). Pathways are color coded to indicate genetic cell type. The prisoner's dilemma as a cancer model 28 2.5 Simulated drug dosing strategies and therapeutic re- sponse Figure 2.7 shows the clear advantage of early stage therapy in our model system. We compare the eect of therapy given at an early stage, mid-stage, and late stages of the Gompertzian growth of the tumor. The dashed black Gompertzian curve is the freely growing cancer cell population. In each of the gures, we depict the eect of a range of drug dose densities, D, where D =ct: (2.12) The dose density is the product of the drug concentration, c, and the time over which the therapy is administered, t, (2.12). Here, the drug concentration value is a weighting (0 c 1) which determines the intensity of the drug treatment. A higher value of c will alter the selection pressure in favor of healthy cells (and to the disadvantage of cancer cells) more dramatically. Figure 2.7 varies the drug dose density by varying the drug concentrations (c = 0:2; 0:4; 0:6; 0:8; 1:0) administered for a constant time (t = 5000 cell divisions, black arrow). The colored curves show the subsequent decline of the cancer cell population under therapeutic pressure. Clearly, to obtain the desired eect of driving the cancer cell population down to man- ageable levels, one needs to (i) use a sucient dose density, and (ii) initiate therapy early enough in the growth history. These gures are meant to paint a broad brush with respect to the simulated advantages of early therapy and to show the capability of the model with The prisoner's dilemma as a cancer model 29 respect to addressing questions of this type in a quantitative way. A detailed investigation is left for a separate publication. (a) (b) (c) Figure 2.7: Eects of varied dose density for early-stage, mid-stage, and late- stage therapies | An average of 25 stochastic simulations of unperturbed tumor growth (N = 10 3 cells, w = 0:5, m = 0:1, no therapy) is plotted (black dashed line). The eect of varied drug dose density (eqn. 2.12), is shown by administering a range of drug concentration values (c = 0:2; 0:4; 0:6; 0:8; 1:0) for constant length of time (t = 5000 cell divisions, black solid arrows). This process is repeated for (a) high growth, early-stage, (b) linear growth, mid-stage, and (c) slowed growth, late-stage. The kill eect is highest for high drug concentration values (i.e. high dose density) and early therapy. An established empirical law which relates drug dose density to its eectiveness in killing o cancer cells is known as the `log-kill' law [9]. The log kill law states that a given dose of chemotherapy kills the same fraction of tumor cells (as opposed to the same number of tumor cells), regardless of the size of the tumor at the time the therapy is administered [9], a consequence of exponential growth with a constant growth rate. This eect is best illustrated on a dose-response curve, plotting the dose density, D, with respect to the probability of tumor cell survival, P S . Thus, the log-kill law states the following: log(P S ) =D: (2.13) As an example, if there are 1000 cancer cells in a tumor population, and the rst therapy dose kills o 90% of them (i.e. = 0:9), then after the rst round of therapy there will The prisoner's dilemma as a cancer model 30 be 100 cancer cells remaining. If a second round of therapy is administered, exactly as the rst round, starting soon enough so that no new cancer cells have formed, then this next round will also kill o 90% of the cells, leaving 10 cells, and so on for each future round of therapy. In a sense, since the rst round killed 900 cells, while the second identical round killed only 90 cells, the population gets increasingly more dicult to kill o using the same treatment on each cycle. The log-kill law, a fundamentally static law (it does not say anything about the relationship of the fraction of cells killed vs. the growth rate of the tumor), is veried in our model system, as shown in the dose response curve in Figure 2.8d. On the x-axis, we increase the dose density D, and we plot the number of surviving cancer cells. The slope of this straight line (verifying the log-kill law) is the rate of regression of the tumor, . Alternatively, can be estimated using an exponential t of i(t) during therapy (i.e. i(t) =i 0 exp((tt 0 )), where i 0 is the initial tumor size and t 0 is the time therapy is initiated). So how is the rate of regression, , related to the growth rate of the tumor, ? This is relevant, since we know from the shape of the Gompertzian curve, the growth rate is highest (exponential) at the beginning stage of tumor development and lowest at the late saturation stage. Figure 2.8a shows therapy is more eective (i.e. a higher rate of regression, ) for earlier stage therapy. These early stage therapies correspond to a higher growth rate, shown in 2.8b. The Norton-Simon hypothesis [7, 88, 88] states that the rate of regression is proportional to the instantaneous growth rate for an untreated tumor of that size at the time therapy is rst administered. A faster growing tumor (early stage) should show a higher rate of regression than a more slowly growing tumor (late stage). This hypothesis is also veried in our model system, and shown clearly in Figure 2.8c. The prisoner's dilemma as a cancer model 31 (a) (b) (c) (d) Figure 2.8: Growth-dependent tumor regression | (a) An average of 25 stochastic simulations of unperturbed tumor growth (N = 10 3 cells, w = 0:5, m = 0:1, no therapy) is plotted (black dashed line) with (b) corresponding instantaneous growth rate, (t), of the unperturbed tumor (red). Tumor regression, , (estimated using an exponential t of i(t) during therapy, shown in legend) during therapy (constant dose density: c = 1:0, t = 2000) is calculated for a high growth, early-stage therapy (purple), linear-growth, mid-stage therapy (green), and late-stage, slowed growth (light blue); (c) This process is repeated for a full range of growth rates (between vertical blue dashed lines). Average values of are plotted with standard deviations. Regression is proportional to growth rate (linear t in red), with higher regression rates associated with high growth rates of early stage tumors. (d) Tumor regression, , can also be calculated as the slope of a dose response curve (red), where therapy is administered for a range of dose densities (0:7c 1:0) for a single timepoint, 8000 cell divisions (i.e. single growth rate). The prisoner's dilemma as a cancer model 32 The reality of this growth-dependent tumor regression rate eect (where early stage faster growing tumors are more vulnerable to therapy than later stage, slowly growing tumors) dramatically reinforces the need to administer drug treatment early in tumor progression when growth rates are high and there are fewer cancer cells to kill o. 2.6 Mathematical modeling and tumor analytics It is important to keep in mind that no mathematical model captures all aspects of reality, so choices must be made which involve prioritizing the features that are most essential in capturing the essence of a complex process and which are not. Most experts now agree that the evolutionary processes in a tumor played out among subpopulations of competing cells are key to understanding aspects of growth and resistance to chemotherapy, which will ultimately lead the way toward a quantitative understanding of tumor growth and cancer progression [10, 49, 53]. The paradigm of the cancer cell subpopulation and the healthy cell subpopulation competing as the defectors and cooperators in a Prisoner's Dilemma evolutionary game has been useful in obtaining a quantitative handle on many of these processes and frames the problem in an intuitive yet predictive way. Nonetheless, the mathematical `taste' of the modeler plays a role in what techniques are se- lected and ultimately where the spotlight shines. This fact makes clinicians uncomfortable and can lead to deep suspicion of the mathematical modeling enterprise as a whole. Aren't the outcomes and predictions of mathematical models a straightforward consequence of the model assumptions? Once those choices are made, isn't the cake already baked? So why should we be surprised if you tell us it tastes good? Why not simply use tried and The prisoner's dilemma as a cancer model 33 true statistical tools like regression methods to curve-t the data directly, with no built in assumptions, and be satised with uncovering correlations and trends? Clinicians (and experimentalists, in general) feel that they are dealing directly with reality, so why mess around with `toy' systems based on possibly `ad hoc' or incorrect assumptions that create articial realities that may or may not be relevant? To a theoretician, calling their as- sumptions ad hoc, as opposed to natural, is as insulting as calling a clinician sloppy and uncaring (try this for yourself at the next conference you go to! But please use the term `somewhat ad hoc' to lessen the blow.) And if you want to deliver an even harsher insult, you could comment that the model seems like an exercise in curve tting. But the usefulness of mathematical models built on simplied assumptions is well estab- lished in the history of the physical sciences, as detailed beautifully in Peter Dear's book, The Intelligibility of Nature: How Science Makes Sense of the World [89]. Bohr's simple model of the structure of the atom was crucial in moving the community forward towards a deeper understanding of cause and eect, and ultimately pushing others to develop more realistic atomic models. The same could be said for many other important, but ultimately discarded models of reality (e.g. the notion of aether used as a vehicle to understand the mysterious notion of action-at-a-distance [89]) now relegated to footnotes in the history of the physical sciences. Lessons from this history highlight the importance of using the principle of Occam's razor (law of parsimony) as a heuristic guide in developing models: (1) keep things simple, but not too simple; (2) see what can be explained by using a given set of assumptions, and try to identify what is either wrong or cannot be explained; (3) add complexity to the model, but do this carefully. Since ultimately, the model will always be wrong (with respect to The prisoner's dilemma as a cancer model 34 some well chosen and specic new question being posed about a system), it is important that it be useful as a vehicle of intelligibility [89] associated with the set of questions surrounding the phenomena it was built to explain. Answers to some new questions will be found using the model as a temporary crutch, and new questions will emerge in the process that had not yet been asked, as their relevance had never previously been realized. A new quantitative language will emerge in which aspects of the model will be associated with the underlying reality it is attempting to describe, predictions will be easier to frame and test, and shortcomings will be exposed. In his famous article [90], Eugene Wigner writes compellingly that `the miracle of the appropriateness of the language of mathematics for the formulation of the laws of physics is a wonderful gift which we neither understand nor deserve. We should be grateful for it and hope that it will remain valid in future research and that it will extend, for better or for worse, to our pleasure, even though perhaps also to our baement, to wide branches of learning.' In general, the more complex the model (as measured, for example, by the number of independent parameters associated with it), the less useful it will be, and the less likely it is to be adopted by the community at large. After all, if the model is as complex as the phenomena it was built to understand, why not stick with reality? Eective models can be thought of as low-dimensional approximations of reality, surrogates that help us bootstrap our way forward. They arise as the outcome of a complex balancing act be- tween simplicity of the ingredients, and complexity of the reality the model is meant to describe. They generally do not arise in a vacuum, but are built in the context of informed and sustained discussions among people with dierent expertise. In the context of medical The prisoner's dilemma as a cancer model 35 oncology, this means physical scientists developing ongoing interactions with clinical oncol- ogists, radiologists, pathologists, molecular and cell biologists and other relevant medical specialists. Appropriate data is a necessary ingredient in developing and testing any successful model, and treasure troves of medical data sit unexamined in patient les and government databases across the country waiting to be put to good use. There is no doubt that they are telling an interesting and important story that we have yet to fully understand. It is not currently possible for the computer to simulate all of the complex, relevant, and systemic ingredients at play to faithfully recreate all aspects of cancer progression and treatment response in patients. It is hard to imagine that a deep and actionable understanding can ever be ob- tained without the combined use of data, models, and computer simulations to help guide us and highlight some of the underlying causal mechanisms of this complex and deadly disease. Chapter 3 An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity driving Gompertzian growth Jerey West, Zaki Hasnain, Paul Macklin, Paul K. Newton An Evolutionary Model of Tumor Cell Kinetics and the Emergence of Molecular Heterogeneity Driving Gompertzian Growth [82] SIAM Review 58:4, 716-736 (2016) 36 An evolutionary model of tumor cell kinetics 37 3.1 Abstract We describe a cell-molecular based evolutionary mathematical model of tu- mor development driven by a stochastic Moran birth-death process. The cells in the tumor carry molecular information in the form of a numerical genome which we represent as a four-digit binary string used to dierentiate cells into 16 molecular types. The binary string is able to undergo stochastic point mu- tations that are passed to a daughter cell after each birth event. The value of the binary string determines the cell tness, with lower t cells (e.g. 0000) dened as healthy phenotypes, and higher t cells (e.g. 1111) dened as malig- nant phenotypes. At each step of the birth-death process, the two phenotypic sub-populations compete in a prisoner's dilemma evolutionary game with the healthy cells playing the role of cooperators, and the cancer cells playing the role of defectors. Fitness, birth-death rates of the cell populations, and overall tumor tness are dened via the prisoner's dilemma payo matrix. Mutation parameters include passenger mutations (mutations conferring no tness ad- vantage) and driver mutations (mutations which increase cell tness). The model is used to explore key emergent features associated with tumor devel- opment, including tumor growth rates as it relates to intratumor molecular heterogeneity. The tumor growth equation states that the growth rate is pro- portional to the logarithm of cellular diversity/heterogeneity. The Shannon entropy from information theory is used as a quantitative measure of hetero- geneity and tumor complexity based on the distribution of the 4-digit binary An evolutionary model of tumor cell kinetics 38 sequences produced by the cell population. To track the development of het- erogeneity from an initial population of healthy cells (0000), we use dynamic phylogenetic trees which show clonal and sub-clonal expansions of cancer cell sub-populations from an initial malignant cell. We show tumor growth rates are not constant throughout tumor development, and are generally much higher in the subclinical range than in later stages of development, which leads to a Gompertzian growth curve. We explain the early exponential growth of the tumor and the later saturation associated with the Gompertzian curve which results from our evolutionary simulations using simple statistical mechanics principles related to the degree of functional coupling of the cell states. We then compare dosing strategies at early stage development, mid-stage (clinical stage), and late stage development of the tumor. If used early during tumor development in the subclinical stage, well before the cancer cell population is selected for growth, therapy is most eective at disrupting key emergent features of tumor development. 3.2 Introduction At the molecular and cellular levels, cancer is an evolutionary process [10, 91, 91, 91] driven by random mutational events [91{94] responsible for genetic diversication which typically arises via waves of clonal and sub-clonal expansions [24, 25], operating over an adaptive tness landscape in which Darwinian selection favors highly proliferative cell phenotypes which in turn drive rapid tumor growth [4, 95, 96]. The tumor environment should be An evolutionary model of tumor cell kinetics 39 viewed as a complex Darwinian adaptive eco-system consisting of cell types which have evolved over many years [10]. As a result, all but the most well designed and tailored therapeutic strategies often deliver disappointing outcomes and potentially introduce a potent new source of selective pressure for the proliferation of variant cells which develop an enhanced ability to resist future therapeutic assaults [56, 97{100]. The prospects for in uencing and controlling such a system are likeliest at the emerging early stages of tumor development when the cell population has not yet been selected for growth and survival, and the tumor size is small. But by the time a typical tumor becomes clinically detectable (often after several years of growth), it already contains O(10 8 ) or more malignant cells (and potentially occupies a volume of 12 mm 3 ), some of which may have entered the blood circulation [95]. Since there is very little human data available in this early subclinical stage of tumor development, computational models can serve as a useful surrogate in this critical developmental stage which clearly in uences and determines many important emergent features of the tumor at later stages. Our goals in this paper are to describe a mathematical model for stochastic cell kinetics in the beginning stages of tumor development (from a single malignant cell) that includes cell reproduction and death, mutations, evolution, and the subsequent emergence of genetic heterogeneity well documented in many soft-tissue tumors [35, 48, 101{107]. The model is a computational one, driven by a stochastic Moran (birth-death) process with a nite cell population, in which birth-death rates are functions of cell tness. The tness is deter- mined by the cell's numerical genome in the form of a four-digit binary string capable of undergoing point mutational dynamics with one digit in the string ipping values stochasti- cally. The corresponding numerical value of the binary string determines whether the cell is An evolutionary model of tumor cell kinetics 40 healthy (low-tness) or cancerous (high tness). These two classes of cells compete against each other at each birth-death event, with tness calculated according to the payo matrix associated with the prisoner's dilemma evolutionary zero-sum game [14, 15, 18, 22]. The healthy cells play the role of cooperators, while the cancer cells play the role of defectors [14, 15]. Our goal is to understand how the model parameters: passenger (m p ) and driver mutation rates (m d ), selection strength (w), birth and death rates, aect tumor growth characteristics, such as tumor growth rates, xation probabilities of malignant and healthy cell types, saturation rates of cancer cells, and the emergence of genetic heterogeneity in a tumor at later stages of development when the tumor is clinically detectable. An important outcome of the model is that growth of the cancer cell population is directly in uenced by the intratumor heterogeneity (represented as the distribution of the 4-digit binary strings throughout the cell population), with high heterogeneity driving more rapid growth. The connection between heterogeneity and growth has been discussed in the liter- ature [83, 85, 107{110]. We quantify heterogeneity in a tumor using tools from information theory [37, 87], as well as quantitative analysis of phylogenetic trees associated with clonal and sub-clonal expansions [48, 111] in the developing tumor. Because our numerical simu- lations are carried out from initial conditions corresponding to a homogeneous population of healthy cells (0000) all the way to a saturated population of cancer cells, we can use the model to test basic dose and scheduling strategies [112, 113] at the very early stages of tumor development in the subclinical range, well before a tumor would be clinically detectable by current technology. Our point of view is that this emerging subclinical tu- mor should be more amenable (and potentially vulnerable) to a well planned therapeutic assault than a more mature tumor comprised (on average) of larger numbers of cells with An evolutionary model of tumor cell kinetics 41 more aggressive proliferative capabilities (having undergone generations of selection), that are potentially in the early stages of migration to other organs. More complex features that might in uence early stage dynamics, like human-immune response [108] and the tumor microenvironment [114] are not included in this model in order to keep things as simple and clear as possible. 3.3 Description of the model The ingredients in our model includes a stochastic birth-death process that is the engine which drives tumor growth, with heritable mutations operating over a tness landscape so that natural selection can play out over many cell division timescales. Genetic mutations (point mutations) are modeled using a four-digit binary string of information that each cell carries with it. 1 This simple sequence divides the cells into 16 dierent \genotypes", ranging from 0000 up to 1111, and this information is passed on to the daughter cell during a birth event. The birth-death replacement process is based on a tness function dened in terms of interactions quantied by the prisoner's dilemma payo matrix which operates on two general classes of cells: healthy (the cooperators), and cancerous (the defectors). Natural selection acts on each generation of the cell population as the computational simulation proceeds on a cell division timescale. In this version of the model we typically simulate up to O(10 11 ) cell divisions, so our mutation rates are chosen to be relatively high to 1 To be clear, the four digit sequence is not meant as a bare-bones representation of the full human genome, but as a simple representation of the relevant dierences in genetic information contained in dier- ent cells, allowing us to course-grain the cells into 16 dierent categories based on their genetic/epigenetic proles. An evolutionary model of tumor cell kinetics 42 accommodate these somewhat modest timescales. See [93] for discussions on mutation rates in cancer. 3.3.1 The Moran birth-death process The stochastic engine [115] that drives tumor growth in our model is a nite cell-based Moran process consisting of a population of N cells, divided into two sub-populations consisting of i cancer cells, and Ni healthy cells. In all of our simulations, N is large enough so that there is not a signicant dierence between the results from our nite- cell model and the (deterministic) replicator equation approach for innite populations, a connection that is discussed in detail in [116]. At each time-step in the simulation, one cell is chosen for reproduction and one cell is chosen for elimination. The cells are chosen randomly, based on their prevalence in the population pool which, in turn, is weighted by the tness function based on a chosen payo matrix. The probability of choosing a cancer cell at any given step isi=N, while the probability of choosing a healthy cell is (Ni)=N. As it unfolds, the process is a stochastic birth-death process where the total population size,N, stays constant and the number of cancer cells in the population,i, is the stochastic state variable. At any given step, the probability of transitioning from i cancer cells to j cancer cells is denotedP ij , withj =i + 1 orj =i 1. These probabilities are determined by the birth/death rates associated with the cancer cell population, which in turn are determined by a cell population tness function. Each cell carries with it a binary string in the form of a four digit binary sequence from 0000 up to 1111. This denes 16 dierent cell types, which are course-grained into two groups: healthy cells (0000 - 1010), and cancer cells (1011-1111). These two sub-populations interact at each birth-death time-step with An evolutionary model of tumor cell kinetics 43 tness dened in terms of the prisoner's dilemma payo matrix. The algorithmic details are shown in the appendix Figure 3.10. To set the stage for more complex simulations, Figure 3.1 shows the result of a stochastic simulation (depicting i) driven by the Moran process alone, with no mutations, and no selection. Figure 3.1 shows three dierent simulations, one leading to the elimination of all cancer cells via random drift (red), another uctuates between a mixed cell population after 10,000 cell divisions (yellow), and a third leading to xation of the cancer cell population (blue) after around 5000 cell divisions. The average of 25 stochastic simulations is also plotted (note that the average will converge to half cancer cells and half healthy cells by the law of large numbers). The mean time to xation of the cancer cell population which starts with `i' cells in this simple setting (no mutations, no selection) is given by k =N 2 4 i X j=1 Ni Nj + N1 X j=i+1 i j 3 5 : (3.1) With no mechanism for natural selection, there is no shape to the growth curves. 3.3.2 The prisoner's dilemma payo matrix To introduce the eect of selection which will regulate cell birth and death rates, we use the prisoner's dilemma evolutionary game in which two players compete against each other for the best payo. Each has to decide whether to cooperate (healthy cell) or defect (cancer cell) and each receives a payo determined from the prisoner's dilemma payo matrix 2 , A: 2 What denes a prisoner's dilemma matrix are the inequalities c > a > d > b. The chosen values in (3.2) are relatively standard, but not unique. More discussion of why the prisoner's dilemma matrix, which models the evolution of defection, is a useful paradigm for cancer can be found in [45] and some of the references therein. An evolutionary model of tumor cell kinetics 44 Figure 3.1: Stochastic Moran birth-death process | Cancer cell population, i(t), during three stochastic simulations of the Moran birth-death process in a population of 100 cells and an initial condition of i = 50 cells. The blue curve leads to xation of the cancer cell population, the red curve leads to elimination of the cancer cell population, and the yellow curve remains uctuating in a mixed population of cells after 10,000 cell divisions. An average of 25 stochastic simulations (black dashed line) is also plotted. A = 0 B B @ a b c d 1 C C A = 0 B B @ 3 0 5 1 1 C C A : (3.2) The essence of the prisoner's dilemma game is the two players compete against each other, and each has to decide what best strategy to adopt in order to maximize their payo. This 2 x 2 matrix assigns the payo (e.g. reward) to each player on each interaction. My options, as a strategy or, equivalently, as a cell type, are listed along the rows, with row 1 associated with my possible choice to cooperate, or equivalently my cell type being healthy, and row 2 associated with my possible choice to defect, or equivalently my cell type being cancerous. Your options are listed down the columns, with column 1 associated An evolutionary model of tumor cell kinetics 45 with your choice to cooperate (or you being a healthy cell), and column 2 associated with your choice to defect (or you being a cancer cell). The analysis of a rational player in a prisoner's dilemma game runs as follows. I do not know what strategy you will choose, but suppose you choose to cooperate (column 1). In that case, I am better o defecting (row 2) since I receive a payo of 5 instead of 3 (if I also cooperate). Suppose instead you choose to defect (column 2). In that case, I am also better o defecting (row 2) since I receive a payo of 1 instead of 0 (if I were to have cooperated). Therefore, no matter what you choose, I am better o (from a pure payo point of view) if I defect. What makes this game such a useful paradigm for strategic interactions ranging from economics, political science, biology [45], and even psychology [15] is the following additional observation. You will analyze the game in exactly the same way I did (just switch the roles of me and you in the previous rational analysis), so you will also decide to defect no matter what I do. The upshot if we both defect is that we will each receive a payo of 1, instead of each receiving a payo of 3 if we had both chosen to cooperate. The defect-defect combination is a Nash equilibrium [22], and yet it is sub-optimal for both players and for the system as a whole. Rational thought rules out the cooperate-cooperate combination which would be better for each player (3 points each) and for both players combined (6 points). In fact, the Nash equilibrium strategy of defect-defect is the worst possible system wide choice, yielding a total payo of 2 points, compared to the cooperate-defect or defect-cooperate combination, which yields a total payo of 5 points, or the best system-wide strategy of cooperate-cooperate yielding a total payo of 6 points. The game becomes even more interesting if it is played repeatedly [14, 15, 18, 22], with each player allowed to decide what strategy to use on each interaction so as to accumulate An evolutionary model of tumor cell kinetics 46 a higher payo than the competition over a sequence of N games. In order to analyze this kind of an evolving set-up, a tness function must be introduced based on the payo matrix A. Let us now switch our terminology so that the relevance to tumor cell kinetics becomes clear. In this case, we randomly select pairs of cells out of the total population at each step, and subject them to a birth-death process, basing our birth rates and death rates on the prisoner's dilemma payo matrix. Thus, in our context, it is not the strategies that evolve, as cells cannot change type based on strategy (only based on mutations), but the prevalence of each cell type in the population is evolving, with the winner identied as the sub-type that rst reaches xation in the population. As the populations evolve, the tness of the two competing sub-populations can be tracked, as well as the overall tness of the combined total population of cells. 3.3.3 The tness landscape Let us start by laying out the various probabilities of pairs of cells interacting and clearly dening payos when there arei cancer cells, andNi healthy cells in the population. The probability that a healthy cell interacts with another healthy cell is given by (Ni1)=(N 1), whereas the probability that a healthy cell interacts with a cancer cell is i=(N 1). The probability that a cancer cell interacts with a healthy cell is (Ni)=(N 1), whereas the probability that a cancer cell interacts with another cancer cell is (i 1)=(N 1). The payos associated with the healthy cells and cancer cells, obtained by weighting the payo matrix values with appropriate probabilities, are given by (following notation in [116]): An evolutionary model of tumor cell kinetics 47 H = 3(Ni 1) + 0i N 1 ; (3.3) C = 5(Ni) + 1(i 1) N 1 : (3.4) This gives rise to the average payo associated with the population of cells: hi = H (Ni) + C (i) N : (3.5) Based on these formulas, we dene the tness of the healthy cells as: f H = 1w H +w H H ; (3.6) and the tness of the cancer cells as: f C = 1w C +w C C : (3.7) Here, (w H ;w C ) are `selection strength' parameters, 0w H 1; 0w C 1, that measure the strength of selection pressure on each of the population of cells. If w H = 0, there is no natural selection acting on the healthy cell population and the dynamics is driven purely by the Moran process. Whenw H = 1, the selection pressure on the healthy cell population is strongest and the prisoner's dilemma payo matrix has maximum eect. Likewise for the An evolutionary model of tumor cell kinetics 48 parameter w C and how it controls selection pressure in the cancer cell population. Since therapy imposes selection pressure on dierent sub-populations of cells, w H and w C are the two parameters we alter to administer simulated therapeutic responses. We discuss this in sectionx3.5. The expected tness of each of the sub-populations are: H = Ni N f H ; (3.8) C = i N f C ; (3.9) with total expected tness: = H i + C : (3.10) From these formulas, we can dene the transition probability of going fromi toi+1 cancer cells on a given step: P i;i+1 = if C if C + (Ni)f H Ni N : (3.11) The rst term represents that probability that a cancer cell is selected for reproduction (weighted by tness), and a healthy cell is selected for death. Likewise, the transition probability of going from i to i 1 cancer cells on a given step is: An evolutionary model of tumor cell kinetics 49 P i;i1 = (Ni)f H if C + (Ni)f H i N : (3.12) Here, the rst term is the probability healthy cell is selected for reproduction (weighted by tness), and a cancer cell is selected for death. The remaining transition probabilities are as follows: P i;i = 1P i;i+1 P i;i1 ; P 0;0 = 1; P N;N = 1: (3.13) It is these simple formulas that drive the subsequent dynamics of the competing populations of cells and determine the emergent features of the forming tumor (cancer cell population). A typical set of simulations of the evolving tness of the healthy cell population, H , the cancer cell population C , and the total tness, , is shown in Figure 3.2 as the selection parameter varies from 0 to 1 (w H = w C w). As the population evolves, the tness of the healthy cell population decreases, the tness of the cancer cell population increases (sometimes reaching a maximum point), while the total population tness decreases. 3.3.4 Passenger and driver mutations Two remaining parameters in our model are the passenger mutation rate, m p and the driver mutation rate, m d [93]. Passenger mutations confer no tness advantage, hence m p controls point mutations that act on the digit strings that dene the 11 levels of healthy cells 0000-1010, and the 5 levels of cancer cells 1011-1111. A mutation diagram is shown in Figure 3.3 depicting all of the possible point mutation transitions at each step. Mutations An evolutionary model of tumor cell kinetics 50 (a) (b) (c) Figure 3.2: Fitness as a function of the selection parameter ww H w C | (a) Monotonically decreasing tness of healthy cell sub-population H ; (b) Fitness of cancer cell sub-population C . Note that C has a maximum ati = N 2 + (N1) 8w , which is between 0 andN forw> 1 4 (1 1=N); (c) Monotonically decreasing tness of the total population, . that stay within either of those two ranges do not alter the cell tness. On the other hand, the driver mutation parameter controls mutations that take a binary string from a healthy cell and, via a point mutation, alter it so that the string becomes a cancer cell 3 . A simple example would be a mutation that alters 1010 (healthy) to 1011 (cancer) by stochastically ipping the rst digit from 0 to 1. The interested reader can consult the ow diagram in Figure 3.10 of the Appendix for more details of the algorithm. The full code is available from the authors upon request. 3.4 Results Gompertzian growth arising from multicellular systems occurs in many settings with dif- ferent physical and biological constraints acting in concert. Hence it appears as if this 3 In our simulations, we assume that driver mutations cannot revert to passenger mutations, i.e. once a cancer cell is born, it stays in that category. We do not know of any evidence in the literature that shows the reversion of a cancer cell to a healthy cell, nor is this particularly a focus of this manuscript. An evolutionary model of tumor cell kinetics 51 Figure 3.3: Markov Point Mutation Diagram | Left: diagram shows 16 genetic cell types based on 4-digit binary string and the eect of a point mutation on each cell type. Blue indicates healthy cell type (0000 | 1010), red indicates cancerous cell type (1011 | 1111). Black arrows indicate passenger mutations (healthy to healthy or cancer to cancer), red arrows indicate driver mutations (healthy to cancer). Top right: 3 scenarios may occur during the reproduction process: no mutation, passenger mutation, or driver mutation. universal growth curve does not depend on specic physical mechanisms (e.g. oxygen dif- fusion, blood supply, tumor microenvironment, etc.) but more on multi-cellularity and the ability for populations of cells to assume a heterogeneous distribution of functional states, as was described most clearly in Kendal's 1985 paper [85] and documented clinically in breast [113] and other tumor types. Alternative bio-mechanistic models of tumor growth at the cellular level have been developed (see [114, 117{121]) although do not generally include molecular information or evolutionary eects. Features of the Gompertzian growth curve dened by eqns (3.14), (3.15) allow us to clearly describe three distinct growth regimes, the earliest being subclinical and the most critical regime in which to in uence future tumor kinetics, the second being the clinical regime where growth measurements An evolutionary model of tumor cell kinetics 52 are typically obtained [95], and the third being the lethal burden phase where growth sat- urates. The growth equation, (3.22), relates tumor heterogeneity to growth rates, and we quantify heterogeneity via the Shannon entropy [37, 87] of the cellular population. One of the main features of our evolutionary simulations is to show how it (i) leads to Gompertzian growth, (ii) how growth is driven by heterogeneity quantitated via Shannon entropy, (iii) how the initiation of heterogeneity and tness can be tracked via dynamic phylogenetic trees, and (iv) how tumor kinetics can be in uenced via therapeutic strategies that target heterogeneity best in earlier growth regimes. In keeping consistent with the notation of the Gompeterzian growth curve, we now represent the tumor growth as the proportion of cancer cells in the population, n G (t). 3.4.1 Gompertzian tumor growth and three growth regimes The basic (top-down) equations giving rise to pure Gompertzian growth [122{125] are the coupled equations: dn G dt = n G ; (3.14) d dt = : (3.15) Here, is the proportion of growing cancer cells in the mixed population, which grows expo- nentially according to (3.14), but with a time-dependent growth rate which is exponentially decaying according to (3.15). It is straightforward to integrate (3.14) to obtain: An evolutionary model of tumor cell kinetics 53 n G (t) =N 0 exp 1 t Z t 0 dt t : (3.16) Then, (3.15) is solved with: (t) = 0 exp(t): (3.17) Plugging (3.17) into (3.16) and integrating yields the Gompertzian curve: n G (t) =N 0 exp h 0 (1 exp(t)) i ; (3.18) where in the long-time limit , the population saturates to the value n 1 =N 0 exp( 0 =); (3.19) which we normalize to one (without loss of generality). The key features of Gompertzian growth are shown in Figure 3.4. As the cancer cell proportionn G (t) increases (Figure 3.4a), there are three distinct growth regimes dened by the in ection point on the n G growth curve (maximum of _ n G and d 2 n G =dt 2 = 0), and the two in ection points on the growth- velocity curve _ n G (maximum/minimum of n G and d 3 n G =dt 3 = 0). As shown in Figure 3.1(a), there are three points that divide the growth curve into four distinct regions. For convenience, and symmetry, we lump the second and third regions together and dene three basic growth regimes: An evolutionary model of tumor cell kinetics 54 Regime 1 (Subclinical): Increasing velocity _ n G , increasing acceleration d 2 n G =dt 2 . Cell population and tumor volume grows at an exponential rate; Regime 2 (Clinical): In this regime, _ n G reaches its maximum value. In the early part of the regime, _ n G is increasing while d 2 n G =dt 2 decreases. In the later part of the regime, _ n G is decreasing andd 2 n G =dt 2 becomes negative (deceleration). Growth rates are clinically typically measured as linear; Regime 3 (Saturation/Lethal): Decreasing tumor velocity _ n G with decreasing de- celeration. Growth rate rapidly slows towards full saturation of the cancer cell pop- ulation. Regime 1, generally speaking, is the subclinical growth regime where the developing tumor has substantially fewer than 10 8 malignant cells with a tumor size smaller than 1 or 2 mm 3 . Typically, the clinically measurable regime is Regime 2, while the lethal stage when the tumor saturates is associated with Regime 3. In reality, the boundaries of these regimes are, of course, not sharp and depend on tumor type and location which in uence detectability. But the clarity of the pure Gompertzian curve gives a useful framework which delineates the three distinct growth regimes based on clear principles associated with growth, velocity, and acceleration. The growth rate curve is shown in Figure 3.4c, with its derivatives shown in Figure 3.4d. It is most instructive to show the average growth rates dened in each of the three regimes, also shown in the Figure 3.4c. The average growth rate in the time interval from t 1 to t 2 is dened as: ave = 1 t 2 t 1 Z t 2 t 1 (t)dt (3.20) An evolutionary model of tumor cell kinetics 55 (a) (b) (c) (d) Figure3.4: Gompertzian equation | Numerical simulation of the Gompertzian equa- tion (3.14), (3.15) with parameters N 0 = 0:001, 0 = 10, and = 0:2895. The three regimes of tumor growth are demarcated by the blue dots in each subgure, representing the maximum and minimum of the second-derivative. (a) Cancer cell proportion, n(t), over time; (b) First- and second-derivatives of the tumor growth curve; (c) Growth rate, (t), over time, with the average growth rate in regimes 1, 2, 3 plotted in red; (d) First derivative of growth rate. The subclinical regime 1 has the highest average growth, whereas regime 2, where tumor growth is typically measured, average growth rates are lower, followed by the lowest average growth in the clinically lethal regime 3. This implies that clinically measured growth rates typically underestimate growth rates that preceded it in the subclinical stage. It also implies that linear extrapolation back from clinically measured growth rates to estimate tumor initiation times (see [4, 95, 96]) will systematically overestimate the amount of time the tumor has been developing before being measured. While this might generally be seen as good news (since the cancer initiation event was more recent than estimated via linear An evolutionary model of tumor cell kinetics 56 extrapolation), it also gives the clinician a shorter window of time in which to act. 3.4.2 Heterogeneity and growth via statistical mechanics Kendal [85] lays out a clear argument of how this growth curve arises from a purely statis- tical mechanics point of view. In a nutshell, his argument can be explained by considering a population of n cells, let the jth cell (j = 1; 2; 3;:::;n) have the potential to assume one ofq j possible states. The number of combinations of states possible within the population, P , can be thought of as a measure of intra-neoplastic diversity: P =q 1 q 2 q 3 :::q n ; (3.21) and is related to the growth rate of a tumor via the equation: dn dt = logP; (3.22) where n(t) is the number of cells capable of proliferation at a given time t and is a parameter that sets the timescale of growth 4 . There are two basic cases to consider. First, suppose the cells have no interaction at all, say in the earliest stages of tumor development, 4 Kendal's formulation [85] assumes a cell population made up of three sub-groups: (1) proliferative cells; (2) nonproliferative and nonclonogentic cells; (3) nonproliferative but clonogenic cells, with an assumption that the neoplasm's growth rate is in uenced by the proportion of proliferating to nonproliferating cells and an expression of each clone's growth potential. The log is chosen based on the fact that heterogeneity is measured as the multiplicative combination of achievable states in the tumor, and the requirement that G(P1 P2) = G(P1) +G(P2) for any two sub-populations P1, P2 and growth function G. The discussion of the relationship between tumor heterogeneity and growth is an ongoing topic in the current literature [48, 101, 102, 104, 106]. An evolutionary model of tumor cell kinetics 57 and let each of the n cells have the ability to assume one of m possible states. Then, P =m n , and the growth equation becomes dn dt =n logm = ( logm)n: (3.23) The solution to this equation is the exponentially growing population: n(t) =N 0 exp(( logm)t): (3.24) Thus, early stage development is characterized by exponential growth (regime 1), with a growth rate proportional to the log of the number of assumable states of the cells comprising the tumor population. This stage is characterized by the Gompertzian curve shown in Figure 3.4a to the left of the rst blue dot, in regime 1. Contrast this with later stages of tumor growth, when the sub-populations of cells communicate and in uence each other's growth characteristics, either via competition, or cooperation (regime 3) within the tumor microenvironment. In eect, this will constrain (reduce) the number of assumable states of each cell, since the population is eectively coupled. In the extreme, suppose P = m n =n n . In other words, supposeP is now inversely related to the total number of possible intercellular interactions. Inserting this into (3.22) yields dn dt = log ( m n ) n =n [logm logn]: (3.25) An evolutionary model of tumor cell kinetics 58 m d m p t Emax t SAT n d n p ave;1 ave;2 ave;3 0.4 0.1 5.50e+5 1.830e+6 1.289e+4 4.68e+4 3.14e-5 3.68e-6 1.448e-7 0.3 0.2 4.88e+5 1.753e+6 1.682e+4 8.26e+4 4.04e-5 4.31e-6 1.677e-7 0.2 0.3 4.85e+5 1.761e+6 1.715e+4 1.230e+4 3.86e-5 4.41e-6 1.729e-7 0.1 0.4 5.40e+5 1.426e+6 1.362e+4 1.836e+4 3.04e-5 3.81e-6 1.658e-7 Table 3.1: m d : driver mutation rate; m p : passenger mutation rate; t Emax : time to maximum entropy;t SAT : time to saturation;n d : number of driver mutations;n p : number of passenger mutations; ave;1 : average growth rate in regime 1; ave;2 : average growth rate in regime 2; ave;3 : average growth rate in regime 3. The solution to this equation is exactly the Gompertzian growth curve (3.18) and accounts for regimes 2 and 3 previously discussed in which tumor growth slows down. The growth equation (3.22) which relates cancer cell population growth to tumor heterogeneity is capa- ble of producing a family of growth curves, depending on details of intercellular coupling, which of course is in uenced by details of the biological and physical constraints in u- encing the tumor microenvironment. Thus, the growth equation (3.22) has the ability to produce dierent detailed shapes based on assumptions associated with intercellular cou- pling. Table 3.1 shows the average growth rates in the three regimes as a function of the key parameters in the model. 3.4.3 Quantitative measures of tumor heterogeneity and growth For our purposes, we measure heterogeneity using the Shannon entropy from information theory [37]: E(t) = N X i=1 p i log 2 p i ; (3.26) An evolutionary model of tumor cell kinetics 59 (here, log is dened as base 2). The probability p i measures the proportion of cells of type i, with i = 1;:::; 16 representing the distribution of binary strings ranging from 0000 to 1111. We then course-grain this distribution further so that cells having strings ranging from 0000 up to 1010 are called \healthy", while those ranging from 1011 to 1111 are \cancerous" 5 . The growth equation (3.22) then becomes dn E dt =E(t): (3.27) It follows from (3.27) that the cancer cell proportion n E (t) can be written in terms of entropy as: n E (t) = Z t 0 E(t)dt: (3.28) The panel in Figure 3.5 shows the results from our cell-based simulations. Figure 3.5a shows the Gompertzian curve associated with the proportion of cancer cells in the population, while Figure 3.5b shows the velocity and accelerations associated with growth, and can be compared with Figure 3.4b. In Figure 3.5c we show the entropy during a typical simulation, marking the maximum entropy point which peaks relatively early in the simulation before the entropy returns back down to zero, re ecting the fact that cancer cells have reached xation and have saturated the population. Figure 3.5d shows the tness of the cancer cell sub-population, healthy cell sub-population, and the overall tumor tness (w H = w C w = 0:5). As a typical simulation proceeds, the cancer cell sub-population tness increases, 5 Our results are relatively insensitive to where we draw the dividing line between healthy and cancerous. An evolutionary model of tumor cell kinetics 60 the healthy cell sub-population tness decreases, while the overall tumor tness decreases. Figure 3.5e, 3.5f shows the Gompertzian growth curves as the selection pressure increases (Figure 3.5e) and as the mutation rate increases (Figure 3.5f). High values for either of these parameters leads to a very steep growth curve, as is expected. Figure 3.6 shows the growth curves linearly extrapolated back to give a prediction of when the rst driver mutation occurred that initiated tumor growth. The growth rates from regime 2 (linear regime) are used to extrapolate back to the initiation event. Since the actual growth rate in regime 1 is much faster than linear, the linear extrapolation extends the event too far back in time as compared to when the event actually occurred. The inset of Figure 3.6 shows histograms of the average growth rates in each of the three regimes as a function of the mutation rate m (here, we take m p =m d =m). A typical stochastic simulation showing the evolution of all 16 possible cell types is shown in Figure 3.7. We also show E(t), where entropy is computed using the most extreme course-grained two-state system comprised of the two sub-populations of healthy cells and cancer cells. We compare in Figure 3.7 the Gompertzian growth curve (eqn. (3.18)) and the corresponding curve obtained from eqn. (3.28) to the stochastic simulation and the agreement is excellent. Likewise, we also show a comparison ofdn=dt with eqn. (3.27) and eqn. (3.14) with E(t) normalized so that limiting values match the stochastic simulation, and the agreement is also excellent. In the beginning, entropy is zero, since the population consists purely of healthy cells, and in the end of the simulation, entropy is again zero as the population consists purely of cancer cells. Entropy peaks somewhere early in the simulation when the mixture of cell types is equally distributed over cancer and healthy An evolutionary model of tumor cell kinetics 61 (a) (b) (c) (d) (e) (f) Figure 3.5: Moran birth-death process with selection | (a) Cancer cell popu- lation, i(t) (w = 0:5, m = 0:2, N = 10 10 ) plotted with a spline curve connecting 200 data points from a single stochastic simulation; (b) First- and second-derivatives of the tumor growth curve in (a) are plotted with maximum and minimum of second-derivative indicated (blue); (c) Entropy of the cell population from eqn. (3.26) as it relates to the growth equation (3.27); (d) Fitness of healthy cell population and cancer cell population and total tness as dened by eqns. (3.8), (3.9), (3.10); (e) Simulations of cancer cell population, i(t), for a range of selection parameter values; (f) Simulations of cancer cell population, i(t), for a range of mutation rate values. types. It is this intermediate but important heterogeneously distributed state that is the key driver of growth, as is clear from eqn. (3.27). An evolutionary model of tumor cell kinetics 62 Figure 3.6: Tumor initiation prediction | Five sample stochastic simulations of tumor growth (N = 10 10 cells, w = 0:5, m = 0:1; 0:2; 0:3; 0:4; 0:5) plotted on a log-linear graph where the model output (i(t), solid lines) is t in the clinical regime (greater than 10 8 cells) using an exponential growth equation and extrapolated backwards in simulation time (dashed lines). The inset bar graph shows the average growth rate in each regime. 3.4.4 Dynamic phylogenetic trees and evolution of tness To track the initiation of cellular heterogeneity from an initially homogeneous state, we follow all of the mutations that take place during the course of a simulation, and orga- nize this in the form of a phylogenetic tree in Figure 3.8 showing the typical size of the genotypic space and the evolution of the genotypic landscape. As the simulation proceeds, the phylogenetic tree dynamically branches out into an increasingly complex structure, with tness characteristics color coded in Figure 3.8a. We also show the bins associated with each of the 16 cell types, the number of cancer cells i(t), and the entropy associated with the sub-population of cell types as a simulation proceeds, in Figure 3.8b. Knowing exactly the types of cells comprising the tumor at any given time allows us to target cell distributions for simulated therapies to test dierent strategies, which we describe next. An evolutionary model of tumor cell kinetics 63 (a) (b) (c) Figure 3.7: Comparison of stochastic Moran birth-death process, Gom- pertzian, and Shannon entropy growth curves | (a) A single stochastic simulation (N = 10 10 cells, m = 0:5, w = 0:5, m p =m d = 0:25) growth curve, n(t), compared with the Gompertzian growth curve, n G (t), eqn. (3.18), and Shannon entropy growth curve, n E (t), eqn. (3.28). Growth curves n G (t) and n E (t) are normalized to equal one in the limit; (b) Comparison of rst-derivatives of n(t), n E (t), n G (t); (c) Comparison of growth rates associated withn(t),n E (t),n G (t), with average growth rates ofn(t) plotted for each regime, eqn. (3.20). An evolutionary model of tumor cell kinetics 64 (a) (b) Figure 3.8: Emergence of genetic heterogeneity | (a) Left: sample dendritic phylogenetic tree tracking point mutations as time extends radially. Right: three snapshots in time of a dendritic tree in a simulation of 30 cells with strong selection (w = 1,m p = 0:1, m d = 0:2). Pathways are color coded to indicate genetic cell type; (b) Linear phylogenetic tree of the same stochastic simulation shown in (a) along with histogram plots of the distribution of genetic cell types and a plot of the cancer cell population i(t) and entropy. An evolutionary model of tumor cell kinetics 65 3.4.5 A comparison of early vs. late therapy In Figure 3.9 we show the results from asking the simple question of how early therapy (administered in regime 1) compares with therapy in the middle stages of tumor devel- opment (regime 2), or in the later stages of development (regime 3). Eqns (3.11), (3.12) are the governing equations controlling birth/death rates of the cancer cell, healthy cell sub-populations as natural selection plays out. Since the proliferation of cancer cells can be thought of as an imbalance of selection pressures on the competing sub-populations in favor of the cancer sub-population, the goal of any therapeutic intervention is to alter this complex imbalance in favor of the healthy cell sub-population. We implement this by adjusting the selection pressure parameters (w H ;w C ) in the formulas (3.6), (3.7). In particular, when therapy is `on', we choosew C = 0, andw H = 1, tilting the selection pres- sure in favor of the healthy cell sub-population. When therapy is `o', the two parameters return to their original baseline values, which here we take as w H = 0:1,w C = 0:1. Figure 3.9 depicts the proportion of cancer cells in the population both in the absence of therapy, and when therapy is administered. As a comparative tool, in each case, we administer the therapy until a xed number of cancer cells remains (in each case, we take this threshold number to be 25 cancer cells), and we compare the amount of time, t, it takes to achieve this low level. The gure clearly shows t 1 < t 2 < t 3 < t 4 , while if therapy is admin- istered too late, as in t 5 , the low threshold is never achieved. The simulations show that a shorter therapeutic time-period is needed if administered earlier to gain the same level of success. The topic of how best to optimize computational therapies is complex, and these simulations are only meant as a conrmation and quantication of how early stage therapy is more eective than late stage therapy. An evolutionary model of tumor cell kinetics 66 Figure 3.9: Simulated therapy | An average of 25 stochastic simulations (N = 10 3 cells, w = 0:5, m = 0:1) where therapy (w H = 1, w C = 0) is administered at dierent time points (t = 6000; 8000; 10000; 12000; 14000 cell divisions) until all cancer cells are eliminated below a small threshold value (25 cells). Time required (t) for tumor elimination increases as the tumor volume increases (i.e. t 1 < t 2 < t 3 < t 4 , blue, red, yellow, purple arrows respectively), until, at later simulation time points, therapy is unable to regress tumor size (t 5 , green arrow). 3.5 Discussion To summarize the main points forming the framework of our model: (i) A tumor is a complex Darwinian ecosystem of competing cells operating on an adaptive tness landscape driven by mutational dynamics and shaped by evolutionary pressures; (ii) The basic competitors in an evolutionary game theory model of tumor development are cell populations with a broad distribution of tness characteristics course grained into two types: healthy cells (cooperators) and cancer cells (defectors). Each of these cell sub-populations attempts to maximize its own tness; An evolutionary model of tumor cell kinetics 67 (iii) Cell tness is associated with reproductive prowess and in this respect, healthy cells are less t than cancer cells; (iv) Primary tumors initiate from a single malignant cell that has undergone the ap- propriate mutational steps and subsequently undergoes clonal and sub-clonal expansion. Polyclonality and heterogeneity are thus seen as emergent features of tumor development; (v) Parameters and distributions measured in the detectable range of tumor growth, such as tumor growth rates and xation probabilities, are emergent features that have developed from a monoclonal state via cell kinetics and evolutionary development taking place in the subclinical regime; (vi) Tumor growth is driven by molecular heterogeneity of the cell population comprising the tumor and re ected in the growth equation (3.22); (vii) Tumor cell populations are more amenable to therapeutic strategies in the early stages of development, before selection for growth and survival have shaped the environment. We believe the simple evolutionary model described in this paper, driven by a Moran process and shaped by heritable mutations with a tness landscape based on the prisoner's dilemma evolutionary game is useful in helping to understand early stage tumor growth and how it is in uenced by the interplay of a few select small number of key parameters. When a malignant tumor cell population has already exceededO(10 8 10 10 ) cells, some of which may have entered the circulation or lymphatics and migrated to other sites, the opportunity to control or even shape future events may be limited. Attacking tumor heterogeneity as soon as it develops seems to be a useful strategy, particularly if heterogeneity is the driver of growth, as in eqn. (3.22). Whether these concepts can be developed in the more general An evolutionary model of tumor cell kinetics 68 context when cell dissemination to other sites is included in the model, and then translated into actionable clinical strategies is a challenge for the future. An evolutionary model of tumor cell kinetics 69 3.6 Appendix Figure 3.10: A ow chart of the Moran process with selection and mutation algorithm | Box 1: mutation rate m (where m = m p +m d ), selection pressure w and the initial state vector x containing N total cells are the inputs for a simulation. Box 2: the prisoner's dilemma game (a = 3; b = 0; c = 5; d = 1) is used to calculate the tness of each healthy and cancer cell type, which is a function of the payo values and the state vector, x. Box 3, 4: a single cell is chosen for death according to the relative proportion of the cell type in the cell population. Simultaneously, a single cell is selected for birth according to the relative proportion, weighted by cell tness. Box 5: During the replication process, the daughter cell inherits a replica of the parent cell's genetic string, with errors occurring at a rate ofm. A single bit of the daughter cell's genetic string may ip during each cell division. The possible mutations can be thought of as a single step random walk on the Markov diagram shown in Figure 3.3. Chapter 4 Chemotherapeutic dose scheduling based on tumor growth rates: the case for low dose metronomic high entropy therapies Jerey West, Paul K. Newton Chemotherapeutic dose scheduling based on tumor growth rates: the case for low dose metronomic high entropy therapies [40] USC Preprint 70 Chemotherapeutic dose scheduling 71 4.1 Abstract We use a stochastic Moran process model of tumor cell kinetics, coupled with a prisoner's dilemma game-theoretic cell-cell interaction model to design chemotherapeutic strategies tailored to dierent tumor growth characteristics. Using the Shannon entropy as a novel tool to quantify dosing strategies, we contrast maximum tolerated dose (MTD) strategies as compared with low dose, high density metronomic strategies (LDM) for tumors with dierent growth rates. Our results show that LDM strategies can outperform MTD strategies. The advantage is magnied for fast growing tumors that thrive on long periods of unhindered growth without chemotherapy drugs present and is not evident after a single cycle of chemotherapy, but grows after each subsequent cycle of repeated chemotherapy. The model supports the concept of designing dierent chemotherapeutic schedules for tumors with dierent growth rates and devel- ops quantitative tools to optimize these schedules for maintaining low volume tumors. Major Findings Model simulations show that metronomic (low dose, high density) therapies can outperform maximum tolerated dose (high dose, low density) therapies. This is due to the fact that tumor cell reduction is more sensitive to changes in dose density than changes in dose concentration, especially for faster growing tumors. This eect is negligible after a single cycle of chemotherapy, but magnied after many cycles. The model also allows for novel Chemotherapeutic dose scheduling 72 chemotherapeutic schedules and quanties their performance according to tumor growth rate. Quick Guide to Equations and Assumptions Assumptions of the model: 1. The model is a computational one, driven by a stochastic Moran (birth-death) process with a nite cell population, in which birth-death rates are functions of cell tness. 2. Two classes of cells (healthy, cancer) compete against each other at each birth-death event, with tness (f H ,f C ) calculated according to the payo matrix associated with the prisoner's dilemma evolutionary zero-sum game. 3. Chemotherapy is a selective agent that alters the selection pressure (w H ,w C ) of each cell population, with two parameters: dose concentration (c) and dose density (d). Key equations: In a Moran nite-population, birth-death process there are i cancer cells in a population of N total cells (where the number of healthy cells is denoted Ni). At each time step in the stochastic evolutionary population dynamics model, a single cell is chosen for birth and a separate single cell is chosen for death. A tumor grows by increasing the number of cancer cells from i to i + 1 in any given time step. The probability that a healthy cell interacts with another healthy cell is given by (Ni 1)=(N 1), whereas the probability that a healthy cell interacts with a cancer cell is i=(N 1). The probability that a cancer Chemotherapeutic dose scheduling 73 cell interacts with a healthy cell is (Ni)=(N 1), whereas the probability that a cancer cell interacts with another cancer cell is (i 1)=(N 1). These probabilities, known as the Moran process, can be extended to include a tness landscape where natural selection can play out over many cell division timescales. The probabilities outlined above are weighted by the \payo" in order to determine the tness function for each subpopulation: healthy (f H ) and cancer (f C ), below. The payo values (a, b, c, d) are associated with the prisoner's dilemma evolutionary game [14, 75]. The prisoner's dilemma is dened by the payo inequalities such that c>a>d>b, but here we assume the relatively standard (but not unique) values of a = 3,b = 0,c = 5, and d = 1. f H = 1w H +w H a(Ni 1) +bi N 1 (4.1) f C = 1w C +w C c(Ni) +d(i 1) N 1 (4.2) Here, (w H , w C ) are \selection strength" parameters, 0 w H 1, 0 w C 1, that measure the strength of selection pressure on each of the population of cells. If w H = 0, there is no natural selection acting on the healthy cell population, and the dynamics is driven purely by the Moran process (i.e. random drift). When w H = 1, the selection pressure on the healthy cell population is strongest, and the prisoner's dilemma payo matrix has maximum eect. Chemotherapeutic dose scheduling 74 From these formulas, we can dene the transition probability of going fromi toi+1 cancer cells on a given step (P i;i+1 ) or from i to i 1 on a given step (P i;i1 ). P i;i+1 = if C if C + (Ni)f H Ni N (4.3) P i;i1 = (Ni)f H if C + (Ni)f H i N (4.4) The rst term in each each equation represents the probability that a cell is selected for reproduction (weighted by tness). The second term represents the probability that a cell is selected for death. The probability of the number of cancer cells remaining the same (P i;i ) is given by the following. There are two absorbing states (P 0;0 , P N;N ). P i;i = 1P i;i+1 P i;i1 ; P 0;0 = 1; P N;N = 1 (4.5) 4.2 Introduction Low dose metronomic chemotherapy (LDM) is the systematic and frequent delivery of chemotherapeutic agents at doses lower than the maximum tolerated dose paradigm (MTD) [38, 39]. It is typically given at a low dose between 1/10th and 1/3rd of the maximum tolerated dose, without a long period of time between subsequent doses, hence it is also associated with higher dose densities [38]. Important features of low dose, high density metronomic chemotherapy include: regular administration of chemotherapy without any interruptions using an optimized dose; preference for oral drugs; low incidence of side Chemotherapeutic dose scheduling 75 eects; low risk of developing resistance; lower cost. In addition, some elderly or frail patients may only be suited for lower dose chemotherapy. Residual toxicity from pre- vious treatment may also reduce consideration for MTD chemotherapy [39]. Metronomic chemotherapy regimens have been associated with lower cost of inexpensive oral drugs such as cyclophosphamide and result in fewer side-eect associated expenditures. Several phase II studies have shown promises of metronomic-like chemotherapy and its excellent safety proles [39]. The lower doses of metronomic chemotherapy regimens are now thought to not only reduce the harmful side eects of toxicity delivered to the patient but perhaps also improve anti-tumor eects [126], by killing endothelial cells in addition to its cyto- toxic eect on cancer cells [127, 128] in an uninterrupted schedule for prolonged treatment periods. Metronomic chemotherapy has been shown to be eective in preclinical trials where cancer cells have developed resistance to the same chemotherapeutics [126]. These LDM regimens are also suited to combination or additive strategies to new targeted and relatively non-toxic anticancer drugs recently developed. While the advantages of LDM chemotherapy may be wide ranging with respect to toxicity, resistance, and anti-angiogenic eects (outside of the scope of our model), the goal of this article is to use an evolutionary mathematical model of cell/tumor growth with the ability to simulate chemotherapeutic scheduling to identify growth regimes where LDM would likely outperform MTD, and to test various scheduling protocols altering dose density and concentration. While there is no simple answer to the question of what types of tumors and growth regimes where LDM would be preferable to MTD, our results show that LDM chemotherapies with an adequate dose can outperform MTD, especially for fast growing tumors that thrive on long periods of drug-free rest with unhindered regrowth. This eect is Chemotherapeutic dose scheduling 76 not evident after a single cycle of chemotherapy, but is magnied after each subsequent cycle of repeated chemotherapy. In the interplay of choosing between high dose chemotherapy (MTD) or low dose, high density chemotherapy (LDM), our results show that increasing dose has diminishing returns, so the higher densities aorded by LDM regimens are an ideal tradeo. These results may have remained hidden even in the advent of helpful theoretical regression laws like Skipper's laws and the Norton-Simon hypothesis because these laws rely on instantaneous rates of regression, rather than the net result of the full chemotherapy cycle operating in an environment with variable growth rates. We explain how our results add to the understanding of these classic growth models and advocate the consideration of tumor growth rates when choosing chemotherapy scheduling. 4.2.1 Administration of metronomic chemotherapy A systematic literature review of the MEDLINE, EMBASE, CENTRAL, and PubMed databases for LDM chemotherapy trials from 2000 to April 2012 performed by Lien et. al. in 2013 [39] revealed a wide variety in dose delivered and dose schedules under the terminology of metronomic chemotherapies. From the 80 studies analyzed, 107 unique treatment regimens were found (including regimens where multiple drugs were used metro- nomically). 38 regimens used LDM only (monotherapy n = 24, doublet LDM therapy n = 14). Of the monotherapy, the relative dose intensity (RDI: measured with respect to the maximum tolerated dose) ranged from 0.27 to 1.58 (median 1.02) and dose density (percentage of days drug is delivered) ranged from 32% to 100% [39]. RDI is calculated by dividing the dose intensity (DI; the sum of the doses given each day of the chemother- apy cycle) for a chemotherapy regimen by the baseline DI value of the conventional MTD Chemotherapeutic dose scheduling 77 schedule. A chemotherapy may deliver a greater overall DI than the MTD (i.e. RDI > 1) if a lower dose is delivered more often, achieving a greater total dose over the course of the full chemotherapy cycle. The lower dose reduces toxicity, allowing for more frequent dosing, a key idea behind the metronomic schedules. For a low dose metronomic chemotherapy, any schedule that administers a lower dose at more frequent intervals (higher dose density) could be classied as \low dose metronomic." But, as seen above, in clinical practice the relative dose intensity delivered and the density of the scheduled are varied without clear consensus. In fact, only one monotherapy treat- ment regimen kept the RDI constant, balancing the lower dose with an equivalent increase in dose density. Of the remaining 23 regimens about half increased RDI (n = 12) while half decreased RDI (n = 11). It is evident that many of the quantitative details of LDM chemotherapy are unresolved including patient selection, choice of drug (or combinations of drugs and treatments), and optimal dose and treatment intervals [39]. With this in mind, the goal of this manuscript is very targeted. We wish to quantify the relationship between dose and dose density delivered using the Shannon entropy index [37] as a quan- titative scheduling and dosage tool. We will rst brie y review the prisoner's dilemma evolutionary game theory model of primary tumor growth that we use to carry out our computational trials [32, 129] as well as the notion of Shannon Entropy as an index to compare chemotherapeutic regimens in order to show that high entropy schedules (with an adequate dose intensity) outperform low entropy schedules. Chemotherapeutic dose scheduling 78 4.2.2 The classic tumor regression laws Benzekry et al.[130] chronicle that, despite a rise in personalized and precision medicine, currently chemotherapeutic agents are often still administered in the maximum tolerated dose paradigm. The author predicts that the forthcoming development of metronomic chemotherapy may pave the way for implementing \computational oncology at bedside, be- cause optimizing metronomic regimen should only be achieved thanks for modeling support." This prediction characterizes a growing eld sometimes referred to as computational or mathematical oncology [100, 131]. First, however, in order to properly understand how al- ternative dosing schedules like the metronomic regimens t into the future of chemotherapy scheduling, it is important to remember the reasons that led to the advent and continued use of MTD paradigms. 4.2.2.1 Skipper's Laws The relationship between dose and tumor cytotoxicity is linear-log (i.e. exponential decay) [6]. Skipper et al. [8] were the rst to develop a set of theoretical laws governing the behavior (and imply the design) of chemotherapy schedules in cancer in the late 1970's. Our understanding of the Gompertzian growth of tumors have made the application of these laws more complex, but the fundamentals of these laws still apply today [3]. In a tumor that grows exponentially (eqn. 4.6 and 4.7) with a constant exponential rate, the rst law states that the tumor volume doubling time is constant over the life of the Chemotherapeutic dose scheduling 79 tumor (d t = log(2)=), _ n = n (4.6) n(t) = n 0 exp(t) (4.7) The second of Skipper's laws is that the percentage of cells killed by a given drug dose, D, is constant, therefore a linear increase in dose causes a log increase in cell kill [100]. As an example, a drug dose, x, that shrinks tumor size from 10 6 to 10 5 cells results in a 90% decrease of tumor population. An identical subsequent drug dose (a total dose of 2x) will further reduce tumor population size according to that same kill constant, to 10 4 . A third dose results in 10 3 cells, a fourth, 10 2 , and so on. The kill law is known as a the `log' kill because the constant fraction is a constant logarithmic amount. Skipper's log-kill law indicates that subsequent dosing has a diminishing return; the last few remaining cells are the most dicult to eliminate. This log-linear relationship can be formulated as follows: logP S =D: (4.8) 4.2.2.2 Norton-Simon Hypothesis One important reason the Skipper-Schabel-Wilcox model is so meaningful is that it con- ceptualizes the tumor growth model (e.g. exponential) and tumor regression (log-kill). Norton and Simon realized the importance of extending these observations to a Gom- pertzian growth model (eqn. 4.11). The log-kill law, a fundamentally static law does not Chemotherapeutic dose scheduling 80 say anything about the relationship between the fraction of cells killed and the growth rate of the tumor, only the relationship between the rate of tumor regression and the dose. In eect, Skipper's second law assumes a constant growth rate, and therefore, a constant regression rate. In Gompertzian growth, the non-constant growth rate results in a range of log-kill rates () corresponding to the instantaneous growth rate ( (t)). Gompertzian growth is given by the following coupled ordinary dierential equations. _ n = n (4.9) _ = (4.10) The Gompertz function reduces to the exponential function when = 0. These coupled ordinary dierential equations may be directly solved, as follows. n(t) =n 0 exp h 0 (1exp(t) i (4.11) The Norton-Simon hypothesis states that tumor regression is positively (linearly) correlated with the instantaneous growth rate just before the treatment of the unperturbed tumor [7, 88]. Generally, smaller tumors are associated with higher growth rates (and therefore, higher regression rates). Mathematically, the Norton-Simon Hypothesis can be formulated, _ n =f(n(t))(1L(t)) (4.12) Chemotherapeutic dose scheduling 81 wheren(t) is the growth rate model of the tumor at timet,f(n(t)) is the growth dynamics associated with the unperturbed tumor (i.e. exponential growth or Gompertzian growth), andL(t) is the loss function of cells resulting from treatment. The growth functionf(n(t)) may be assumed to be exponential, (eqn. 4.7) or Gompertzian, (eqn. 4.9 and 4.10). Remembering that Skipper's second law states that cell kill follows rst-order kinetics, we may assume for the time being that L(t) const:, or that the rate of cell removal due to treatment is constant. In other words, each dose of chemotherapy is associated with some value of L. The goal is to nd the optimal dose concentration and dose density that maximize the loss rate of cell kill, L. 4.2.3 The Implications of the Norton-Simon Hypothesis Norton and Simon hypothesized that chemotherapy will only be eective in targeting cells that are in active proliferation (and as such are directly contributing the growth of the tu- mor in equation 4.12). Their model demonstrated ability to t data preclinical experiments [112] and predict future tumor growth and regression after a few initial measurements and data from clinical trials in breast cancer [7]. The model has several key implications. First, the model predicts a higher regression for higher dose delivered. The highest dose tolerable to the patient should be chosen. Second, tumor regrowth during rest periods of chemotherapy necessitates a shorter rest period and subsequently, a shorter time of tumor regrowth. The next round of a dose dense chemotherapy will attack a smaller tumor (with higher growth rate) and lead to higher regression. Both implications give rise to the invention of the MTD paradigm to attack the Chemotherapeutic dose scheduling 82 tumor with the highest dose, coupled with shortest rest. These predictions were conrmed by clinical trials in which chemotherapy schedules were densied from 21 to 14 days [132]. The hypothesis also predicts that tumors with an identical tumor burden may have varied responses. The growth rate of the tumor determines the response to chemotherapy. As such, early administration is important, implying a better response when the tumor is in initial stages of high growth. Similar models using the ratio of tumor volume to the host-in uenced tumor carrying capacity (which corresponds inversely to the instantaneous growth rate of the tumor) has been shown to inversely predict radiotherapy response [133]. Fundamentally however, the Norton-Simon hypothesis provides no predictions for the eect of dose and dose density on regression. The Norton-Simon hypothesis (equation 4.12) conceals the fact that the rate of cell-kill, L(t) will be dependent on two factors: drug concentration and the number of days the drug is administered. The goal of this manuscript is to extend the classical and well-accepted predictions of Norton-Simon hypothesis from instantaneous regression rates (i.e. the derivative) to the cumulative eect (i.e. the integral) over one (or many) cycles of chemotherapy. Chemotherapy \strategies," or schedules are quantied using the Shannon entropy [37] by their total cell reduction (TCR) over the course of the full schedule, rather than the initial regression rate (). The evolutionary model we introduce in this paper is compared with regression data from murine models (see gure 4.3) and shown to be in good agreement. Chemotherapeutic dose scheduling 83 4.3 Materials and methods 4.3.1 Chemotherapeutic agents alter the tness landscape It is now well established that cancer is an evolutionary and ecological process [10, 11]. Studying cancer as a disease of clonal evolution has major implications on tumor progres- sion, prevention and therapy [24, 25]. The evolutionary forces at play inside the tumor such as genetic drift with heritable mutations and natural selection operating on a t- ness landscape are in uenced by tumor microenvironment and the interactions between competing cell types. Increased selection will in uence the rates of proliferation and sur- vival, which cause the population of cells within a tumor to progress toward more invasive, metastatic, therapeutic resistant cell types. The role of chemotheraputic agents is to kill proliferating cancer cells, which can be modelled as altering the evolutionary trajectory of the tumor: the relative tness of the surrounding healthy cell subpopulation is greater than the targeted proliferating cancer cells. In order to model these complex evolutionary forces in cancer, many theoretical biolo- gists have used an evolutionary game theory (EGT) framework, pioneered by Nowak, to study cancer progression (see [14, 15, 17, 19, 134]). Evolutionary game theory provides a quantitative framework for analyzing contests (competition) between various species in a population (via the association of `strategies' with birth/death rates and relative sub-clonal populations) and provides mathematical tools to predict the prevalence of each species over time based on the strategies [17, 64, 66, 72]. More specically, the framework of EGT allows the modeler to track the relative frequencies of competing subpopulations with dierent Chemotherapeutic dose scheduling 84 traits within a bigger population by dening mutual payos among pairs within the group. From this, one can then dene a tness landscape over which the subpopulations evolve. 4.3.2 The model The model presented in [32, 129] and used in this paper to carry out our computational trials is a framework of primary tumor growth used to test the eect of various chemotherapuetic regimens, including MTD and LDM. The model is a stochastic Moran (nite-population birth-death) process [115] that drives tumor growth, with heritable mutations [50] operating over a tness landscape so that natural selection can play out over many cell division timescales (described in more detail in [32, 129]). The birth-death replacement process is based on a tness landscape function dened in terms of stochastic interactions with payos determined by the prisoner's dilemma game [14, 75]. This game incorporates two general classes of cells: healthy (the cooperators) and cancerous (the defectors) [135, 136]. During tumor progression, each cell is binned into one of two tness levels, corresponding to their proliferative potential: healthy (low tness) and cancer (high tness). In our model, we can think of a cancer cell as a formerly cooperating healthy cell that has defected and begins to compete against the surrounding population of healthy cells for resources and reproductive prowess. The model demonstrates several simulated emergent `cancer-like' features: Gompertzian tumor growth driven by heterogeneity [36, 48, 85], the log-kill law which (linearly) relates therapeutic dose density to the (log) probability of cancer cell survival, and the Norton-Simon hypothesis which (linearly) relates tumor regression rates to tumor growth rates, and intratumor molecular heterogeneity as a driver of tumor growth [129]. Chemotherapeutic dose scheduling 85 (a) (b) 0 0.2 0.4 0.6 0.8 1 i/N 0 0.5 1 1.5 2 2.5 3 fitness (c) 0 0.2 0.4 0.6 0.8 1 i/N 0 0.5 1 1.5 2 2.5 3 fitness (d) 0 0.2 0.4 0.6 0.8 1 i/N 0 0.5 1 1.5 2 2.5 3 fitness (e) Figure 4.1: Chemotherapy is a selective agent that alters the tness landscape of cells | (a) The dose strength parameter, c; (0c 1), alters the selection pressure parameter, w; (0 w 1), in favor of the healthy cell population (w H > w) and to the disadvantage of the cancer cell population (w C <w). (b) Total dose density delivered in the one chemotherapeutic cycle, D, is the product of the dose strength (c; 0c 1) and dose interval (d; 0d 1) such thatD =ct (eqn. 4.13 (0D 1). (c,d,e) Plots showing the tness of the healthy cell subpopulation (f H , blue) and the cancer cell subpopulation (f C , red) for no therapy, low dose therapy, and high dose therapy. Others have presented mathematical models to study evolutionary dynamics of tumor response to targeted therapy [113] in either combination or sequential therapy [34, 137], and optimal drug dosing schedules to prevent or delay the emergence of resistance or optimize tumor response [138{140]. We are interested testing \strategies," or drug schedules that control the number of cancer cells, i, in a population of N cells comprising the simulated tissue region. (Note: the size of the tumor population,i, is variable and changing according to the tness landscape, detailed in equations 4.1 through 4.5. The carrying capacity,N, is a parameter in the model, but all plots shown here are normalized byN, so the proportion of cancer cells, i=N, is used to track tumor growth, without loss of generality.) The model Chemotherapeutic dose scheduling 86 presented here uses a parameter, w, to control the eect of selection pressure. A value of w = 0 corresponds to neutral drift (no selection) and a value of w = 1 corresponds to strong selection. We break w into two separate parameters, w H , the selection pressure on the healthy population, andw C , the selection pressure on the cancer population (see gure 4.1a). Each dose of chemotherapy is associated with a dose concentration, c, which alters the selection pressure as indicated in gure 4.1a. Here, we assume drug concentration will be measured as as a fraction of the conventional maximum tolerated dose (MTD) dosages, hence 0c 1 (see Figure 4.1b). As c increases, the selection pressure is altered in favor of the healthy cells (w H > w) and to the detriment of cancer cells (w C < w) as shown in Figure 4.1a. Before therapy, the tness landscape of an untreated tumor is that of a prisoner's dilemma (Figure 4.1c), where the tness of the cancer subpopulation is greater than the healthy population for the entire range of cancer proportion, i=N. The change in tness landscape for a moderate value of c (c = 0:4) is shown in Figure 4.1d, which gives the healthy cell population a tness advantage over the cancer population. The advantage is lessened as the tumor size (i=N) increases (which contributes to the emergence of the Norton-Simon model in Figure 4.2, explained in detail below). For a strong dose of therapy (such as c = 0:8, shown in Figure 4.1e), the eect on the tness landscape is exaggerated. Thus, a higher dose leads to a higher kill rate of cancer cells. In literature, two models have been proposed to model loss functions due to a drug: 1) non-cycle specic (where the loss function is linear with tumor size) [141] and 2) cycle- specic (where loss function is linear with tumor growth rate) [7, 8]. Cycle-specic drugs are considered here, and thus a model of regression that is linear with tumor growth rate is chosen. The loss function of the Norton Simon hypothesis in equation 4.12 shows an Chemotherapeutic dose scheduling 87 (a) (b) Figure 4.2: Classical Tumor Regression Laws | (a) The Norton-Simon hypothesis states that tumor regression is proportional to the growth rate of an unperturbed tumor of that size. Unperturbed tumor growth, n U (t) (black) in a representative population of N = 10 3 cells, and growth rate, (t) (red) is shown. Therapy is administered at various timepoints in the growth of the tumor and then regression, n T (t), is plotted (blue). Rate of regression, , is the best-t slope on the log-plot. (b) The average regression rate was calculated for 25 stochastic simulations, and plotted as a function of at the time of therapy with error bars indicating the standard deviation of values. A linear best t (predicted to be linear by the Norton-Simon hypothesis) is calculated to be (t) = 3:0865 +5:2359e 05. example of cycle-specic drug modeling. The instantaneous growth rate (i.e. )of a stochastic Moran process model (see equations 4.1 through 4.5) is proportional to the selection pressure, w. This indicates that varying w linearly with dose concentration c (shown in gure 4.1a) is comparable to the previous cycle-specic drug models. This linearity is conrmed by plotting instantaneous regression rates throughout the life of a tumor in gure 4.2. Identical, continuous chemotherapy is administered at dierent time points in the life of the tumor, corresponding to dierent instantaneous growth rates (gure 4.2a). A linear relationship between the instantaneous growth rate ( ) and instantaneous regression rate ()(gure 4.2b) emerges from the model, comparable to the predictions of the Norton-Simon hypothesis. A separate justication of the linear model of the eect of drug concentration on selection Chemotherapeutic dose scheduling 88 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.1 CM.43 T.1 (a) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.1 CM.43 T.1 (b) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.2 CM.43 T.2 (c) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.2 CM.43 T.2 (d) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.3 CM.43 T.3 (e) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) 0 20 40 60 80 100 120 t i m e 10 -7 10 -6 10 -5 10 -4 10 -3 n ( t ) CM.41 T.3 CM.43 T.3 (f) Figure 4.3: Response of murine tumors to 5-Fluorouracil (5-FU) treatment with model best-t | Data (reproduced from [142]) from two treated mice: CM.41 (left panel) and CM.43 (right panel), receiving doses of 50mg/kg or 100mg/kg, respectively, on 2 days out of 7 (biweekly). Biweekly measurements of tumor volume are recorded for untreated (black circles) until 3-4mm in size and treated volumes (black x's) are measured until tumor reaches 1cm size. A Gompertzian function is best-t (dashed line) and the Prisoner's dilemma model is t using w and c as parameters (solid line). The model t performs well for the wide range of tumor growth rates found in six tumors (w = [0:18; 0:08; 0:21; 0:08; 0:35; 0:12] and c = [0:30; 0:49; 0:34; 0:34; 0:36; 0:32] left to right, top to bottom, respectively). Chemotherapeutic dose scheduling 89 pressure is shown in gure 4.3. A best-t was performed to nd the optimal parameters of w andc to t data reproduced from mouse models quantifying inter-mouse and intra-mouse variability and response to 5-Fluorouracil (5-FU) in two treatment groups: 50mg/kg (gure 4.3 left panel) and 100mg/kg (gure 4.3 right panel) [142]. Tumor size measurements were taken from visibility until a mouse tumor reached 3 to 4 mm in size, and drug treatment was administered weekly until 1cm in size. The prisoner's dilemma model (black solid lines) appears to accurately capture both the growth dynamics (solid black circles) and the treatment dynamics (black x's). The dashed lines are Gompertzian best-t functions of the unperturbed pre-treament data (black circles), showing good aggreement with the prisoner's dilemma (black solid lines). Previously, we have reported the model's success in capturing current unperturbed growth models (i.e. Gompertzian growth) as an emergent phenomena of this evolutionary model [129]. 4.3.3 Dose concentration versus dose density Despite a growing trend toward personalized and precision medicine, treatment goals have shifted from complete cure to an optimization of long term management of the disease; rather than trying to nd the silver bullet, we might utilize the advances in mathematical models to optimize existing therapeutic options [130, 143]. For this reason, we have decided to test the merit of various chemotherapeutic regimens by comparing the total tumor cell reduction (TCR). Presumably, a therapy regimen with a higher value of TCR will provide a greater level of tumor control, a longer time to relapse, and better prognosis. Chemotherapeutic dose scheduling 90 A drug dose,D, (equation 4.13) is generally measured in units of mg/m 2 /week (here, aver- age body surface area assumed to be 1.8m 2 ). Yet, doseD consists of two components: dose concentration (parameterc in our model) and dose time factor (parametert in our model). The time factor, called the dose density when normalized by the intercycle time, represents the percentage of days a dose is administered. In order to compare the importance of each term on tumor cell reduction, we hold one term constant and vary the other in Figure 4.4. D =ct: (4.13) Clearly seen in Figure 4.4a, there is a diminishing return on increasing the dose strength of a given chemotherapy regimen. Although there is a positive relationship (an increase in dose leads to a higher regression) that relationship lessens as the dose is increased further. However, in Figure 4.4b, the relationship between dose density and regression is linear, showing no signs of diminishing returns of increasing density. The point has an important subtlety: the dose cannot be continually lowered in favor of density. The dose must be sucient to overcome the growth rate of the tumor; some doses are not adequate for tumor regression regardless of the density. This is seen for values below the dotted line in Figure 4.4a and 4.4b. Chemotherapeutic dose scheduling 91 (a) (b) Figure 4.4: Diminishing returns of dose escalation compared to linear rela- tionship of dose density | (a) Dose Escalation: The percent regression of a tumor for a range of dose strength (constant dose interval: t = 10 days,T = 14 days) are shown for a range of selection pressure: w = 0:1 (blue), w = 0:2 (red), and w = 0:3 (yellow). For each subsequent increase in dose strength, the dose escalation approach to chemotherapy shows diminishing returns in percent tumor regression. (b) Dose Density: The percent regression of a tumor for a range of dose interval (constant dose strength: c = 1:0) are shown for a range of selection pressure: w = 0:1 (blue), w = 0:2 (red), and w = 0:3 (yellow). Dose density shows a linear relationship between densifying chemotherapy and percent tumor regression. 4.4 Results 4.4.1 Quantifying chemotherapeutic strategies via entropy metric In clinical practice today, there are three common chemotherapy regimens in use considered here: Maximum Tolerated Dose (MTD), Low Dose Metronomic weekly (LDMw) and Low Dose Metronomic daily (LDMd). These three chemotherapy strategies are shown in Figure 4.5, left. Each regimen consists of identical cycles that are repeated until the tumor is eradicated. The MTD (left, top) regimen delivers the maximum dose on a single day, repeated once every 2 weeks. The LDMw (left, middle) regimen lowers the dose, but Chemotherapeutic dose scheduling 92 doubles the dose density from 1 to 2 days out of 14. The LDMd (left, bottom) regimen has the highest density (there is a dose administered on 100% of the days), but the lowest dose. Figure 4.5: Shannon entropy as an index to compare treatment strategies | Left: 3 common chemotherapy schedules are shown for one cycle (N = 14 days). Maximum Tolerated Dose (left, top) is a high dose (administered once at the beginning of every 2 week cycle) and low dose density (d = 0:071, see equation 4.16) regimen. Low Dose Metronomic Weekly (left, middle) is a lower dose, higher density (d = 0:143) regimen, while Low Dose Metronomic Daily is the lowest dose, highest density (d = 1:00). Right: Similarly, chemotherapy regimens can be simulated for a range of dose, density, and entropy values. Pictured from top to bottom are a range of representative regimens from low entropy (i.e. high dose, low density) to high entropy (i.e. low dose, high density) for a cycle of N = 4 days. On each ith day, treatment of dose c i is administered. The treatment strategy's Shannon Entropy, E, is calculated according to equation 4.14 and the total dose delivered is calculated according to equation 4.15. All treatment strategies are front loaded (monotonically decreasing) regimens. It should be noted that LDM-like regimens correspond to a high entropy value (bottom, left and right). There are hundreds of such choices of chemotherapy regimens when considering varying doses across many days or weeks (Figure 4.5, right), each varying the total dose delivered, Chemotherapeutic dose scheduling 93 D, and the density,d. We propose using a Shannon Entropy index,E, of a given chemother- apy schedule as a measure that can quantify and synthesize information of both the dose on a given day and the distribution of unique, daily doses across the entire chemotherapy regimen into a single metric. The entropy is calculated as follows, where c i is the dose strength (often simply referred to as `dose') on day i. E = N X i2c i >0 c i logc i (4.14) The assumption in equation 4.13 that an identical dose is delivered every day can be relaxed, and the total dose delivered is found by summing the dose on each ith day (c i ) multiplied by the length of the dose in days (t i ). We assume that the smallest resolution of discrete times between doses, t i is a single day, or t i = 1 for all i. N is the number of days between cycles, also known as the intercycle time. D = N X i=1 c i t i = N X i=1 c i (4.15) The dose density of a regimen can be found by summing the number of days where a non- zero dose is delivered, and dividing by the intercycle time in days, N. Thus, the density will be a non-dimensionalized parameter such that (0d 1). d = N X i2c i >0 t i N (4.16) Chemotherapeutic dose scheduling 94 The Shannon entropy metric, E, is an ideal metric for comparing chemotherapy regimens because it separates the existing cases already in clinical practice today: MTD (low entropy, characterized by high doses with long periods of rest), metronomic regimens (high entropy, characterized by low doses with short or no periods of rest), as well as any arbitrary strategy of varied doses administered in a cycle of arbitrary length of days. All of the simulated therapy regimens were assumed to be frontloaded (non-increasing, with the highest dose on day 1 and equal or lower subsequent doses). Backloaded regimens give similar but slightly disadvantageous results, because backloaded regimens often start with a period of rest, giving the tumor time to grow to a larger tumor, which is associated with a lower growth rate (and therefore lower regression). 4.4.2 LDM versus MTD chemotherapies Computational simulations of 1000 unique chemotherapy schedules were run with identical initial conditions (N = 1e6 cells; i=N = 1e3). Mean values of tumor cell regression percentage for 50 simulations were calculated and plotted in a pictorial histogram according to regression percentage (Figure 4.6). Both slow growing tumors (w = 0:1) and fast growing tumors (w = 0:2) were simulated. Each block represents a chemotherapy regimen, which has an associated Shannon entropy index (eqn. 4.14). The background color of the blocks of the chemotherapy regimens are color-coded from red (low entropy) to blue (high entropy). The smaller white squares within each block indicate the strength of the therapy dose for each day (c i ). Pictured are 1000 combinations of N = 4 day chemotherapy schedules, but similar trends are seen for Chemotherapeutic dose scheduling 95 chemotherapy schedules of longer length of days. All regimens are equivalent total dose (D = 0:3), non-increasing, and are repeated for 8 cycles of chemotherapy and the tumor cell regression (TCR) is recorded. The histograms clearly show a color-shift from red toward blue for low TCR toward high TCR. This indicates that high entropy (blue) therapies outperform low entropy therapies and consistently lead to higher tumor cell reduction. These high entropy regimens are low dose, more dose-dense chemotherapies, characteristic of LDM chemotherapy. In Figure 4.7, the analysis is repeated for varied tumor growth rates (i.e. varied selection pressure) for w = 0:1, (blue) w = 0:2, (red) and w = 0:3 (yellow). The dierence in reduction is shown for 1 cycle, 8 cycles, and 16 cycles. Fast growing tumors have a high slope on a least-squares linear t approximation of the entropy-TCR plot, which means that high entropy therapies (LDM) are more eective for fast growing tumors than for slow growing tumors. By contrast, slow growing tumors have a lower slope on the entropy-regression plot, which means that all regimens have relatively similar performance outcomes. Fast growing tumors, therefore, have a higher likelihood of beneting from a more LDM-like chemotherapy, provided the dose is adequate to lead to tumor regression. The eect is almost negligible after a single cycle (Figure 4.7a). The appeal of the impli- cation of Norton-Simon toward an MTD approach to chemotherapy lies in the high initial response of tumors to a high dose. The metronomic chemotherapies take more cycles to overtake the initial quick response of the MTD, but after the 8 cycles (Fig. 4.7b) and 16 cycles (Fig. 4.7c), the cumulative eect is evident and metronomic chemotherapies outper- form MTD therapies. For each growth rate (i.e. selection pressure), there is a correspond- ing optimal chemotherapy schedule. In each case, the optimal solution corresponds to the Chemotherapeutic dose scheduling 96 Figure 4.6: High entropy, LDM-like chemotherapies outperform low entropy MTD-like chemotherapies | Two pictorial histograms are plotted where each block (color-coded from red: low entropy to blue: high entropy) represents a chemotherapy regimen. The top histogram is for a slow-growing tumor (w = 0:1) and the bottom histogram is faster growth (w = 0:2). All regimens are equivalent total dose (D = 0:3), monotonically decreasing, and are repeated for 8 cycles of chemotherapy and the tumor cell reduction (TCR) is recorded. The dose density, d, and dose concentration, c i , are varied between regimens. The histogram clearly shows a color-shift from red toward blue for low TCR, ineective therapies toward high TCR, eective therapies. High entropy (blue) therapies outperform low entropy therapies. The data was t to a Weibull distribution (shown in upper left panel; top: k = 14:251; = 65:882, bottom: k = 6:647; = 46:758), overlaid in green. Chemotherapeutic dose scheduling 97 highest entropy (which corresponds to the low-dose metronomic chemotherapy schedule). 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 8 cycles 16 cycles 1 cycle y=3.72x+2.23 y=3.62x 0.29 y=3.33x 2.20 y =19.4x+19.0 y =22.6x+0.95 y = 24.0x 15.4 y =24.0x+37.9 y =34.9x+6.32 y = 44.9x = 28.8 tumor cell reduction (TCR) entropy (E) tumor cell reduction (TCR) tumor cell reduction (TCR) entropy (E) entropy (E) Figure 7 (a) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 8 cycles 16 cycles 1 cycle y=3.72x+2.23 y=3.62x 0.29 y=3.33x 2.20 y =19.4x+19.0 y =22.6x+0.95 y = 24.0x 15.4 y =24.0x+37.9 y =34.9x+6.32 y = 44.9x = 28.8 tumor cell reduction (TCR) entropy (E) tumor cell reduction (TCR) tumor cell reduction (TCR) entropy (E) entropy (E) Figure 7 (b) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -20 -10 0 10 20 30 40 50 60 70 80 w = 0.1 w = 0.2 w = 0.3 8 cycles 16 cycles 1 cycle y=3.72x+2.23 y=3.62x 0.29 y=3.33x 2.20 y =19.4x+19.0 y =22.6x+0.95 y = 24.0x 15.4 y =24.0x+37.9 y =34.9x+6.32 y = 44.9x = 28.8 tumor cell reduction (TCR) entropy (E) tumor cell reduction (TCR) tumor cell reduction (TCR) entropy (E) entropy (E) Figure 7 (c) Figure 4.7: High entropy strategies lead to an increase in tumor regression | The relationship between tumor cell reduction (TCR) and entropy (H) is shown for a single cycle of chemotherapy (a), 8 cycles (b), and 16 cycles (c). The simulations (averages of 25 stochastic simulations for total dose delivered D = 0:3) are repeated for slow (w = 0:1, blue), medium (w = 0:2, red), and fast growing tumors (w = 0:3, yellow). The low slope value in (a) indicates negligible advantage of high entropy strategies after only a single cycle. After many cycles, the advantage of high entropy strategies is apparent (b,c). Also note that the slope associated with faster growing tumors (yellow) is higher than those of slower growing tumors (blue). This indicates that at high entropies, TCR for the fast growing tumors is closer to those for slow growing tumors, as compared with low entropies. 4.5 Discussion We use a stochastic Moran process model coupled with a prisoner's dilemma evolutionary game (cellular interactions) to contrast LDM and MTD chemotherapies with respect to their eect on tumor growth. The Shannon entropy was identied as a useful metric to compare chemotherapy strategies. The metric is useful in quantifying LDM strategies (which correspond to high entropy values), MTD strategies (low entropy), as well as novel strategies with intermediate entropy values. Our results show that high dose chemotherapy strategies outperform low dose, although there are some subtleties associated with the growth rates of the tumors. Dosing consists of Chemotherapeutic dose scheduling 98 a product of concentration and density and our results show that an increase in density is more eective than the same percentage increase in concentration. In other words, higher dose concentrations shown diminishing returns. The eectiveness of density in leading to a higher tumor cell reduction allows the LDM chemotherapies (which are more dose dense) to outperform MTD strategies. This eect is magnied for fast growing tumors that thrive on long periods of unhindered growth without chemotherapy drugs present. This eect is not evident after a single cycle of chemotherapy, but is magnied after each subsequent cycle of repeated chemotherapy. We could ask if there is any evidence of this eect in the literature on clinical trials already performed. We rst point to a paper comparing dierent chemotherapeutic schedules for prostate tumors [144] (relatively slow-growth rates). In this phase 3 study, docetaxel dosing given every three weeks was compared to dosing every week. The mean survival was only slightly higher for the rst group (three weeks) compared with the second (weekly), showing no obvious benet to a low-dose high density treatment. By contrast, a phase 2 trial for small cell lung cancer (SCLC) [145] was performed, a tumor with typically higher growth rates than prostate tumors. For this group, the drug topotecan was administered on a higher dose weekly basis with disappointing results, pointing out the advantages of the LDM therapies for this fast-growing tumor type. Thus our model points to the benets of choosing dosing strategies based on tumor growth rates, something not currently done in medical practice. The concept of choosing dosing schedules based on tumor growth rates could well be a fruitful avenue to test further in clinical trials focused on this question. Others have attempted to estimate prospective patient-specic tumor growth rates to make clinical decisions about treatment scheduling and fractionization, using measurements at diagnosis and rst day of treatment [133, 146]. Chemotherapeutic dose scheduling 99 Furthermore, the promise of LDM chemotherapy on mitigating the risk of resistance [126] and metastasis [130] could be a separate line of future investigation. Chapter 5 Harnessing the evolutionary cost of chemotherapeutic resistance by shaping the tness landscape of a tumor Jerey West, Yongqian Ma, Paul K. Newton Harnessing the evolutionary cost of chemotherapeutic resistance by shaping the tness landscape of a tumor [41] USC Preprint 100 Harnessing the evolutionary cost of resistance 101 5.1 Abstract Pre-existing resistant clones present in a tumor at the start of treatment re- mains a major problem in cancer therapeutics today. Tumor relapse and the development of chemotherapeutic resistance is now thought largely to be a consequence of the mechanism of \competitive release" of pre-existing resis- tant cells in the tumor which are selected for growth after chemotherapeutic agents attack the sub-population of chemo-sensitive cells which had previously dominated the collection of competing sub-clones in the tumor. We present an evolutionary game theory model based on the replicator equation to de- scribe the clonal competition according to a tness landscape, where the less t resistant clone is released from competition after continuous chemotherapy, leading to a resistant tumor upon relapse. We explain the important param- eters (cost of resistance and initial fraction of resistance) in anticipating the evolutionary adaptations of the tumor in order to design therapies that exploit or mitigate the harmful eects of potential future adaptions. As the tumor is growing according to Darwinian principles, so must our treatment schedules be designed with these same Darwinian principles in mind. Tumor adaptations only adapt to local and current selection forces and never anticipate the future. In contrast, treatment protocols must be dynamic and models must be evolu- tionary in scope, so we present a dynamic adaptive therapy control paradigm to indirectly control the resistant population to shape the tness landscape by carefully choosing the timing of \on" and \o" chemotherapy periods. Harnessing the evolutionary cost of resistance 102 5.2 Introduction In his now classic 1961 study of competition for space between two species of barnacles in the intertidal zone o the Scottish coast, Joseph Connell [147] discovered something interesting. The blue barnacles Balanus normally occupied the intertidal zone, while the brown barnacles Chthamalus occupied the coast above high tide. Despite the commonly held belief that each occupied their own niche because of dierent adaptations to local micro-conditions, Connell hypothesized that the colonization of the intertidal zone by Bal- anus was actually preventing Chthamalus from inhabiting this region. To test this, he removed the blue barnacles from the intertidal zone and tracked the subsequent penetra- tion of Chthamalus into this region. He concluded that Chthamalus had undergone relief from competition with Balanus which allowed it to ourish where previously it could not. The point, he emphasized, was there was nothing inherent about the micro-environment of the intertidal zone that was preventing Chthamalus from occupying this region. It was simply the competition against a more dominant species that was holding it back. With- out the presence of that species, Chthamalus happily claimed both zones as fundamental niches. Thus, the important notion of competitive release was formulated (see Grant [148]). When two (or more) sub-species compete for the same resources, with one species domi- nating the other, if the dominant species is removed, this can provide the needed release from competition that can allow the less dominant species to ourish. The mirror image of competitive release is the related notion of character displacement developed by Brown and Wilson [149] in which competition can serve to displace one or more morphological, ecological, behavioral, or physiological characteristics of two closely related species that Harnessing the evolutionary cost of resistance 103 develop in close proximity. These concepts are now well established as part of the overall framework of co-evolutionary ecology theory. Since co-evolution among competing sub-clones is now a well established [10, 11, 24, 25] phylogenetic process in malignant tumors, the mechanism of competitive release should be expected to play a role and aect the chemotherapeutic strategies one might choose to eliminate or control tumor growth [150]. Indeed, tumor relapse and the development of chemo-therapeutic resistance is now thought largely to be a consequence of the mechanism of competitive release of pre-existing resistant cells in the tumor which are selected for growth after chemotherapeutic agents attack the sub-population of chemo-sensitive cells which had previously dominated the collection of competing sub-clones. Anticancer thera- pies strongly target sensitive cells in a tumor, selecting for resistance cell types and, if total eradication of all cancer cells is not accomplished, the tumor will recur as derived from resistant cells that survived initial therapy [3, ch. 1]. Subsequent application of identical therapies will have a diminished eect. A recent (2012) systematic literature analysis of cancer relapse and therapeutic research showed that while evolutionary terms rarely ap- peared in papers studying therapeutic relapse since 1980 (< 1%), the language usage has steadily increased more recently, due to a huge potential benet of studying therapeutic relapse from an evolutionary perspective [151]. Our goal in this paper is to describe an evolutionary mathematical model of competi- tive release in a tumor in order to better quantify and understand what we feel is a key mechanism responsible for the evolution of chemotherapeutic resistance with the hope that understanding it could ultimately prove crucial for controlling and harnessing the evolu- tionary engine that drives its growth. We also describe quantitative tools to shape the Harnessing the evolutionary cost of resistance 104 tness landscape of a tumor using chemotherapeutic strategies that are as dynamic as the tumor. 5.2.1 Pre-existing resistance Cancer therapies have shown success in reducing tumor burden for signicant time peri- ods, but eventual relapse and resistance have led many to use evolutionary principles and mathematical modeling to answer whether resistance arises de novo during therapy or is pre-existing before therapy. Pre-existing resistant sub-clones should generally be present in all patients with late-stage metastatic disease (for single point mutations which confer resistance) as told by probabilistic models [152] and from samples taken prior to treatment [153, 154] which have been reported for melanoma [155], prostate cancer [156], colorec- tal cancer [157, 158], ovarian cancer [159], and medulloblastoma [160]. This implies that treatment failure is not due to evolution of resistance but rather pre-treatment presence of resistant phenotypes that are relatively sheltered from the toxic eects of therapy [161]. Pre-existing resistance has important therapeutic implications. If we assume no pre- existing resistance then most models predict maximum dose-density therapy will reduce the probability of resistance (because this treatment minimizes the number of cell-divisions, thereby minimizing the risk of de novo resistance) [161]. To be clear, when curative therapy is possible: treatment strategy must be designed to achieve that result. But in pre-existing resistance scenarios, the maximum dose-density therapy strategy lends itself to competitive release due to the evolutionary nature of tumor progression. Most preclinical eorts that aim to maximize the short-term eect of the drug on sensitive cells does not signicantly Harnessing the evolutionary cost of resistance 105 aect the long-term control of cancer [152]. This is because the phenomenon of competitive release can occur via the harsh selective pressure imposed by the tumor microenvironment after cancer therapies diminish the presence of the dominant (i.e. the chemo-sensitive) clone, or the process of metastasis may allow a resistant subclone in the primary tumor to emerge to higher frequencies [150]. These pre-existing mutations that are responsible for conferring resistance may be dele- terious and may be associated with a phenotypic cost, or a reduced tness, compared to sensitive cells. Even with this tness cost, deleterious mutations are still expected to be present in late-stage metastatic cancers [162]. This cost can come in many ways, such as an increased rate of DNA repair or an active pumping out of the drug across cell mem- brane. All these strategies use up energy that would otherwise be available for invasion into non-cancerous tissues or proliferation. Tumors that have not yet undergone treatment may possess resistant cells in small numbers because a tness disadvantage allows the sen- sitive population to suppress the resistant population to a smaller population. During drug therapy, the rapid removal of chemo-sensitive cells permits unopposed proliferation of the resistant cell population, an evolutionary phenomenon known as competitive release. A graphical representation of clonal competition within a growing tumor known as a \sh- plot" (sometimes known as a Muller plot) was rst utilized in cancer to compare modes of clonal evolution in acute myeloid leukemia (see [163]). A shplot shows the tumor burden (vertical axis) over time (horizontal axis) and the clonal lineages (a subclone is encased inside of the founding parent clone in the graph). A sample shplot, seen in gure 5.1, displays a schematic of unhindered tumor growth after the rst driver mutation (gure 5.1, left). Before diagnosis, the tumor grows exponentially, during which time a resistant Harnessing the evolutionary cost of resistance 106 Figure 0 total tumor burden sensitive subpopulation resistant subpopulation first driver mutation exponential growth of chemo-sensitive population first resistant mutation begin continuous chemotherapy treatment good initial response to chemotherapy competitive release of resistant population Chemo-sensitive Chemo-resistant time tumor burden Figure5.1: Clonal evolution of competitive release | A shplot (sometimes known as a Muller plot), showing the tumor size (vertical axis) and composition (sensitive: red; resistant: green) over time (horizontal axis, left to right) with important events annotated. After rst driver mutation (left), initial exponential growth of sensitive population occurs until diagnosis (dashed line). Continous therapy targeting the chemo-sensitive population responds well with a decrease in tumor burden. In the absence of sensitive cells, the resistant population (existing in small numbers before the start of therapy) grow to become the dominant clone at relapse, albeit with lower exponential growth rate. mutation occurs (gure 5.1, middle). After diagnosis, a regimen of continuous chemother- apy shows initial good response and tumor regression the resistant population grows back (although at a slower growth rate) unhindered by competition, leading to relapse (gure 5.1, right). 5.2.2 Using evolutionary principles to model chemotherapy There is a need for a change in strategy in the war on cancer: eradicating most disseminated cancers may be impossible, undermining the typical goal of cancer treatment of killing as many tumor cells as possible [143]. The underlying assumption of this approach has been that a maximum cell-kill will either lead to a cure or, at worst, the maximum life extension. Taking cues from agriculturists who have long abandoned the goal of complete eradication of pests in favor of applying insecticides only when infestations exceeds a threshold in Harnessing the evolutionary cost of resistance 107 the name of \control" over \cure," so we must also change from the cure paradigm in cancer treatments to a control paradigm [143]. The rst step in this paradigm shift is viewing tumor progression from an evolutionary lens. As such, any therapeutic methods should take the following parameters into account: the pre-existing fraction of the resistant population in the tumor before therapy and the relative growth rates (i.e. the tness cost) of resistant subclones. Several treatment strategies have been proposed to exploit or predict the evolutionary trajectory of tumor growth and adaptations, such as targeting the trunk driver events that are present in every tumor cell, target parallel evolutionary events, forcing the tumor down a specic evolutionary path, resulting in acquired sensitivity (sequential therapy), and dynamic therapies that maintain a stable population of treatment-sensitive cells [150]. Some have proposed modelling tumorigenesis as a process by which the homeostasis that characterizes healthy tissue is disrupted, explaining how the order of treatments can take advantage of an evolutionary double bind [26, 29, 139], thereby predicting tumor adapta- tions and exploiting that prediction using fundamental evolutionary principles. Regaining homoeostasis might not mean tumour eradication but instead may represent a new state where we live with cancer as a controllable, yet chronic disease. Treatments can be syner- gized such that evolving resistance to a single drug will increase susceptibility to a dierent drug. Others are modeling and planning \evolutionary enlightened" therapies, known as \adaptive therapies" that respond to the tumor's adaptations in order to make future treatment decisions. A theoretical framework for these adaptive therapies rst developed by Gatenby [42], leverages that pre-existing resistance is typically present only in small population numbers due to a cost of resistance. This less t phenotype is suppressed in Harnessing the evolutionary cost of resistance 108 the Darwinian environment of the untreated tumor but treatments that are designed to kill maximum numbers of cells remove the competition for the resistant population and ulti- mately select for that population during tumor relapse. (It's important to note that both high-dose, maximum tolerated dose schedules and low-dose, metronomic dose schedules have this cumulative goal of achieving maximum cell-kill over the course of many cycles of treatment.) In contrast, the goal of an adaptive therapy is to maintain a stable tumor burden that permits a signicant population of chemo-sensitive cells for the purpose of suppressing the less t but chemo-resistant populations. It takes an evolutionary strategy to combat an evolutionary tumor. An optimal treatment will indeed maintain a stable population of chemo-sensitive cells that suppress resistant populations when therapy is absent. In other words, any treatment with either the explicitly or implicitly stated goal of eliminating the maximum amount of tumor cells will likely never achieve cure or the maximum survival time. These ideas were tested experimentally using mouse models to optimize adaptive strate- gies designed to maintain a stable, controllable tumor volume [164, 165]. The two-phase adaptive therapy involved an initial high-dose phase to treat the exponential growth of the tumor and a second phase designed for stable tumor control using a variety of strategies (such as decreasing doses or skipping doses when stability is achieved). Findings suggest that adaptive therapies based on evolutionary treatment strategies that maintain a resid- ual population of chemo-sensitive cells may be clinically viable, and is currently extended to an on-going clinical trial for metastatic castration resistant prostate cancer patients (NCT02415621). Harnessing the evolutionary cost of resistance 109 Healthy Sensitive Resistant 3 species model of competitive release therapy Tumor is dominated by sensitive cells which out- compete both healthy cells and resistant cells. Pre-treatment After chemotherapy the sensitive cell population is reduced, leaving resistant cells free to re-populate without competition. Resistant cells eventually repopulate a larger proportion of the tumor, rendering the subsequent rounds of chemotherapy less effective. re-growth Post-treatment Tumor re-growth Figure 5.2: Schematic of competitive release in a tumor | Prior to treatment (left), a tumor consists of a large population of sensitive cells (red) and a small population of less t resistant cells (green) competing for resources with the surrounding healthy cells (blue). Chemotherapy targets the sensitive population (middle), selecting for the less t resistant population that thrives in the absence of competition from the sensitive population. Upon regrowth, the tumor composition has larger numbers of resistant cells, rendering the subsequent rounds of treatment less eective. The goal of this manuscript is to introduce an evolutionary framework to model the im- portant parameters of competitive release (cost, initial fraction of pre-existing resistance) and use that framework to better understand therapeutic implications of the tumor evolu- tion, then present a few sample simulations of adaptive therapy techniques and ideas. A schematic of the three compartment model of competitive release is shown in gure 5.2, where the tumor consisting of sensitive and resistant cells is competing with the surround- ing normal/healthy tissue. At diagnosis (see gure 5.2, left), the tumor is dominated by sensitive cells (red) which outcompetes the surrounding healthy population (blue) during unhindered tumor progression. A small portion of resistant cells (green) remains in small Harnessing the evolutionary cost of resistance 110 numbers, suppressed by the larger sensitive population. After several rounds of chemother- apy, the tumor shrinks, leaving the resistant population largely unaected (gure 5.2, mid- dle). Inevitably, the tumor relapses due to the small number of cancer cells remaining after therapy (gure 5.2, right). In the absence of competition from the dominant sensitive pop- ulation, the resistant cells grow unhindered, rendering subsequent rounds of chemotherapy less eective. In the next section, we will introduce the details of the three compartment model using the replicator equation dynamics and introduce the idea of shaping the tness of landscape of the evolving tumor in order to combat the competitive release phenomena. 5.3 Materials and Methods Previously, a model predicting and quantifying competitive release has been proposed to track the relative tumor volume, v(t), after treatment as a function of the exponential death rate of the sensitive cells,d, the exponential growth rate of the resistant cells,g, and the initial fraction of resistant cells, f [152]. The model can be written as follows. v(t) = (1f)e dt +fe gt (5.1) This model, shown to be a good description of the changing tumor size during therapy for colorectal, prostate, and multiple myeloma cancers, identies the important parameters in competitive release: initial fractional resistance (f), and birth/death rates (g,d) for each resistant and sensitive populations, respectively. The regrowth rate of the resistant population (g) aects the eectiveness of a continuous therapy (see gure 5.3a). A greater Harnessing the evolutionary cost of resistance 111 \cost" of resistance (re ecting by a lower regrowth rate, g) leads to a longer time to relapse. Although the tumor might recur with a slower growth rate, subsequent treatment is ineective due to resistant population in large numbers (again, see gure 5.1). The initial fraction of resistant cells present at the time of treatment, f, also aects treatment eectiveness (gure 5.3b). An increase in initial fraction leads to a shorter time to relapse. A small dierence in resistant regrowth rate (two-fold decrease fromg = 0:2 to 0:1) leads to an earlier relapse time compared to a large dierence in initial fractional resistant (hundred- fold increase from f = 10 6 to 10 4 ), indicating that the regrowth rate parameter has a greater eect on the eectiveness of a continuous therapy (gure 5.3c). However, the model provides no evolutionary information or concepts, a keystone principle behind competitive release. The next section will outline a model that encompasses the same characteristics of competitive release, but using evolutionary modeling concepts in order to test adaptive, \evolutionary enlightened" chemotherapy strategies in the nal section. 5.3.1 The Model Evolutionary game theory (EGT, see [14, 15, 17, 19, 134]), which has been used extensively to study cancer progression, provides a quantitative framework for analyzing contests (com- petition) between various species in a population (via the association of `strategies' with birth/death rates and relative sub-clonal populations) and provides mathematical tools to predict the prevalence of each species over time based on the strategies [17, 64, 66, 72, 166]. The framework denes mutual payos among pairs within the subclonal populations com- peting for resources in the tumor. From this, one can dene a tness landscape (related to population growth rates) over which the subclonal populations evolve in time. The aim Harnessing the evolutionary cost of resistance 112 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 Figure 0: EXPONENTIAL LINEAR COMBO MODEL Initial fraction of resistance, f Regrowth rate vs initial fraction Regrowth rate, g 0 10 20 30 40 50 10 -4 10 -3 10 -2 10 -1 10 0 f = 10 -6 , g = 0.2 f = 10 -4 , g = 0.1 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 g = 0.10 g = 0.15 g = 0.20 g = 0.25 g = 0.30 (a) 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 Figure 0: EXPONENTIAL LINEAR COMBO MODEL Initial fraction of resistance, f Regrowth rate vs initial fraction Regrowth rate, g 0 10 20 30 40 50 10 -4 10 -3 10 -2 10 -1 10 0 f = 10 -6 , g = 0.2 f = 10 -4 , g = 0.1 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 g = 0.10 g = 0.15 g = 0.20 g = 0.25 g = 0.30 (b) 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 Figure 0: EXPONENTIAL LINEAR COMBO MODEL Initial fraction of resistance, f Regrowth rate vs initial fraction Regrowth rate, g 0 10 20 30 40 50 10 -4 10 -3 10 -2 10 -1 10 0 f = 10 -6 , g = 0.2 f = 10 -4 , g = 0.1 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 g = 0.10 g = 0.15 g = 0.20 g = 0.25 g = 0.30 (c) Figure 5.3: Dynamics of competitive release under continuous therapy | The dynamics described by the simple model of fraction, f, of resistant cells at the start of continuous therapy (equation 5.1). During therapy, sensitive cell population decreases at rate d and resistant population grows at rate g. Left: An increase in the resistant regrowth rate (g = [0:1; 0:15; 0:2; 0:25; 0:3] in equation 5.1) leads to decreased eectiveness of therapy (f = 1e 5;d = 0:3;). Middle: an increased initial fraction of resistance (f = [1e 6; 1e 5; 1e 4; 1e 3; 1e 2]) also leads to decreased therapy eectiveness (g =d = 0:3;). Right: a small dierence in resistant regrowth rate (two-fold decrease from g = 0:2 to 0:1) compared to a large dierence in initial fractional resistance (hundred-fold increase fromf = 10 6 to 10 4 ). The regrowth rate parameter has a greater eect on the eectiveness of a therapy. of this research is to develop an evolutionary model based on the replicator equation dy- namical framework (described below) of tumor progression and competitive release due to continuous chemotherapy, and use that model to identify and test treatment implications. This work builds on previous work from the authors (see [32, 129]) used to carry out computational trials testing chemotherapy drug scheduling strategies (see [40]) where two cell types (a healthy and a cancerous cell type) compete according to a tness landscape described by the evolutionary prisoner's dilemma game (see [14, 75]). This game incorpo- rates two general classes of cells: healthy (the cooperators) and cancerous (the defectors) [135, 136]. Natural selection acts on each generation of the cell population as the com- putational simulation proceeds, selecting (on average) for cells with higher tness [129], Harnessing the evolutionary cost of resistance 113 leading to tumor growth and eventual saturation of the cancerous population. The work is extended to include a third player: the resistant cell type. The tness landscape of interactions between all three players (healthy: blue; chemo-sensitive cancer cells: red; and chemo-resistant cancer cells: green) is depicted in gure 5.4 before therapy (top) and dur- ing therapy (bottom). In our model, we can think of a cancer cell (either a chemo-sensitive or chemo-resistant cancer cell) as a formerly cooperating healthy cell that has defected (via a driver mutation) and begins to compete against the surrounding population of healthy cells for resources and reproductive prowess. This is indicated in gure 5.4, showing the chemo-sensitive cancer cell as a higher tness (red) than a healthy cell (blue). In the pro- cess of acquiring mutations during early tumor growth, a resistant cell may arise (green). As described previously, the development of a resistant-conferring mutation is modelled as a tness cost; thus the resistant cell is placed further down on the tness scale (gure 5.4, top: green) yet still a higher tness than the healthy cell type. In this way, the pre-treatment dynamics are described by 3 separate prisoner's dilemmas. The interactions between healthy and sensitive are described by the prisoner's dilemma (tness of sensitive is greater than healthy, as the cancerous sensitive cells have a higher proliferative potential than the surrounding normal healthy cells). The interactions be- tween healthy and resistant are also described by the prisoner's dilemma (tness of re- sistant is greater than healthy). Finally, the sensitive population tness is greater than resistant tness (due to the cost of resistance): a third prisoner's dilemma. During therapy, the tness of chemo-sensitive cancer cells dramatically decreases and the healthy population tness increases (see gure 5.4, bottom). The resistant cell type is un- aected by the chemotherapy treatment, leaving the tness the same as prior to treatment Harnessing the evolutionary cost of resistance 114 (gure 5.4, bottom: green). Figure 2: Fitness paradigm Fitness advantage of driver mutation (prisoner’s dilemma) no therapy Fitness Sensitive Cancer Cell Resistant Cancer Cell Healthy Cell Fitness cost of resistance Fitness of resistant subclones remain constant before, during therapy Fitness increase due to altered selection pressure during therapy Low Fitness High Fitness Fitness decrease due to altered selection pressure Figure 5.4: Fitness landscape before and during therapy | A schematic of the tness of each subpopulation before therapy (top) and during therapy (bottom). A driver mutation leads to a tness advantage of the cancer cell (red), determined by the prisoner's Dilemma. A subsequent resistant-conferring mutation comes at a tness cost (green). The tness of the resistant population is unaected by therapy's selective pressure, but the healthy population is given an advantage over the chemo-sensitive population. 5.3.2 Replicator equation dynamics The dynamics of the tness landscape of the three competing cell types depicted in gure 5.4 are described by the replicator equation (see [116]). Each ith cell type (i = 1; 2; 3) competes according to the replicator equation (eqn. 5.2, below), where x 1 , x 2 , x 3 are the corresponding frequency of healthy, sensitive and resistant cells, respectively, such that Harnessing the evolutionary cost of resistance 115 P i x i = 1. The change in frequency of each cell type, _ x i , is a function of the proportion (x i ) and tness of that cell type (f i ). The tness is a function of the selection pressure, w i ; 0 w 1 1, and the payo matrix, A. A value of w = 0 corresponds to neutral drift (no selection) and a value of w = 1 corresponds to strong selection. _ x i = (f i hfi)x i (5.2) f i = 1w i +w i (A~ x) i (5.3) Here,~ x is the vector~ x = (x 1 ;x 2 ;x 3 ) T and (Ax) i is theith element of vectorAx. The preva- lence of each sub-population, x i , changes over time according to the changing population tness, f i , and the average tness of all populationshfi =f 1 x 1 +f 2 x 2 +f 3 x 3 . Before therapy, the selection pressure is constant across all cell type (i.e. w i =w). During therapy, the values ofA remain constant (described in more detail below) and the selection pressure values are altered as shown in equations 5.4, 5.5, and 5.6 (see gure 5.4 for explanation of changing tness landscape). Therapy can be administered at dierent doses (i.e. values of the drug concentration: c; 0 c 1). A higher value of c indicates a stronger dose of chemotherapy drug (described in more detail in [40]). The selection pressure parameters are linear functions of the drug concentration, comparable to other cycle-specic drug therapy models [7, 8], explained in greater detail in our previous work here [40]. Harnessing the evolutionary cost of resistance 116 w 1 = (1 +c)w (healthy) (5.4) w 2 = (1c)w (sensitive) (5.5) w 3 = w (resistant) (5.6) 5.3.3 The pairwise prisoner's dilemmas The tness landscape is described by the entries of the payo matrix, where each pair- wise interaction is described by the row and column values, which are parameters in the population tness, equation 5.3. H S R H S R 0 B B B B B B @ a b c d e f m n g 1 C C C C C C A = 0 B B B B B B @ 1:2 1 1 1:2 +p 2 1:01 1:03 1:2 +p 3 1:01p 1 1:02 1 C C C C C C A =A (5.7) As described previously, there are three pairwise games between healthy (H), sensitive (S), and resistant (R) populations: 1. (H,S), 2. (H,R), and 3. (R,S) which are all described as prisoner's dilemma (cooperators, defectors). The three pairwise prisoner's dilemmas necessitates the following inequalities: d>a>e>b, m>a>g >c, and g >f >n>e. More discussion of why the prisoner's dilemma matrix, which models the evolution of defection, is a useful paradigm for cancer can be found in our previous work [32, 129]. Values are chosen as small perturbations from 1, consistent with the necessary inequalities. Harnessing the evolutionary cost of resistance 117 Thus the 9 parameters of the payo matrix, A, are reduced to 3 parameters: p 1 ,p 2 ,p 3 . A new set of parameters ( 2 , 3 , see equation 5.8) is dened such that i =w i p i , where i is the initial growth rate of the sensitive ( 2 ) and resistant ( 3 ) populations. 8 > > < > > : _ x 2 j x 3 =0;x 2 !0 w (da)x 2 =wp 1 x 2 = 1 x 2 _ x 3 j x 2 =0;x 3 !0 w (ma)x 3 =wp 2 x 3 = 2 x 3 (5.8) An important feature of the model is the cost of resistance. Strictly speaking, a slower growth rate of resistant cell population (i.e. 3 < 2 ) implies a cost of resistance. Im- portantly, a second aspect of the cost must also be considered: the prisoner's dilemma between sensitive and resistant cells (i.e. p 1 > 0) must be present in order for the sensitive population to suppress the resistant cell population in the absence of therapy. The rest of the work assumes a constant value of selection pressure (w = 0:1) and a constant value for suppression (p 1 = 0:7) in order to investigate the cost parameter (cost = 3 2 ). 5.4 Results The model of competitive release has two important parameters: the phenotypic cost of resistance (the dierence between the unperturbed growth rates of the sensitive and resistant populations: 2 3 ), and the initial fraction of resistance (f). The eect of both parameters is shown in gure 5.5. An increased cost of resistance (gure 5.5a) decreases the regrowth rate and leads to a longer relapse time under continuous therapy. The initial regression rate is identical at the beginning of therapy, but the regrowth rate is diminished by the cost of resistance. An increase in initial fraction (gure 5.5b) of resistant at the Harnessing the evolutionary cost of resistance 118 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 Figure 5: Growth rates, fractional resistance, and fitness cost Phenotypic cost of resistance, ! 1 -! 2 Initial fraction of resistance, f 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 0 10 20 30 40 50 60 10 -2 10 -1 10 0 Fraction vs Cost f = 10 -6 , f = 10 -2 , 2 3 = 0.08 2 3 = 0.06 2 3 = 0.04 2 3 = 0.02 2 3=0 2 3=0 2 3 = 0.1 (a) 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 Figure 5: Growth rates, fractional resistance, and fitness cost Phenotypic cost of resistance, ! 1 -! 2 Initial fraction of resistance, f 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 0 10 20 30 40 50 60 10 -2 10 -1 10 0 Fraction vs Cost f = 10 -6 , f = 10 -2 , 2 3 = 0.08 2 3 = 0.06 2 3 = 0.04 2 3 = 0.02 2 3=0 2 3=0 2 3 = 0.1 (b) 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 Figure 5: Growth rates, fractional resistance, and fitness cost Phenotypic cost of resistance, ! 1 -! 2 Initial fraction of resistance, f 0 20 40 60 80 10 -3 10 -2 10 -1 10 0 f = 10 -6 f = 10 -5 f = 10 -4 f = 10 -3 f = 10 -2 0 10 20 30 40 50 60 10 -2 10 -1 10 0 Fraction vs Cost f = 10 -6 , f = 10 -2 , 2 3 = 0.08 2 3 = 0.06 2 3 = 0.04 2 3 = 0.02 2 3=0 2 3=0 2 3 = 0.1 (c) Figure 5.5: The eect of growth rate, cost, and fractional resistance under continuous therapy | (a) Increasing the phenotypic cost of resistance ( 2 3 ), results in an identical initial regression rate but extended relapse time ( 2 = 0:02,f = 10 3 ). The slope of the regrowing tumor decreases with cost. (b) Increasing the fractional resistance results in an identical initial regression rate but extended relapse time ( 2 = 3 = 0:02,f = 10 3 ). (c) Patient 1 (blue) has a low initial fractional resistance (f = 10 6 ) with no cost to resistance. Patient 2 (red) has four orders of magnitude greater fractional resistance with a relatively small cost 2 3 = 0:1 and yet has a shorter time to relapse. This implies that the cost of resistance (i.e. the diminished regrowth rate of the tumor) is more important to relapse than the initial fraction of resistance. start of therapy leads to shorter relapse times. The initial regression rate and the regrowth rate are identical for all simulations but a small fractional resistance leads to a lower minimum tumor burden achieved and longer relapse times. Immediately, it is important to ask which of these parameters has the most signicant eect on time to relapse. In gure 5.5c, the relapse times for two patients are compared. Patient 1 (blue) has a low initial fractional resistance (f = 10 6 ) with no cost to resistance. Patient 2 (red) has four orders of magnitude greater fractional resistance (f = 10 2 ) with a relatively small cost 2 3 = 0:1 and yet has a shorter time to relapse. This implies that the cost of resistance (i.e. the diminished regrowth rate of the tumor) is more important to relapse than the initial fraction of resistance. Even a small cost of resistance leads to a signicant change in tumor relapse dynamics. This is consistent with previous models, shown in gure 5.3c. Harnessing the evolutionary cost of resistance 119 Evolutionary game theory models allow testing two dierent treatment paradigms: 1) treatments that target the strategies (or, subpopulations of cells) and 2) treatments that target the game, by anticipating or driving the evolutionary trajectory to more amenable dynamical spaces [167]. Traditional treatment methods (i.e. maximum tolerated dose schedules or low dose metronomic schedules) have the goal of complete eradication of all tumor cells (or, at least the eradication of all chemo-sensitive cells), encourage the competitive release phenomenon in the presence of pre-existing resistance. Figure 5.5 and 5.6 (below) detail the important outcomes of the rst paradigm: targeting the strategies (i.e. indiscriminately targeting chemo-sensitive cells, without a view to tumor evolution). In the nal section, gure 5.8 is a proof-of-concept adaptive strategy aimed at anticipating the evolutionary trajectory (i.e. \targeting the game") to improve outcomes. 5.4.1 Treating the strategy: competitive release Figure 5.6 details the relationship between dose concentration and two important measures of therapy eectiveness: progression free survival (PFS) and time to relapse. Measuring the eectiveness of a chemotherapy based on the killing rate or progression free survival alone are not sucient predictive measures of long-term cancer control [152]. As seen in gure 5.6, left, an increased dose (3 therapies are simulated on identical initial conditions: 2 = 0:020, 3 = 0:018, f = 10 3 ) corresponds to a slightly shorter PFS (an increase from blue to red to yellow), but an increased time to relapse to the initial tumor volume. However, despite the increase relapse times, none of these doses optimize tumor control, as seen in the shplots (gure 5.6, right). At the point of relapse to the initial tumor volume, the tumor is dominated by the presence of resistant clones (green), rendering Harnessing the evolutionary cost of resistance 120 future treatments ineective. Oftentimes, the eectiveness of a new chemotherapy drug may determined by PFS times when drugs that have high killing rates of sensitive cells may have shorter times to progression and lower total tumor burden at all times, everything else equal. The gure clearly shows that all treatments have similar progression free times but with a greater range of relapse times (even though continuous treatment always eventually leads to relapse). Figure 4: Time to relapse and progression free survival 0 20 40 60 80 100 10 -4 10 -3 10 -2 10 -1 10 0 time to relapse progression free survival (PFS) Therapy 1 (c = 0.5) Therapy 2 (c = 0.6) Therapy 3 (c = 0.7) Therapy 1 Therapy 2 Therapy 3 Chemo-sensitive Chemo-resistant t = 0 20 40 60 80 Figure 5.6: The eect of dose on tumor relapse and progression free survival under continuous therapy | Left: 3 therapies are simulated on identical initial condi- tions ( 2 = 0:020, 3 = 0:018, f = 10 3 , w = 0:1). Time to relapse signicantly increases with increasing dose while the progression free survival shows marginal, but decreasing, dierence. Right: the same 3 therapies are shown in a sh plot, where the simulation is stopped at the point of relapse to initial tumor size (now consisting entirely of chemo- resistant population). In order to better understand the dynamics of competitive release to anticipate the evolu- tionary trajectory of the tumor, dynamical state space diagrams are shown in Figure 5.7 for no therapy (c = 0; 5.7a) and with therapy (c = 0:6; top right triangle). A schematic ex- plaining the tri-linear coordinates (sometimes referred to as barycentric coordinates where Harnessing the evolutionary cost of resistance 121 P i x i = 1) of the triangular simplex (a representation of the phase space for every possible value of ~ x) is shown in gure 5.7c. The corners represent saturation of a single cell type (e.g. the top corner represents ~ x = [1; 0; 0], or all healthy cells). First, let us consider the no treatment dynamics: the top left triangle in gure 5.7a. The healthy (H; top corner), sensitive (S; bottom left corner), and resistant (R; bottom right cor- ner) populations compete according to equation 5.2 and follow trajectories shown (black). Instantaneous relative velocity is indicated by background color gradient (red to blue). All internal trajectories (pre-therapy) lead to tumor growth and eventual saturation of the sensitive population (bottom left corner of triangle). The resistant population nullcline (line of zero growth; _ x R = 0) is plotted (dashed green line). After an initial slow positive growth, the resistant population is suppressed by the faster growing sensitive population; all trajectories lead to the bottom left corner: saturation of the sensitive population. The healthy population nullcline (dashed blue line) is also plotted. The healthy nullcline is the dividing line of positive tumor growth and negative tumor growth (with the tumor consisting of both the sensitive and resistant populations). The sensitive nullcline (dashed red line) is on the boundary of the triangle, indicated there is no internal point for which there is negative sensitive growth. Chemotherapy alters the selection pressure to the disadvantage of chemo-sensitive cancer population and advantage of the healthy population (see equations 5.4, 5.5, and 5.6), altering dynamical outcomes: shown in the top right triangle in gure 5.7b. Internal solution trajectories (black) show initial trajectory toward healthy saturation (triangle top) but eventual relapse toward resistant population (bottom right of triangle) upon passing the resistant nullcline (dashed green line). Interestingly, the trajectories pass the resistant Harnessing the evolutionary cost of resistance 122 nullcline before the healthy nullcline (dashed blue line), indicating the resistant population may have positive growth ( _ x R > 0) even though the tumor has positive regression ( _ x H < 0), indicated as region 3. The sensitive nullcline is also plotted (dashed red line), showing initial conditions where tumor burden is too large and the dose concentration is inadequate to achieve positive regression of the sensitive population (regions 1, 5, and 6 such that _ x S > 0). The nullcline information will be used in the next section to determine timing schedules of adaptive therapy (see gure 5.8). 5.4.2 Treating the game: adaptive therapy In contrast to targeting a specic player in the game (i.e. heavily targeting the chemo- sensitive population), adaptive therapies take advantage of the important recognition that tness is contextual and changes during therapy or on drug holidays. The mechanism for control is also contextual: the suppression of the growth of resistant cell population occurs during periods of rest or weaker doses of therapy (drug-sensitive cells have a tness advantage in these conditions); suppression of the growth of the sensitive cell population occurs during treatment. Figure 5.8 shows how the idea of contextual tness can be applied to therapeutic strategies. A simple control paradigm is proposed to indirectly control the resistant population. Ther- apy targets only the chemo-sensitive cells, but the resistant population can be controlled by systematically choosing when to administer therapy and when to give drug holidays. Therapy \on" is for the purpose of killing sensitive cells. Therapy \o" is for the purpose of allowing a stable amount of sensitive cells to remain, in order to suppress the resistant Harnessing the evolutionary cost of resistance 123 Barycentric Coordinates Schematic length = 1 S H R 1 xR 1 xS 1 xH No therapy: c = 0 3 2 1 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ; ; ; 1 2 3 ; ; ; S H R With therapy: c = 0.6 4 3 2 1 6 5 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH > 0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xS < 0 ˙ xS < 0 ˙ xS < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 ˙ xR > 0 ; ; ; ; ; ; 1 2 3 4 5 6 ; ; ; ; ; ; S H R therapy ON; therapy OFF therapy ON; therapy OFF therapy ON; therapy OFF Resistant nullcline (therapy ON) Adaptive Therapy Schematic 1 2 3 ˙ xR < 0 ˙ xR > 0 Resistant nullcline (therapy OFF)˙ xR =0 ˙ xR =0 3 2 1 Initial condition at start of therapy Therapy ON trajectory (c > 0) Continuous therapy trajectory ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 Therapy OFF trajectory (c = 0) S H R low high low high (a) Barycentric Coordinates Schematic length = 1 S H R 1 xR 1 xS 1 xH No therapy: c = 0 3 2 1 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ; ; ; 1 2 3 ; ; ; S H R With therapy: c = 0.6 4 3 2 1 6 5 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH > 0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xS < 0 ˙ xS < 0 ˙ xS < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 ˙ xR > 0 ; ; ; ; ; ; 1 2 3 4 5 6 ; ; ; ; ; ; S H R therapy ON; therapy OFF therapy ON; therapy OFF therapy ON; therapy OFF Resistant nullcline (therapy ON) Adaptive Therapy Schematic 1 2 3 ˙ xR < 0 ˙ xR > 0 Resistant nullcline (therapy OFF)˙ xR =0 ˙ xR =0 3 2 1 Initial condition at start of therapy Therapy ON trajectory (c > 0) Continuous therapy trajectory ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 Therapy OFF trajectory (c = 0) S H R low high low high (b) length = 1 S H R Trilinear Coordinates Schematic x R x H x S (c) Barycentric Coordinates Schematic length = 1 S H R 1 xR 1 xS 1 xH No therapy: c = 0 3 2 1 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ; ; ; 1 2 3 ; ; ; S H R With therapy: c = 0.6 4 3 2 1 6 5 ˙ xR =0 ˙ xS =0 ˙ xH =0 ˙ xH > 0 ˙ xH > 0 ˙ xH > 0 ˙ xH < 0 ˙ xH < 0 ˙ xH < 0 ˙ xS > 0 ˙ xS > 0 ˙ xS > 0 ˙ xS < 0 ˙ xS < 0 ˙ xS < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 ˙ xR > 0 ; ; ; ; ; ; 1 2 3 4 5 6 ; ; ; ; ; ; S H R therapy ON; therapy OFF therapy ON; therapy OFF therapy ON; therapy OFF Resistant nullcline (therapy ON) Adaptive Therapy Schematic 1 2 3 ˙ xR < 0 ˙ xR > 0 Resistant nullcline (therapy OFF)˙ xR =0 ˙ xR =0 3 2 1 Initial condition at start of therapy Therapy ON trajectory (c > 0) Continuous therapy trajectory ˙ xR < 0 ˙ xR < 0 ˙ xR > 0 ˙ xR > 0 Therapy OFF trajectory (c = 0) S H R low high low high (d) Figure 5.7: Dynamic phase portraits before and during chemotherapy | Top left: before chemotherapy, the healthy (H), sensitive (S), and resistant (R) populations compete on a dynamical tness landscape, with several solution trajectories shown (black) and the instantaneous relative velocity indicated by background color gradient (red to blue). All internal trajectories lead to tumor growth and eventual saturation of the sen- sitive population (bottom left corner). Each population nullcline (line of zero growth: _ x i = 0) is plotted: healthy (dashed blue), sensitive (dashed red), and resistant (dashed green). Top right: chemotherapy alters the selection pressure to the disadvantage of chemo-sensitive cancer population and advantage of the healthy population (shown for c = 0:6: 2 = 0:020, 3 = 0:018, w = 0:1). Solution trajectories (black) show initial trajectory toward healthy saturation (triangle top) but eventual relapse toward resistant population (bottom right of triangle) upon passing the resistant nullcline. Bottom Left: a schematic of proposed adaptive therapy concept using the resistant nullclines to determine therapy \on" and \o" times in order to trap the tumor in the controllable region 2, and reach approximate cycle that repeats back on itself in red. The continuous therapy is also plotted in dashed blue, for comparison. Harnessing the evolutionary cost of resistance 124 population. The control paradigm is as follows: a continuous dose of therapy is adminis- tered until the nullcline ( _ x R = 0) is reached (see gure 5.7b). This is the starting point of positive growth for the resistant population (further therapy would result in _ x R > 0). At this point, a drug holiday (no therapy administered) is imposed until the second nullcline is reached (see gure 5.7a). The sensitive population is allowed to regrow until it is large enough to suppress the resistant population once again (and when _ x R = 0). Therapy is administered to allow the tumor to cycle back and forth between the two nullclines. Target- ing the resistant population for control allows an extension of relapse times by shaping the tness landscape at the right time in order to return to an amenable, controllable tumor burden and tumor composition. A schematic of this control paradigm is shown in gure 5.7d. The control paradigm is tested in gure 5.8 for identical initial conditions using a range of drug concentration dose values (low dose: blue; medium dose: red; high dose: yellow). The low and medium dose adaptive therapy strategies adequently outperform the continuous, constant dose (dashed lines). As seen in the shplots on the right, the resistant population (green) is suppressed during the \o" times of drug holidays, leading to an extended time without relapse for low and medium dose concentration (therapy 1 and 2). However, a higher dose (during therapy \on") leads to less tumor control and eventual relapse. The low- and medium-dose adaptive therapy strategies are successful for two reasons. First, the drug holidays allow an adequate sensitive population size to suppress the growth of the lower-tness resistant population. Second, the resistant population is never allowed to reach a positive growth ( _ x R > 0) during therapy \on" time periods. This positive growth region for the resistant population is the \point of no return" where the tumor will no Harnessing the evolutionary cost of resistance 125 0 50 100 150 t i m e ( m o n t h s ) 10 -3 10 -2 10 -1 10 0 v ( t ) 0 50 100 150 t i m e ( m o n t h s ) 10 -3 10 -2 10 -1 10 0 v ( t ) 0 50 100 150 t i m e ( m o n t h s ) 10 -3 10 -2 10 -1 10 0 v ( t ) Therapy 1: low dose (c = 0.4) Therapy 2: medium dose (c = 0.6) Therapy 3: high dose (c = 0.8) Therapy 1: low dose (c = 0.4) Therapy 2: medium dose (c = 0.6) Therapy 3: high dose (c = 0.8) t = 0 25 50 75 100 125 150 Sensitive Resistant Therapy 1 Therapy 2 Therapy 3 Figure 5.8: Adaptive therapy strategy to control resistant population | Left: Low dose (blue; top: c = 0:4), medium dose (red; middle: c = 0:6) and high dose (yellow; bottom: c = 0:8) are shown for continuous therapy (dashed line) and adaptive therapy (solid line). Right: the respective shplot is shown to the right of each therapy (low, medium, high) dose level. The adaptive therapy control paradigm administers therapy (ON) until the nullcline is reached (x R = 0; see gure 5.5, right). At this point, further therapy will result in the growth of the resistant population (i.e. x R > 0). Therapy is then turned OFF (a drug \holiday") until the second nullcline is reached (see gure 5.5, left), allowing the sensitive population to regain a stable and adequate size useful to suppress the resistant population. This control paradigm works well for low and medium dose, but stable control is not achieved for high dose. longer be controllable via the same drug. 5.5 Discussion Previous models of competitive release (e.g. equation 5.1) have been shown to be good de- scriptions of changing tumor size during therapy for many types of cancer [152], yet include no evolutionary concepts. As these models cannot be used to test adaptive chemotherapy drug schedules that are designed to anticipate and correct for the evolutionary trajectory of Harnessing the evolutionary cost of resistance 126 a tumor, a model was proposed here using a pairwise prisoner's dilemma in the replicator equation. Key parameters of the model include the initial growth rates of the chemo- sensitive and chemo-resistant population ( 2, 3), the cost of resistance ( 2 3), and the initial fraction of resistance. The most important parameter was determined to be the cost of resistance, with a small change exhibiting a large dierence in time to relapse. The key concept of the tness landscape in the replicator equation is that the drug, designed to kill proliferating chemo-sensitive cells, can be modeled as a change in the tness landscape via the selection pressure parameter, w. Resistant subclones, by denition, resist any change in tness due to the drug. Traditionally, chemotherapy drugs are tested by maximizing the kill rate, or maximizing the progression free survival time. Time to relapse is shown to be a more meaningful parameter to judge continuous chemotherapy drugs, even though all simulated tumors eventually relapse if there is pre-existing resistance at the start of therapy. Finally, an adaptive therapy concept was introduced to indirectly control the resistant population by leaving a stable tumor burden designed to suppress the resistant population. The adaptive therapy concept was informed by the resistant nullclines in the dynamics phase plane for the therapy \on" and \o" state space. The concept was shown to have good success at extending relapse times for low- or medium-dose therapies. Tumor adaptations only adapt to local and current selection forces and never anticipate the future. Here, oncologists and modelers may have an advantage if we can anticipate the evolutionary adaptations of the tumor and design therapies that exploit or mitigate the harmful eects of potential future adaptions. As the tumor is growing according to Darwinian principles, so must our treatment schedules be designed with these same Dar- winian principles in mind. Treatment must be dynamic and models must be evolutionary Harnessing the evolutionary cost of resistance 127 in scope. The goal here isn't to develop a new drug, but to optimize scheduling for existing drugs and better understand the dynamics of competitive release with the goal of control, and therefore, life extension. This life extension may give researchers the necessary time to improve existing drugs or develop new targeted drugs. Chapter 6 High Performance Computing Techniques The Center for High-Performance Computing (HPC) at the University of Southern Cali- fornia claims the 8th fastest supercomputer cluster in the United States, a resource made available to researchers at the university. The HPC is comprised of 2,225-node, 4-core, 6-core, and 12-core dual-processors manufactured by Dell, Oracle Sun, Hewlett-Packard and IBM connected on a 10-gigaabit Myrinet backbone. This computing cluster is ideal for the development of massive computational research jobs that require high-speed com- muncations among parallel computational elements [168]. In order to fully utilize the resources on the HPC for this research project, algorithms were implemented using the standard message passing system Message Passing Interface (MPI) and the shared memory multiprocessing API OpenMP. The purpose of using these 128 High Performance Computing Techniques 129 high performance computing resources is to scale simulations to large cell numbers, increase the speed of single stochastic simulations, and to scale the number of stochastic simulations output by the previously described models. 6.1 Hybrid MPI + OpenMP Parallelization In order to scale the number of cells in the single site, primary tumor ecology models, a hybrid message-passing and shared-memory (MPI + OpenMP) system was implemented. Figure 6.1 shows a schematic diagram of the process. During a single simulation, n p nodes are allocated, and n td threads are spawned on each of those nodes. The communication between nodes is handled by MPI, and the thread communication is handled by OpenMP. In Algorithm 1, the single thread, single node process for the primary tumor ecology model is shown (a moran process with mutations and selection). The procedure calculates the probability vectors for birth and death (p birth , p death ), chooses a single cell for birth and death each time step, and checks if the daughter cell incurred a mutation during the birth event during each time step, or cell division. The birth, death, and mutation processes are then parallelized on n p nodes (i.e. an MPI task) each with n td threads, as seen Algorithm 2. The single thread algorithm runs simul- taneously on each thread. The MPI Send command sends the p birth and p death vectors to each node. Next, each node forks into n td threads. Each thread simulates nThreadLoops High Performance Computing Techniques 130 Figure6.1: Schematic diagram of the Hybrid MPI + OpenMP Parallelization technique. In order to scale primary tumor growth simulations to large numbers of cells, a double Fork/Join process was used. MPI was used to fork to n p nodes and OpenMP was used to fork to n td threads. A join process for all nodes and threads is used in order to keep the possible error in the p birth and p death vectors to a minimum. birth, death and mutation processes without updating the p birth and p death vectors. Each thread keeps a local version of the change in state vector, x. 6.2 Error in Hybrid MPI + OpenMP Parallelization During the fork/join process, error in the p birth and p death vectors accumulates. This occurs because each node (and each subsequent thread) operates independently during the duration of the fork process. The state of the system (i.e. the total number of cancer cells) is tracked independently for each of the forks (a local system state), and the global state is estimated by the local system state on each fork. Not until the fork comes to an end at the High Performance Computing Techniques 131 Figure6.2: Strong scalability of Hybrid Parallelization | A sample primary tumor ecology model simulation using Algorithm 2 (Hybrid MPI + OpenMP) for N = 10 6 cells, mutation rate m = 0:2, and selection pressure w = 0:5. The simulations were performed on Dual Quadcore nodes with AMD Opteron 2.3 GHz processors with 16GB Memory. The simulation was duplicated with a single node (n p = 1, blue line) and a 2 node (n p = 1, red line) setup with varying thread count (n td = 1; 2; 4; 8). Left: total wall time required to nish the simulation (2:5e7 timesteps) for each value of n p and n td . Middle: speedup calculated from 6.1. Right: eciency of scaling, calculated from 6.2. Algorithm 1 Moran Process using Single-Thread Computation 1: procedure moran 2: calculate p birth and p death . birth and death probability vectors 3: while t<totalTime do 4: update p birth and p death 5: choose cell type for birth 6: check for mutation during birth process 7: choose cell type for death 8: update state vector, ~ x 9: t =t + 1 10: end while 11: end procedure MPI Allreduce step is the global state known to the each fork. Meanwhile, the local state is used in the calculation of the birth and death probability vectors. Between updates to the global state, the number of timesteps executed will ben total =n p n td nThreadLoops. In general, the longer the fork, the higher error accumulated. The error is also increased High Performance Computing Techniques 132 Algorithm 2 Moran Process using Mutli-Thread Computation 1: procedure moran 2: calculate p birth and p death . birth and death probability vectors 3: while t<totalTime do 4: update p birth and p death 5: MPI SEND p birth and p death to nodes 0 - n p 6: FORK into n td threads on each node 7: loop . each thread performs nThreadLoops 8: choose cell type for birth 9: check for mutation during birth process 10: choose cell type for death 11: update local ~ x 12: end loop 13: JOIN 14: MPI ALL REDUCE each node's ~ x local state vector to update the global state vector~ x 15: t =t + (n p n td nThreadLoops) 16: end while 17: end procedure for lower number of cells and higher number of threads or nodes. It is useful to understand how the error in the probability of birth and death for cancer cells changes during the fork. In Figure 6.3, this error is documented. The maximum error in the birth and death vectors can be easily computed by considering the case that maximizes p birth and p death in the n total timesteps. High Performance Computing Techniques 133 Figure 6.3: Error in the Hybrid MPI + OpenMP primary tumory growth simulation. Left: The probability of a cancer birth (blue) and cancer cell death (red) for w = 0:5 and N = 10 6 cells. A Fork/Join process for n p = 2, n td = 8 and nThreadLoops = 1000 was employed. Right: the maximum p birth and maximum p death possible during the fork process. High Performance Computing Techniques 134 6.3 Speedup and Eciency In order to measure the eectiveness of the hybrid parallelization of the primary tumor growth model, the speedup and eciency, speedup and eciency was plotted for xed- problem sized scaling. The speedup for n p MPI processes and n td OpenMP threads is shown in equation 6.1 and eciency is shown in equation 6.2. In Figure 6.2, the speedup and eciency is plotted. As the number of nodes or processes increase, the speedup departs from ideal scaling, with a maximum of 3 times speedup for n p = 2 and n td = 8. S p = T (n p ;n td = 1) T (n p ;n td =p) (6.1) E p = S p p (6.2) 6.4 Performance Tunability The performance of a Hybrid MPI + OpenMP algorithm can be tuned to nd the best combination of threads and nodes. Each simulation is run with xed-size problem (N = 10 6 cells), keeping the total number of processesP =n p n td constant. The results are shown in Table 6.1. The optimal combination is 1 node with 8 threads. It is important to note that most simulations with large number of cells (> 10 9 ) will require more than a combination of 8 threads and nodes, but the underlying principle here is that a maximum number of threads should be used, over a maximum of nodes. High Performance Computing Techniques 135 Table 6.1: Performance Tuning Number of OpenMP threads, n td Number of MPI processes, n p Wall Time, W T 1 8 1.4709 2 4 2.5882 4 2 4.5893 8 1 0.5955 Chapter 7 Future Work Mathematical models of the type developed in this thesis will play an increasingly important role in identifying and implementing novel chemotherapeutic strategies and will join the arsenal of tools that clinicians will turn to as they continue to develop strategies to ght cancer. Future progress will increasingly require the interactions of an interdisciplinary array of scientists, including applied mathematicians, to develop mathematical models of tumorigenesis and therapy response as a complex process [169]. Active collaboration among experimentalists, clinicians, and modelers will help produce an integrative and comprehensive view of the disease. With this in mind, we outline brie y some extensions and new developments we plan to pursue related to this thesis. 136 Future Work 137 7.1 Chemotherapy scheduling and metastasis As the role of specic genetic events that enable the tolerance and propagation of genome instability is elucidated, more data will become available to predict evolutionary trajectory and subsequent metastatic initiation and evolution [29]. Mathematical modeling will help uncover the role of competition as well as the role of mutations in driving or supporting metastasis. Current evolutionary models (such as those developed in this thesis) can be extended to include cell tracking through circulatory and lymphatic systems, integrating metastatic data [87, 170]. The successful development of targeted cancer drugs has been limited by drug resistant subclones, especially in a metastatic setting. Full cell-tracking models can address questions about the tracking of clones as well as the eect of drug treatment on circulating tumor cells. 7.2 Control theory methods to minimize stable tumor bur- den The replicator equation models competitive release [41, 116] as an evolutionary trade o of developing resistance mechanisms at the cost of the growth rates of resistant populations. The model is an ideal one for testing chemotherapeutic strategies designed to exploit this tness cost by allowing a stable burden of sensitive population that suppresses the growth of the resistant population. One such example of an adaptive therapy (see chapter 5) uses the nullclines to determine the time boundaries of bang-bang (on-o) control. Future Work 138 Further research can explore more complicated methods of control including varying the dose concentration as well as the length of drug holidays. Importantly, other groups have developed models of competitive release using spatial evolutionary models, which allows for drug penetration to be modelled [171]. Microenvironmental eects are shown to be important to the development of resistance, and can be considered in future models [172, 173]. The current model (chapter 5) assumes perfect information about the fraction of resistant population and the underlying phase portrait dynamics, an impractical assumption in clin- ical settings. Further research can integrate probabilistic models predicting the occurrence of resistance based on known or estimated mutation rates (see [56, 152, 174{176] Others have looked at optimal combinations or sequences of three or more drugs [34, 138, 177], or the evolutionary relationship of multiple treatments (see [139]). The current model (chapter 5) can be extended to include multiple non-cross resistant drugs [137, 166]. 7.3 Heterogeneity drives tumor growth There are a wide array of classical models used to describe tumor growth [96]. The underly- ing biological mechanisms of such growth models, often remain unclear [4, 178]. Kendal has developed a description of the role of heterogeneity in driving tumor growth using methods borrowed from statistical mechanics [85], which has also been reviewed in the context of models developed for this thesis (see [32]). These methods can be extended to modern growth models to better understand the role of heterogeneity in more sophisticated model paradigms. These statistical methods could be used to explore dierences between neutral Future Work 139 models of evolution [179], Darwinian models of evolution (linear and branched models) [150] and models of clonal interactions [29]. 7.4 Integrating preclinical and clinical trial data A changing paradigm shift has led to a change in the strategy of treating cancer [143] from \cure" to \control" paradigms. Evolutionary models have had success in predicting existence of an evolutionary double bind [26, 29, 139] in the order of administering multiple treatments. Other strategies target trunk mutation driver events that are present in every tumor cell, target parallel evolutionary events, forcing the tumor down a specic evolu- tionary path, resulting in acquired sensitivity (sequential therapy), and dynamic therapies that maintain a stable population of treatment-sensitive cells [150]. The challenge will be in integrating molecular, genomic, phylogenetic, and serum marker (e.g. PSA counts in prostate cancer) data to inform existing models and aid the de- velopment of new evolutionary and ecological models of tumor initiation, with a goal of predicting tumor adaptations and exploiting that prediction using fundamental evolution- ary principles. Data describing the evolutionary eect of cancer treatments are becoming increasingly available, in the form of preclinical trials [164, 165] and is currently extended to an on-going clinical trial (NCT02415621). Bibliography [1] Douglas Hanahan and Robert A Weinberg. The hallmarks of cancer. Cell, 100(1): 57{70, 2000. [2] Philip Gerlee. The model muddle: in search of tumor growth laws. Cancer Research, 73(8):2407{2411, 2013. [3] Michael Clinton Perry. The Chemotherapy Source Book. Lippincott Williams & Wilkins, 2008. [4] Anna Kane Laird. Dynamics of tumour growth. British journal of cancer, 18(3):490, 1964. [5] Benjamin Gompertz. On the nature of the function expressive of the law of hu- man mortality, and on a new mode of determining the value of life contingencies. Philosophical transactions of the Royal Society of London, 115:513{583, 1825. [6] Emil Frei and George P Canellos. Dose: a critical factor in cancer chemotherapy. The American Journal of Medicine, 69(4):585{594, 1980. [7] L Norton and R Simon. Tumor size, sensitivity to therapy, and design of treatment schedules. Cancer Treatment Reports, 61(7):1307, 1977. 140 Bibliography 141 [8] Howard E Skipper, Frank M Schabel Jr, and William S Wilcox. Experimental evalu- ation of potential anticancer agents. xxi. on the criteria and kinetics associated with curability of experimental leukemia. Cancer Chemotherapy Reports. Part 1, 35:1, 1964. [9] Howard E Skipper. Laboratory models: some historical perspective. Cancer Treat- ment Reports, 70(1):3{7, 1986. [10] Lauren MF Merlo, John W Pepper, Brian J Reid, and Carlo C Maley. Cancer as an evolutionary and ecological process. Nature Reviews Cancer, 6(12):924{935, 2006. [11] Camille Stephan-Otto Attolini and Franziska Michor. Evolutionary theory of cancer. Annals of the New York Academy of Sciences, 1168(1):23{51, 2009. [12] Steven A Frank. Dynamics of Cancer: Incidence, Inheritance, and Evolution. Prince- ton University Press, 2007. [13] Mel Greaves. Cancer: The Evolutionary Legacy. Oxford University Press, 2001. [14] Robert Axelrod, David E Axelrod, and Kenneth J Pienta. Evolution of cooperation among tumor cells. Proceedings of the National Academy of Sciences, 103(36):13474{ 13479, 2006. [15] Martin A Nowak. Evolutionary Dynamics. Harvard University Press, 2006. [16] IPM Tomlinson. Game-theory models of interactions between tumour cells. European Journal of Cancer, 33(9):1495{1500, 1997. Bibliography 142 [17] Sabine Hummert, Katrin Bohl, David Basanta, Andreas Deutsch, Sarah Werner, G unter Theien, Anja Schroeter, and Stefan Schuster. Evolutionary game theory: cells as players. Molecular BioSystems, 10(12):3044{3065, 2014. [18] Thomas L Vincent. Carcinogenesis as an evolutionary game. Advances in Complex Systems, 9(04):369{382, 2006. [19] Franziska Michor, Yoh Iwasa, and Martin A Nowak. Dynamics of cancer progression. Nature Reviews Cancer, 4(3):197{205, 2004. [20] Martin A Nowak and Karl Sigmund. Evolutionary dynamics of biological games. Science, 303(5659):793{799, 2004. [21] John von Neumann. Sur la th eorie des jeux. Comptes Rendus de l'Acad emie des Sciences, 186(25):1689{91, 1928. [22] John Maynard Smith. Evolution and the Theory of Games. Cambridge University Press, 1982. [23] William Poundstone. Prisoner's Dilemma. Anchor, 2011. [24] Peter C Nowell. The clonal evolution of tumor cell populations. Science, 194(4260): 23{28, 1976. [25] Mel Greaves and Carlo C Maley. Clonal evolution in cancer. Nature, 481(7381): 306{313, 2012. [26] David Basanta and Alexander RA Anderson. Exploiting ecological principles to better understand cancer progression and treatment. Interface focus, 3(4):20130020, 2013. Bibliography 143 [27] Robert A Gatenby, Robert J Gillies, and Joel S Brown. Of cancer and cave sh. Nature Reviews Cancer, 11(4):237{238, 2011. [28] Thomas S Deisboeck, Zhihui Wang, Paul Macklin, and Vittorio Cristini. Multiscale cancer modeling. Annual review of biomedical engineering, 13:127{155, 2011. [29] David Basanta and Alexander Anderson. Homeostasis back and forth: An eco- evolutionary perspective of cancer. bioRxiv, page 092023, 2016. [30] S Venkatesan and C Swanton. Tumor evolutionary principles: How intratumor het- erogeneity in uences cancer treatment and outcome. In American Society of Clini- cal Oncology educational book. American Society of Clinical Oncology. Meeting, vol- ume 35, pages e141{9, 2016. [31] Robert J Gillies, Daniel Verduzco, and Robert A Gatenby. Evolutionary dynamics of carcinogenesis and why targeted therapy does not work. Nature Reviews Cancer, 12(7):487{493, 2012. [32] J West, Z Hasnain, P Macklin, J Mason, and PK Newton. An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity and gompertzian growth. SIAM Review, 58(4):716{736, 2016. [33] Robert A Gatenby and Thomas L Vincent. Application of quantitative models from population biology and evolutionary game theory to tumor therapeutic strategies. Molecular cancer therapeutics, 2(9):919{927, 2003. [34] Ivana Bozic, Johannes G Reiter, Benjamin Allen, Tibor Antal, Krishnendu Chat- terjee, Preya Shah, Yo Sup Moon, Amin Yaqubie, Nicole Kelly, Dung T Le, et al. Bibliography 144 Evolutionary dynamics of cancer in response to targeted combination therapy. Elife, 2:e00747, 2013. [35] Rick Durrett, Jasmine Foo, Kevin Leder, John Mayberry, and Franziska Michor. Intratumor heterogeneity in evolutionary models of tumor progression. Genetics, 188(2):461{477, 2011. [36] Rebecca A Burrell, Nicholas McGranahan, Jiri Bartek, and Charles Swanton. The causes and consequences of genetic heterogeneity in cancer evolution. Nature, 501 (7467):338{345, 2013. [37] Thomas M Cover and Joy A Thomas. Elements of Information Theory. John Wiley & Sons, 2012. [38] Douglas Hanahan, Gabriele Bergers, and Emily Bergsland. Less is more, regularly: metronomic dosing of cytotoxic drugs can target tumor angiogenesis in mice. The Journal of Clinical Investigation, 105(8):1045{1047, 2000. [39] K Lien, S Georgsdottir, L Sivanathan, K Chan, and U Emmenegger. Low-dose metronomic chemotherapy: a systematic literature analysis. European Journal of Cancer, 49(16):3387{3395, 2013. [40] Jerey West and P.K. Newton. Chemotherapeutic dose scheduling based on tu- mor growth rates: The case for low dose metronomic high entropy therapies. USC Preprint, 2017. [41] Jerey West, Yongqian Ma, and Paul K Newton. Jerey west, yongqian ma, paul k. newton Bibliography 145 harnessing the evolutionary cost of chemotherapeutic resistance by shaping the tness landscape of a tumor. USC Preprint, 2017. [42] Robert A Gatenby, Ariosto S Silva, Robert J Gillies, and B Roy Frieden. Adaptive therapy. Cancer Research, 69(11):4894{4903, 2009. [43] Marc Crommelinck, Bernard Feltz, and Philippe Goujon. Self-organization and emer- gence in life sciences. Springer, 2006. [44] AR Kansal, S Torquato, EA Chiocca, and TS Deisboeck. Emergence of a subpop- ulation in a computational model of tumor growth. Journal of Theoretical Biology, 207(3):431{441, 2000. [45] Jerey West, Zaki Hasnain, Jeremy Mason, and Paul K Newton. The prisoner's dilemma as a cancer model. Convergent Science Physical Oncology, 2(3):035002, 2016. [46] Bernard Crespi and Kyle Summers. Evolutionary biology of cancer. Trends in Ecology & Evolution, 20(10):545{552, 2005. [47] Steven A Frank. Dynamics of cancer: incidence, inheritance, and evolution. Prince- ton University Press, 2007. [48] Charles Swanton. Intratumor heterogeneity: evolution through space and time. Can- cer Research, 72(19):4875{4882, 2012. [49] Robert Weinberg. The biology of cancer. Garland science, 2013. [50] I~ nigo Martincorena and Peter J Campbell. Somatic mutation in cancer and normal cells. Science, 349(6255):1483{1489, 2015. Bibliography 146 [51] Lucy R Yates and Peter J Campbell. Evolution of the cancer genome. Nature Reviews Genetics, 13(11):795{806, 2012. [52] A Farrell, E Hutchinsin, B Marte, and N McCarthy. Nature milestones: cancer. Nat, 440:S7{S23, 2006. [53] Robert Allan Weinberg. One renegade cell: how cancer begins. Basic Books, 1998. [54] John S Spratt and John A Spratt. What is breast cancer doing before we can detect it? Journal of surgical oncology, 30(3):156{160, 1985. [55] James H Goldie and Andrew J Coldman. Drug resistance in cancer: mechanisms and models. Cambridge University Press, 2009. [56] Yoh Iwasa, Martin A Nowak, and Franziska Michor. Evolution of resistance during clonal expansion. Genetics, 172(4):2557{2566, 2006. [57] John Von Neumann and Oskar Morgenstern. Theory of games and economic behavior. Princeton university press, 2007. [58] J Maynard Smith. Evolutionary game theory. In Vito Volterra symposium on math- ematical models in biology, pages 73{81. Springer, 1980. [59] M Nowak. Evolutionarydynamics: Exploringtheequationso ife. Cambridge: Harvar- dUniversityPress, 2006. [60] Martin A Nowak and Karl Sigmund. Evolutionary dynamics of biological games. science, 303(5659):793{799, 2004. Bibliography 147 [61] Stefan Schuster, Jan-Ulrich Kreft, Anja Schroeter, and Thomas Pfeier. Use of game- theoretical methods in biochemistry and biophysics. Journal of biological physics, 34 (1-2):1{17, 2008. [62] John F Nash et al. Equilibrium points in n-person games. Proceedings of the national academy of sciences, 36(1):48{49, 1950. [63] John Nash. Non-cooperative games. Annals of mathematics, pages 286{295, 1951. [64] Josef Hofbauer and Karl Sigmund. Evolutionary Games and Population Dynamics. Cambridge University Press, 1998. [65] Robert M Axelrod. The evolution of cooperation. Basic books, 2006. [66] David Basanta and Andreas Deutsch. A game theoretical perspective on the so- matic evolution of cancer. Modeling and Simulation in Science, Engineering and Technology, page 97, 2008. [67] Michael Doebeli and Christoph Hauert. Models of cooperation based on the prisoner's dilemma and the snowdrift game. Ecology letters, 8(7):748{766, 2005. [68] Michael Doebeli, Christoph Hauert, and Timothy Killingback. The evolutionary origin of cooperators and defectors. science, 306(5697):859{862, 2004. [69] Robert A Gatenby and Thomas L Vincent. An evolutionary model of carcinogenesis. Cancer Research, 63(19):6212{6220, 2003. [70] Robert A Gatenby and Thomas L Vincent. Application of quantitative models from population biology and evolutionary game theory to tumor therapeutic strategies. Molecular cancer therapeutics, 2(9):919{927, 2003. Bibliography 148 [71] C Hauert. Evolution from cellular to social scales. Evolutionary Dynamics, pages 11{44, 2008. [72] Christoph Hauert and Gy orgy Szab o. Game theory and physics. American Journal of Physics, 73(5):405{414, 2005. [73] Katrin Bohl, Sabine Hummert, Sarah Werner, David Basanta, Andreas Deutsch, Stefan Schuster, G unter Theien, and Anja Schroeter. Evolutionary game theory: molecules as players. Molecular BioSystems, 10(12):3066{3074, 2014. [74] Irina Kareva. Prisoner's dilemma in cancer metabolism. PloS one, 6(12):e28576, 2011. [75] Martin Nowak. Stochastic strategies in the prisoner's dilemma. Theoretical Popula- tion Biology, 38(1):93{112, 1990. [76] Martin A Nowak and Robert M May. Evolutionary games and spatial chaos. Nature, 359(6398):826{829, 1992. [77] Arne Traulsen and Christoph Hauert. Stochastic evolutionary game dynamics. Re- views of nonlinear dynamics and complexity, 2:25{61, 2009. [78] J orgen W Weibull. Evolutionary game theory. MIT press, 1997. [79] Richard Dawkins. The selsh gene. Oxford university press, 2016. [80] Zhen Wang, Satoshi Kokubo, Marko Jusup, and Jun Tanimoto. Universal scaling for the dilemma strength in evolutionary games. Physics of life reviews, 14:1{30, 2015. [81] Erez Lieberman, Christoph Hauert, and Martin A Nowak. Evolutionary dynamics on graphs. Nature, 433(7023):312{316, 2005. Bibliography 149 [82] Jerey West, Zaki Hasnain, Paul Macklin, and Paul K Newton. An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity driving gompertzian growth. SIAM Review, 58(4):716{736, 2016. [83] Anna Kane Laird. Dynamics of tumour growth: comparison of growth rates and extrapolation of growth curve to one cell. British Journal of Cancer, 19(2):278, 1965. [84] Steven A Frank and Martin A Nowak. Problems of somatic mutation and cancer. Bioessays, 26(3):291{299, 2004. [85] WS Kendal. Gompertzian growth as a consequence of tumor heterogeneity. Mathe- matical Biosciences, 73(1):103{107, 1985. [86] Thomas M Cover and Joy A Thomas. Elements of information theory. John Wiley & Sons, 2012. [87] Paul K Newton, Jeremy Mason, Brian Hurt, Kelly Bethel, Lyudmila Bazhenova, Jorge Nieva, and Peter Kuhn. Entropy, complexity, and markov diagrams for random walk cancer models. Scientic reports, 4, 2014. [88] Larry Norton and Richard Simon. Growth curve of an experimental solid tumor following radiotherapy. Journal of National Cancer Institute, 58(6), 1977. [89] Peter Dear. The intelligibility of nature: How science makes sense of the world. University of Chicago Press, 2008. [90] EP Wigner. The unreasonable eectiveness of mathematics in the natural sciences. In Philosophical Re ections and Syntheses, pages 534{549. Springer, 1995. Bibliography 150 [91] CO Nordling. A new theory on the cancer-inducing mechanism. British journal of cancer, 7(1):68, 1953. [92] Alfred G Knudson. Mutation and cancer: statistical study of retinoblastoma. Pro- ceedings of the National Academy of Sciences, 68(4):820{823, 1971. [93] Ian PM Tomlinson, MR Novelli, and WF Bodmer. The mutation rate and cancer. Proceedings of the National Academy of Sciences, 93(25):14800{14803, 1996. [94] Charles M Rubin. The genetic basis of human cancer. Annals of Internal Medicine, 129(9):759{759, 1998. [95] Sten Friberg and Stefan Mattson. On the growth rates of human malignant tumors: implications for medical decision making. Journal of surgical oncology, 65(4):284{297, 1997. [96] S ebastien Benzekry, Clare Lamont, Afshin Beheshti, Amanda Tracz, John ML Ebos, Lynn Hlatky, and Philip Hahnfeldt. Classical mathematical models for description and prediction of experimental tumor growth. PLoS Comput Biol, 10(8):e1003800, 2014. [97] Richard Durrett and Stephen Moseley. Evolution of resistance and progression to disease during clonal expansion of cancer. Theoretical population biology, 77(1):42{ 48, 2010. [98] Natalia L Komarova and Dominik Wodarz. Drug resistance in cancer: principles of emergence and prevention. Proceedings of the National Academy of Sciences of the United States of America, 102(27):9714{9719, 2005. Bibliography 151 [99] Natalia Komarova. Stochastic modeling of drug resistance in cancer. Journal of theoretical biology, 239(3):351{366, 2006. [100] James H Goldie and Andrew J Coldman. Drug resistance in cancer: mechanisms and models. Cambridge University Press, 2009. [101] Nicholas Navin, Alexander Krasnitz, Linda Rodgers, Kerry Cook, Jennifer Meth, Jude Kendall, Michael Riggs, Yvonne Eberling, Jennifer Troge, Vladimir Grubor, et al. Inferring tumor progression from genomic heterogeneity. Genome research, 20 (1):68{80, 2010. [102] Andriy Marusyk and Kornelia Polyak. Tumor heterogeneity: causes and conse- quences. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer, 1805(1):105{117, 2010. [103] Robert G Abbott, Stephanie Forrest, and Kenneth J Pienta. Simulating the hall- marks of cancer. Articial Life, 12(4):617{634, 2006. [104] Gloria H Heppner. Tumor heterogeneity. Cancer research, 44(6):2259{2265, 1984. [105] Michail Shipitsin, Lauren L Campbell, Pedram Argani, Stanislawa Weremowicz, Noga Bloushtain-Qimron, Jun Yao, Tatiana Nikolskaya, Tatiana Serebryiskaya, Rameen Beroukhim, Min Hu, et al. Molecular denition of breast tumor hetero- geneity. Cancer cell, 11(3):259{273, 2007. [106] Marco Gerlinger, Andrew J Rowan, Stuart Horswell, James Larkin, David Endes- felder, Eva Gronroos, Pierre Martinez, Nicholas Matthews, Aengus Stewart, Patrick Bibliography 152 Tarpey, et al. Intratumor heterogeneity and branched evolution revealed by multire- gion sequencing. N Engl j Med, 2012(366):883{892, 2012. [107] Gloria H Heppner and Bonnie E Miller. Tumor heterogeneity: biological implications and therapeutic consequences. Cancer and Metastasis Reviews, 2(1):5{23, 1983. [108] John Adam and Nicola Bellomo. A survey of models for tumor-immune system dynamics. Springer Science & Business Media, 2012. [109] Georey B West, James H Brown, and Brian J Enquist. A general model for onto- genetic growth. Nature, 413(6856):628{631, 2001. [110] Robyn P Araujo and DL Sean McElwain. A history of the study of solid tumour growth: the contribution of mathematical modelling. Bulletin of mathematical biol- ogy, 66(5):1039, 2004. [111] Randall T Schuh. Biological systematics: principles and applications. Cornell Uni- versity Press, 2000. [112] Tiany A Traina, Ute Dugan, Brian Higgins, Kenneth Kolinsky, Maria Theodoulou, Cliord A Hudis, and Larry Norton. Optimizing chemotherapy dose and schedule by norton-simon mathematical modeling. Breast Disease, 31(1):7{18, 2010. [113] LH Abbott and F Michor. Mathematical models of targeted cancer therapy. British Journal of Cancer, 95(9):1136{1141, 2006. [114] Paul Macklin and John Lowengrub. Nonlinear simulation of the eect of microenvi- ronment on tumor growth. Journal of theoretical biology, 245(4):677{704, 2007. Bibliography 153 [115] Sheldon M Ross. Introduction to Probability Models: Solutions Manual for Introduc- tion to Probability Models. Solu, 4. Acad. Press, 1989. [116] Arne Traulsen, Jens Christian Claussen, and Christoph Hauert. Coevolutionary dynamics: from nite to innite populations. Physical review letters, 95(23):238701, 2005. [117] Katarzyna A Rejniak and Alexander RA Anderson. Hybrid models of tumor growth. Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 3(1):115{125, 2011. [118] Vittorio Cristini, John Lowengrub, and Qing Nie. Nonlinear simulation of tumor growth. Journal of mathematical biology, 46(3):191{224, 2003. [119] Steven M Wise, John S Lowengrub, Hermann B Frieboes, and Vittorio Cristini. Three-dimensional multispecies nonlinear tumor growth|i: model and numerical method. Journal of theoretical biology, 253(3):524{543, 2008. [120] Hermann B Frieboes, Fang Jin, Yao-Li Chuang, Steven M Wise, John S Lowengrub, and Vittorio Cristini. Three-dimensional multispecies nonlinear tumor growth|ii: tumor invasion and angiogenesis. Journal of theoretical biology, 264(4):1254{1278, 2010. [121] Yangjin Kim, Magdalena A Stolarska, and Hans G Othmer. The role of the mi- croenvironment in tumor growth and invasion. Progress in biophysics and molecular biology, 106(2):353{379, 2011. [122] CL Frenzen and JD Murray. A cell kinetics justication for gompertz'equation. SIAM Journal on Applied Mathematics, 46(4):614{629, 1986. Bibliography 154 [123] Frank Kozusko and Zeljko Bajzer. Combining gompertzian growth and cell popula- tion dynamics. Mathematical Biosciences, 185(2):153{167, 2003. [124] P Tracqui. Biophysical models of tumour growth. Reports on Progress in Physics, 72(5):056701, 2009. [125] Thomas Robert Malthus. An essay on the principle of population: or, A view of its past and present eects on human happiness. Reeves & Turner, 1888. [126] Robert S Kerbel and Barton A Kamen. The anti-angiogenic basis of metronomic chemotherapy. Nature Reviews Cancer, 4(6):423{436, 2004. [127] Timothy Browder, Catherine E Buttereld, Birgit M Kr aling, Bin Shi, Blair Marshall, Michael S O'Reilly, and Judah Folkman. Antiangiogenic scheduling of chemotherapy improves ecacy against experimental drug-resistant cancer. Cancer Research, 60 (7):1878{1886, 2000. [128] Giannoula Klement, Sylvain Baruchel, Janusz Rak, Shan Man, Katherine Clark, Daniel J Hicklin, Peter Bohlen, and Robert S Kerbel. Continuous low-dose therapy with vinblastine and vegf receptor-2 antibody induces sustained tumor regression without overt toxicity. The Journal of Clinical Investigation, 105(8):R15{R24, 2000. [129] J West, Z Hasnain, J Mason, and PK Newton. The prisoner's dilemma as a cancer model. Converg. Sci. Phys. Oncol., 2(3), 2016. [130] S ebastien Benzekry, Eddy Pasquier, Dominique Barbolosi, Bruno Lacarelle, Fabrice Barl esi, Nicolas Andr e, and Joseph Ciccolini. Metronomic reloaded: Theoretical Bibliography 155 models bringing chemotherapy into the era of precision medicine. Seminars in Cancer Biology, 35:53{61, 2015. [131] Dominik Wodarz and Natalia L Komarova. Dynamics of Cancer: Mathematical Foundations of Oncology. World Scientic, 2014. [132] Gerhard Held, Joerg Schubert, Marcel Reiser, Michael Pfreundschuh, German High- Grade Non-Hodgkin-Lymphoma, and Study Group. Dose-intensied treatment of advanced-stage diuse large b-cell lymphomas. Seminars in Hematology, 43(4):221{ 229, 2006. [133] Sotiris Prokopiou, Eduardo G Moros, Jan Poleszczuk, Jimmy Caudell, Javier F Torres-Roca, Kujtim Lati, Jae K Lee, Robert Myerson, Louis B Harrison, and Heiko Enderling. A proliferation saturation index to predict radiation response and personalize radiotherapy fractionation. Radiation Oncology, 10(1):1, 2015. [134] Philip Gerlee and ARA Anderson. An evolutionary hybrid cellular automaton model of solid tumor growth. J. Theor. Bio., 246(4):583{603, 2007. [135] Michael Doebeli and Christoph Hauert. Models of cooperation based on the prisoner's dilemma and the snowdrift game. Ecology Letters, 8(7):748{766, 2005. [136] Michael Doebeli, Christoph Hauert, and Timothy Killingback. The evolutionary origin of cooperators and defectors. Science, 306(5697):859{862, 2004. [137] Robert A Beckman, Gunter S Schemmann, and Chen-Hsiang Yeang. Impact of genetic dynamics and single-cell heterogeneity on development of nonstandard per- sonalized medicine strategies for cancer. Proceedings of the National Academy of Bibliography 156 Sciences, 109(36):14586{14591, 2012. [138] Jasmine Foo and Franziska Michor. Evolution of resistance to anti-cancer therapy during general dosing schedules. Journal of Theoretical Biology, 263(2):179{188, 2010. [139] David Basanta, Robert A Gatenby, and Alexander RA Anderson. Exploiting evolu- tion to treat drug resistance: combination therapy and the double bind. Molecular Pharmaceutics, 9(4):914{921, 2012. [140] David Liao, Luis Est evez-Salmer on, and Thea D Tlsty. Conceptualizing a tool to optimize therapy based on dynamic heterogeneity. Physical Biology, 9(6):065005{18, November 2012. [141] RB Martin. Optimal control drug scheduling of cancer chemotherapy. Automatica, 28(6):1113{1123, 1992. [142] Charalambos Loizides, Demetris Iacovides, Marios M Hadjiandreou, Gizem Rizki, Achilleas Achilleos, Katerina Strati, and Georgios D Mitsis. Model-based tumor growth dynamics and therapy response in a mouse model of de novo carcinogenesis. PloS one, 10(12):e0143840, 2015. [143] Robert A Gatenby. A change of strategy in the war on cancer. Nature, 459(7246): 508{509, 2009. [144] Ian F Tannock et al. Docetaxel plus prednisone or mitoxantrone plus prednisone for advanced prostate cancer. The New Eng. J. of Med., 351(15):1502{1512, 2004. Bibliography 157 [145] David R Spigel et al. A phase 2 study of higher dose weekly topotecan in relapsed small-cell lung cancer. Clinical Lung Cancer, 12(3):187{191, 2011. [146] Maxwell Lewis Neal, Andrew D Trister, Tyler Cloke, Rita Sodt, Sunyoung Ahn, Anne L Baldock, Carly A Bridge, Albert Lai, Timothy F Cloughesy, Maciej M Mru- gala, et al. Discriminating survival outcomes in patients with glioblastoma using a simulation-based, patient-specic response metric. PLoS One, 8(1):e51951, 2013. [147] Joseph H Connell. The in uence of interspecic competition and other factors on the distribution of the barnacle chthamalus stellatus. Ecology, 42(4):710{723, 1961. [148] Peter R Grant. Convergent and divergent character displacement. Biological journal of the Linnean Society, 4(1):39{68, 1972. [149] William L Brown and Edward O Wilson. Character displacement. Systematic Zool- ogy, 5(2):49{64, 1956. [150] Charles Swanton. Tumor evolutionary principles: How intratumor heterogeneity in uences cancer treatment and outcome. American Society of Clinical Oncology, 2016. [151] C Athena Aktipis, Virginia SY Kwan, Kathryn A Johnson, Steven L Neuberg, and Carlo C Maley. Overlooking evolution: a systematic analysis of cancer relapse and therapeutic resistance research. PloS One, 6(11):e26100, 2011. [152] Ivana Bozic and Martin A Nowak. Resisting resistance. Annual Review of Cancer Biology, 2017. Bibliography 158 [153] K-A Kreuzer, P Le Coutre, O Landt, I-K Na, M Schwarz, K Schultheis, A Hochhaus, and B D orken. Preexistence and evolution of imatinib mesylate-resistant clones in chronic myelogenous leukemia detected by a pna-based pcr clamping technique. An- nals of Hematology, 82(5):284{289, 2003. [154] Catherine Roche-Lestienne and Claude Preudhomme. Mutations in the abl kinase domain pre-exist the onset of imatinib treatment. In Seminars in Hematology, vol- ume 40, pages 80{82. Elsevier, 2003. [155] Kristel Kemper, Oscar Krijgsman, Paulien Cornelissen-Steijger, Aida Shahrabi, Fleur Weeber, Ji-Ying Song, Thomas Kuilman, Daniel J Vis, Lodewyk F Wessels, Emile E Voest, et al. Intra-and inter-tumor heterogeneity in a vemurafenib-resistant melanoma patient and derived xenografts. EMBO Molecular Medicine, 7(9):1104{ 1118, 2015. [156] Alessandro Romanel, Delila Gasi Tandefelt, Vincenza Conteduca, Anuradha Ja- yaram, Nicola Casiraghi, Daniel Wetterskog, Samanta Salvi, Dino Amadori, Zafeiris Zafeiriou, Pasquale Rescigno, et al. Plasma ar and abiraterone-resistant prostate cancer. Science Translational Medicine, 7(312):312re10{312re10, 2015. [157] Luis A Diaz Jr, Richard T Williams, Jian Wu, Isaac Kinde, J Randolph Hecht, Jordan Berlin, Benjamin Allen, Ivana Bozic, Johannes G Reiter, Martin A Nowak, et al. The molecular evolution of acquired resistance to targeted egfr blockade in colorectal cancers. Nature, 486(7404):537{540, 2012. [158] Pierre Laurent-Puig, Deniz Pekin, Corinne Normand, Steve K Kotsopoulos, Philippe Nizard, Karla Perez-Toralla, Rachel Rowell, Je Olson, Preethi Srinivasan, Delphine Bibliography 159 Le Corre, et al. Clinical relevance of kras-mutated subclones detected with pico- droplet digital pcr in advanced colorectal cancer treated with anti-egfr therapy. Clin- ical Cancer Research, 2014. [159] Roland F Schwarz, Charlotte KY Ng, Susanna L Cooke, Scott Newman, Jillian Temple, Anna M Piskorz, Davina Gale, Karen Sayal, Muhammed Murtaza, Peter J Baldwin, et al. Spatial and temporal heterogeneity in high-grade serous ovarian cancer: A phylogenetic analysis. PLoS Med, 12(2):e1001789, 2015. [160] A Sorana Morrissy, Livia Garzia, David JH Shih, Scott Zuyderduyn, Xi Huang, Patryk Skowron, Marc Remke, Florence MG Cavalli, Vijay Ramaswamy, Patricia E Lindsay, et al. Divergent clonal selection dominates medulloblastoma at recurrence. Nature, 529(7586):351{357, 2016. [161] Pedro M Enriquez-Navas, Jonathan W Wojtkowiak, and Robert A Gatenby. Applica- tion of evolutionary principles to cancer therapy. Cancer research, 75(22):4675{4680, 2015. [162] Ivana Bozic and Martin A Nowak. Timing and heterogeneity of mutations associated with drug resistance in metastatic cancers. Proceedings of the National Academy of Sciences, 111(45):15964{15968, 2014. [163] Li Ding, Timothy J Ley, David E Larson, Christopher A Miller, Daniel C Koboldt, John S Welch, Julie K Ritchey, Margaret A Young, Tamara Lamprecht, Michael D McLellan, et al. Clonal evolution in relapsed acute myeloid leukaemia revealed by whole-genome sequencing. Nature, 481(7382):506{510, 2012. Bibliography 160 [164] Pedro M Enriquez-Navas, Yoonseok Kam, Tuhin Das, Sabrina Hassan, Ariosto Silva, Parastou Foroutan, Epifanio Ruiz, Gary Martinez, Susan Minton, Robert J Gillies, et al. Exploiting evolutionary principles to prolong tumor control in preclinical models of breast cancer. Science Translational Medicine, 8(327):327ra24{327ra24, 2016. [165] Sarah Seton-Rogers. Chemotherapy: Preventing competitive release. Nature Reviews Cancer, 16(4):199{199, 2016. [166] Robert A Beckman and Chen-Hsiang Yeang. Nonstandard personalized medicine strategies for cancer may lead to improved patient outcomes. Personalized Medicine, 11(7):705{719, 2014. [167] Artem Kaznatcheev, Robert Vander Velde, Jacob G Scott, and David Basanta. Can- cer treatment scheduling and dynamic heterogeneity in social dilemmas of tumour acidity and vasculature. British Journal of Cancer, 116(6):785{792, 2017. [168] University of Southern California. Hpc infrastruction | hpc | usc, 2015. URL https://hpcc.usc.edu/about/hpcc-infrastructure/. [Online; accessed 6-July- 2015]. [169] Alexander RA Anderson and Vito Quaranta. Integrative mathematical oncology. Nature Reviews Cancer, 8(3):227{234, 2008. [170] Paul K Newton, Jeremy Mason, Kelly Bethel, Lyudmila A Bazhenova, Jorge Nieva, and Peter Kuhn. A stochastic markov chain model to describe lung cancer growth and metastasis. PloS one, 7(4):e34637, 2012. Bibliography 161 [171] Jill A Gallaher, Pedro M Enriquez-Navas, Kimberly A Luddy, Robert A Gatenby, and Alexander RA Anderson. Adaptive therapy for heterogeneous cancer: Exploiting space and trade-os in drug scheduling. bioRxiv, page 128959, 2017. [172] Olivier Tr edan, Carlos M Galmarini, Krupa Patel, and Ian F Tannock. Drug re- sistance and the solid tumor microenvironment. Journal of the National Cancer Institute, 99(19):1441{1454, 2007. [173] Mark B Meads, Robert A Gatenby, and William S Dalton. Environment-mediated drug resistance: a major contributor to minimal residual disease. Nature reviews cancer, 9(9):665{674, 2009. [174] Cristian Tomasetti and Doron Levy. An elementary approach to modeling drug resistance in cancer. Mathematical biosciences and engineering: MBE, 7(4):905, 2010. [175] Franziska Michor, Timothy P Hughes, Yoh Iwasa, Susan Branford, Neil P Shah, Charles L Sawyers, and Martin A Nowak. Dynamics of chronic myeloid leukaemia. Nature, 435(7046):1267{1270, 2005. [176] Natalia L Komarova, Lin Wu, and Pierre Baldi. The xed-size luria{delbruck model with a nonzero death rate. Mathematical biosciences, 210(1):253{290, 2007. [177] Jasmine Foo and Franziska Michor. Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies. PLoS Compututa- tional Biology, 5(11):e1000557, 2009. Bibliography 162 [178] Peter W Sullivan and Sydney E Salmon. Kinetics of tumor growth and regression in igg multiple myeloma. Journal of Clinical Investigation, 51(7):1697, 1972. [179] Andrea Sottoriva, Haeyoun Kang, Zhicheng Ma, Trevor A Graham, Matthew P Sa- lomon, Junsong Zhao, Paul Marjoram, Kimberly Siegmund, Michael F Press, Darryl Shibata, et al. A big bang model of human colorectal tumor growth. Nature genetics, 47(3):209{216, 2015.
Abstract (if available)
Abstract
Tumor development is an evolutionary process in which a heterogeneous population of cells with different growth capabilities compete for resources in order to gain a proliferative advantage. We describe a cell-molecular based evolutionary mathematical model of tumor development driven by a stochastic Moran birth-death process, in which cancer cells and healthy cells compete for dominance in the population. Each are assigned payoffs according to a prisoner’s dilemma evolutionary game where the healthy cells are the cooperators and the cancer cells are the defectors. The cells in the tumor carry molecular information in the form of a numerical genome which we represent as a four-digit binary string used to differentiate cells into 16 molecular types. The binary string is able to undergo stochastic point mutations that are passed to a daughter cell after each birth event. The value of the binary string determines the cell fitness, with lower fit cells (e.g. 0000) defined as healthy phenotypes, and higher fit cells (e.g. 1111) defined as malignant phenotypes. The model is used to explore key emergent features associated with tumor development, including tumor growth rates (e.g. Gompertzian growth) as it relates to intratumor molecular heterogeneity (e.g. heterogeneity is a driver of tumor growth). ❧ Next, the model is extended to include tumor regression due to chemotherapy, allowing for the design of chemotherapeutic “strategies,” or schedules, tailored to different tumor growth characteristics. Using the Shannon entropy as a novel tool to quantify dosing strategies, we contrast maximum tolerated dose (MTD) strategies as compared with low dose, high density metronomic strategies (LDM) for tumors with different growth rates, concluding that LDM strategies can outperform MTD strategies, especially for fast growing tumors that thrive on long periods of unhindered growth without chemotherapy. The model shows good agreement with tumor growth data for unperturbed tumor growth and regression under drug therapy. ❧ The final section of this thesis is focused on the pre-existing resistant clones present in a tumor at the start of treatment, which remains a major problem in cancer therapeutics today. Tumor relapse and the development of chemotherapeutic resistance is now thought largely to be a consequence of the mechanism of “competitive release” of pre-existing resistant cells in the tumor which are selected for growth after chemotherapeutic agents attack the sub-population of chemo-sensitive cells which had previously dominated the collection of competing sub-clones in the tumor. We explain the important parameters (cost of resistance and initial fraction of resistance) in anticipating the evolutionary adaptations of the tumor in order to design therapies that exploit or mitigate the harmful effects of potential future adaptions. We introduce an “adaptive therapy” concept of utilizing information in the dynamical phase plane to target the growth of the resistant population indirectly, by introducing periods of drug holidays.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance
PDF
Stochastic multidrug adaptive chemotherapy control of competitive release in tumors
PDF
Ancestral inference and cancer stem cell dynamics in colorectal tumors
PDF
Reinforcement learning based design of chemotherapy schedules for avoiding chemo-resistance
PDF
A stochastic Markov chain model to describe cancer metastasis
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
PDF
Molecular signature of aggressive disease and clonal diversity revealed by single-cell copy number analysis of prostate cancer cells across multiple disease states
PDF
Tumor-immune agent-based modeling: drawing insights from learned spaces
PDF
Engineering immunotoxin and viral vectors for cancer therapy
PDF
Stochastic and multiscale models for urban and natural ecology
PDF
Integrative mathematical and computational approaches to investigate metastatic cancers
PDF
Heterogeneous graphs versus multimodal content: modeling, mining, and analysis of social network data
PDF
Insights from a computational approach to natural killer cell activation
PDF
A parallel computation framework for EONS synaptic modeling platform for parameter optimization and drug discovery
PDF
Efficient connectivity assessment of heterogeneous porous media using graph theory
PDF
Understanding anti-angiogenic signaling and treatment for cancer through mechanistic modeling
PDF
Combination therapy for solid tumor
PDF
Probabilistic medical image imputation via Wasserstein and conditional deep adversarial learning
PDF
Feature and model based biomedical system characterization of cancer
PDF
Investigating the complexity of the tumor microenvironment's role in drug response
Asset Metadata
Creator
West, Jeffrey Benjamin
(author)
Core Title
Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
09/05/2017
Defense Date
07/26/2017
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
chemotherapy modeling,mathematical oncology,OAI-PMH Harvest,resistance,tumor evolution,tumor progression
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Newton, Paul K. (
committee chair
), Nakano, Aiichiro (
committee member
), Pahlevan, Niema (
committee member
)
Creator Email
runningwest@gmail.com,westjb@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c40-426878
Unique identifier
UC11265787
Identifier
etd-WestJeffre-5707.pdf (filename),usctheses-c40-426878 (legacy record id)
Legacy Identifier
etd-WestJeffre-5707.pdf
Dmrecord
426878
Document Type
Dissertation
Rights
West, Jeffrey Benjamin
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
chemotherapy modeling
mathematical oncology
tumor evolution
tumor progression