# Optimal treatment and stochastic modeling of heterogeneous tumors

- Hamidreza Badri
^{1}and - Kevin Leder
^{1}Email author

**Received: **21 March 2016

**Accepted: **7 July 2016

**Published: **23 August 2016

## Abstract

In this work we review past articles that have mathematically studied cancer heterogeneity and the impact of this heterogeneity on the structure of optimal therapy. We look at past works on modeling how heterogeneous tumors respond to radiotherapy, and take a particularly close look at how the optimal radiotherapy schedule is modified by the presence of heterogeneity. In addition, we review past works on the study of optimal chemotherapy when dealing with heterogeneous tumors.

**Reviewers**: This article was reviewed by Thomas McDonald, David Axelrod, and Leonid Hanin.

### Keywords

Tumor heterogeneity Radiotherapy Stochastic modeling## Background

In recent years there have been many exciting studies that have observed the high levels of diversity present within tumors, (e.g., [42, 74, 82]). In addition to this genomic diversity it is possible for intra-tumor diversity to show up through cell cycle asynchrony or variability in microenviroment. This intra-tumor diversity has the potential to alter the evolutionary trajectory of the tumor cell population under therapy. An important question this raises is how we design optimal treatment strategies when dealing with heterogeneous populations. For example, if we have multiple therapies available there might be tumor subpopulations that respond better to certain therapies. The question then becomes how we optimally administer the various therapies. Addressing this question requires the use of mathematical models to understand the heterogeneity present, and in addition the development of optimization techniques to treat heterogeneous tumors.

In this work, we review past literature that has studied the question of optimal treatment of heterogeneous tumors, as well as stochastic modeling of heterogeneous tumors. The primary focus of the review is on the structure of optimal radiotherapy fractionation schedules when incorporating intra-tumor heterogeneity. A reason for focusing on the radiotherapy setting is that, in simple models of radiotherapy there are well established results for the structure of optimal radiotherapy schedules; see e.g. Badri et al. [5]. Therefore, it is possible to investigate the changes in the optimal schedule as a result of incorporating tumor heterogeneity.

We also review past literature on stochastic modeling and stochastic optimization for the treatment of heterogeneous tumors with chemotherapy or targeted therapy. In this section we look at works that considered stochastic models of heterogeneity in response to therapy, looking at works on both stochastic optimization and stochastic analysis.

## Modeling and optimization in radiotherapy

In this section we will review previous works that studied mathematical modeling and optimization of radiotherapy for heterogeneous tumor cell populations using Linear-Quadratic (LQ) model and its various extensions based on timing effects, cell cycle, hypoxia and cancer stem cell.

### Background on the linear-quadratic model

The LQ equation is widely used to describe the effects of ionizing radiation on normal and neoplastic tissue (For a review see [73]). The basic model states that the fraction of cells that survives a radiation dose of *d* Gy is given by exp ( − *αd* − *βd*
^{2}) where the radiosensitivity parameters, *α* and *β*, account for non-repairable lesions to DNA and the lethal mis-repair events occurring in the repair process of DNA double strand breaks (DSB), respectively [49, 71]. The initial model has been extended to include the four ‘Rs’ of radiobiology, repopulation of the tumor cells during the treatment period by surviving tumor cells, reoxygenation of hypoxic cells, repair of radiation-induced damage between fractions and redistribution of cells in the cell cycle [96]. These four phenomena are often extended by a fifth ‘R’, which is intrinsic radiosensitivity, defined as the considerable variability between different cell types [85]. These are important determinants of local tumor control after fractionated irradiation, and significantly change the optimal fractionation schemes. In this section, we review several studies that model tumor heterogeneity in radiation fractionation problem and discuss how the conventional optimal fractionation protocols change when considering intra-tumor heterogeneity.

Despite a significant history of predicting doses response curves by the LQ model [13], there is a significant amount of debate as to whether the LQ is appropriate for measuring high dose per fraction effects in stereotactic high-dose radiotherapy (e.g., see [53, 81]). The application of the LQ model is thought to underestimate tumor control at high doses (larger than 10 Gy). Several models have been proposed for improving the prediction of high dose survival curves, e.g. see the models developed by Hanin [51, 52] and Hanin and Zaider [53] or the review by Brown et al. [15] which discusses the validity of LQ model to high dose irradiation of tumors in detail. Since the LQ model is the most widely used model for quantitative predictions of dose/fractionation dependencies in radiotherapy and most models for heterogeneous tumors have been developed based on the same principal structure of the LQ model, we will mainly focus on the LQ model and its extensions in this study.

There are two widely used approaches for delivering radiotherapy: fractionated and continuous radiation. Assuming sufficiently large inter-fraction time, in fractionated radiation, the damage induced in a cell by an acute dose of radiation either causes cell death or complete repair of the cell before the next exposure. Therefore this model leads to memoryless kinetic that can be captured using Markov processes. However this is not the case for continuous irradiation where a longer biological memory of the irradiated cells is stored. See the work of Hanin et al. [57] and experimental studies cited therein for a more detail discussion on how the processes of damage repair/misrepair, cell proliferation and cycling can be modeled by a non-Markovian model. The remainder of this work will largely focus on models of fractionated radiotherapy.

*N*treatment fractions in which radiation dose

*d*

_{ i }is administered in fraction

*i*(

*i*= 1,..,

*N*) is given by

where [*α*/*β*] is a tissue-specific radio-sensitivity parameter. The normal tissues toxicity constraints in radiotherapy fractionation problem are mathematically modeled by insisting that BED levels for various OAR stay within prescribed levels [5] or keeping the total number of functional proliferating normal cells more than the required threshold [54]. These constraints can be satisfied by keeping the total dose, fractional dose or dose rate in continuous irradiation within some acceptable levels.

Two possible solutions to the fractionation problem are hyper-fractionated and hypo-fractionated schedules. In hyper-fractionated schedules small fraction sizes are delivered over a long period of time whereas in hypo-fractionated schedules, large fraction sizes are administrated during a short period of radiation delivery. If we maximize tumor control probability (TCP) at the conclusion of treatment, it has been observed that whether hyper or hypo-fractionation is optimal depends on the radio-sensitivity parameters of the normal and cancerous tissue [5, 70, 92]. More specifically if tumor *α*/*β* ratio is smaller than effective *α*/*β* ratio for all normal tissues (defined as (*α*
_{
i
}/*β*
_{
i
})/*γ*
_{
i
}, where *α*
_{
i
}/*β*
_{
i
} and *γ*
_{
i
}, denote the radio-sensitivity parameter and sparing factor, respectively, in *i*th OAR), then a single-dosage solution (hypo-fractionated schedule) is optimal, whereas a multiple-dosage solution with equal doses (uniform schedule) is optimal otherwise (hyper-fractionated schedule) [5, 6].

These results are based on the assumption that irradiated cell survival curves are explained by the LQ model, therefore TCP is invariant under rearrangement of fractional doses. However considering more complicated models [55] or different objectives such as minimizing metastatic risk [7, 8] instead of maximizing TCP may result in the optimality of non-standard fractional schedules. These schedules are formed using the front loading principle: administering the maximum possible dose as soon as possible. Moreover, as a result of these emerging alternative models and objectives, other factors such as the time point at which the performance criteria is evaluated may play an important role in the structure of optimal schedules, e.g. see Zaider and Hanin [102] and Badri et al. [7].

### Intra-tumor heterogeneity

The uncertainties in radiotherapy treatment can be categorized into two groups: inter-patient variability and intra-tumor heterogeneity. Inter-patient variability stems from heterogeneity in patient-specific variables such as the sensitivity of their normal tissues and tumor to radiation (*α*/*β* ratio), the growth rate of their tumor or the healing kinetics of normal tissues. Several studies addressed these uncertainties using different techniques. Badri et al. [6] proposed a stochastic optimization formulation to incorporate inter-patient variability in tumor and normal tissue radiosensitivity parameters (*α* and *β*) and sparing factor of the OAR into the scheduling optimization problem. Hanin and Zaider [54] developed a mechanistic approach that models post-irradiation normal tissue toxicity when considering inter-patient variation of kinetic parameters. On the other hand, to improve the efficacy of radiation therapy, it is necessary to study the role of intra-tumor heterogeneity, since it significantly changes the tumor response curves [34, 58]. The range of cell sensitivity comes from inherent genetic and epigenetic differences among the tumor cells and from temporal variations arising from the asynchronous cell cycle phases and variable micro environmental conditions during therapy. The focus of the present work is to review studies that model intra-tumor heterogeneity and present where possible novel optimization problems that arise from these models.

In [56], Hanin et al. studied the role of radiosensitivity variation amongst cancer cells on optimal radiotherapy fractionation schemes. They used a new criterion developed by Rachev and Yakovlev that considers the difference between weighted survival probabilities for normal and neoplastic cells, where tumor cell radiosensitivity is considered as a random variable with a known distribution function [76]. For several special cases, the exact solution of optimal fractionation is obtained and an iterative approximation methodology is designed when it is not possible to compute the exact optimal fractionation schedule.

