- Research
- Open Access

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

- Andrew Dhawan
^{1, 2}, - Kamran Kaveh
^{1}, - Mohammad Kohandel
^{1, 3}Email author and - Sivabal Sivaloganathan
^{1, 3}

**11**:49

https://doi.org/10.1186/1742-4682-11-49

© Dhawan et al.; licensee BioMed Central Ltd. 2014

**Received: **29 August 2014

**Accepted: **21 October 2014

**Published: **22 November 2014

## Abstract

### 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.

## Keywords

## Introduction

*α*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:

*λ*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,

*ζ*=

*ζ*

_{0}+

*h*(

*t*), where

*h*(

*t*) is known as the hazard function and is related to the radio-sensitivity 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:

*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 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 two-compartment 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

*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:

*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 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,

*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}}$,

*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.

_{ 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:

*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:

*a*

^{∗}and

*q*

^{∗}are:

*c*

_{1},

*c*

_{2}, and

*c*

_{3}are defined as

*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

_{ ∞ }):

_{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 fixed-point (

*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 Figure1. 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 Figure1 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, 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.

_{ ∞ }= 0,1 for different values of the division and conversion rates

*μ*and

*γ*in Figure2. 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

*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. Then, by the chain rule:

*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:

*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)).

*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 curves${\stackrel{\u0304}{x}}_{1},{\stackrel{\u0304}{x}}_{2},\cdots \phantom{\rule{0.3em}{0ex}},{\stackrel{\u0304}{x}}_{n}$. Moreover, based on the above discussion, we observe that

The values${\stackrel{\u0304}{x}}_{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 become${\stackrel{\u0304}{x}}_{i}({t}^{\ast})$. In this case, notice that the computation of the function *g*(*t*) can be *vastly* simplified if the functions *f*_{
i
}(*x*_{1}(*t*),*x*_{2}(*t*),⋯,*x*_{
n
}(*t*),*t*) do not depend explicitly on *t*. That is, if${f}_{i}({x}_{1}(t),{x}_{2}(t),\cdots \phantom{\rule{0.3em}{0ex}},{x}_{n}(t),t)={\widehat{f}}_{i}({x}_{1},{x}_{2},\cdots \phantom{\rule{0.3em}{0ex}},{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 that${\stackrel{\u0304}{x}}_{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

^{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:

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}=\sum _{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 *d* *T* = - ln(*p*_{1})/*T*_{
L
}. (5) Update time variable by adding time step computed *t* = *t* + *d* *T*. (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.

*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 Figure3. 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.

*f*. As discussed in Section ‘Stochastic two-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 one-compartment and two-compartment models. As shown in Figure4, this is in fact the case and the TCP(

*t*;

*f*) plots are almost indistinguishable.

*t*;

*f*) to reveal itself. One example is plotted in Figure5, 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 Figure1 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 have analytically derived all regimes of the phase diagram of TCP = 0,1 in steady-state, 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 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.

## Declarations

### Acknowledgements

This work was financially supported by the NSERC/CIHR Collaborative Health Research Grant (to MK and SS).

## Authors’ Affiliations

## References

- Sinclair W: The shape of radiation survival curves of mammalian cells cultured in vitro. Biophys Aspects Radiat Qual. 1966, 58 (1): 21-43.Google Scholar
- Munro T, Gilbert C: The relation between tumour lethal doses and the radiosensitivity of tumour cells. Br J Radiol. 1961, 34 (400): 246-251. 10.1259/0007-1285-34-400-246.View ArticlePubMedGoogle Scholar
- Källman P, Ågren A, Brahme A: Tumour and normal tissue responses to fractionated non-uniform dose delivery. Int J Radiat Biol. 1992, 62 (2): 249-262. 10.1080/09553009214552071.View ArticlePubMedGoogle Scholar
- Horas JA, Olguín OR, Rizzotto MG: Examining the validity of poissonian models against the birth and death tcp model for various radiotherapy fractionation schemes. Int J Radiat Biol. 2010, 86 (8): 711-717. 10.3109/09553001003734618.View ArticlePubMedGoogle Scholar
- Hanin LG: Iterated birth and death process as a model of radiation cell survival. Math Biosci. 2013, 169 (1): 89-107.View ArticleGoogle Scholar
- Zaider M, Minerbo G: Tumour control probability: a formulation applicable to any temporal protocol of dose delivery. Phys Med Biol. 2000, 45 (2): 279-10.1088/0031-9155/45/2/303.View ArticlePubMedGoogle Scholar
- Dawson A, Hillen T: Derivation of the tumour control probability (tcp) from a cell cycle model. Comput Math Methods Med. 2006, 7 (2–3): 121-141.View ArticleGoogle Scholar
- Gong J, Dos Santos MM, Finlay C, Hillen T: Are more complicated tumour control probability models better?. Math Med Biol. 2013, 30 (1): 1-19. 10.1093/imammb/dqr023.View ArticlePubMedGoogle Scholar
- Leith JT, Quaranto L, Padfield G, Michelson S, Hercbergs A: Radiobiological studies of pc-3 and du-145 human prostate cancer cells: X-ray sensitivity in vitro and hypoxic fractions of xenografted tumors in vivo. Int J Radiat Oncol Biol Phys. 1993, 25 (2): 283-287. 10.1016/0360-3016(93)90350-5.View ArticlePubMedGoogle Scholar
- Maler A, Lutscher F: Cell-cycle times and the tumour control probability. Math Med Biol. 2010, 27 (4): 313-342. 10.1093/imammb/dqp024.View ArticlePubMedGoogle Scholar
- Hillen T, De Vries G, Gong J, Finlay C: From cell population models to tumor control probability: including cell cycle effects. Acta Oncologica. 2010, 49 (8): 1315-1323. 10.3109/02841861003631487.View ArticlePubMedGoogle Scholar
- Swanson KR, True LD, Lin DW, Buhler KR, Vessella R, Murray JD: A quantitative model for the dynamics of serum prostate-specific antigen as a marker for cancerous growth: an explanation for a medical anomaly. Am J Pathol. 2001, 158 (6): 2195-2199. 10.1016/S0002-9440(10)64691-3.PubMed CentralView ArticlePubMedGoogle Scholar
- Basse B, Baguley BC, Marshall ES, Joseph WR, van Brunt B, Wake G, Wall DJ: A mathematical model for analysis of the cell cycle in cell lines derived from human tumors. J Math Biol. 2003, 47 (4): 295-312. 10.1007/s00285-003-0203-0.View ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the CreativeCommons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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.