Stochastic model for tumor control probability: effects of cell cycle and (a)symmetric proliferation

Background Estimating the required dose in radiotherapy is of crucial importance since the administrated dose should be sufficient to eradicate the tumor and at the same time should inflict minimal damage on normal cells. The probability that a given dose and schedule of ionizing radiation eradicates all the tumor cells in a given tissue is called the tumor control probability (TCP), and is often used to compare various treatment strategies used in radiation therapy. Method In this paper, we aim to investigate the effects of including cell-cycle phase on the TCP by analyzing a stochastic model of a tumor comprised of actively dividing cells and quiescent cells with different radiation sensitivities. Moreover, we use a novel numerical approach based on the method of characteristics for partial differential equations, validated by the Gillespie algorithm, to compute the TCP as a function of time. Results We derive an exact phase-diagram for the steady-state TCP of the model and show that at high, clinically-relevant doses of radiation, the distinction between active and quiescent tumor cells (i.e. accounting for cell-cycle effects) becomes of negligible importance in terms of its effect on the TCP curve. However, for very low doses of radiation, these proportions become significant determinants of the TCP. We also present the results of TCP as a function of time for different values of asymmetric division factor. Conclusion We observe that our results differ from the results in the literature using similar existing models, even though similar parameters values are used, and the reasons for this are discussed.


Introduction
External beam radiotherapy remains one of the most common treatment options for various cancers. However, the dose distribution of radiation must be optimized to reduce the risk of side effects of radiotoxicity and damage to healthy tissues surrounding the tumor volume. A widely used model for radiation treatment is the linear-quadratic (LQ) model [1,2]. This model estimates the surviving fraction of cancer cells after each treatment based on the total dose, and has the form: where α and β are sensitivity parameters (which depend on the tissue and the type of the applied beam) and D is the total dose delivered during the radiation treatment. To include stochastic effects, a binomial or Poisson model has been used to describe the random variable representing the number of surviving cells after a treatment, centered upon a mean value determined by the linear-quadratic model of cell survival (see, for example, [3,4]). An iterated birth and death process has been also suggested as a model of radiation cell survival [5]. A related quantity of interest is the tumor control probability (TCP) which is the extinction probability of the clonogenic cell population after radiation therapy. A model for the TCP accounting for cell proliferation dynamics was suggested by [6]. Their model is a birth-death process for the probability distribution function of the tumor cells, p n (t), and the corresponding master equation of such a birth-death model is: dp where λ and ζ are the birth and death rates, respectively, and n is the population of tumor cells. The effect of radiation is reflected as a time-dependent part in the death rate, where h(t) is known as the hazard function and is related to the radiosensitivity parameters α and β through the LQ model (Eq. 1). From Eq. 2, Zaider and Minerbo were able to calculate the extinction probability, p 0 (t), as a function of time and dose fractions (which is encoded in the form of h(t)). Thus, in their model, the TCP is given by: where n 0 is the initial number of tumor cells and S(t) is the exponential of the integral of the hazard function: with D(t) being the dose in Gy delivered until time t and its time derivative representing dose rate (Gy/day). Paper [7] expanded this approach to include the effect of cell cycle sensitivity in the TCP. They considered a two-compartment model for the active (M, G 1 , S, and G 2 phases) and the quiescent (G 0 phase) cells (see also [8]). The radio-sensitivity of resting cells and active cells are significantly different; the radio-sensitivity is typically much higher for actively proliferating cells [9]. This model was discussed both deterministically and stochastically in [7], but the stochastic master equation is solved under the assumption that the joint probability distribution function of two populations, p n a ,n q , can be written in a factorized form as if the two random variables n a and n q are independent. However, this is clearly not true for small tumor populations, as pointed out by [10]. Small tumor populations can arise from a number of possible clinically relevant scenarios; for example, this would be the case for adjuvant radiation applied after surgery or chemotherapy, irradiation of micrometastases, as well as at the final stages of radiation therapy, when the tumor has shrunk to a few milimeters in size. Thus, as one approaches the limit of small tumor cell http://www.tbiomed.com/content/11/1/49 populations, a proper stochastic approach is needed to estimate the extinction probability, i.e. the TCP. Moreover, in previous cell cycle models of the TCP [7,10], it is assumed that the proliferation is such that upon each cell division the daughter cells go into the G 0 (quiescent) state soon thereafter. In the following, we consider a more general situation where there is a probability f , such that one of the daughter cells goes into the resting phase upon division [11]; the master equation is again solved with the same assumption of independent random variables for the subpopulations of cells which breaks down in the key limit of small cell populations. In the following, we investigate thoroughly the TCP for such a model throughout the range of pertinent parameter values and plot a phase diagram of the model using a generating function method (see Section 'Stochastic twocompartment model with (a)symmetric proliferation'). In Section 'Numerical solutions: final-value method', we solve the differential equation for a probability generating function for the number of tumor cells using a novel final-value method of characteristics and in Section 'Gillespie solution' we validate this with a Gillespie algorithm solution of the master equation.

