Repopulation of interacting tumor cells during fractionated radiotherapy: Stochastic modeling of the tumor control probability

Purpose: Optimal treatment planning for fractionated external beam radiation therapy requires inputs from radiobiology based on recent thinking about the “ﬁve Rs” (repopulation, radiosensitivity, reoxygenation, redistribution, and repair). The need is especially acute for the newer, often individualized, protocols made feasible by progress in image guided radiation therapy and dose conformity. Current stochastic tumor control probability (TCP) models incorporating tumor repopulation effects consider “stem-like cancer cells” (SLCC) to be independent, but the authors here propose that SLCC-SLCC interactions may be signiﬁcant. The authors present a new stochastic TCP model for repopulating SLCC interacting within microenvironmental niches. Our approach is meant mainly for comparing similar protocols. It aims at practical generalizations of previous mathematical models. Methods: The authors consider protocols with complete sublethal damage repair between fractions. The authors use customized open-source software and recent mathematical approaches from stochastic process theory for calculating the time-dependent SLCC number and thereby estimating SLCC eradication probabilities. As speciﬁc numerical examples, the authors consider predicted TCP results for a 2 Gy per fraction, 60 Gy protocol compared to 64 Gy protocols involving early or late boosts in a limited volume to some fractions. Results: In sample calculations with linear quadratic parameters α = 0.3 per Gy, α / β = 10 Gy, boosting is predicted to raise TCP from a dismal 14.5% observed in some older protocols for advanced NSCLC to above 70%. This prediction is robust as regards: (a) the assumed values of parameters other than α and (b) the choice of models for intraniche SLCC-SLCC interactions. However, α = 0.03 per Gy leads to a prediction of almost no improvement when boosting. Conclusions: The predicted efﬁcacy of moderate boosts depends sensitively on α . Presumably, the larger values of α are the ones appropriate for individualized treatment protocols, with the smaller values relevant only to protocols for a heterogeneous patient population. On that assumption, boosting is predicted to be highly effective. Front boosting, apart from practical advantages and a possible advan-tage as regards iatrogenic second cancers, also probably gives a slightly higher TCP than back boosting. If the total number of SLCC at the start of treatment can be measured even roughly, it will provide a highly sensitive way of discriminating between various models and parameter choices. Updated mathematical methods for calculating repopulation allow credible generalizations of earlier results.


1.A. Current treatment planning
Modern technological developments for fractionated external beam radiotherapy include IMRT (intensity modulated radiation therapy), IGRT (image guided radiation therapy), and SBRT (stereotactic body radiation therapy). 1 Such techniques allow planning and efficient delivery of very conformal treatments as well as escalated doses to appropriate volumes, and therefore offer new versatility in treatment strategies. Nonhomogeneous dose distributions and nonuniform dose fractionations are becoming more common. Due to the wide spectrum of emerging protocols and to the substantial timelag between treatment and the availability of data on outcomes, radiotherapy treatment planning needs ongoing new inputs from radiobiology. Inputs on tumor control probability (TCP) are based on modern versions (reviewed in Ref. 2) of the "five Rs" of radiobiology, including repopulation.

1.B. TCP and stochastic tumor repopulation during radiotherapy
It has been known since the seminal research of Tucker and co-workers 3-5 that tumor repopulation during treatment must often be considered in quantitative estimates of TCP, may have to be estimated stochastically, and should then take into account non-Poisson statistics of malignant cell numbers. Stochastic estimates can take eradication of cells into account more systematically than can deterministic ones. Prospective estimates of repopulation are especially important for novel protocols. Reviews of the mathematical, computational, and translational aspects of stochastic and deterministic TCP estimates taking into account repopulation include. [6][7][8] Repopulation models indirectly incorporate some aspects of reoxygenation (reviewed in Ref. 9).
Stochastic estimates involve a mathematical/ computational approach to repopulation that uses "stochastic processes" (also known as "random processes"); roughly speaking a stochastic process is a time-dependent variable that is subject to random fluctuations; 10 in stochastic TCP analyses it can be specified by a time-dependent probability distribution function for malignant cell number. A stochastic-exponential process applicable to proliferating, noninteracting cells (the Feller-Arley birth-death process with time dependent, state-independent per-cell rates, as characterized, e.g., in Ref. 11, Example 4.1) has been studied for more than 70 years. 12 It has been applied, as an extension of the previously standard deterministic-exponential equation, to TCP estimates. 3,6 It leads to a front-loading theorem: 13 other things being equal, giving dose earlier increases the TCP.