### Timing and 4R’s effects on tumor heterogeneity

The fraction of surviving cells after a dose of radiation not only depends on dose and tumor radio-sensitivity parameters, but also it typically depends on the time-course of dose delivery [88]. Timing affects cell killing due to several reasons such as DNA repair and misrepair, tumor repopulation, redistribution and reoxygenation [48, 49]. The basic LQ model typically assumes that tumor radio-sensitivity parameters (*α* and *β*) and repopulation are constant over the time course of radiotherapy. This implies the failure of the simple version of the LQ equation exp ( − *αd* − *βd*
^{2}) to capture the dynamics of reoxygenation and repopulation throughout the course of treatment. Mathematical models for ionizing radiation therapy, applied to multicellular populations whose cells have time-dependent radio sensitivity have been studied widely [17, 60]. However in some cases such as heterogeneity associated with cell sensitivity and proliferation rate when fractionated irradiation with sufficiently many fractions or protracted continuous radiation is implemented, it is possible to only consider the homogeneous subpopulations of the most resistant and/or fastest growing cells. This is due to the fact that usually slowly growing tumor cells and sensitive subpopulations die out after commencement of therapy, and therefore it is sufficient to design the therapy to target the fast growing tumor cells and resistant population. As an example see the mathematical model developed by [54] to model the number of proliferating as well as non-proliferating normal cells as a function of time post treatment when incorporating the selection of the fastest growing subpopulation to capture the tissue damage at the conclusion of therapy and of the subsequent healing kinetics.

*α*at time

*t*, i.e.

*n*(

*α*,

*t*), we can write the equation explaining the fluctuating diversity of a population with fixed size using a Kolmogorov forward equation as (see [60] for more details)

*u*denotes the average number of DSB per cell, \( \frac{1}{2}\kappa\ {u}^2 \) shows the average rate at which binary misreapirs removes DSB by lethal rearrangements,

*k*displays the rate at which cells change their radiation sensitivity, and

*α*

_{0}and

*σ*

^{2}represent the mean and variance of random variable

*α*, respectively. Note that in the case of a homogeneous tumor,

*σ*= 0, Eq. (1) becomes the deterministic model developed by Sachs et al. in [80] which adds the enzymatic modification of the immediate damage through a Markov process to the basic LQ model. Considering tumor population in the long term, it was shown that the solution to the Eq. (1) gives the surprising simple result of

*N*(

*t*) shows the total population at time

*t*,

*D*is total radiation dose delivered for period (0,T), and

*G*is the Lea-Catcheside function [60]. Equation (2) can be considered as the elementary LQ model with

*α*being replaced by its average

*α*

_{0}, and

*β*being replaced by its modified value. Results of their analysis support the hypothesis that the therapeutic paradigm of low dose rate or fractionated radiation can help conquer radioresistance in hypoxic tumors [91, 97]. This is due to the fact that a large fractionation interval (parameter

*T*in (2)) allows the tumor population to complete the reoxygenation process and thereby the tumor population radio-resistance due to oxygenation status will be minimized. This phenomenon is supported by a smaller coefficient for

*D*

^{2}in Eq. (2). One year later, Brenner et al. developed a parsimonious model to include the resensitization effect into the LQ model. In the extended model, designated LQR, survival is written as a function of dose

*d*as

*α*[14]. The cell survival values based on Brenner et al. model (Eq. (3)) are plotted in Fig. 2 for values of

*σ*

^{2}= 0, 0.01 and 0.09 for cell population without, low and high diversity, respectively. By comparison of cellular diversity effect for tumors with different values of \( \frac{\alpha }{\beta }, \) we observe a more significant effect for tumors with large values of \( \frac{\alpha }{\beta } \), e.g. 10 for prostate cancer (Fig. 2b), compared to tumors with a small value of \( \frac{\alpha }{\beta } \), e.g. 3 for head and neck cancer (Fig. 2a).

Optimization of radiotherapy treatment within the Hlatky model which includes time dependence of sublethal damage repair has been studied by Yang and Xing in [99]. It has been observed that incorporating these effects into the LQ model may give rise to optimal non-uniform fractionation schedules where fractional doses at the beginning and end of each irradiated week become significantly greater than others. Furthermore it was observed that the hyper-fractionation schedule gives an insignificant advantage over hypo-fractionation or a standard regimen.

### Cell cycle

*G*

_{0}-phase of the cell cycle, quiescent cells, possess a lower level of radio-sensitivity than proliferating cells that are in the

*G*

_{1},

*S*,

*G*

_{2},

*M*-phases [23, 73]. Therefore for tumors with asynchronous cells, increasing radiation delivery time,

*T*, increases tumor radiosensitivity. This makes sense because at first the radiation kills the cells in more sensitive phases, and then radioresistant cells, e.g. those are in

*G*

_{0}-phase, have time to reach more sensitive phases. Also due to cell arrest in the most sensitive phases of cell cycle, protracted radiation promotes synchronization. Chen et al. studied the effect of cell cycle redistribution on the population resensitization when ignoring the quadratic misrepair of radiation damage,

*β*[17]. They used a Mc-Kendrick-von Foerster equation adjusted for the first track radiation cell kill to model the age dependent cell dynamics as

where *n*(*a*, *t*)*da* shows the density number of cells in the age range (*a*, *a* + *da*) at time *t* and *α*(*a*) shows the tumor radiosensitivity at age *a*. They observed that the tumor population resensitization effect occurs as the duration *T* of irradiation is increased from essentially zero times to short, and sufficiently small finite times. They concluded that population resensitization is proportional to *T*
^{2} and \( \exp\ \left(-\alpha (a)D\right)\ {\left(\frac{Dd\alpha }{da}\right)}^2 \) and the resensitization happens when *T* is small and the cell population is in a stable age-distribution phase before irradiation, which in this case happens regardless of how the radiation cell kill, function *α*(*a*), depends on age. Hahnfeldt and Hlatky generalized the model proposed by Chen et al. beyond constant-dose-rate irradiation and small *T* in more explicit terms [48]. They have used the same equation described in (4) and have shown mathematically that variation with the time of resensitisation due to redistribution is not monotonic but damped oscillatory. They found that spreading a dose of *d* Gy over a longer period of time in any way is more desirable and results in higher TCP than delivering an acute dose of equal magnitude. They proved that this result continues to apply regardless of age-dependent sensitivity and mitosis rate functions chosen.

*G*

_{1},

*S*,

*G*

_{2},

*M*) and a quiescent phase(

*G*

_{0}). If the clonogenic cells do not enter a

*G*

_{0}phase, which is modeled with considering the transition between both compartments during radiotherapy, then the model equally applies for a splitting into

*S*,

*G*

_{2},

*M*and

*G*

_{1}phases. The key assumption is that actively proliferating cancer cells are much more susceptible to radiation damage than quiescent cells. The basic model states that the expected number of cells in active,

*N*

_{ A }, and quiescent compartments,

*N*

_{ Q }, satisfy a system of differential equations as

where *μ* is the rate of active cell division, *ν* describes the transition from quiescent compartment into the cell cycle, *λ*(*t*) shows the death rate of different types of cells at time *t* and *h*(*t*) explains the radiation induced death rate in different compartments. Note that since active cells are more radiosensitive, we have*h*
_{
A
}(*t*) > *h*
_{
Q
}(*t*). The original model of Dawson and Hillen have been taken and extended to describe more complex systems or models with more compartments [25, 32, 59, 68]. Analysis of Dawson and Hillen active-quiescent radiation model and its comparison to LQ model confirms that a larger *α*/*β* ratio relates to a fast cell cycle and indicates the presence of a significant quiescent compartment, while a smaller ratio is associated with a slow cell cycle [23]. These comparisons were performed under the LQ model assumptions which allowed the authors to construct a relationship between proliferation and transition rates, *μ* and *ν* in Eq. (5), respectively, in their model with *α* and *β* parameters in the LQ model. Therefore we can conclude that for the tumor population with a substantial quiescent compartment, which indicates a large value for *α*/*β* ratio, hyper-fractionated schedules provide a better TCP than the hypo-fractionated schedules (see [70] or Badri et al. [5]). These types of analysis are indeed the future direction of the cell cycle modeling in TCP, i.e. the inclusion of cell cycle and diversity of the cellular radiosensitivity of a tumor in optimization of radiation dosing schedules.

### Hypoxia

*α*and

*β*, and the tumor net repopulation rate,

*γ*, depend upon the cell radial distance,

*r*, from the center of the tumor, and on the current tumor radius,

*R*[94]. Then if we assume all of these parameters take on a fixed well-oxygenated level at the tumor surface (i.e.

*α*=

*α*

_{0},

*β*=

*β*

_{0}and

*γ*=

*γ*

_{0}at

*r*=

*R*) and decrease linearly as

*r*decreases, we can compute the radio-sensitivity parameters and tumor growth rate as a function of

*r*∊ [0,

*R*] for

*R*<

*r*

_{0}as (see Fig. 3 and [94] for more details)

*r*∊ (

*R*−

*r*

_{0},

*R*] and

*R*>

*r*

_{0}as