Stochastic two-compartment model with (a)symmetric proliferation
Here we consider a two compartment model of active cells (A) and quiescent cells (Q), with the following dynamics: active cells can divide into either: (1) two quiescent cells or (2) one quiescent and one active, or (3) two active cells; assuming each active offspring is born with probability f and each quiescent with probability 1 − f while the proliferation rate for active cells is μ. Note also that quiescent cells may, after a certain time, move from the G 0 to the G 1 phase of the cell cycle, and thereby become active. We assume this happens at a constant rate γ . Death rates for the cells in the active and quiescent compartments are denoted by a and q , respectively: The deterministic ordinary differential equations (ODEs) for the above dynamics are given by: where n a and n q are the population of the active and quiescent compartments. The death rates of active and quiescent cells, a (t) and q (t), are dose-dependent through the LQ formula (Eq. 1) and the given radiation protocol. Similarly, we can determine the stochastic dynamics of the model Eq. 5 as follows. Denoting the joint probability distribution of http://www.tbiomed.com/content/11/1/49 having a population of n a active cells and n q of quiescent cells at time t by p n a ,n q (t), the master equation then reads, dp n a ,n q (t) + a (n a + 1)p n a +1,n q (t) + q (n q + 1)p n a ,n q +1 (t) − ( a + μ)n a p n a ,n q (t) − ( q + γ )n q p n a ,n q (t).
The model in [7] and [10] corresponds to f = 0 in Eq. 7, while the Zaider and Minerbo model [6] corresponds to f = 1. We define the probability generating function for the joint probability distribution, p n a ,n q , p n a ,n q (t)a n a q n q .
Using Eq. 7 and Eq. 8, we obtain the following partial differential equation (PDE) for V (a, q, t): Taking n a,0 and n q,0 to be the initial numbers of active and quiescent cells, respectively, we have the initial condition V (a, q, 0) = a n a,0 q n q,0 and the boundary condition V (1, 1, t) = 1, where the boundary condition comes from the definition of the generating function.
The TCP is defined as the extinction probability of the tumor cells after radiation therapy, i.e. TCP(t) = p n a =0,n q =0 (t) = V (0, 0, t). In the steady state, TCP(∞) = lim t→∞ V (0, 0, t). In the case of a constant radiation dose, we can find an analytical solution that relates TCP to all the parameters appearing in the model, especially the values of the death rates, a and q . In the limit of a large -but finite -total number of cells N, we expect the steady state of the system to have two absorbing states of either zero population of either active or quiescent cells or both populations together reaching their maximum limits, N a and N q (N a,q 1). This means that in the steady state, the form of the generating function V (a, q, t → ∞) is: The first term indicates that there is a non-zero probability for either population to become extinct and the second term is indicative of the possibility that eventually one or both populations reach large population limits -details to be determined by the values of N a and N q . The coefficients A and B are the extinction and survival probabilities of the dynamical system, respectively, while q * and a * are the fixed points of Eq. 9, which satisfy the following relations: (11) http://www.tbiomed.com/content/11/1/49 One fixed point is (1, 1), and the non-trivial solutions for a * and q * are: where the coefficients c 1 , c 2 , and c 3 are defined as Using the initial and boundary conditions mentioned above, we can obtain the values of A and B. We are interested in the value of A, which is the extinction probability in the long run. This is the TCP in the steady state (TCP ∞ ): In parameter space, the phase boundary can be defined in the parameter space in terms of the model parameters such as a,q , γ , and μ, when these parameters are such that (a * , q * ) = (1, 1). For the region of the phase diagram where TCP ∞ = 0, the (a * , q * ) fixed point is attractive while the (1, 1) fixed point is a saddle-point. As the parameters such as death rates a,q increase, one moves into the TCP ∞ = 1 regime where now the fixedpoint (a * , q * ) vanishes and the only fixed point is (1, 1) which is globally attractive. The phase boundary for variable death rates is plotted in Figure 1. To provide a comparison between the results of [7], we use identical parameter values, namely a constant radiation dose rate of 2.75 Gy/day, and the division rate μ and the conversion rate γ were taken to be 0.065 day −1 [12] and 0.047 day −1 [13], respectively. Death rates, which were effectively derived from a limit of the LQ model, are given by: q = 0.4/Day and a = 1.5/Day. The death rates were derived by using the dose-dependent survival fraction given by the LQ model, creating a hazard function from that, and substituting in values for the radiosensitivity parameters α a = 0.487Gy −1 , α q = 0.155Gy −1 , β a = β q = 0.055Gy −2 taken from [9], where the subscript a or q indicates active cells or quiescent cells, respsectively. Also, note that a constant radiation dose is not necessarily a clinical possibility for treatment, but is used in order to facilitate direct comparison with the results of [7].
Our plots in Figure 1 for the phase boundary between TCP = 0 and TCP = 1.0 regimes show the interesting evolution of the two regimes of the one-compartment model of [6] into the two-compartment model of [7]. It can be noted that the two ends of the phase boundary at the a -axis and q -axis are in fact μ and γ for the fully two-compartment model ( f = 0), i.e. the values for the cutoff death rates are determined by the proliferation and conversion potentials μ and γ . For values of a,q 's in these regions one expects to get an unsuccessful therapy or TCP ∞ = 0. The implication of this is the fact that the values of a,q estimated for real irradiation protocols lie deep inside the TCP ∞ = 1 phase for all the values of the asymmetric proliferation factor, f , and thus the division of the population into different compartments based on the cell-cycle has a negligible effect on the TCP, http://www.tbiomed.com/content/11/1/49 given that the single and two compartment models utilize identical parameters. That is, given a real treatment schedule, the effect of f on the TCP curve itself becomes negligible.
We have also plotted the phase boundary for TCP ∞ = 0, 1 for different values of the division and conversion rates μ and γ in Figure 2. A similar evolution between a one-compartment and two-compartment model can be observed in this case. The phase boundaries for the TCP = 0 and TCP = 1 regimes can be used to determine a crude cutoff dose below which treatments will not work, and above which treatments will work in finite time. However, we note that for clinical treatments, parameter values must be deep inside the TCP = 1 regime to succeed within a reasonable timescale. In the next two sections we will focus on the time-dependence of the TCP via two different approaches.