1.C. Interactions of stem-like cancer cells (SLCC)
For TCP estimates the relevant repopulating tumor cells in a heterogeneous tumor are those that can individually regrow the entire tumor population when in the appropriate microenvironment. Such cells are sometimes called stem-like cancer cells; 14 we use SLCC for short in what follows. However, in this context no additional "stem" properties are implied 15 and various other terms have been used (e.g., "clonogens" 8,16 ).
Previous stochastic TCP models have assumed independent, i.e., noninteracting, SLCC. However, cell-cell interactions are ubiquitous in vivo, 17 and there is increasing evidence that tumors include microenvironmental niches within which important SLCC-SLCC interactions occur. 18 A niche can be an actual biological structure, such as a crypt in the case of colon cancer; 19 or be a locale for radio-resistant clones; 20 or it may merely approximate low cell mobility with dominant cell-cell interactions localized.

1.D. Preview
This paper develops a new radiobiologically based formalism for analyzing repopulation of interacting SLCC. The formalism emphasizes subpopulations of SLCC, each within its own niche, 21  niches are visualized as comparatively distant from each other. Generalizing the stochastic-exponential model described above and its deterministic counterpart, all kinds of LSCC-LSCC interactions are allowed, though only among SLCC within the same niche (Fig. 1).
Interactions within one niche are comparatively easy to model because the number of SLCC per niche is comparatively small. Interactions among SLCC in different niches are assumed negligible.
The following ten comments preview motivations, characteristic aspects, and limitations of our approach.   Birth rate when a niche contains j SLCC Calculated from m, σ , and C Preliminary comments above, d(j) "Endogenous" death rate a when a niche contains j SLCC a Applies between fractions, and also after the last fraction before recovery ends.
6. The number of SLCC within one niche is not assumed to follow a Poisson distribution, but the number of SLCC-occupied niches per patient just before treatment starts is. A corresponding mathematical technique is called "poissonization" (reviewed in Ref. 26) and has become common in probability theory because it simplifies many calculations and/or adds to the realism. The Poisson assumption on initial niches per patient is at least as reasonable as an alternative assumption, often made (explicitly or implicitly) that, after stratifying for age, gender, nationality, and cancer stage, all grouped patients have exactly the same number of SLCC before treatment starts. 7. Niche time evolution obeys "IID" (independent, identically distributed) assumptions. For example, it is assumed that the probability distribution function for the number of SLCC per niche just before radiation starts is the same for all niches. The IID assumption does not imply that all the SLCC niches initially have the same number of SLCC, just the same probabilities. Also, two niches that happen to have the same number of SLCC at some given time need not always continue to have the same number, but the SLCC in both niches have the same birth and death rates at that particular time. The IID assumption for niches is more general than the corresponding assumption usually made for individual SLCC because SLCC within a niche need not be independent. 8. When giving illustrative numerical results, we will, as is often done, take TCP as the probability that all SLCC have been eradicated (due to either radiation or stochastic fluctuations) at a specified time; gener-alizations are mentioned in item 10. For our stochastic niche models we will take the time as 23 days (a weekend plus 3 weeks) after the end of therapy. Here the specific number 23 days is chosen for three main reasons: (a) it is short enough that long-term processes 8,27 such as microscopic tumor dormancy 28,29 have not yet had time to set in; (b) it is long enough that for our stochastic sigmoid niche calculations it gives the same numerical TCP (to within less than 0.1%) as does artificially extending the calculations to periods too long to be within their domain of validity; and (c) it corresponds roughly to an actual recovery period observed clinically. 9. Our approach is meant mainly for comparing similar but somewhat different protocols. 10. It uses a number of intentional oversimplifications and has a number of limitations, though not, we would argue, more than similar approaches. Examples of limitations include the following: (a) the TCP criterion just described above glosses over various possibilities, e.g., complete remission without SLCC eradication, and it circumvents modeling of long-term processes; (b) we do not here attempt to give novel mathematical results on the other four Rs, on early or late deterministic normal tissue complications, or on iatrogenic second cancers. These aspects will be briefly discussed, but without quantifications. Table I The 60 Gy protocol is one that was used 30,31 in NSCLC treatments that gave a (dismal) TCP of 14.5%. The LQ parameters are consensus values for an individual patient (e.g., Ref. 32). Smaller values of α are used in calculations intended to model the shape of the TCP vs dose curve near TCP 50 for a heterogeneous patient population 33 so two of our auxiliary examples will use α = 0.03 per Gy instead. The default value of the (one-cell) Malthusian rate m corresponds to a commonly adopted minimum doubling time of 3 days (e.g., Ref. 34). The binomial distribution for J(t) before treatment starts, the value J = 10, and the value C = 20 of the carrying capacity are merely illustrative; several auxiliary models will be used to estimate the sensitivity of the conclusions to varying these values. Finally, the biological motivation for assuming a substantial default endogenous death rate (σ = 0.5) is the observed high cell loss in pulmonary tumors (e.g., Refs. 35 and 36).