As discussed by the authors, the linearity assumptions in Eqs. (6) and (7) may not be compatible with the physics of oxygen diffusion and were chosen for their parsimony and computational feasibility. The actual situation in vitro and in vivo is significantly more complex, e.g. the oxygen enhancement ratio depends on the fraction size [49], therefore a more complicated model is required to explain the tumor radiosensitivity as a function of radial location. Also Eqs. (6) and (7) are based on the assumptions that the net rate of spontaneous cell death decreases as the tumor shrinks, which is applicable for most types of tumors (e.g., well-differentiated squamous cell cancers) and is consistent with experimental results [87, 89].

*t*by

*R*

_{ t }and

*n*

_{ t }, respectively, and we assume the density of cells per unit volume in the spherical tumor to be

*θ*, then we have \( {n}_t=\frac{4}{3}\theta \pi\ {R}_t^3 \). Using LQ formulation adjusted for exponential tumor growth [49], the expected change in number of tumor cells after a dose of size

*d*is (see [94] for more details)

where *μ* is tumor repair rate. Substituting Eqs. (6) and (7) and using equation \( {\overset{.}{n}}_t=4\theta \pi\ {R}_t^2\ {\overset{.}{R}}_t \), we can write Eq. (8) in terms of *Ṙ*
_{
t
} and *R*
_{
t
} and forms the basis of the optimal control problem. Wein et al. proposes a dynamic programming approach to numerically solve this problem. The resulting optimal protocols suggest a non-standard time varying schedules with irregular time intervals between fractions, administering larger fractions before longer breaks, such as afternoon sessions or Fridays, and shorter fractions before shorter breaks, such as morning sessions [94]. Wein et al. proposed two main reasons for this phenomenon. First the large fractions make up for tumor repopulation during overnight or weekend breaks. Second the tumor size is smaller at the end of the week, i.e. Fridays, and smaller tumors are more sensitive to radiation. They also observed that as the tumor shrinks during therapy, it is optimal to increase the doses on Friday afternoons. Based on their model, as the tumor shrinks, *α*(*R*)/*β*(*R*) becomes smaller which leads to the optimality of hypo-fractionated schedules.

### Cancer stem cell

The existence of cellular heterogeneity in solid tumors may originate from a number of sources, including hypoxia, cell cycle asynchrony, infiltration of normal cells, vascular structures and stroma into the tumor and the hierarchical structure of the cell populations from which cancers arise. The cancer stem cell (CSC) model of tumorigenesis has received significant attention in recent years. CSC refers to a subset of tumor cells that has the ability to self-renew and generate differentiated progeny which make up the bulk of a tumor [77]. Existence of CSCs has been identified in different cancers such as acute myeloid leukemia [26] breast cancer [1] and brain tumors [84]. The definition of CSC implies that an anticancer therapy can control a tumor, i.e. permanent local tumor control, only if all CSCs are eradicated. Therefore it is possible that removal of CSCs is the crucial determinant in curing cancer and eradicating tumor cells [10].

*α*and

*β*). The radiation response model is constructed as

*S*(

*d*) represents the fraction of surviving cells after delivering an acute dose of radiation,

*F*represents the fraction of CSCs out of all cells, and (

*α*

_{ s },

*β*

_{ s }) and (

*α*

_{ d },

*β*

_{ d }) show the radiosensitivity parameters in CSC and DCC, respectively. The interplay between CSCs and DCCs can be modeled by using the ODE introduced in Hillen et al. [59] as (10)

where *N*
_{
s
}(*t*) and *N*
_{
d
}(*t*) are the volume fractions of CSCs and DCCs, respectively. The function *N*(*t*) is the total volume of tumor normalized between 0 and 1 which is equal to *N*
_{
s
}(*t*) + *N*
_{
d
}(*t*), *p* is the probability of symmetric CSC division, and *μ*
_{
s
}, *μ*
_{
d
}and *a*
_{
v
} define the CSC growth, DCC growth and DCC apoptosis rate, respectively. *k*(*N*(*t*)) is a constraint defined as max {1 − *N*(*t*)^{
σ
}, 0} for a *σ* ≥ 1 and keeps the total volume fraction less than 1. In [4], Bachman and Hillen used the ODE Eqs. in (10) and showed that the differentiation therapy proposed by Youssefpour et al. [100], which is defined as the combination of radiotherapy and chemotherapy where the chemotherapeutic agent pushes CSCs into the differentiation stage, can have large beneficial effects in head and neck cancer, brain cancers and breast cancer for the patient increasing treatment success and reducing side effects.

*ρ*∊ (0, 1], i.e.

*α*

_{ s }=

*ρα*

_{ d }and

*β*

_{ s }=

*ρβ*

_{ d }. This simplifying assumption enabled the authors to characterize the sensitivity of CSCs to radiation by a single parameter,

*ρ*. The model is described in Fig. 4. The model stipulates that

*t*hours after the previous dose of radiation, the fraction of DCCs capable of reversion to CSCs is given by \( \gamma (t)={\gamma}_0{e}^{-{\left(t-{a}_0\right)}^2/{a}_1^2} \) (note that

*γ*(

*t*) =

*γ*

_{0}for the first dose of radiation), for some constants

*γ*

_{0},

*a*

_{0}and

*a*

_{1}and the fraction of surviving cells can be computed based on the LQ model. They predicted several optimal radiation strategies that substantially enhanced survival in experimental studies using a mouse model of glioblastoma. The resulting optimized schedules recommend a non-uniform schedule delivering larger fractions at the beginning and toward the end of the therapy. In a follow up work, Badri et al. used the Leder model to consider fractionated schedules that have optimal survival while, maintaining acceptable levels of toxicity in early- and late-responding tissues [5]. They derived the closed form solution to the problem and proved that the optimization problem can be split into two separate optimization tasks that can be tackled independently. The first model involves optimization of dose per fraction and the optimal total dose, and the second model optimizes inter-fraction intervals between radiation doses. It was observed that normal tissues sparing factors and radiosensitivities, and the magnitude of the

*α*/

*β*ratio for tumor are determinant factors defining the optimal radiation scheme, i.e. for low (high) values of tumor

*α*/

*β*ratio, the hypo-fractionated (hyper-fractionated) schedule is optimal. For the time-dependent model, the optimal inter-fraction intervals only depend on the time dynamics of the dedifferentiation process and treatment duration. In particular it was observed that optimal inter-fraction intervals are equal to the dose spacing that leads to the maximal amount of cell reversion to the stem-like state, i.e.

*a*

_{0}.

Several stochastic and cellular automata models have been used in more complicated simulation based studies of complete tumor cell kinetics during radiation therapy. Gao et al. used an integrated experimental and cellular Potts model to simulate glioblastoma population growth and response to irradiation [41]. They found that in order to maintain the tumor population following radiotherapy, surviving glioma CSCs in vitro increase their rate of self-renewal, i.e. the fraction of CSCs in the populations is increased after radiation. By comparing acute and fractionated irradiation response, the authors observed that the relative increase in fraction of CSCs in tumor population after fractionated treatment cannot be explained merely by radioresistance of CSCs. This simulation based model suggests that repeated exposure to radiation might increase the symmetric division rate of CSCs, which eventually may lead to accelerated repopulation of CSCs. A series of in vivo 4D simulation models for GBM explore the tumor growth dynamics and response to radiation, considering vasculature, oxygen supply and radiosensitivity [2, 27, 28]. These works clustered cells into dynamic classes based on the mean cell cycle phase durations over the various cell cycle phases and used a linear quadratic model to describe the number of cells killed. They associated p53 mutations with increased radioresistance and inefficient clinical outcome for patients with GBM, as suggested by Haas-Kogan et al. [46]. Evaluating the response to treatment for different fractionation regimens revealed that hyper-fractionated schedules may lead to an improvement in local tumor control compared to standard schedules.

## Chemotherapy

In this section we will review previous works that studied stochastic modeling and optimization of chemotherapy for heterogeneous tumor cell populations.

### Optimization models

*x*is the cancer cell population size and

*u*is the drug concentration level. Also

*f*and

*g*are arbitrary functions that represent density dependence and drug induced cell kill respectively. Then a cost function is specified, e.g.

and the goal is then to use optimal control methodology to numerically identify the optimal drug concentration profile *u*. There is a wide range of works on models of this kind and we refer the reader to the reviews Shi et al. [83], Swan [86] and Kimmel and Swiernak [64] for further examples.

Given the large amount of literature on this topic, we focus in the remainder of this section on works related to optimization of stochastic models of the treatment process for heterogeneous tumors with resistant subtypes.

### Optimization of stochastic models

*m*possible cell types, and all cells of a given type behave in a statistically identical fashion independently of all other cells present. In particular a cell of type-

*i*well wait an exponentially distributed amount of time with mean 1/

*a*

_{ i }before a birth/death/mutation event. During this event the type-

*i*cell produces offspring of type (

*j*

_{1}, …,

*j*

_{ m }) with probability

*p*

^{(i)}(

*j*

_{1}, …,

*j*

_{ m }), where

*j*

_{1}+ … +

*j*