Numerical solutions: final-value method
In the previous section, we discussed the steady-state behavior and the fixed points of Eq. 9. In this section, we derive the time dependence of the TCP as it approaches unity for a given radiation protocol. Solving (Eq. 9), i.e. the PDE for the generating function, with a combination of initial and boundary conditions is a difficult task. We approach the problem by a novel application of the method of characteristics. Consider a PDE of the form: Recall that the method of characteristics relies upon finding a set of characteristic curves t(s), x 1 (s), · · · , x n (s) such that f (s) = V (x 1 (s), x 2 (s), · · · , x n (s), t(s)) is a constant.  By comparing the form of this differential equation with the form of the equation we wish to solve, we observe that to find these characteristic curves, the following set of ordinary differential equations must be solved: . . .
Note that we constrain t(0) = 0, so that the initial conditions of the system can be used in the calculation of f (0). The last equation in the system, given the initial condition t(0) = 0, can be solved. Thus we obtain the following system: We also notice that for this particular set of characteristic curves, We define f 0 as the function relating the initial values of the characteristic functions to the initial conditions for the PDE. From the given initial condition for our PDE, we have This gives f (s) = f 0 (x 1 (0), x 2 (0), · · · , x n (0)).
Recall that we are only interested in the function g(t) = V (0, · · · , 0, t), and not the entire solution to the PDE since g(t) represents the extinction probability of the tumor at the time t, which is exactly the TCP. Thus, to compute g at a fixed t = t * , the only characteristic curve that needs to be considered is such as x 1 (t * ) = x 2 (t * ) = · · · = x n (t * ) = 0. We denote these characteristic curvesx 1 ,x 2 , · · · ,x n . Moreover, based on the above discussion, we observe that The valuesx i (0) are determined by the set of ODEs in (17), with the final value condition that x 1 (t * ) = x 2 (t * ) = · · · = x n (t * ) = 0. Thus, to obtain g(t), at any set of time points, the final value problem must be solved independently to obtain the initial values of the characteristic curve, which must then be substituted into the initial condition for the PDE.
We note that taking t → t * − t will transform the aforementioned final value problem into an initial value problem, where the desired values becomex i (t * ). In this case, notice that the computation of the function g(t) can be vastly simplified if the functions x n ), then observe that for every t * , the set of ODEs that must be solved is the same, and all have the same initial condition thatx i = 0. Thus, in this case, computation of the function g(t) can be done for all t in a given interval, by solving the set of coupled ODEs once. If this simplification cannot be made, then the method will still solve the PDE, but for each time point, the set of ODEs that must be solved will be different.