1.E. Explanations of and motivations for the default parameter values in
The parameter σ is a quantity with no real analogue in deterministic models. It can be interpreted as follows: 37 if net endogenous birth rate b-d is fixed, σ = 1 corresponds to a highly stochastic situation such as neutral drift for SLCC population number between fractions, while σ = 0 usually gives results much closer to a corresponding deterministic model with b-d as growth rate. Four of our auxiliary examples will explore the effects of changing σ compared to its default value 0.5.

2.A. Fractionated radiotherapy
We consider throughout a multiweek protocol of K dose fractions, administered at times T 1 = 0, T 2 , . . . , T K so far apart that complete sublethal damage repair occurs between fractions. The respective doses, which can be different for different fractions, are denoted by d k and the total dose is denoted by D = K k=1 d k . For example, Fig. 2 shows a protocol with front boosting of dose.
During therapy, SLCC are killed by radiation but also repopulate. On average the SLCC number per niche decreases markedly, as indicated. Most SLCC niches will be cleared of SLCC during treatment, indicated here by omitting such "eradicated" niches. If all niches are eradicated, the tumor is controlled. Otherwise, remaining SLCC niches may die out accidentally after therapy during a recovery period. Dur- ing recovery there could also be an increase in the number of SLCC per niche, as shown. If some SLCC remain, the outcome during the following months and years (indicated as "long term") could be recurrence, but could also be different, e.g., formation of tumors that are dormant, microscopic, and thus asymptomatic. 28,29 Long-term processes are not analyzed in this paper. We define S(k) as the LQ cell survival probability for the kth fraction: If there are j SLCC present just before the kth fraction the probability of having i SLCC present just after is thus We can and shall regard this expression as a square matrix with rows labeled by i and columns labeled by j. It is a "stochastic" and "triangular" matrix (i.e., each row sum is 1 and the (i,j) entry is 0 if i > j).

2.B. Customized computer programs
Subsection 2.C describes calculations for our main example. Those for the 14 auxiliary examples were similar. Computer calculations used the R language. Section S1 of the supplementary material gives open source R scripts which can be used to reproduce the main calculations. 44

