|Home | About | Journals | Submit | Contact Us | Français|
Mathematical models of radiation carcinogenesis are important for understanding mechanisms and for interpreting or extrapolating risk. There are two classes of such models: (1) long-term formalisms that track pre-malignant cell numbers throughout an entire lifetime but treat initial radiation dose–response simplistically and (2) short-term formalisms that provide a detailed initial dose–response even for complicated radiation protocols, but address its modulation during the subsequent cancer latency period only indirectly. We argue that integrating short- and long-term models is needed. As an example of this novel approach, we integrate a stochastic short-term initiation/inactivation/repopulation model with a deterministic two-stage long-term model. Within this new formalism, the following assumptions are implemented: radiation initiates, promotes, or kills pre-malignant cells; a pre-malignant cell generates a clone, which, if it survives, quickly reaches a size limitation; the clone subsequently grows more slowly and can eventually generate a malignant cell; the carcinogenic potential of pre-malignant cells decreases with age.
Biologically motivated mathematical modeling of background and ionizing radiation-induced carcinogenesis has a history spanning more than 50 years (Nordling 1953). Many of the models can be characterized as short-term, in that they focus on processes occurring during and shortly after irradiation (Hahnfeldt and Hlatky 1998; Radivoyevitch et al. 2001; Mebust et al. 2002; Schollnberger et al. 2002; Sachs and Brenner 2005; Hofmann et al. 2006; Shuryak et al. 2006; Little 2007; Sachs et al. 2007; Schneider and Walsh 2008). The main advantage of such models is that they provide a detailed initial dose–response relation for short-term endpoints, which are used as surrogates for carcinogenesis. The main disadvantage is that the possibly substantial modulations of the magnitude and shape of this initial dose–response during the lengthy period (multiple years-decades) between irradiation and manifestation of typical solid tumors are not considered mechanistically.
In contrast, another class of biologically motivated models can be characterized as long-term, in the sense that they track carcinogenesis rates throughout the entire life span, e.g., the Armitage-Doll model (Armitage and Doll 1954; Armitage 1985), the Moolgavkar-Venzon-Knudson two-stage clonal expansion (TSCE) model (Moolgavkar 1978, 1980; Moolgavkar and Knudson 1981), the two-stage logistic model (Sachs et al. 2005), and many others (Yakovlev and Polig 1996; Pierce and Mendelsohn 1999; Wheldon et al. 2000; Little and Wright 2003; Pierce and Vaeth 2003; Ritter et al. 2003; Little and Li 2007). The main advantage of long-term models is the more detailed treatment of slow carcinogenesis processes, including the modulation of the radiation dose–response during the long latency period. The main disadvantage is that the initial dose–response is typically treated in a non-mechanistic, phenomenological manner.
The lack of detailed treatment of radiation-specific effects typically limits risk predictions from long-term models to exposure conditions where a simple dose–response relationship holds. Exposures where this relationship is more complex, such as high fractionated radiotherapeutic doses that can lead to treatment-induced cancers in nearby organs (Travis et al. 1996, 1997, 2002, 2003; van Leeuwen et al. 2003; Travis et al. 2005) are difficult to describe with current long-term models. Conversely, the more detailed dose–responses produced by sophisticated short-term models can be scaled to cancer risks only by considering the effects of factors such as background risks, age at exposure, and time since exposure, which are not explicitly taken into account by the short-term formalisms. A new, unified approach of integrating short- and long-term techniques is needed, where a detailed initial dose–response for pre-malignant cell numbers is produced over a wide range of doses, and changes to the shape of this dose–response during the latency period before the development of cancer are also analyzed in detail. A schematic representation of model unification is provided in Fig. 1.
Here, we provide a specific example of integrating short- and long-term radiation carcinogenesis models. In that a variety of short- and long-term models have been proposed, our approach should be considered as illustrative, intended to investigate the practicality of integrating the two model types. The goal was to produce a novel formalism, which can describe typical patterns of background and radiogenic carcinogenesis with the smallest possible number of adjustable parameters. This goal required multiple simplifying assumptions about the complex multi-step carcinogenesis process. We suggest that for the purposes of generating a preliminary integrated model, the consequent reduction in biological realism is compensated by the reduction in model complexity and number of parameters. More complicated examples of unifying long- and short-term models are certainly possible, e.g., models analyzing multiple pre-malignant cell stages, analyzing genomic instability, and/or analyzing stochastic effects at both time scales, rather than just at short-time scales as we shall do.
The short-term part of our model (Fig. 1a) is the more mathematically intensive part of the integrated formalism. It is based on existing initiation, inactivation and repopulation (iir) models originally designed for analyzing second cancers induced by radiotherapy (Lindsay et al. 2001; Sachs and Brenner 2005; Shuryak et al. 2006; Little 2007; Sachs et al. 2007). This short-term part analyzes normal and pre-malignant stem cells during complex radiation exposure regimens and during the following weeks of tissue recovery. It tracks individual pre-malignant cell clones rather than just the total number of pre-malignant cells. The probability that a pre-malignant clone becomes extinct during radiotherapy can be substantial, so a fully stochastic formalism is used for the population dynamics of pre-malignant cells.
The long-term part of our model (Fig. 1b) approximates carcinogenesis by a process where normal target (e.g., stem) cells can be initiated by spontaneous or radiation-induced mutations to become pre-malignant cells, which can clonally expand and perhaps produce by mutation a fully malignant cell, which then gives rise to cancer after some lag period. The basic assumptions are somewhat similar to those of several long-term models cited above, including the TSCE model. Specifically, our long-term formalism uses the following main assumptions:
The integration of long- and short-term models will have two features typical of multi-timescale modeling (Engquist and Runborg 2005): (1) information is passed in both directions between the two components; and (2) a formally infinite time interval in the short-term model represents a short time interval, here typically several months, in the long-term model (Fig. 1).
The deterministic long-term equations provide the mean number of niches filled with pre-malignant stem cells and the mean number of pre-malignant cells per niche just before irradiation. The stochastic short-term equations then provide the number of these niches that are eradicated by radiation as well as the number of pre-malignant clones that are induced by, and survive, the radiation exposure. Each of these clones is assumed to be independent (i.e., located in a different part of the organ), and to be capable of quickly filling a niche with pre-malignant cells. Therefore, the total number of pre-malignant niches soon after irradiation can be calculated by the stochastic short-term model. The mean of this number is the initial condition for the deterministic long-term equations, which are applied from this point onwards until old age (Fig. 1). Cancer incidence is assumed to be proportional to the total number of pre-malignant cells in all niches, shifted by a lag time. A mathematical description of the model is provided below.
The mathematical techniques for implementing the assumptions for our specific example of integrating long- and short-term modeling are discussed next. A schematic representation of the concepts is provided in Fig. 2. The notation and interpretations for model parameters are listed in Table 1. Data on both background and radiation-induced cancers are used to estimate the parameters.
The long-term model was intended mainly to place the results of the stochastic short-term calculations in an appropriate context, enabling estimation of the effects of age at exposure and time since exposure on predicted cancer risk. Simplicity and parsimony were emphasized. The long-term model in the absence of radiation consists just of a few rather simple deterministic equations, as follows.
Each pre-malignant stem cell in any niche has a certain small probability per unit time (q units = cells−1 × time−1) of transforming into a fully malignant cell, and, after a fixed lag time (L), into clinical cancer. These are common assumptions made in many long-term carcinogenesis models, e.g., in the TSCE model. The average number of new fully malignant cells per unit time (A units = time−1) is the product of q and three other variables: the number of stem cell niches filled with pre-malignant cells (M units = niches), the average number of pre-malignant stem cells per niche (ρ units = cells × niche−1), and the probability that the pre-malignant cells remain non-senescent and capable of producing cancer (P). The number of pre-malignant stem cells per niche (ρ) is assumed to be homeostatically regulated (in a qualitatively similar manner to the regulation of normal stem cells), so that it always tends towards a constant number—a carrying capacity Z (units = cells × niche−1). For convenience, the number of pre-malignant stem cells per niche can be redefined as a dimensionless normalized fraction C = ρ/Z. The proportionality constant qZ can be removed from the expression for A by defining N = qZM (units = time−1), where N is proportional to the number of pre-malignant niches. Consequently, A can be written in the following simplified form:
Equation (1) is an approximation for the incidence hazard function: it gives the mean expected number of new malignant cells (i.e., eventual cancers) per individual per year, whereas the hazard (H) given in (2) is the yearly probability that a malignant cell occurs in a previously healthy individual at age t:
The functions A and H are numerically very similar for realistic parameter values; they diverge substantially only if cancer incidence is high. We use the exact hazard H for data fitting in the companion paper, but use the simpler expression for A in the equations below, keeping in mind its interpretation and limitations.
The function N is described by the following differential equation (3), where t is patient age, the constant a (units = time−2) is proportional to the spontaneous stem cell initiation rate, and the constant b (units = time−1) is the pre-malignant niche replication rate:
Equation (3) can be adjusted to accommodate an assumption that more than one mutation is needed to initiate a stem cell by replacing the constant a with the composite term atκ−1, where κ is the necessary number of mutations. However, doing so does not substantially improve the fit of the model to the data sets analyzed, and so was not used as the default because it introduces an extra adjustable parameter κ. The results for fitting multi-stage extensions of our model (i.e., where κ > 1) to SEER data for female breast and male lung cancers are shown in Fig. 3.
The function C is regulated by the following logistic differential equation, where δ (units = time−1) is a rate constant representing the strength of homeostatic control of the number of pre-malignant cells per niche:
For background carcinogenesis the solution C = 1 is used.
The probability that the pre-malignant cells remain non-senescent and capable of producing cancer (P) is described by a Gaussian function with an adjustable parameter c (units = time−2), which cannot become negative even as age approaches infinity, contrary to the expression used by Pompei et al. (2001) and Pompei and Wilson (2002):
For convenience later, when radiation-induced excess relative risk (ERR) will be calculated, it is useful to express patient age (t) as the sum of age at exposure (Tx) and time since exposure (Ty). The background cancer risk function without radiation (Abac(Tx, Ty)) follows from (1) and (3–5), assuming that no pre-malignant cells are present at birth (i.e., that N(0) = 0):
Here, the duration of radiation exposure is ignored, because it is considered to be short compared with the multi-year time scale of the long-term carcinogenesis processes.
The exposure scenario analyzed here is radiotherapy for an existing malignancy, where there are K daily dose-fractions, all of equal size d in some nearby organ, with treatment gaps during the weekends. Straightforward generalizations to variable doses per fraction, to other temporal patterns, and/or to cases where one must consider dose-volume histograms are omitted for brevity. A single acute dose exposure is a simple special case, where the number of fractions is K = 1.
The short-term part of the model considers initiation, inactivation, and repopulation (iir) of normal and pre-malignant stem cells during the radiation regimen and a recovery period of about a month (Fig. 1). Niche takeover by pre-malignant cells that survive the irradiation and modulation of niche size by radiogenic promotion (discussed below) are here also considered short-term processes. They are assumed to be essentially completed by a few months following exposure, before the long-term processes of spontaneous initiation and niche replication, which operate on the time scale of multiple years, start to have an appreciable effect (Fig. 1).
Normal stem cells are treated deterministically because their number per organ is assumed to be large. In contrast, the number of pre-malignant stem cells per clone is much smaller, and extinction of some pre-malignant clones is a real possibility. Consequently, a stochastic approach described below is used to estimate the average number of live pre-malignant niches when post-irradiation long-term processes begin.
For simplicity, both normal and pre-malignant stem cells are assumed to be equally radiosensitive. The cell surviving fraction (S) is described by the standard linear-quadratic (LQ) function with parameters α and β for dose-fractions of size d:
The normal stem cell number (n) in the entire organ just before the kth dose-fraction is denoted by n−(k), and the number just after the kth dose-fraction by n+(k). The reduction in n due to initiation of a few normal stem cells is neglected. Inactivation is calculated as follows:
Repopulation of normal stem cells during radiotherapy is assumed to be regulated by homeostatic mechanisms—if some stem cells are killed by radiation, surviving stem cells are induced to proliferate (Dörr and Kummermehr 1990; Pabst et al. 2004) with the goal of restoring the total number of normal stem cells in the organ (n) to the number present before irradiation (ν). It should be noted that ν is conceptually and numerically distinct from the carrying capacity Z introduced above for pre-malignant stem cells per pre-malignant niche. The homeostatic regulation of cell proliferation during the time gaps between dose-fractions and after the last fraction is described by the Logistic differential equation, where λ (with units = time−1) is the maximum net proliferation rate:
Equations (7–9) allow the calculation of the normal stem cell number n(t) in the entire organ at all relevant times t throughout the radiation exposure and recovery periods. Explicit results for n(t) are provided in the Appendix.
We assume that the long-term growth advantage of pre-malignant cells manifests itself only on the scale of years and decades, and is negligible on the much shorter time scale of a few weeks of radiotherapy. Consequently, the maximum net proliferation rates for normal and pre-malignant stem cells are assumed to be equal, described by the parameter λ.
As discussed above, radiation initiates some normal stem cells, making them pre-malignant. The number of cells initiated in the kth dose-fraction is assumed to be Poisson distributed, with average aXI(k), where X is an adjustable parameter (units = time × dose−1) and the function I(k) (units = dose) is given by:
Our assumptions allow the birth, death, and initiation rates for normal stem cells calculated deterministically by (7–10) to be used as parameters for a stochastic formalism for pre-malignant cell clones. Each such clone is initiated as a single cell in some random niche by the kth dose-fraction, and can fluctuate in cell number during subsequent radiotherapy due to the opposing effects of inactivation and repopulation. We count the number of clones that contain at least one viable cell when radiotherapy ends, because only these surviving clones are capable of eventually taking over their niches. Their number is determined by the probability F(k) that a live stem cell initiated by the kth dose-fraction produces a clonal lineage, which survives all subsequent dose-fractions. Using analytic results on stochastic birth–death processes with variable rates (Tan 2002; Hanin 2004) the Appendix derives an equation, (22), for F(k), which is repeated here:
Here, F(0) is the probability that a pre-malignant cell that was present before irradiation began produces a lineage which survives all dose-fractions.
Using the facts that a random thinning of a Poisson distribution is Poisson and the sum of independent Poisson distributions is Poisson, one can show that the total number of surviving pre-malignant clones at the end of the short-term period is Poisson distributed. The above arguments give the mean value as Ninit = aXISf(D), where D is the total radiation dose (i.e., the sum of all doses per fraction, dK) and ISf(D) (units = dose) represents a net outcome of initiation, inactivation, and cell repopulation during exposure, The probability that a pre-malignant niche that was present before exposure is not inactivated by radiation (i.e., that at least one pre-malignant stem cell in the niche survives) is given by:
The above short-term processes are regarded as effectively instantaneous relative to the long-term ones. The stochastic results for the short-term exposure period, in the form of the functions ISf(D) and Sf(Z, D), are inserted into the deterministic equations for long-term carcinogenesis processes. On the long-term time scale, niches already pre-malignant at the end of irradiation and recovery then increase in number by replication; we will refer to all the resulting niches as “old” niches. Other, “new” niches are formed by spontaneous initiation after the end of the exposure period. The contributions of niches in these two categories to the cancer risk are called NradE(Tx, Ty) and NradN(Tx, Ty), respectively. At a given time (Ty) after irradiation they are given by the following solutions for (3):
Radiation is commonly assumed to promote hyper-proliferation of pre-malignant cells (Heidenreich et al. 2007). Here, we interpret this effect as an increase in the number of pre-malignant stem cells per surviving niche above the niche carrying capacity Z, i.e., C becomes >1. For simplicity, the initial excess of C is assumed to be linearly dependent on radiation dose with the coefficient Y (units = dose−1). Promotion may be eventually reversed, because the number of pre-malignant cells per pre-malignant niche may gradually return to pre-irradiation carrying capacity Z. This process may occur concurrently with extinction of some radiation-induced niches (Zhang et al. 2001; Ullrich 1986). Since only the product of the number of niches and the number of cells per niche is relevant for cancer risk (1), extinction of some niches and/or shrinkage of niche size do not need to be modeled separately. The net effect—i.e., a gradual reversal of promotion—can be modeled using only one adjustable parameter (δ), as done here. According to these assumptions, at any given time (Ty) after exposure to radiation the average normalized number of pre-malignant cells per surviving niche (Crad(Tx, Ty)) can be calculated by solving (4):
The approximation for the absolute cancer risk after radiation (Arad(Tx, Ty)) can now be calculated at any age at exposure and time since exposure by using the equations above:
Note that the senescence parameter (c) cancels out of this ERR expression; due to the way we have defined our radiation initiation parameter (X), the spontaneous initiation parameter (a) also cancels out. The term Q1 can be interpreted as the normalized size of old pre-malignant stem cell niches. Q2 is proportional to the number of such niches. Q3 is proportional to the number of new pre-malignant niches, and Q4 is proportional to the total number of pre-malignant niches under background conditions.
To display the properties of the model, we produced Figs. 4–7 using generic parameter values guided by the best-fit values for specific cancer sites found in the accompanying paper (Shuryak et al. 2009).
Figure 4 shows an example of a model-generated age-dependent background incidence curve for an adult-onset solid cancer, using the four relevant parameters a, b, c, and L. The shape of this curve agrees well with the data for many cancers (e.g., SEER database, http://seer.cancer.gov), including the downturn in incidence at old age, which is not described as well by standard models.
Figure 5 shows a typical radiation dose–response, which is determined by the balance of cell initiation, inactivation (killing), and repopulation. Fractionation of the dose increases cancer ERR because repopulation of both normal and pre-malignant cells during the inter-fraction intervals compensates for much of the cell killing.
Figures 6 and and77 show the effects of age at exposure and time since exposure. The components of cancer ERR produced by radiation-induced initiation and radiation-induced promotion exhibit very different dependences on age at exposure. Radiation is assumed to initiate the same number of cells per unit dose independent of age, whereas the background number of pre-malignant cells, which is essentially the denominator of ERR, grows with age. Consequently, the initiation-driven component of ERR decreases with age at exposure. This process dominates the ERR for ages <20 years in Fig. 6. In contrast, promotion is assumed to be a multiplicative amplification of the background number of pre-malignant cells per niche. Consequently, the promotion-driven component of ERR is approximately constant over most ages at exposure. This process dominates the ERR for older ages in Fig. 6.
Promotion-driven ERR can also be modulated by time after exposure. This occurs due to the assumption that the number of pre-malignant cells per niche is homeostatically regulated (by parameter δ), so that radiation-induced hyper-proliferation of cells within their niches can be reversed, as pre-irradiation cell birth/death rates, and hence niche sizes, are restored. If δ > 0, ERR due to promotion will decrease over time following exposure. This effect is seen in Figs. 6, ,7:7: If the risk is measured at some constant age (e.g., 70 years), which is a sum of age at exposure and time since exposure, a decrease in ERR with time since exposure due to a δ > 0 will appear as an increase in the ERR with age at exposure (Fig. 6). A decrease in promotion-driven ERR with time since exposure also can have a conceptually important effect on the dose–response—at longer times after exposure, not only the magnitude, but also the shape, of the dose–response can change (Fig. 7).
The formalism presented here is the first comprehensive attempt to unify short- and long-term modeling approaches. The short-term part of the model belongs to the previously discussed iir category. It tracks the numbers of pre-malignant cells throughout irradiation stochastically. The long-term part of the model builds on the concepts developed in previous two-stage formalisms by adding an analysis of some aspects of tissue architecture (i.e., stem cell niches/compartments) and aging of pre-malignant stem cells. The particular short and long-term models that we have chosen to use are not crucial—the real issue is the integration. Certainly other long-term models could perfectly well be integrated into this short-/long-term framework—and we hope they will be.
The unified formalism has a number of advantages. The short-term part can generate reasonable predictions even at high doses, such as those in cancer radiotherapy. The long-term part analyzes the entire lifetime of the individual, putting the short-term predictions in an appropriate context by estimating the effects of age at exposure and time since exposure. The combined approach therefore allows the dose–response for the number of pre-malignant cells to be examined at any time point, from the start of irradiation until development of cancer years to decades later, which is not possible using either short- or long-term models alone.
Our model can be used for estimating risks of second malignancies induced by radiotherapy. This issue is growing in importance (Brenner et al. 2000; Ron 2006) as patients are treated earlier in life and the number of cancer survivors increases (Editorial 2004); the lifetime risk of radiation-induced second cancers is not negligible (Brenner et al. 2007). Direct measurement of second cancer incidence requires decades of follow-up because the latency period for radiogenic solid tumors is long (Tokunaga et al. 1979; Brenner et al. 2000; Ivanov et al. 2004, 2009). Meanwhile, radiotherapy protocols are rapidly changing. Our model, calibrated using data from older protocols, can produce risk predictions for any modern or prospective radiotherapeutic protocol. This application of the model is discussed in the accompanying paper (Shuryak et al. 2009).
Research supported by National Cancer Institute grant 5T32-CA009529 (IS), Department Of Energy DE-FG02-03ER63668 and National Institutes of Health CA78496 (PH), National Aeronautics and Space Administration 03-OBPR-07-0059-0065(LH), National Aeronautics and Space Administration NSCOR NNJ04HJ12G/NNJ06HA28G (RKS), National Institutes of Health grants P41 EB002033-09 and P01 CA-49062 (DJB).
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
This appendix gives some details on the equations for the short-term calculation and on their derivation. We first analyze repopulation effects for normal stem cells deterministically (16–18); then we analyze survival probabilities for pre-malignant clones stochastically. The equations derived are extensions of results given previously (Sachs et al. 2007).
Denote the time of the kth dose-fraction by T(k). We will derive recursive equations valid for k = 1,…, K, where the number of fractions is K as in the main text. We formally define T(K + 1) ∞; this infinite value represents the end of the recovery period (Fig. 1). It reflects a standard procedure in multi-timescale analyses (Engquist and Runborg 2005), an infinite time interval in a short-timescale model represents a short time interval in the next larger timescale model (here several months in our long-timescale model). We correspondingly set n−(K + 1) = ν, the set point value attained at the end of the recovery period.
Solving (9) of the main text gives n(t) for any time t after the kth dose-fraction and before the (k + 1)st:
In particular, for n−(k + 1), 16 gives:
To analyze pre-malignant clones stochastically, consider a clone that starts with a single live stem cell initiated by the kth dose-fraction, and is followed in time. The time evolution of the clone is modeled as a time inhomogeneous birth–death process with birth rate b(t) and death rate r(t) (Tan 2002). The appropriate death rate, taking the spontaneous death rate as zero to minimize the number of adjustable parameters, is the death rate due to the dose-fractions, namely
where δ [t] is the Dirac delta function.
Equation (19) corresponds to the statements that on average the surviving cell fraction for the kth dose-fraction is given by (18) and that pre-malignant stem cells have the same radiosensitivity to inactivation as do normal stem cells. By our assumption that, during the comparatively short radiotherapy and recovery periods, normal and pre-malignant cells have effectively the same proliferation rate, the appropriate birth rate is
In (20) n(t) is a known function of time, determined as discussed above.
It is well known (Tan 2002, pp. 169–171) that by integrating an appropriate partial differential equation for the probability generating function one can deduce the following expression for the probability F(k) that the pre-malignant clone survives all subsequent dose-fractions:
Here, T+(k) denotes the time just after the kth fraction. Because of the way a Dirac δ function behaves, using T+ rather than T is important in the expressions for ψ and ξ.
Equation (22) is valid for k = 0, 1,…, K, with k = 0 referring to pre-malignant cells present before radiation starts and F(K) = 1. It is the primary mathematical result needed for the data analysis discussed in the main text.
We have some extensions, not needed in the present paper. Generalizing to situations where the spontaneous death rate is non-zero and/or one needs the probability that a clone has a given number of cells at the final time, can readily be done by using results given by Tan (2002). In addition, it can be shown that (22) holds even if the recovery equation for normal cells is different from the logistic equation used in this paper, e.g., is Gompertzian.
An erratum to this article is available at http://dx.doi.org/10.1007/s00411-011-0378-5.
Igor Shuryak, Email: ishuryak/at/gmail.com.
David J. Brenner, Phone: +1-212-3059930, Fax: +1-212-3053229, Email: djb3/at/columbia.edu.