Gillespie solution
In order to simulate the stochastic process representing the cellular dynamics within the model framework, Gillespie's algorithm for stochastic simulation was implemented. This algorithm simulates one realization of the time evolution of the system by first computing propensities for the events that can occur at any time step (i.e. the set of cell births/deaths in the above model). Subsequently, the time before the next event occurs is computed via an exponential distribution, and the event that occurs at this time step is chosen by a distribution weighted by the total propensity of all events (i.e. the likelihood that any reaction would occur). Thus, the events occur individually, with a likelihood proportional to their individual propensity, and the times between the individual events is based on an exponential distribution of waiting times, weighted by the total propensity of all events. Each simulation describes one specific time course for the system. This is then repeated a large number of times, typically 10 5 in our simulations, and for each, an indicator function known as the treatment success indicator is defined: TS i (t) = 1 if at time t, the tumor is controlled (i.e. there are zero cells remaining), and 0 otherwise. Then, after M such simulations, the TCP function is defined to be: (21) http://www.tbiomed.com/content/11/1/49 The process to calculate TS i (t) is: (1) Compute likelihood of each cellular reaction occurring (L i for reaction R i ). (2) Sum together all likelihoods into quantity T L = i L i .
(3) Compute uniformly distributed random numbers p 1 and p 2 in the interval (0, 1). (4) Compute the next time step of a likelihood reaction, assuming exponentially distributed times dT = − ln(p 1 )/T L . (5) Update time variable by adding time step computed t = t + dT. (6) Determine which reaction to carry out: if L i−1 /T L ≤ p 1 ≤ L i /T L , carry out reaction R i . (7) Update the cellular population variables, assuming reaction R i was carried out. (8) If number of stem cells is zero, treatment success is one and terminate program, else treatment success is zero and repeat step 1. (9) If time is greater than the max simulation time, treatment success is zero and terminate program.
We illustrate the effectiveness of the numerical method presented in solving for the TCP for the active quiescent model that was outlined previously. To do this, we compare the TCP as computed by a high number of Gillespie simulation runs with the TCP as computed by the output of the numerical method.
To obtain a proper stochastic limit, we use a small number of each type of cell, letting a 0 = 10 2 = q 0 . Using these and the rest of the parameter values mentioned in Section 'Stochastic two-compartment model with (a)symmetric proliferation', we obtain the TCP plot depicted in Figure 3. In this plot, both the numerical solution, computed by an implementation of the method presented above, as well as the Gillespie solution are plotted, to highlight the high degree of similarity between the curves. In order to quantify the degree to which these curves agree, we sample both curves at the nine time points corresponding to t = 0, 3, 6, ·, 24 and compute a root-mean-square distance between the two vectors representing the TCP values of the Gillespie and numerical solutions to obtain 0.022, which is indeed very small.
Next, to check the relevance of the two-compartment model, we plot the TCP vs. time for different values of asymmetric division, f . As discussed in Section 'Stochastic two- Figure 3 Numerical results. A plot of the TCP computed by the numerical method outlined above, and by Gillepsie's algorithm. http://www.tbiomed.com/content/11/1/49 compartment model with (a)symmetric proliferation', we do not expect any difference as the physical parameters estimated from clinical data indicate a high-death rate for both the active and quiescent cells which lie deep inside the overlap region of the onecompartment and two-compartment models. As shown in Figure 4, this is in fact the case and the TCP(t; f ) plots are almost indistinguishable.
However, if one decreases both death rates, from the values in the phase diagram in Figure 1, one should expect any difference between the TCP(t; f ) to reveal itself. One example is plotted in Figure 5, with death rates a = 0.08/Gy and a = 0.1/Gy. Using the phase diagram, we can see that these values correspond to a point in the a − q plane very close to the f = 1 phase-boundary. This explains why the TCP graph for f = 1 in Figure 1 appears to approach unity on a much longer time scale than the other graphs. Similarly, one can expect the characteristic saturation time of the TCP (i.e. the time to reach unity) to tend to infinity as we choose death rates (by varying the dose of radiation) that cross the phase boundary corresponding to that asymmetric proliferation factor f .