_{ m }∊ {0, 2} (see Fig. 5). The multi-type branching process is specified by the vector \( \overrightarrow{a} = \left[{a}_1\cdots\ {a}_m\right] \) and the vector valued mapping

The long term behavior of a multi-type branching process is easily deduced from this information. In particular, one forms a mean matrix *M* = (*m*
_{
ij
}), where *m*
_{
ij
} is the expected number of type-*j* offspring a type-*i* offspring will produce. If the maximum eigenvalue of the matrix *M* is less than or equal to one then the branching process is guaranteed to go extinct, while if it is greater than 1 then the branching process can either go extinct or its size diverge to positive infinity. Therefore understanding the long term behavior of the branching process is straightforward. When studying the problem of drug resistance in cancer one is often interested in the behavior of the process over a long (but finite) time interval, and therefore it is not sufficient to simply look at the maximal eigenvalue of *M*. For an example of other techniques that can be used see e.g. Durrett and Moseley [30], Iwasa et al. [61], Haeno et al. [47], or Durrett et al. [31].

When modeling drug resistance in chemotherapy, a standard approach would be to assume that initially most cells are type-1, which is assumed to be sensitive to some first line therapy. Thus during treatment with this first line therapy the type-1 cells will begin to decrease; however, these cells may mutate to a different type of cell that can grow under the first line therapy. This type of cell may decline under a second line therapy; however it may now mutate to a type of cell resistant to both types of therapy. In this model then the question becomes, how does one administer the various therapies so that the risk of total treatment failure (no more viable drugs) is minimized.

Seminal work was done in this field by Coldman and Goldie in several papers, e.g., Goldie and Coldman [43], Goldie et al. [44] and Coldman and Goldie [19]. We will focus on Coldman and Goldie [19], since it generalizes the previous works. It is assumed that there are *n* treatments available *T*
_{1}, …, *T*
_{
n
}, and 2^{
n
} different cell types present, each type specified by which subset of therapies the constituent cells are resistant to. Specifically, \( {R}_{i_1,\dots,\ {i}_m}(t) \) is the number of cells at time *t* that are resistant to the therapies \( {T}_{i_1}\dots,\ {T}_{i_m} \) and sensitive to all other therapies. The cell type *R*
_{0} is sensitive to all therapies. In the absence of therapy it is assumed that all cells behave according to a pure birth process with birth rate *λ* per cell. During cell division events, mutations may occur and cells can acquire resistance to new types of drugs. Chemotherapy is modeled as an instantaneous probabilistic reduction in population of all sensitive cells according to a log cell kill rule. The authors then derive formulas for the probability of evolution of cells resistant to therapies within a finite time horizon. Coldman and Goldie, consider the case of two therapies and three distinct resistant cells in depth. In particular, let *P*
_{12}(*t*) be the probability that no cells with resistance to both therapies evolve by time*t*. Under symmetry assumptions on the efficacy of the two therapies and the behavior of the two singly resistant mutants, Coldman and Goldie [19] establishes that alternating therapies maximizes*P*
_{12}(*t*). Day computationally investigated relaxation of the symmetry assumptions and found that some non-alternating schedules could outperform the alternating schedule in that scenario [24]. In particular Day proposed a ‘worst drug first’ rule, this rule was investigated in further depth by Katouli and Komarova who considered a wide range of possible cyclic therapies [62]. In later works Murray and Coldman [72] and Coldman and Murray [20] extended the original model of Coldman and Goldie [19] to allow for toxicity constraints on normal tissues, simultaneous administration of multiple drugs, and included the possibility of inter-patient heterogeneity. In [18], Chen et al. further investigated the effects of asymmetry in the efficacy of the two possible therapies, and derived general conditions for the identification of optimal drug administration sequences. One potential shortfall of the Goldie and Coldman model is that the tumor cell populations grow exponentially ignoring possible effects of resource depletion. Chapter 9 of the monograph Martin and Teo [69] develops a deterministic model that allows for logistic and Gompertz growth in the tumor population. In this model they have four types of cells and two therapies, the authors develop an algorithm that searches for the schedule of therapies that leads to the maximal time until treatment failure. Note that this algorithm only identifies local optima though.

Despite the large amount of work done in this field there are still significant challenges remaining. In particular, previous works have looked at optimal schedules with only a small number of resistant types and potential therapies. Going forward, an important extension will be to develop methodologies that allow for the optimization of administration schedules for larger number of therapies and a larger number of resistant types. Another possible extension is to study minimization of probability of resistance in more complex stochastic models. Until now the stochastic models have all had essentially exponential growth properties, which are known to be inconsistent with tumor growth curves. An exciting challenge for the future is to minimize resistance probabilities in stochastic models that include density dependence.

### Stochastic analysis

*t*by

*Z*

_{0}(

*t*) and the drug resistant cell population by

*Z*

_{1}(

*t*), with

*Z*

_{0}(0) =

*n*and

*Z*

_{1}(0) = 0. The sensitive cell population is modeled as a subcritical binary branching process, that produces resistant cells at rate

*μ*and each resistant cell initiates a super-critical branching process with random net growth rate. In these works the properties of the cancer cell population is investigated at the ‘crossover-time’:

In particular, Foo and Leder [35] study the relationship between *ξ* and the extinction time of the sensitive cell process. While Foo et al. [38] study the diversity properties of the resistant cell population at the time *ξ*. There are several standard metrics for diversity of a population, e.g. the number of distinct species present, the Simpson’s Index (probability two randomly chosen cells are genomically identical), and Shannon’s Index (related to Shannon’s Entropy, see e.g. [22]). In Foo et al. [38] they consider all three of these diversity measures. Lastly the work of Foo et al. [39] establishes a central limit theorem for *ξ* in the limit as the initial population *Z*
_{0}(0) goes to infinity, and identifies the effect of the random fitness distribution on the large *n* behavior of the crossover time *ξ*.

There are lots of open problems remaining in the topic of stochastic models of cancer cells undergoing therapy. An interesting extension would be to investigate the treatment process when spatially explicit models (such as [95]) are used.

## Discussion

Viewing tumors as an evolving population of cells has proven to be a useful tool in the study of cancer. Anti-cancer therapy clearly has the potential to impact the evolutionary trajectory of the tumor cell population. The behavior of this evolution is extremely interesting in the context of diverse tumor cell populations. For example, one might expect that therapy will select for cells with therapy resistance, thus leaving a more difficult to treat tumor. In order to achieve the best possible therapeutic results it is thus seems necessary to create treatment strategies that take into account the diversity present within a tumor and the evolutionary changes the tumor might undergo during therapy.

There has clearly been a significant amount of work done in the field of cancer therapy optimization. However, there are still lots of exciting problems remaining to be investigated. For example, there are few theoretical results about the structure of optimal radiotherapy schedules when studying heterogeneous populations. In the chemotherapy setting there are no suitable optimization methods for dealing with large amounts of heterogeneity present, i.e., large numbers of distinct cell types. There are several interesting open problems in the stochastic modeling and optimization framework. In particular, more work needs to be done in this area that incorporates cellular competition.

Perhaps the biggest challenge in the field of designing optimal cancer therapies, is bringing these optimized therapeutic schedules into the clinic. While there have been successes in the laboratory setting, e.g., Leder et al. [67], Gao et al. [41], successes in a clinical setting are quite rare.

## Abbreviations

CSC, cancer stem cell; DCC, differentiated cancer cells; DLQ, dual-compartment linear-quadratic; DSB, double strand breaks; LQ, linear-quadratic; OAR, organs-at-risk; TCP, tumor control probability

## Reviewers’ comments

### Reviewer’s report 1 Thomas McDonald, Biostatistics and Computational Biology, Dana-Farber Cancer Institute

Reviewer comments:

**Summary:**

The review provides a general overview of modeling therapy of tumors. It separates into radiotherapy and chemotherapy discussing historical and more recent models of each and the impact of heterogeneity that affect tumor response. The authors do a good job discussing radiotherapy beginning with the Linear-Quadratic model before moving into the various extensions too account for the four ‘Rs’. The section on Chemotherapy covers a wide range of work from the Coldman and Goldie models up to modern methods used and include a discussion of the next necessary steps and issues to tackle. Ultimately, this review provides a useful recap of the work done in mathematical modeling of radiotherapy and chemotherapy.

**Reviewer recommendations to authors:**

Major: The main suggestion is to include a few more pictures of some of the processes mentioned. The radiotherapy models could be illustrated with curves and example plots of tumor response curves showing the impact of heterogeneity as modeled in some of the articles cited.

The second part on chemotherapy seems lacking in the detail that the radiotherapy section got, and it may deserve a little more time or mathematical description of the work. The focus of the work is clearly radiotherapy, but explaining some of the chemotherapy models in a little more depth or describing a quick background of branching processes and their use would make the review more complete. A more careful proofreading is necessary. There are minor grammatical errors scattered throughout. An incomplete list is given below.

The first section on radiotherapy may be separated into subsections since the authors jumped between models abruptly.

Author’s response: *Thank you for your careful reading of the manuscript and helpful suggestions, we have addressed these comments*.