2.C. Time evolution of SLCC
For repopulation between fractions or after treatment stops, the assumptions stated and motivated in the Introduction imply the following for the main example, i.e., for the stochastic-logistic process with the default parameter values in Table I.
(b) The probability distribution of SLCC number in one niche just before treatment starts is (c) The specific numerical values are C = 20, m = ln(2)/3 per day, and σ = 0.5.
Using Eqs. (3) and (4) we integrated 20 ordinary differential equations (ODE) numerically 20 times to construct a 21 × 21 repopulation matrix, which in the jth row gives the probability distribution for ending one interfraction day with i SLCC (i = 0, 1, . . . , 20) given that the day started with j SLCC. We then used this repopulation matrix, the survival matrix given by Eq. (2), and a day-by-day iteration to find the probability P[J(62) = 0] that a niche has no SLCC at the end of the recovery period (Fig. 2). This matrix approach was spot-checked with Monte Carlo simulations.
These calculations and poissonization allowed calculation of (a) the average initial number of SLCC per patient (n 0 in Table I) from data on the default protocol and (b) the TCP for boosted protocols, given this calculated n 0 .

3.A. Results on examples
Recall the following: (a) our "default" protocol is 30 weekday doses of 2 Gy each and by adjusting parameters is assigned the 14.5% TCP found in some NSCLC studies; 30,31 (b) our "illustrative TCP criterion" is that there are no SLCC at day 62 (Fig. 2); and (c) our "main" example of niche dynamics during periods of no radiation is a stochasticlogistic process with the default parameters in Table I. In our main example, the basic Poisson color theorem, 10 that a random thinning of a Poisson process is a Poisson process, implied that under our assumptions the number of nonempty malignant niches N(t) per patient remains Poisson distributed at all times, in contrast to the number of SLCC per niche. The mean of N(t) is the initial mean N times the probability that J(t) = 0. Using our criterion that TCP is given by the Poisson zero class thus gave where q is the probability that a niche has at least 1 SLCC at time T K+1 . Equation (5) applied to the default protocol gave N ∼ 5.2 × 10 5 . Using this N in 64 Gy front or back loaded protocols gave row 1 in Table II

3.B. More detailed results on time dependence
In all 13 stochastic examples, the probability distribution of SLCC per niche differed strongly from a binomial or Poisson distribution even after <10 dose fractions-for example, the surviving fraction was much smaller than for a Poisson distribution of the same mean.
It was found that the niche model can give fully detailed time-dependent results. As one specific example, consider the following procedure for estimating with any of the stochastic examples the probability that, say, a patient has exactly M niches that contain exactly j SLCC at a particular time <62 days. By the Poisson color theorem this probability is Here p is the probability that a niche has exactly j SLCC at the given time, calculated from the niche dynamics as described in Sec. 2. It was also found that the special case where each niche starts out with exactly one cell and the niche-dynamics is stochastic-exponential gives stochastic-exponential dynam- ics, corresponding to noninteracting SLCC, for the SLCC population as whole. Section S3 of the supplementary materials analyzes some previously known TCP results 3, 6, 13 obtained assuming such dynamics, thereby supplying some additional mathematical and intuitive explanations of the results shown in Table II.

4.A. Review
We argued that biologically based prospective TCP estimates, though difficult, are needed. A mathematical niche model was given for estimating TCP. The calculations emphasized: (a)

4.B. Effectiveness of moderate boosting
For α = 0.3 per Gy, Table II shows that modest boosts are predicted to be very effective, with an extra 4 Gy improving TCP from 14.5% to >70% irrespective of the particular model or the other parameters. We suggest that the contrary result predicted by assuming α = 0.03 per Gy (rows 12 and 13 in Table II) is not applicable to individually designed protocols for two reasons. First, these small values of α are estimated not from the average of heterogeneous patient α values but as a formal consequence of the effect of heterogeneity in flattening the observed TCP 50 slope. 33 Second, the small values of initial SLCC number predicted, 20 for row 12 and 34 for row 13, are not plausible. For example, if, as has been suggested, 38 >1% of all NSCLC cells are SLCC, then the initial number of SLCC must be >1 000 000, rather than 34. Table II   Table II shows that if n 0 can be estimated from separate, biological data it becomes a highly sensitive way to select preferred models and/or appropriate parameters; for example, when the assumed value of α is decreased by one log as above, the resulting estimate of initial SLCC number decreases by more than five logs. In the present context, n 0 in Table II serves as an indicator of how effectively SLCC are eradicated: higher effectiveness means more initial SLCC must be present to lead to 14.5% survival for the default protocol. Proliferation is predicted to have a marked effect on loss of control, as exemplified, e.g., by the greater n 0 in row 11 vs row 1. Thus unplanned interruptions are, as usual, deprecated.