Discussion
In this work, we have investigated a two-compartment stochastic model for the tumor control probability by including the asymmetric nature of division of active cells into either quiescent cells or active cells. We argue that the method suggested by [11] does not properly address the coupled nature of the joint probability distribution of the active and quiescent populations and have presented an alternative consistent approach. We  [7]. The value for the dose delivery rate, (2.75 Gy/day), is so high that the differences between different TCP plots are indistinguishable. http://www.tbiomed.com/content/11/1/49 have analytically derived all regimes of the phase diagram of TCP = 0, 1 in steadystate, for variable division and conversion rates and also separately the phase diagram of TCP = 0, 1 for varying death rates. From the phase diagram, we may conclude that the two-compartment model diminishes the effects of the original birth-death model of [6] while the significantly lower death rates (dose delivery rates and radio-sensitivities) can be addressed with a two-compartment model which includes cell cycle effects. The phase boundaries obtained for the TCP = 0 and TCP = 1 regimes can be used to crudely determine a dose cutoff suitable for tumor control for tumors comprised of different populations of active and quiescent cells, when death rates are low enough between treatments being compared so that parameters such as f become significant. We note that the time to achieve tumor control depends on the distance from the phase boundary, and those parameters within the TCP = 1 regime but very close to the phase boundary may not be able to achieve tumor control in a realistic time frame.
We also note now that there is a significant difference in the results computed via the method presented here and the results presented in [10] using similar parameter values. In [10], the computed TCP curve shows that the time to cure for a population of 1000 cells in total is approximately 20 hours, which is much less than the 20 days predicted by the model presented here (for a smaller population of 100 cells).
To complete the analysis we have presented a comprehensive numerical approach to compute the TCP as a function of time. The numerical method (which we call the Final Value Method), when implemented to solve the TCP problem for the above case and http://www.tbiomed.com/content/11/1/49 parameter set, can be seen to solve the PDE, producing nearly identical solutions to that of the Gillespie algorithm, which is a good approximation to the true solution. Based on the work presented here, we may conclude that the final value method is a new way to numerically solve any PDE with an initial condition that is of a form appropriate for the method of characteristics. In the case presented above, this method has been utilized to solve the real-world problem of computing the TCP for a model based on incorporating cell-cycle effects into radiotherapy treatment planning, by using a two-compartment model for the active and quiescent cells.
To compare the analytical and numerical results, we used a constant dose of radiation, which might be useful (for example) for in vitro experiments. However, our numerical methods are robust enough to account for the time dependent cases, and thus can be applied to compare the TCP for various fractionation protocols. There are several assumptions that can be improved in future works. For instance, it is experimentally known that radiation may induce cell quiescence. Thus, an extension of our model is to include a transition from survived active cells to a quiescent state, as a function of radiation administered. In addition, we considered a minimal two-compartment model to take into account the effects of the cell cycle. Thus, a significant generalization of the model is to include the details of cell cycle phases and duration. Finally, one should note that the death rates described in this paper are dose-dependent death rates for radiotherapy, but could easily be interpreted as death rates from chemotherapy for instance. In fact, it is well-known that the cytotoxic effects of chemotherapy primarily impact cells actively proliferating within the cell cycle, so here the division between active and quiescent cell populations become important. Thus, one may anticipate that the framework presented in this paper can be extended to study the effects of other treatments for tumor control, such as chemotherapy.