### Reviewer’s report 2 David Axelrod, Rutgers University

Reviewer comments:

**Summary:** Recommendation status: Endorse publication as a Review. Reviewers report: Summary of some mathematical modeling to optimize radiotherapy and chemotherapy, with brief mention of open problems, but little indication of whether or not the modeling has had a clinical impact, and if not why not. Not comprehensive or original, although useful as an entrance to the literature.

Author’s response: *Thank you for your careful reading of the manuscript and helpful suggestions, we have addressed these comments.*

### Reviewer’s report 3 (Author’s response included in italics) Leonid Hanin, Idaho State University

Reviewer comments:

**Summary:** The authors attempted to review a huge research field (mathematical models of radiation therapy/chemotherapy and stochastic models of tumor cell populations) through a prism of optimal cancer treatment schedules and intra-tumor heterogeneity. From among many hundreds of relevant articles and dozens of books published on these subjects the authors selected a relatively small fraction for their discussion. The review is somewhat sketchy and oftentimes superficial. In my opinion, it does not delve deep enough into biological, clinical and mathematical issues. The overall picture of the field does not come through clear enough as well. I believe the review is missing some general guiding ideas that would make the discussion of the subject coherent and captivating from methodological and historical standpoints.

**Reviewer recommendations to authors:** The article provides a brief overview of the following areas of biomathematical research: cancer radiotherapy, chemotherapy and stochastic modeling of cancer cell populations. The umbrella topic that gives a common thread to the reviewed modeling approaches and results is the effects of heterogeneity of cancer cell populations on optimal radiotherapy or chemotherapy schedules. I believe the review is too sketchy, incomplete and lacking technical details to do justice to these extensive and important areas of research. Specifically, the following important questions have not been addressed in sufficient depth and detail:

(1) What are the biological and mathematical assumptions underlying the quoted studies?

(2) What are the sources of heterogeneity (spatial, oxygenation level, radiosensitivity, cell cycle phases, variation in kinetic parameters, inter-patient variation, etc.) accounted for and disregarded in any particular study? Without this, results of the cited research can be neither fully appreciated nor compared.

(3) What are the types of radiation involved (fractionated, continuous with constant dose rate, brachytherapy, etc.)?

(4) Are the results and conclusions theoretical or numerical and, in the latter case, how were model parameters selected?

(5) What is the basis for various equations discussed in the text?

Author’s response: *Thank you so much for your comments, we have included more discussion on the underlying assumptions on the models, equations, various types of radiation, and classification of the heterogeneity sources in radiotherapy*

General Comments

1. Due to selection effects some of the heterogeneity issues seem irrelevant in the case of fractionated irradiation with sufficiently many fractions, protracted continuous radiation and chemotherapy. For example, this is the case for heterogeneity with respect to cell sensitivity and proliferation rate, for sensitive and slowly growing tumor subpopulation will disappear soon after the start of treatment, so it seems feasible to deal from the outset with the homogeneous subpopulations of the most resistant and/or fastest growing cells. See e.g. [1] in the list of references below for the discussion of selection of the fastest growing subpopulation. This fairly obvious but important consideration provides the missing evolutionary context to the discussion of heterogeneity.

Author’s response: *Thank you so much for you valuable comment. We have added a paragraph in the paper to explain this matter.*

2. The article disregarded a profound difference between fractionated and continuous radiation. While the former leads to memoryless kinetic models that can be described using Markov processes, the latter brings about long biological memory (due to the arrest of irradiated cells in the most radiosensitive phases of the cell cycle and non-markovian kinetics of radiation damage accumulation and repair/misrepair), see e.g. [2] and experimental studies quoted therein.

Author’s response: *Thank you so much for pointing out this shortcoming. We have added a paragraph to address this issue.*

3. As it was briefly mentioned in the article, clinically relevant approaches to radiotherapy optimization must involve constraints accounting for damage to normal tissue. However, no details were provided and no results reviewed. Modeling normal tissue complication probability (NTCP) leads to many mathematical and biomedical challenges including heterogeneity issues [1].

Author’s response: *We appreciate reviewer comment. The focus of the present work is to review the studies that properly model the intra-tumor heterogeneity. However to add some discussion about this important topic, we have added a paragraph that explains the difference between inter-patient and intra-tumor heterogeneity. We also cited reference [1] to provide some additional sources for covering this important topic briefly.*

4. The article deals with the linear-quadratic (LQ) model of irradiated cell survival and its extensions. This model is based on a fairly sophisticated mechanistic description of the kinetics of sublesion generation, repair/misrepair and pairwise interaction that produces lethal lesions (typically chromosomal aberrations). However, converting this formalism into cell survival probability is based on a highly unrealistic assumption that the distribution of the number of lesions is Poisson. Although this and other critical flaws of the LQ model have been uncovered about four decades ago, see [3] for a more recent discussion, LQ model is still in business. This is especially surprising given that alternatives have been proposed, see e.g. [3] where a parsimonious model based on rigorous microdosimetric analysis and overcoming many flaws of the LQ model was introduced. Another fundamental problem of the LQ model is that it is inapplicable to large doses (>10 Gy) [4]. For example, it was shown in [3] that for such doses LQ model underestimates cell survival (compared to the more realistic model developed in [3]) by several orders of magnitude! Insisting on the LQ model confines researchers to a mathematical abstraction that in many cases has little to do with clinical reality.

Author’s response: *We appreciate reviewer comments on this shortcoming. We have added a paragraph in the paper to explain this shortcoming and present our reasoning for using LQ model.*

5. Discussion of optimal radiation schedules is overly ad hoc. Addressing the question as to whether some general principles are true in a given biological/modeling setting would bring much needed structure and clarity. For example, is TCP invariant under rearrangement of fractional doses? Does it satisfy the front loading principle (i.e. “hit the tumor as hard and as early as possible”) true? Is the uniform radiation schedule optimal? For a discussion of these and other general principles, see [5-7].

Author’s response: *Thank you so much, we have added a few sentences to explain this topic briefly.*

6. The article says nothing about the time point at which TCP is evaluated. As discussed in [8], its selection is consequential.

Author’s response: *We have added a few sentences to explain its importance.*

7. Among the many Rs of radiation biology repopulation is probably the most important, yet its discussion in the article is scarce thus missing many aspects of the subject at hand. In the case of fractionated radiation, the TCP in the repopulation setting was computed in closed form in [9, 10] under arbitrary time-dependent birth and spontaneous death rates, arbitrary time post-treatment, arbitrary radiation schedule and arbitrary dose–response function. Moreover, computed in these two works was not only TCP but the entire distribution of the number of surviving clonogenic cancer cells.

Author’s response: *We have added these two papers to our review study.*

7. Discussion of optimization problems never mentioned constraints on the total dose, fractional doses and dose rates. What are they and where did they come from?

Author’s response: *Thank you. We have added a paragraph to explain this important matter in detail.*

Technical Comments

1. P. 4, lines 22-27. The basic LQ model contains also the time-dependent Lea-Catcheside dose-protraction factor that accounts for the temporal pattern of radiation dose delivery and depends on the rate of sublethal damage repair. Models accounting for repopulation and reoxygenation do not have to be of LQ type.

Author’s response: *Thank you. We agree that Lea-Catcheside model is time dependent, however we were referring to the basic equation of LQ model* exp(−*αd* − *βd*
^{2}) . *We have added this equation to the sentence to clarify it.*

2. What does variable n is Eqs (1) and (4) represent: absolute or relative number of cells?

Author’s response: *Thank you so much for bringing this subtle point. We modified the definition of n for both equations, which is absolute in Eq (1) and density in Eq (4).*

3. In Eq. (1), is the distribution of radiosensitivity fixed or changes in time?

Author’s response: *Eq (1) represents a standard Ornstein-Uhlenbeck process for a cell population of fixed size undergoing “convection” and “diffusion” in a “radiation sensitivity space” parametrized by α and centered on α*
_{0}
*. We added few words to point it out in the text.*

4. It follows from Eq. (2) that for sufficiently large sigma the number of cancer cells will eventually exceed N(0). How could this happen in the absence of cell proliferation? Also, for sigma = 0 the formula does not coincide with the LQ model. Finally, what is the meaning of T? The two questions about large sigma and sigma = 0 relate to Eq. (3) as well.

Author’s response: *Thank you so much for your comment, there were two typos in these two equations which are fixed. Also we defined the parameter T in the text.*

5. P. 5, lines 24-26. This observation is unclear.

Author’s response: *We have added two references to support this statement and also few sentences to clarify it.*

6. P. 5, line 60 and p. 6, line 4. Due to cell arrest in the most sensitive phases of cell cycle, protracted radiation promotes synchronization.

Author’s response: *Thank you, we have added a sentence to the article about this topic*.

7. P. 6, Eq. (4). What is g(a)? Also, n(a, t) on line 22 should be n(a, t) da.

Author’s response: *Function g(a) is the function modeling mitosis rate at age a. We have added to the text and changed n(a,t) to n(a,t)da.*

8. P. 7, line 24. Shouldn’t t be an argument of function h rather than A and Q?

