- Open Access
Stochastic model for tumor control probability: effects of cell cycle and (a)symmetric proliferation
© Dhawan et al.; licensee BioMed Central Ltd. 2014
- Received: 29 August 2014
- Accepted: 21 October 2014
- Published: 22 November 2014
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.
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.
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.
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.
- Tumor control probability
- Cell cycle
- Mathematical modeling
- Stochastic birth-death process
- Method of characteristics
- Gillespie algorithm
with D(t) being the dose in Gy delivered until time t and its time derivative representing dose rate (Gy/day).
Paper expanded this approach to include the effect of cell cycle sensitivity in the TCP. They considered a two-compartment model for the active (M,G1,S, and G2 phases) and the quiescent (G0 phase) cells (see also). The radio-sensitivity of resting cells and active cells are significantly different; the radio-sensitivity is typically much higher for actively proliferating cells. This model was discussed both deterministically and stochastically in, but the stochastic master equation is solved under the assumption that the joint probability distribution function of two populations,, 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. 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 G0 (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; 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.
Taking na,0 and nq,0 to be the initial numbers of active and quiescent cells, respectively, we have the initial condition and the boundary condition V(1,1,t) = 1, where the boundary condition comes from the definition of the generating function.
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 into the two-compartment model of. 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.
This gives f(s) = f0(x1(0),x2(0),⋯,x n (0)).
The values are determined by the set of ODEs in (17), with the final value condition that x1(t∗) = x2(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. In this case, notice that the computation of the function g(t) can be vastly simplified if the functions f i (x1(t),x2(t),⋯,x n (t),t) do not depend explicitly on t. That is, if, 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. 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.
The process to calculate TSi(t) is: (1) Compute likelihood of each cellular reaction occurring (L i for reaction R i ). (2) Sum together all likelihoods into quantity. (3) Compute uniformly distributed random numbers p1 and p2 in the interval (0,1). (4) Compute the next time step of a likelihood reaction, assuming exponentially distributed times d T = - ln(p1)/T L . (5) Update time variable by adding time step computed t = t + d T. (6) Determine which reaction to carry out: if Li-1/T L ≤ p1 ≤ 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.
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 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 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 using similar parameter values. In, 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.
This work was financially supported by the NSERC/CIHR Collaborative Health Research Grant (to MK and SS).
- 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
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.