4.C. Additional comments on
Front boosting is slightly favored over back boosting in most models. A different consideration that favors front boosting is that often treatment is preceded by a positron emission tomography scan and the location of the gross tumor volume is more accurately known at the start of the treatment than during later fractions.

4.D. Parsimony
Equations (1), (3) and (4), and (5) imply that our main model uses the following adjustable parameters: α, β, m, σ , C, and n 0 = JN. For comparison, the only previously used stochastic model, the stochastic-exponential (Feller-Arley) model, uses all of these except C. The deterministicexponential equation for the time evolution of the average SLCC number, which has long been the mainstay of TCP estimates, uses all but C and σ . Thus this paper introduces one new parameter, C, to take intercellular interactions into account, and, as in previous stochastic calculations, uses σ . In some cases, such as colon cancer, 39,40 niches are identifiable and independent estimates for C are then available; σ is observable, e.g., by observing net endogenous birth rate and cell death. In Table II the TCP estimates change by ∼5% or less when σ or C are varied, in contrast to the more sensitive dependence of TCP on α and n 0 .
The insensitivity of our conclusion about boost effectiveness to many parameter alterations (Table II, rows 3-6 and 8-11) means that for this purpose the niche model is no less parsimonious than the earlier models that assume no SLCC-SLCC interactions. In one way it may actually be more parsimonious. The calculations that showed almost all surviving niches in the sigmoid models have been filled up to the carrying capacity by the end of a 23 day recovery period indicate that perhaps long-term primarily the TCP and the carrying capacity are relevant, not additional details on the SLCC per niche probabilities.

4.E. Some further limitations of the present approach
Results on TCP must always be considered in the context of other relevant factors, such as the probability of late normal tissue complications and the influence of other "Rs." In this paper this broader context was not quantified. For example, redistribution, radiosensitivity, and especially reoxygenation might be modified by radiation and therefore affect TCP dependence on front vs back boosting. Such factors can depend in a complex manner on dose timing, tumor type, and stage. Scenarios increasing TCP and scenarios decreasing TCP can both be envisaged.
In the related problem of radiogenic second cancers, 27 increasing attention is now being paid to the long-term tumor dynamics. 21,41 For TCP the corresponding question concerns events during the years after treatment (Fig. 2), which determine whether surviving SLCC, if any, lead to recurrence or remain controlled at sub-clinical levels, and also determine the timing of any recurrences that do recur. 8 The assumption made here (and in most other TCP analyses) that every SLCC must be eradicated by the end of the recovery period for control to occur is an oversimplification; 8,21,41 this particular assumption tends to underestimate TCP. We did not attempt postrecovery period analyses here because finding parsimonious, well-motivated models for the long-time behavior is, in our opinion, a difficult, separate task. 27 For example, some long-term processes increase niche number 21 and then the basic Poisson color theorem cannot be used in corresponding long-term calculations.
Another complication is that, apart from deterministic late normal tissue effects, early boosts might modify the probability of radiation-induced, independent secondary cancers, where malignant or premalignant cell repopulation during treatment can be important. 27,42 Stochastic calculations suggest that early boosts are less dangerous in this respect than late boosts. 21,43 The intuitive reason is that new malignant clones initiated by early boosts can be eradicated by later fractions, whereas those initiated by late boosts are more likely to survive the entire treatment.
Finally, we here considered only individualized protocols, not protocols designed for a cohort, where interpatient variations in basic radiobiological parameters are known to be important to the TCP. For the latter the predictions of drastically improved TCP for moderate boosts are overoptimistic, as discussed above.

CONCLUSIONS
The application of our results to full treatment planning protocols, which attempt to consider all relevant major effects, thus remains problematical. But we would argue: r Prospective, biologically motivated calculations that produce numerical TCP estimates are important, especially for designing model-assisted, physician controlled, individualized protocols. r The niche model can give such estimates using a more general, at least equally practical, and at least equally credible approach compared to older models.