Author’s response: *Thank you so much, we have addressed this issue.*

9. P. 7, lines 29-41. This whole paragraph is obscure. What are the assumptions here and how does the argument work?

Author’s response: *We have added a few sentences to explain the comparison of the active-quiescent model by Dawson and Hillen and the LQ model.*

10. P. 8, Eqs (6) and (7). Are the linearity assumptions compatible with the physics of oxygen diffusion? Also, does this imply that the rate of spontaneous cell death is decreasing with time as found in [11, 12]?

Author’s response: *Thank you so much for your comment. We agree with your point about linearity assumption. We have added a few sentences after equation (6) and (7) to explain these shortcomings. Also as it is mentioned in the original article [wein et al. 2000], based on equation (6) and (7), the tumor’s net repopulation rate increases as the tumor shrinks, which is consistent with the experimental evidence showing that the growth fraction and sensitivity in solid tumors decrease as the distance from the nutrient supply increases, therefore tumor death rate decreases with time based on the equation (6) and (7).*

11. P. 9, line 16. Studied in Hanin et al. 1993 were the effects of radiosensitivity variation among cancer cells without any spatial considerations. A more detailed discussion was presented in the book [13].

Author’s response: *Thank you so much for pointing this out. The spatial term has been removed.*

12. P. 9, line 39. The paper Dick 1997 deals with acute myeloid leukemia that does not form tumors.

Author’s response: *Thank you so much for pointing this out. We removed this paper in that sentence.*

13. P. 10, line 17. s(d) should be S(d).

Author’s response: *Corrected.*

14. P. 10, line 34. “…total volume of tumor with respect to a desired volume…” What do you mean?

Author’s response: *The function N*(*t*) *is the total volume of tumor normalized between 0 and 1 which is equal to N*
_{
s
}(*t*) + *N*
_{
d
}(*t*)*. We modified it in the text.*

15. P. 10, lines 55-56. Was such dedifferentiation observed and what is its mechanism?

Author’s response: *This topic is discussed in the referenced manuscript*.

16. P. 10, line 57. Beta depends on many kinetic parameters accounting for damage production, repair, misrepair and pairwise interaction, see [4]. Therefore, the stated proportionality does not seem likely and, in any case, requires discussion of the underlying assumptions.

Author’s response: *This is simplifying assumption to enable authors to characterize the radio-sensitivity of CSC by a single parameter. We have added few sentences to clarify this.*

17. P. 15, line 34. What are the Simpson’s and Shannon’s indices?

Author’s response: *Definitions were added.*

References

1. Hanin L and Zaider M (2013). A mechanistic description of radiation-induced damage to normal tissue and its healing kinetics. Phys Med Biol 58: 825-839.

2. Hanin LG, Hyrien O, Bedford J and Yakovlev AY (2006). A comprehensive stochastic model of irradiated cell populations in culture. J Theor Biol 239(4): 401-416.

3. Hanin L and Zaider M (2010). Cell-survival probability at large doses: an alternative to the linear-quadratic model. Phys Med Biol 55: 4687-4702.

4. Sachs RK, Hahnfeld P and Brenner DJ (1997). The link between low-LET dose-response relations and the underlying kinetics of damage production/ repair/misrepair. Int J Radiat Biol 72: 351–74.

5. Sachs RK, Heidenreich WF, Brenner DJ (1996). Dose timing in tumor radiotherapy: considerations of cell number stochasticity. Math Biosci 138: 131-146.

6. Fakir H, Hlatky L, Li H and Sachs R (2013). Repopulation of interacting tumor cells during fractionated radiotherapy: stochastic modeling of the tumor control probability. Med Phys 40(12):121716.

7. Hanin L and Zaider M (2014). Optimal schedules of fractionated radiation therapy by way of the greedy principle: biologically-based adaptive boosting, Phys Med Biol 59: 4085-4098.

8. Zaider M and Hanin L (2011). Tumor Control Probability in radiation treatment. Med Phys 38 (2): 574-583.

9. Hanin LG (2001). Iterated birth and death process as a model of radiation cell survival. Math Biosci 169(1): 89-107.

10. Hanin LG (2004). A stochastic model of tumor response to fractionated radiation: limit theorems and rate of convergence. Math Biosci 191: 1–17.

11. Jones B and Dale RG (1995). Cell loss factors and the linear-quadratic model. Radiother Oncol 37:136-139.

12. Ang KK, Thames HD, Jones SD, Jiang G-L, Milas L and Peters LJ (1988). Proliferation kinetics of a murine fibrosarcoma during fractionated irradiation. Radiat Res 116: 327-336.

13. Hanin LG, Pavlova LV and Yakovlev AY (1994). Biomathematical Problems in Optimization of Cancer Radiotherapy. CRC Press, Boca Raton, FL.

## Declarations

### Acknowledgments

Both HB and KL were supported by NSF grant CMMI- 1362236.

### Authors’ contributions

HB and KL wrote the manuscript. Both authors read and approved the manuscript.

### Competing interests

The authors declare that they have no competing interests.

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Al-Hajj M, Wicha MS, Benito-Hernandez A, Morrison SJ, Clarke MF. Prospective identication of tumorigenic breast cancer cells. Proc Natl Acad Sci. 2003;100(7):3983–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Antipas VP, Stamatakos GS, Uzunoglu NK, Dionysiou DD, Dale RG. A spatio-temporal simulation model of the response of solid tumors to radiotherapy in vivo: parametric validation concerning oxygen enhancement ratio and cell cycle duration. Phys Med Biol. 2004;49(8):1485.View ArticleGoogle Scholar
- Athreya K, Ney P. Branching processes. Dover Books on Mathematics Series. Mineola: Dover Publications; 2004.Google Scholar
- Bachman JW, Hillen T. Mathematical optimization of the combination of radiation and differentiation therapies for cancer. Front Oncol. 2013;3:52.View ArticlePubMedPubMed CentralGoogle Scholar
- Badri H, Pitter K, Holland E, Michor F, Leder K. Optimization of radiation dosing schedules for proneural glioblastoma. Journal of Mathematical Biology. 2016;72(5):1–36Google Scholar
- Badri H, Watanabe Y, Leder K. Optimal radiotherapy dose schedules under parametric uncertainty. Phys Med Biol. 2015;61(1):338.View ArticlePubMedGoogle Scholar
- Badri H, Ramakrishnan J, Leder K. Minimizing metastatic risk in radiotherapy fractionation schedules. Phys Med Biol. 2015;60(22):N405.View ArticlePubMedGoogle Scholar
- Badri H, Salari E, Watanabe Y, Leder K. Optimizing chemoradiotherapy to target multi-site metastatic disease and tumor growth. 2016. http://arxiv.org/pdf/1603.00349.pdf. Accessed June 2016
- Bao S, Wu Q, McLendon RE, Hao Y, Shi Q, Hjelmeland AB, Dewhirst MW, Bigner DD, Rich JN. Glioma stem cells promote radioresistance by preferential activation of the DNA damage response. Nature. 2006;444(7120):756–60.View ArticlePubMedGoogle Scholar
- Baumann M, Krause M, Thames H, Trott K, Zips D. Cancer stem cells and radiotherapy. Int J Radiat Biol. 2009;85(5):391–402.View ArticlePubMedGoogle Scholar
- Bernhard EJ, Maity A, Muschel RJ, McKenna WG. Effects of ionizing radiation on cell cycle progression. Radiat Environ Biophys. 1995;34(2):79–83.View ArticlePubMedGoogle Scholar
- Bouchat V, Nuttens VE, Michiels C, Masereel B, Feron O, Gallez B, Vander Borght T, Lucas S. Radioimmunotherapy with radioactive nanoparticles: biological doses and treatment efficiency for vascularized tumors with or without a central hypoxic area. Med Phys. 2010;37(4):1826–39.View ArticlePubMedGoogle Scholar
- Brenner D. The linear-quadratic model is an appropriate methodology for determining isoeffective doses at large doses per fraction. Seminars Radiation Oncology. 2008;18:234–9.Google Scholar
- Brenner DJ, Hlatky LR, Hahnfeldt PJ, Hall EJ, Sachs RK. A convenient extension of the linear-quadratic model to include redistribution and reoxygenation. Int J Radiat Oncol Biol Phys. 1995;32(2):379–90.View ArticlePubMedGoogle Scholar
- Brown JM, Carlson DJ, Brenner DJ. The tumor radiobiology of SRS and SBRT: are more than the 5 Rs involved? International Journal of Radiation Oncology* Biology* Physics. 2014;88(2):254–262.Google Scholar
- Buffa FM, West C, Byrne K, Moore JV, Nahum AE. Radiation response and cure rate of human colon adenocarcinoma spheroids of different size: the significance of hypoxia on tumor control modelling. Int J Radiat Oncol Biol Phys. 2001;49(4):1109–18.View ArticlePubMedGoogle Scholar
- Chen PL, Brenner DJ, Sachs RK. Ionizing radiation damage to cells: effects of cell cycle redistribution. Math Biosci. 1995;126(2):147–70.View ArticlePubMedGoogle Scholar
- Chen JH, Kuo YH, Luh HP. Optimal policies of non-cross-resistant chemotherapy on Goldie and Coldmans cancer model. Math Biosci. 2013;245:282–98.View ArticlePubMedGoogle Scholar
- Coldman AJ, Goldie JH. A model for the resistance of tumor cells to cancer chemotherapeutic agents. Math Biosci. 1983;65:291–307.View ArticleGoogle Scholar
- Coldman AJ, Murray JM. Optimal control for a stochastic model of cancer chemotherapy. Math Biosci. 2000;168:187–200.View ArticlePubMedGoogle Scholar
- Conger AD, Ziskin MC. Growth of mammalian multicellular tumor spheroids. Cancer Res. 1983;43(2):556–60.PubMedGoogle Scholar
- Cover TM, TM, Thomas JA. Elements of information theory. New York, USA: John Wiley & Sons; 2012Google Scholar
- Dawson A, Hillen T. Derivation of the tumor control probability (TCP) from a cell cycle model. Computational and Mathematical Methods in Medicine. 2006;7(2-3):121–41.View ArticleGoogle Scholar
- Day RS. Treatment sequencing, asymmetry, and uncertainty: protocol strategies for combination chemotherapy. Cancer Res. 1986;46:3876–85.PubMedGoogle Scholar
- A. Dhawan, K. Kaveh, M. Kohandel, S. Sivaloganathan. Stochastic model for tumor control probability: effects of cell cycle and (a) symmetric proliferation. arXiv preprint arXiv:1312.7556, 2013Google Scholar
- Dick D. Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell. Nature Med. 1997;3:730–7.View ArticlePubMedGoogle Scholar
- Dionysiou DD, Stamatakos GS. Applying a 4d multiscale in vivo tumor growth model to the exploration of radiotherapy scheduling: the effects of weekend treatment gaps and p53 gene status on the response of fast growing solid tumors. Cancer Informat. 2006;2:113.Google Scholar
- Dionysiou DD, Stamatakos GS, Uzunoglu NK, Nikita KS, Marioli A. A four-dimensional simulation model of tumor response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation. J Theor Biol. 2004;230(1):1–20.View ArticlePubMedGoogle Scholar
- Durrett R. Branching process models of cancer. Cham, Switzerland: Springer; 2015.Google Scholar
- Durrett R, Moseley S. Evolution of resistance and progression to disease during clonal expansion of cancer. Theor Popul Biol. 2010;77(1):42–8.View ArticlePubMedGoogle Scholar
- Durrett R, Foo J, Leder K, Mayberry J, Michor F. Evolutionary dynamics of tumor progression with random fitness values. Theor Popul Biol. 2010;78:54–66.View ArticlePubMedPubMed CentralGoogle Scholar
- Enderling H, Chaplain MA, Hahnfeldt P. Quantitative modeling of tumor dynamics and radiotherapy. Acta Biotheor. 2010;58(4):341–53.View ArticlePubMedGoogle Scholar
- Fla T, Rupp F, Woywod C. Deterministic and stochastic dynamics of chronic myelogenous leukaemia stem cells subject to hill-function-like signaling. In: Recent Trends in Dynamical Systems, pages 221-263. Basel, Switzerland: Springer; 2013.Google Scholar
- Fletcher GH. Textbook of radiotherapy. Philadelphia, USA: Lea & Febiger; 1973Google Scholar
- Foo J, Leder K. Dynamics of cancer recurrence. Ann Appl Probab. 2013;23:1437–68.View ArticleGoogle Scholar
- Foo J, Michor F. Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies. PLoS Comput Biol. 2009;5(11):e1000557.View ArticlePubMedPubMed CentralGoogle Scholar
- Foo J, Michor F. Evolution of resistance to anti-cancer therapy during general dosing schedules. J Theor Biol. 2010;263(2):179–88.View ArticlePubMedGoogle Scholar
- Foo J, Leder K, Mumenthaler S. Cancer as a moving target: understanding the composition and rebound growth kinetics of recurrent tumors. Evol Appl. 2013;6:54–69.View ArticlePubMedGoogle Scholar
- Foo J, Leder K, Zhu J. Escape times for branching processes with random mutational fitness effects. Stochastic Processes and Their Applications. 2014;124:3661–97.View ArticleGoogle Scholar
- Fowler JF. The phantom of tumor treatment-continually rapid proliferation unmasked. Radiother Oncol. 1991;22(3):156–8.View ArticlePubMedGoogle Scholar
- Gao X, McDonald JT, Hlatky L, Enderling H. Acute and fractionated irradiation differentially modulate glioma stem cell division kinetics. Cancer Res. 2013;73(5):1481–90.View ArticlePubMedGoogle Scholar
- Gerlinger M, Horswell S, Larkin J, Rowan AJ, Salm MP, Varela I, Fisher R, McGranahan N, Matthews N, Santos CR, et al. Genomic architecture and evolution of clear cell renal cell carcinomas defined by multi-region sequencing. Nat Genet. 2014;46(3):225–33.View ArticlePubMedPubMed CentralGoogle Scholar
- Goldie JH, Coldman AJ. A mathematical model relating the drug sensitivity of tumors to their spontaneous mutation rate. Cancer Treat Rep. 1979;63:1727–33.PubMedGoogle Scholar
- Goldie JH, Coldman AJ, Gudauskas GA. A rationale for the use of alternating noncross resistant chemotherapy. Cancer Treat Rep. 1982;66:439–49.PubMedGoogle Scholar
- Gray LH, Conger AD, Ebert M, Hornsey S, Scott O. The concentration of oxygen dissolved in tissues at the time of irradiation as a factor in radiotherapy. Br J Radiol. 1953;26(312):638–48.View ArticlePubMedGoogle Scholar
- Haas-Kogan DA, Yount G, Haas M, Levi D, Kogan SS, Hu L, Vidair C, Deen DF, Dewey WC, Israel MA. p53-dependent G1 arrest and p53-independent apoptosis in influence the radiobiologic response of glioblastoma. Int J Radiat Oncol Biol Phys. 1996;36(1):95–103.View ArticlePubMedGoogle Scholar
- Haeno H, Iwasa Y, Michor F. The evolution of two mutations during clonal expansion. Genetics. 2007;177(4):2209–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Hahnfeldt P, Hlatky L. Resensitization due to redistribution of cells in the phases of the cell cycle during arbitrary radiation protocols. Radiat Res. 1996;145(2):134–43.View ArticlePubMedGoogle Scholar
- Hall E, Giaccia A. Radiobiology for the radiologist. Philadelphia, USA: Wolters Kluwer Health; 2006Google Scholar
- Hambardzumyan D, Becher OJ, Rosenblum M, Pandol PP, Manova-Todorova K, Holland EC. Pi3k pathway regulates survival of cancer stem cells residing in the perivascular niche following radiation in medulloblastoma in vivo. Genes Dev. 2008;22(4):436–48.View ArticlePubMedPubMed CentralGoogle Scholar
- Hanin LG. Iterated birth and death process as a model of radiation cell survival. Math Biosci. 2001;169(1):89–107.View ArticlePubMedGoogle Scholar
- Hanin LG. A stochastic model of tumor response to fractionated radiation: limit theorems and rate of convergence. Math Biosci. 2004;191:1–17.View ArticlePubMedGoogle Scholar
- Hanin LG, Zaider M. Cell-survival probability at large doses: an alternative to the linear-quadratic model. Phys Med Biol. 2010;55:4687–702.View ArticlePubMedGoogle Scholar
- Hanin LG, Zaider M. A mechanistic description of radiation-induced damage to normal tissue and its healing kinetics. Phys Med Biol. 2013;58(4):825–39.View ArticlePubMedGoogle Scholar
- Hanin LG, Zaider M. Optimal schedules of fractionated radiation therapy by way of the greedy principle: biologically-based adaptive boosting. Phys Med Biol. 2014;59:4085–98.View ArticlePubMedGoogle Scholar
- Hanin L, Rachev S, Yakovlev AY. On the Optimal Control of Cancer Radiotherapy for Non-Homogeneous Cell Populations. Advances in Applied Probability. 1993;25(1):1–23. doi: 10.2307/1427493.
- Hanin LG, Hyrien O, Bedford J, Yakovlev AY. A comprehensive stochastic model of irradiated cell populations in culture. J Theor Biol. 2006;239(4):401–16.View ArticlePubMedGoogle Scholar
- Hendry JH, Moore JV. Is the steepness of dose-incidence curves for tumor control or complications due to variation before, or as a result of, irradiation? Br J Radiol. 1984;57(683):1045–6.View ArticlePubMedGoogle Scholar
- Hillen T, VrIeS GD, Gong J, Finlay C. From cell population models to tumor control probability: including cell cycle effects. Acta Oncol. 2010;49(8):1315–23.View ArticlePubMedGoogle Scholar
- Hlatky LR, Hahnfeldt P, Sachs RK. Influence of time-dependent stochastic heterogeneity on the radiation response of a cell population. Math Biosci. 1994;122(2):201–20.View ArticlePubMedGoogle Scholar
- Iwasa Y, Nowak MA, Michor F. Evolution of resistance during clonal expansion. Genetics. 2006;172(4):2557–66.View ArticlePubMedPubMed CentralGoogle Scholar
- Katouli AA, Komarova NL. The worst drug rule revisited: mathematical modeling cyclic cancer treatments. Bull Math Biol. 2011;73:549–84.View ArticlePubMedGoogle Scholar
- Kimmel M, Axelrod D. Branching processes in biology. 2nd ed. Springer-Verlag; 2015Google Scholar
- Kimmel M, Swiernak A. Control theory approach to cancer chemotherapy: benefiting from phase dependence and overcoming drug resistance. Lect Notes Math. 2006;1872:185–221.View ArticleGoogle Scholar
- Komarova N. Stochastic modeling of drug resistance in cancer. J Theor Biol. 2006;239:351–66.View ArticlePubMedGoogle Scholar
- Komarova N, Wodarz D. Drug resistance in cancer: principles of emergence and prevention. Proc Natl Acad Sci U S A. 2005;102(27):9714–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Leder K, Pitter K, LaPlant Q, Hambardzumyan D, Ross B, Chan T, Holland E, Michor F. Mathematical modeling of PDGF-driven glioblastoma reveals optimized radiation dosing schedules. Cell. 2014;156(3):603–16.View ArticlePubMedPubMed CentralGoogle Scholar
- Maler A, Lutscher F. Cell-cycle times and the tumor control probability. Mathematical Medicine and Biology, 2009. doi:10.1093/imammb/dqp024.
- Martin R, Teo KL. Optimal control of drug administration in cancer chemotherapy. New Jersey, USA: World Scientific; 1994Google Scholar
- Mizuta M, Takao S, Date H, Kishimotoi N, Sutherland K, Onimaru R, Shirato H. A mathematical study to select fractionation regimen based on physical dose distribution and the linear-quadratic model. Int J Radiat Oncol Biol Phys. 2012;84(3):829–33.View ArticlePubMedGoogle Scholar
- Munro TR, Gilbert CW. The relation between tumor lethal doses and the radiosensitivity of tumor cells. Br J Radiol. 1961;34(400):246–51.View ArticlePubMedGoogle Scholar
- Murray JM, Coldman AJ. The effect of heterogeneity on optimal regimens in cancer chemotherapy. Math Biosci. 2003;185:73–87.View ArticlePubMedGoogle Scholar
- ORourke SFC, McAneney H, Hillen T. Linear quadratic and tumor control probability modelling in external beam radiotherapy. J Math Biol. 2009;58(4-5):799–817.View ArticleGoogle Scholar
- Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344(6190):1396–401.View ArticlePubMedPubMed CentralGoogle Scholar
- Phillips TM, McBride WH, Pajonk F. The response of CD24-/low/CD44+ breast cancer-initiating cells to radiation. J Natl Cancer Inst. 2006;98(24):1777–85.View ArticlePubMedGoogle Scholar
- Rachev ST, Yakovlev AY. Theoretical bounds for the tumor treatment efficacy. Syst Anal Model Simul. 1988;5(1):37–42.Google Scholar
- Reya T, Morrison SJ, Clarke MF, Weissman IL. Stem cells, cancer, and cancer stem cells. Nature. 2001;414(6859):105–11.View ArticlePubMedGoogle Scholar
- Rockwell S, Dobrucki IT, Kim EY, Marrison ST, Vu VT. Hypoxia and radiation therapy: past history, ongoing research, and future promise. Current Molecular Medicine. 2009;9(4):442.View ArticlePubMedPubMed CentralGoogle Scholar
- Ruggieri R. Hypofractionation in non-small cell lung cancer (NSCLC): suggestions from modelling both acute and chronic hypoxia. Phys Med Biol. 2004;49(20):4811.View ArticlePubMedGoogle Scholar
- Sachs RK, Chen PL, Hahnfeldt PJ, Hlatky LR. DNA damage caused by ionizing radiation. Math Biosci. 1992;112(2):271–303.View ArticlePubMedGoogle Scholar
- Sachs RK, Hahnfeldt PJ, Brenner DJ. The link between low-LET dose-response relations and the underlying kinetics of damage production/ repair/misrepair. Int J Radiat Biol. 1997;72:351–74.View ArticlePubMedGoogle Scholar
- Shah SP, Morin RD, Khattra J, Prentice L, Pugh T, Burleigh A, Delaney A, Gelmon K, Guliany R, Senz J, et al. Mutational evolution in a lobular breast tumor profiled at single nucleotide resolution. Nature. 2009;461(7265):809–13.View ArticlePubMedGoogle Scholar
- Shi J, Alagoz O, Erenay F, Su Q. A survey of optimization models on cancer chemotherapy treatment planning. Ann Oper Res. 2014;221:331–56.View ArticleGoogle Scholar
- Singh SK, Clarke ID, Terasaki M, Bonn VE, Hawkins C, Squire J, Dirks PB. Identification of a cancer stem cell in human brain tumors. Cancer Res. 2003;63(18):5821–8.PubMedGoogle Scholar
- Steel GG, McMillan TJ, Peacock JH. The 5Rs of radiobiology. Int J Radiat Biol. 1989;56(6):1045–8.View ArticlePubMedGoogle Scholar
- Swan GW. Role of optimal control theory in cancer chemotherapy. Math Biosci. 1990;101:237–84.View ArticlePubMedGoogle Scholar
- Tannock IF. The relation between cell proliferation and the vascular system in a transplanted mouse mammary tumor. Br J Cancer. 1968;22(2):258.View ArticlePubMedPubMed CentralGoogle Scholar
- Thames HD, Hendry JH. Fractionation in radiotherapy. London: Taylor and Francis; 1987.Google Scholar
- Thomlinson RH, Gray LH. The histological structure of some human lung cancers and the possible implications for radiotherapy. Br J Cancer. 1955;9(4):539.View ArticlePubMedPubMed CentralGoogle Scholar
- Tucker SL, Thames HD. The effect of patient-to-patient variability on the accuracy of predictive assays of tumor response to radiotherapy: a theoretical evaluation. Int J Radiat Oncol Biol Phys. 1989;17(1):145–57.View ArticlePubMedGoogle Scholar
- Turesson I. Radiobiological aspects of continuous low dose-rate irradiation and fractionated high dose-rate irradiation. Radiother Oncol. 1990;19(1):1–15.View ArticlePubMedGoogle Scholar
- Unkelbach J, Craft D, Salari E, Ramakrishnan J, Bortfeld T. The dependence of optimal fractionation schemes on the spatial dose distribution. Phys Med Biol. 2013;58(1):159.View ArticlePubMedGoogle Scholar
- Victoria YY, Nguyen D, Pajonk F, Kupelian P, Kaprealian T, Selch M, Low DA, Sheng K. Incorporating cancer stem cells in radiation therapy treatment response modeling and the implication in glioblastoma multiforme treatment resistance. Int J Radiat Oncol Biol Phys. 2015;91(4):866–75.View ArticleGoogle Scholar
- Wein LM, Cohen JE, Wu JT. Dynamic optimization of a linear-quadratic model with incomplete repair and volume-dependent sensitivity and repopulation. Int J Radiat Oncol Biol Phys. 2000;47(4):1073–83.View ArticlePubMedGoogle Scholar
- Williams T, Bjerknes R. Stochastic model for abnormal clone spread through epithelial basal layer. Nature. 1972;236:19–21.View ArticlePubMedGoogle Scholar
- Withers H. Four R’s of radiotherapy. Adv Biol. 1975;5:241–7.Google Scholar
- Withers HR. Biological basis of radiation therapy for cancer. Lancet. 1992;339(8786):156–9.View ArticlePubMedGoogle Scholar
- Woodward WA, Chen MS, Behbod F, Alfaro MP, Buchholz TA, Rosen JM. Wnt/β-catenin mediates radiation resistance of mouse mammary progenitor cells. Proc Natl Acad Sci. 2007;104(2):618–23.View ArticlePubMedPubMed CentralGoogle Scholar
- Yang Y, Xing L. Optimization of radiotherapy dose-time fractionation with consideration of tumor specific biology. Med Phys. 2005;32:3666.View ArticlePubMedGoogle Scholar
- Youssefpour H, Li X, Lander AD, Lowengrub JS. Multispecies model of cell lineages and feedback control in solid tumors. J Theor Biol. 2012;304:39–59.View ArticlePubMedPubMed CentralGoogle Scholar
- Zagars GK, Schultheiss TE, Peters LJ. Inter-tumor heterogeneity and radiation dose-control curves. Radiother Oncol. 1987;8(4):353–61.View ArticlePubMedGoogle Scholar
- Zaider M, Hanin LG. Tumor Control Probability in radiation treatment. Med Phys. 2011;38(2):574–83.View ArticlePubMedGoogle Scholar
- Zaider M, Minerbo GN. Tumor control probability: a formulation applicable to any temporal protocol of dose delivery. Phys Med Biol. 2000;45(2):279.View ArticlePubMedGoogle Scholar