Beyond Two-Stage Models for Lung Carcinogenesis in the Mayak Workers: Implications for Plutonium Risk

Mechanistic multi-stage models are used to analyze lung-cancer mortality after Plutonium exposure in the Mayak-workers cohort, with follow-up until 2008. Besides the established two-stage model with clonal expansion, models with three mutation stages as well as a model with two distinct pathways to cancer are studied. The results suggest that three-stage models offer an improved description of the data. The best-fitting models point to a mechanism where radiation increases the rate of clonal expansion. This is interpreted in terms of changes in cell-cycle control mediated by bystander signaling or repopulation following cell killing. No statistical evidence for a two-pathway model is found. To elucidate the implications of the different models for radiation risk, several exposure scenarios are studied. Models with a radiation effect at an early stage show a delayed response and a pronounced drop-off with older ages at exposure. Moreover, the dose-response relationship is strongly nonlinear for all three-stage models, revealing a marked increase above a critical dose.


Introduction
Cancer is a genetic disease. In the widely held theory of somatic evolution [1], a cell's path toward the malignant state is portrayed as a series of mutations or epigenetic events, lending it successive selective advantages. These advantages, as summarized in the "hallmarks of cancer" [2], essentially amount to an increasingly uncontrolled proliferation.
Those essential features-mutations accompanied by proliferation-have long been identified as key ingredients in modeling carcinogenesis. Beginning with the seminal multi-step models by Armitage/Doll and Nordling [3,4], this eventually led to the stochastic two-stage model with clonal expansion due to Moolgavkar, Venzon, and Knudson [5,6], which has by now become an established tool to understand and predict cancer risk [7][8][9].
What fundamentally distinguishes such mechanistic models from conventional epidemiological ones is that they do not directly model the endpoint-say, the cancer mortality ratebut rather the process thought to lead to it, parametrized via mutation and proliferation rates. This may prove useful especially when the mortality rate is highly convoluted by an exposure to carcinogens such as ionizing radiation, as is the case in this study.
Such a mechanistic approach can be readily generalized so as to build in known biological effects, such as multiple genetic pathways [7,10] or a more realistic number of stages [11,12]. Indeed, for colorectal cancer, where the understanding of the cellular mechanisms is comparatively advanced [13,14], a number of extended models have been put forward to account for the role of genomic instability [10,[15][16][17], the rather large number of premalignant stages [12,18], as well as the intricate dynamics during progression [19].
By contrast, far less is known regarding other cancer types. For lung cancer, mechanistic modeling studies are abundant but have focused almost exclusively on the two-stage model. These indicate that for the two main risk factors, smoking [20][21][22] and α-particle radiation (most relevant being Radon decay products) [23][24][25][26][27][28][29][30][31][32], the best description is afforded by an enhancement of the proliferation rate of premalignant cells. However, for radiation, this conclusion has been disputed [33] because it lacks a conclusive biological mechanism, in contrast to the accepted mutagenic effect of radiation. This debate has been further sparked by the single analysis to date going beyond the two-stage model [34]. Comparing fits to the Colorado-miners data using the two-stage model with those from a subclass of three-stage models, the authors suggested that a proliferation effect were confined to the two-stage model, whereas a better fit quality was achieved within a three-stage framework with a mutational radiation action.
The objective of this paper is to perform a comparative analysis of lung-cancer risk associated with α-radiation using different multi-stage models-specifically two-and three-stage models as well as a two-pathway model for (radiation-induced) genomic instability. Our goal is to identify the mechanisms of radiation action suggested by those models, as well as to lay out to what extent their predicted radiation risks are consistent. To this end, we apply these models to the Mayak-workers cohort [35,36]. These workers, employed at the formerly Soviet Plutonium-production plant, have been exposed to substantial doses of 239 Pu via inhalation and exhibit a large number of lung-cancer deaths, 895 in total [37]. A notable feature of Plutonium exposure is its strong protraction, which might facilitate the assessment of risk on the long time scales relevant to indoor Radon. Furthermore, information on the strongest risk factor, smoking, is available for most workers. We will show that certain three-stage models give an improved description of the data, and we elucidate how these lead to predictions for the risk that are qualitatively different from both two-stage mechanistic and standard descriptive models.

Materials and Methods
Mayak-workers cohort dominant cancer-mortality endpoint is lung cancer (defined here as ICD-9 code 162), with a total of 895 mortality cases.
Cohort definition. To obtain a sufficiently homogeneous data set amenable to mechanistic modeling, several selection criteria have been applied, similar to previous studies [31,32]. Our reduced (sub)cohort excludes females, since these make up less than 25% of the whole cohort and exhibit very different mortality rates. Moreover, full information is required on smoking/alcohol status and annual internal dose(rate)s-i.e., 239 Pu doses must be measured or assumed to vanish (for workers outside the radiochemical and Pu plants). Although missing risk-factor information may be taken into account via additional fit categories, we found that this does not reduce the parameter uncertainties, presumably due to the noise introduced this way.
Finally, the follow-up period is restricted as follows. If a Pu measurement has been performed at time t Pu , then the entry date is set to two years after the measurement date, t Pu + 2a. This is done to avoid selection bias, specifically healthy-survivor effects (due to extended follow-up periods for persons surviving until t Pu ) and diagnostic bias (in case the measurement has been caused by imminent health problems). To ensure complete follow-up, the exit date is cut off at the end of 2008 or, in the case of migrants, 2003.
The reduced cohort includes 8,604 persons and 388 lung-cancer deaths. Risk factors. Let us briefly highlight some aspects of the major risk factors (see Ref. [37] for details). The main interest here is in internal radiation, with measured nonzero doses available for 3,667 persons. Due to slow degradation in the lungs and the long half-life of 239 Pu, exposures are highly protracted: First exposure peaks around age 20, and typically continues until the end of follow-up. Cumulative lung doses are well described by a log-normal distribution. Among those measured, it is peaked at 4mGy, much below the mean dose D ¼ 0:12Gy. Restricted to lung-cancer cases, the dose distribution is shifted to higher values, with D ¼ 0:44Gy and lower/upper 5% quantiles of 7.6mGy and 2.3Gy.
The overall smoking fraction is about 3/4. Alcohol status, although not known to be a risk factor for lung cancer, may serve as an indicator for smoking habits because the fraction of smokers increases with alcohol consumption. We group heavy/chronic drinkers as one category, a = 1, and otherwise set a = 0.
As with any cancer, age is a crucial intrinsic risk factor. The ages of cohort members range broadly between about 18 and 81 years, as defined by the lower/upper 5% quantiles of entry/ exit age, with an average of 27 years spent in the cohort. By contrast, 90% of lung-cancer cases are found only between ages 49-78, with a mean cancer age of 65 years.
Ethics statement. The study of the Mayak-workers cohort has been reviewed and approved by the Southern Urals Biophysics Institute's Review Board for issues related to privacy and personal data protection. All patient records were anonymized and de-identified prior to analysis.

Statistical analysis
Different classes of multi-stage models with clonal expansion are applied to the Mayak data. Specifically, these are two-and three-stage models as well as a model with two distinct pathways. In the following, we will sketch the basic assumptions underlying these mechanistic models and some properties of their solutions. We then discuss how external risk factors, such as radiation, are included in this framework, before laying out the procedure for model fitting and selection.
Multi-stage models. The common rationale of all mechanistic models studied here is to radically reduce the complexity of carcinogenesis to essentially two key processes: (i) mutations -or, generally, a series of (epi-)genetic transitions from healthy via pre-malignant to malignant stem-like cells-and (ii) proliferation(i.e., symmetric division; cell inactivation or death) of pre-malignant cells with a selective advantage [40]. Note that this simplified single-cell picture does not explicitly include cellular interactions. In particular, non-stem cells are not taken into account.
Mathematically, this is modeled as continuous-time Markov processes for the (stochastic) numbers of cells, X i (t), at the different stages (i = 0, . . ., k). Specifically, at age t = 0, one starts with X 0 N healthy cells, which can make a transition (modeled as a Poisson process) with rate μ 0 to a first, premalignant stage, X 1 . These cells can then undergo a birth/death process with rates α 1 /β 1 , leading to a net proliferation rate γ 1 % α 1 − β 1 . Further transitions eventually lead to malignant cells, X k . The occurrence of the first malignant cell is assumed to lead to cancer after an effective lag time t lag , typically on the order of a few years. A cartoon depiction is shown in Fig 1; for k = 2, this corresponds to the standard two-stage model with clonal expansion (TSCE).
The mathematical model above can be solved for the survival function S(t) and, equivalently, the hazard (here: lung-cancer mortality rate), h ¼ À d dt lnS, using the method of characteristics [41]. For the two-stage model, assuming rate parameters to be piecewise age-independent, an exact closed-form solution can be attained [42]. In the general case, an extension of this solution is valid approximately if age bins are small enough for intermediate-cell numbers to change slowly.
As an illustration, let us highlight some generic features shared by all such multistage models. At earlier ages, the hazard is well described by a deterministic model, h ' m kÀ1 X kÀ1 , in terms of the mean numbers of cells, _ X j ¼ m jÀ1 X jÀ1 þ g j X j [41]. This leads to an initially polynomial growth, h(t) ' + Nμ 0 Á Á Áμ k − 1 t k − 1 /(k − 1)!, followed by a rapid proliferation-driven phase, h(t) = O(e γt ), where γ denotes the maximum growth rate. However, this deterministic approximation fails to account for the effective reduction in available premalignant cells as new malignancies arise: Whenever a person reaches the cancer endpoint, those cells can no longer lead to further cancer cases. At older ages, approximately around the mean cancer age, a steady state is reached between growth and effective "loss" of premalignant cells [41], and the hazard levels off to a constant limit, h 1 * Nμ 0 γ 1 /α 1 .
These borderline cases also illustrate a more general point: Not all biological rates can be determined from fits to the cancer data alone. Generically, the hazard only depends on certain parameter combinations [43]. The combinations chosen for the fits here are shown in Tables 3  and 4. For example, in the two-stage case, Nμ 0 μ 1 may be interpreted as an overall scale factor for the hazard; γ α 1 − β 1 − μ 1 is essentially the net proliferation rate and yields the system's (inverse) time scale, whereas α 1 μ 1 effectively reduces the mean cancer age, at which saturation of the hazard sets in.
In the discussion so far, we have tacitly assumed a linear series of transitions for simplicity. However, we also test a model where mutations may occur along two different pathways, as sketched in Fig 1(bottom). Such a model has been introduced by Little et al. [10] to account for genomic instability, inspired by models for colon cancer [14]. The underlying idea is that a second path, activated via transition rates σ j , corresponds to the loss of a gene involved in maintaining genomic integrity. This may lead to mutation rates much larger than for the genomically stable (upper) path, m GI j ) m j [44]. Despite the more complex topology, the model equations are constructed and solved using exactly the same principles as outlined above.
Risk modeling. In the framework of multi-stage models, the hazard is fully determined by the (generally time-dependent) mutation and growth rates. These parameters are essentially assumed to have an age-independent background value which may be modified by external risk factors such as radiation dose rate, d(t), and smoking, s(t), here assumed to start from 18 years on. In practice, both risk factors are allowed to independently increase any of the rate parameters ϑ l 2 {μ 0 , μ 1 , . . ., γ 1 , . . .} in a suitable parametrization, e.g. ϑ l (d); from all possible combinations {ϑ l (d), ϑ l 0(s)}, the best-fitting model is selected. Depending on which of these rates show a radiation effect, the radiation risk varies in a characteristic fashion with age [42], as well as with modifiers such as duration of, or age at, exposure. This is in marked contrast to conventional descriptive models employed in radiation epidemiology [45], which we also use to benchmark our results. In a descriptive model, the hazard function is modeled directly-rather than its underlying mechanism. Here we use the conventional parametrization [37] hðtÞ where the baseline, h bsl ðtÞ ¼ e cðtÞþc rf , is factored into terms for age-dependent background, cðtÞ P 2 j¼0 c j ln j ð t 60a Þ, and other risk factors, ψ rf (birth year, smoking and alcohol, etc.). The excess relative risk, ERR, is factored into dose-response shape-typically a function of cumulative dose, D(t − t lag )-and time-dependent modifiers such as attained age or age at exposure; see Appendix. Generally, to test specific covariables, we first use categorical fits so as to explore the qualitative dependence on this covariable, before constructing adequate analytic parametrizations.
Fit procedure. All model parameters ϑ are estimated by maximizing their likelihood, L(ϑ) = ∏ i ℓ i (ϑ), constructed using the individual likelihoods of all cohort members [46]. These are the probabilities for survival throughout each member's follow-up period, multiplied-for cancer cases-by the probability for cancer occurring during the final year. Equivalently to maximizing L, we minimize the deviance, D = −2 ln L > 0, using the Minuit function-minimization library [47]. For model selection, we rely on the likelihood-ratio test for nested models so as to retain only parameters significant at a 95% confidence level. The same level is also adopted throughout this paper for confidence intervals. To rank non-nested models, the entropy-based Akaike index is used [48], AIC = D + 2n, which effectively penalizes overfitting for models with larger number of parameters n. Without loss of information, we will denote only the difference ΔAIC with respect to the benchmark descriptive model. Table 1 provides an overview of the best-fitting models. Before explaining in detail their mechanisms and the implications for radiation risk, let us anticipate some general patterns. All highest-ranked multi-stage models share a Plutonium-induced enhancement of proliferation rates. More specifically, the fits suggest that 3-stage models with a radiation effect on an early stage of proliferation (models A, B) may yield an improved description of the Mayak data, compared with an effect on the penultimate stage (C). Although model A exhibits by far the lowest AIC, we will present the radiation risks for several three-stage models so as to give an impression of the range of model predictions, as discussed in the Appendix.

Results
No evidence is found for a radiation-induced second pathway. We note that all dose rates here refer to internal Pu exposure; external radiation is not significant in the two-stage and descriptive models and thus not considered in the following.

Two-stage model
The two-stage model has been applied to previous follow-ups of the Mayak workers [31,32]. We shall therefore discuss it here as a benchmark, but also to illustrate some mechanisms inherent in any multistage model. We find that the effects of both Pu dose rate as well as smoking status are by far best described as additively enhancing the net proliferation rate, γ(s, d) = γ (0) (s) + δγ(d). (An additional radiation effect on the initiating rate, μ 0 (d), is highly insignificant when fitted to the data and thus not included in our model. This may be simply due to lack of data power and does not imply absence of radiation-induced initiating events.) A categorical fit of the dose dependence strongly suggests a saturation at larger dose rates (Fig 2), and we find it best modeled by an exponentially leveling function Here r * 5/Gy (see Table 3) governs the linear, low-dose response, δγ * r d, and γ 1 * 0.3/a denotes the rate approached as d ) dÃ γ 1 /r; here dÃ * 0.06Gy/a. This is qualitatively in line with previous analyses [31,32] but also several Radon-risk studies [23,25,27,29,30], see Discussion. It is worth noting that the data fit equally well to a response in terms of the accumulated dose D-i.e., δγ(D) in Eq (2)-which may relate to the long protracted exposures of Mayak workers. It should be stressed that, even though the main risk factors, smoking and Plutonium, enter the growth rate additively, the actual risk will exhibit an interaction between them. To illustrate this, in Fig 3(a) we display the age-dependent hazard, h(t), for a representative exposure scenario at a constant dose rate from age t 1 = 25a to 60a, with dose D = 0.2Gy. For a wide age range, coinciding with the phase of exponential proliferation, the relative risks of radiation and smoking approximately multiply (note the log scale). It is only at larger ages (≳ 60a) that the combined risk drops below the multiplicative value. Such sub-multiplicity agrees with trends glimpsed in a descriptive analysis [37]. Here, it follows naturally because the hazard of Pu-exposed smokers levels off much earlier, reflecting an earlier onset of cancer (see Methods). To better link these findings to those of descriptive models (see Appendix), we will from now on consider the excess relative risk, ERR h/h bsl − 1, defined relative to the zero-exposure baseline risk. From this angle, multiplicity of risk implies an equal ERR for smokers and nonsmokers-i.e., a ratio ERR (s = 1) /ERR (s = 0) = 1. Fig 3(b) shows that, under the scenario above, the risk is indeed multiplicative until older ages (t≲ 60a). It then becomes sub-multiplicative, and the ERR ratio drops markedly below unity after (time-lagged) exposure ends.
Since most cases are related to smoking, we will now concentrate on the risk for smokers. Moreover, to separate age and dose dependencies, we scale the ERR by the accumulated dose, ERR(D;t)/D(t − t lag ), with the lag time t lag = 5a. Fig 4 displays the age-dependent ERR/D for the scenario above but at D = 0.5Gy. For the two-stage model, it reveals a characteristic increase during exposure (owing to proliferation), followed by a marked drop-off after exposure has ended. The latter is similar to the (non-significant) attained-age trend found in the descriptive model, cf. Appendix.
The age dependence just discussed pertains to one specific scenario. Let us now indicate how the risk is modified by different exposure patterns. The dependence upon age at first exposure, t 1 , is shown in Fig 5(a). We have chosen a scenario with D = 0.2Gy and a typical duration Δt = 20a; the ERR is recorded at t 1 + Δt + t lag . Clearly, for this two-stage model, the variation is rather mild for all but very early exposures (not encountered at Mayak) and very late ones. In the former case, virtually no premalignant cells are available for proliferation. At older ages, in turn, they are increasingly lost to new malignancies; thus the risk is attenuated. In the descriptive model, a weak (non-significant) trend also suggests a slight decrease of ERR with older ages at exposure (not shown). Fig 5(b) reveals a characteristic influence of exposure duration, Δt (where D = 0.2Gy, t 1 = 20a). For very short exposures, Δt ( D/dÃ * 3a, the ERR is strongly suppressed: The short duration cannot be compensated by an increased dose rate, since the growth rate saturates at dose rates larger than dÃ. This inverse dose-rate effect is thus inherently connected to a saturating radiation response as in Eq (2). It has been observed also in mechanistic Radon studies [23,25,30], although the Mayak data are not powerful enough to support this at the descriptive level. At sufficiently long exposures, Δt ) D/dÃ, the inverse dose-rate effect disappears and eventually gives way to a slight direct effect. This is related to the leveling of the hazard upon reaching malignancy.
To wrap up this discussion, let us consider the dose-response relationship implied by this proliferation-based model. Fig 6 illustrates that, in contrast to the linear dose response typically assumed in descriptive modeling (ERR(D) = cD, here c % 4.7/Gy), this model exhibits a nonlinear response. Although the quantitative values depend on the exposure scenario (here: exposure during t = 25 − 60a), some general features apply to any two-stage model with a proliferation enhancement similar to Eq (2): The dose response is characterized by a linear low-dose regime, ERR ' r × D for D ( 1/r, and leveling at sufficiently large dose(rate)s and/or ages. Intermittently, typically an exponential increase is seen, reflecting the exponential growth of premalignant cells due to proliferation. However, this effect tends to be washed out by the leveling at older ages. In a descriptive model fit, where essentially an averaging occurs over all exposure histories in the cohort, this pronounced nonlinearity will be even harder to resolve. Still, it is noteworthy that a TSCE-inspired descriptive model, ERR = e f(D) − 1 with f ðDÞ c 1 ð1 À e Àc lin D=c 1 Þ, yields an improved fit (but similar AIC), with a linear response c lin * 3.2/Gy and leveling at DÃ c 1 /c lin * 0.8Gy.

Three-stage models
As mentioned earlier, the highest-ranked 3-stage models fall into two categories.
The best two models (A, B; see Table 4) show a radiation effect on the earlier stage of proliferation, γ 1 (d), differing only in their background parameters-specifically the smoking Beyond Two-Stage Models for Lung Carcinogenesis response on μ 1 (A) or γ 2 (B). This early impact leads to a substantially delayed radiation response compared to both the 2-stage model and 3-stage model C, as is seen from Fig 4: The ERR/D is initially zero and, once it has peaked, tends to drop off more mildly. (In some scenarios, the ERR/D may even increase when exposure stops.) This lag is intuitive, as the insulted cells need to pass through an additional mutation stage before becoming malignant. In stark contrast, model C-with an effect on the penultimate stage-predicts a much higher risk, ERR/ D ' r * 10/Gy, right after exposure. In the scenario depicted in Fig 4, it further displays an almost monotonic drop with attained age, not unlike the trend seen in the descriptive model risk.
However, it must be stressed that this specific scenario conceals an underlying complexity not present in the 2-stage model. The reason is that the 3-stage model is determined by two competing growth rates, γ 1 and γ 2 . Let us consider models A/B: For a small enough dose rate such that g ð0Þ 1 þ rd < g ð0Þ 2 , the clonal dynamics is governed by the largest baseline rate, g ð0Þ 2 . In other words, below the critical dose rate d crit ðg ð0Þ 2 À g ð0Þ 1 Þ=r $ 0:01Gy=a, both baseline and excess risk grow with the same exponential rate-hence the ERR would level off with age. It is only above that critical dose rate that an exponential increase is seen in the ERR. Likewise, for model C, the critical point, ðg ð0Þ 1 À g ð0Þ 2 Þ=r, marks the dividing line between exponential increase and leveling with age. Notice that the value in Fig 4 is just at the borderline.
The dose-rate dependence just described is reflected in the dose response (Fig 6). At low doses, a linear increase is seen, just as for the two-stage model. However, that low-dose response is typically much lower than what we saw for two-stage and descriptive models. It is only at the critical dose, here d crit Δt * 0.35Gy, that a rapid exponential increase sets in. Compared to the two-stage case, this exponential increase is much sharper owing to the higher response coefficients, r * (10 − 17)/Gy. As before, the response levels off once the corresponding dose rates exceed dÃ * (0.03 − 0.05)Gy/a.
Another consequence of an early-stage radiation effect is a notably suppressed risk for older ages at exposure, t 1 [Fig 5(a)]. To understand this, recall that the radiation effect consists in multiplying the available pool of stage-1 cells, X 1 (t 1 ), prior to their making the transition to Beyond Two-Stage Models for Lung Carcinogenesis stage 2. Since that number grows only very slowly at a rate g ð0Þ 1 < g ð0Þ 2 , the head start given to already existing stage-2 cells, X 2 (t 1 ), becomes overwhelming the later in life exposure starts. Thus the ERR is suppressed exponentially as e Àðg ð0Þ 2 Àg ð0Þ 1 Þt 1 for large t 1 . By contrast, the mechanism for model C is fairly similar to that of the 2-stage model: Both radiation and baseline risk are governed by the growth of existing stage-2 cells. This is why virtually no dependence on t 1 is seen other than a mild drop-off for older ages at exposure, due to the onset of malignancies.
Differences with respect to the two-stage model may also be observed in the exposure-duration dependence [Fig 5(b)]. Like the 2-stage models, all 3-stage models exhibit an inverse doserate effect for durations Δt ≲ D/dÃ. As explained previously, this merely hinges on the saturation of the proliferation rate for higher dose rates. It may be noted that the early-stage models A/B lead to a stronger risk suppression at such short durations, due to the extra mutational step to be passed. However, a more qualitative difference is found for longer durations: Here a marked direct effect occurs, the ERR falling off as Δt −1 . This may be viewed as a linear doserate modification for perturbatively small dose rates D/Δt in the limit of long durations.
Finally, let us comment on the interaction between radiation and smoking risks. We saw earlier that the two-stage model predicts a largely multiplicative interaction or, equivalently, a near-unit ratio of smokers' and non-smokers' ERR [Fig 3(b)]. The situation is less clear cut for the three-stage models. Model C, with a later-stage radiation effect, is most comparable to the 2-stage case: Initially, it also exhibits near-multiplicity, which reflects that radiation simply leads to multiplication of existing stage-2 cells. However, the ERR ratio then drops very rapidly. This is essentially because for smokers, the dose rate is below the critical value, ðg ð0Þ 1 À g ð0Þ 2 Þ=r $ 0:01Gy=a, as smoking leads to a very large growth rate γ 1 (s = 1) % 0.16a −1 . Thus the ERR levels off, in contrast to the exponentially growing ERR for non-smokers, where there is no threshold. An analogous mechanism is at work in the model B, which is submultiplicative throughout.
By contrast, for model A, the risk in the scenario shown in Fig 3(b) is strikingly super-multiplicative except for ages t ≳ 70a. At first glance, this deviation from multiplicity might seem surprising: For constant rates, the hazard is typically proportional to μ 1 (s), and the smoking dependence should thus drop out of the relative risk. However, smoking is assumed to start at age 18, only a few years prior to irradiation. Hence the baseline risk-initially proportional to the number of existing stage-2 cells-mostly stems from those cells created before smoking started, and is thus comparable for non-/smokers. Thus, the large ERR for smokers reflects the much higher excess risk due to freshly mutated stage-1 cells. It is only long after the start of exposure that the smokers' baseline risk increases sufficiently to compensate for this.

Two-pathway models
We have tested a family of multi-path models of genomic instability (GI). Let us briefly outline the key assumptions made here to reduce the number of parameters-in total, 6 mutation rates and 6 cell-division/death rates, as sketched in Fig 1(bottom). Most premises are motivated by the biological mechanisms thought to underlie GI [14,44]. First, the destabilizing transition rates are set equal, σ j σ, as they pertain to the inactivation of the same genomic-integrity gene. Since, presumably, GI per se does not lead to any growth advantage, we set the birth/ death rates at stage 0 GI equal, a GI 0 ¼ b GI 0 , their precise value being marginal. Mutation rates following GI, m GI j , are supposed to be larger, or at least equal, compared to their genomically stable (upper-branch) counterparts.
Although GI may well be present as a sporadic mechanism, the baseline data are not thought to provide sufficient structure to distinguish complex background models (see Discussion; cf. also Ref. [15]). We have thus focused on the case of radiation-induced GI-i.e., an activation of the otherwise silent GI path (σ) by radiation. This may be modeled as σ = r GI × d (with vanishing background rate). In addition, radiation may affect regular mutation or growth rates.
None of the tested two-pathway models has led to any significant, numerically stable improvement beyond the benchmark two-stage model. It is stressed that relaxing the assumptions of linearity or vanishing background rate do not yield a qualitative improvement.

Mechanism of radiation action
A robust result of our analysis is an α-radiation-induced enhancement of proliferation rates of premalignant cells. This corroborates a consistent finding in many studies based on fitting two-stage models to lung-tumor data, both from other epidemiological cohorts (Radonexposed miners [23,25,27,29,30]) and animal experiments (see Ref. [49] for an overview). It also alleviates concerns that a proliferation effect might be confined to the two-stage model [34].
However, it has long been criticized [33] that no radiobiological evidence exists for such a premalignant-growth-enhancing effect. Although experimental evidence is still sparse [50,51], let us briefly discuss some theoretical models put forward to explain how radiation might lead to enhanced proliferation of premalignant cells.
The most elaborate model is based on the following idea [52]: Radiation kills cells, which in turn triggers division of neighboring stem cells so as to replenish the lost tissue. This may lead to net proliferation of premalignant cells if these have a selective advantage, that is, a slightly higher rate for cell division than needed for homeostasis. (Strictly speaking, cells need not necessarily be killed; it might be sufficient for them to have a proliferative disadvantage-e.g., by effected cell-cycle arrest-such that these are subsequently repelled by premalignant cells.) That hypothesis has been shown [53][54][55] to lead to a dose-rate response, γ(d), which is qualitatively similar to that found in cohort studies (see, e.g. , Fig 2), albeit with quantitatively modest agreement. Saturation of the growth rate much higher than a characteristic dose rate, d ) d s , occurs if more cells are killed than can be substituted for by premalignant cells.
To obtain a very rough estimate for that critical dose rate, note that α-particle hits of a cell nucleus are independent and rare. Thus the number of hits, N, is Poisson-distributed, P N ¼ e À N N N =N!, which means the fraction of cells not hit is P N¼0 ¼ e À N . Assuming (i) a linear dose response, N ¼ n D, with n * 4/Gy [56], and (ii) delivery of the dose D d τ over a characteristic time of order the interval between cell cycles (τ * 1a for basal stem cells [57]), we have P N = 0 (d) = e −ndτ . At the characteristic dose rate, d s , about one out of, say, six nearest neighbors would be hit, such that P 0 (d s ) = 1 − p, p 1/6, yielding a characteristic dose rate of This is on the same order of magnitude as the value found in this study, dÃ * (0.03 − 0.06)Gy/a. In this light, the model results presented here and the repopulation mechanism may be interpreted to be compatible. However, it is important bearing in mind that these estimates are naturally crude. Matters are further complicated by the spatially inhomogeneous energy deposition within different spots in the lungs, an effect not reflected in whole-lung doses used here [58].
As an alternative mechanism, a radiation-induced disturbance of cell communication has been suggested [9,59]. This may lead to, e.g., up-regulated growth signals or a reduction of apoptosis [60,61], with a higher effect on intermediate cells because those tend to evade homeostatic control. It has even been proposed that a proliferation enhancement, mediated by such a bystander signaling, might be the generic mechanism for the response to densely ionizing radiation [62]. However, no mechanistic model has been put forward explaining in detail how this might lead to a dose-rate response, γ(d).
Even so, should the radiation response indeed be governed by the bystander effect, then a similar behavior ought to be expected as for bystander-mediated mutation induction [63]. In microbeam experiments [64,65], it has been found that for low doses-corresponding to less than * 10% of cells being hit-the mutagenic yield was strongly amplified as bystander cells also received signals. (A similar pattern has been found for intercellular induction of apoptosis [66].) For much higher doses, in turn, the response essentially saturated. Along the lines leading to Eq (3), we can estimate the characteristic dose rate d s by assuming the crossover to occur at a fraction p = 0.1 of cells being hit. This yields d s % p/nτ * 0.025Gy/a, under the same caveats as mentioned above. From this standpoint, both bystander signaling and the repopulation hypothesis appear compatible with the dose-rate response found in this study.

Comparison with previous studies
As discussed earlier on, for mechanistic models, the risk is essentially determined by its structure, particularly, the radiation response. A dose-rate dependent proliferation rate, γ(d), saturating for larger dose rates (as in Eq 2) is found in many studies applying the two-stage model to α-particle-induced lung cancer. In particular, our response quantitatively agrees with that of the preferred two-stage model by Jacob et al. for the Mayak cohort, both for Plutonium and smoking [32]. Concordantly, their risk estimates are similar to those of the present TSCE model-such as a cohort-averaged excess risk of ERR(t = 60a)/D * 4/Gy, a nonlinear dose dependence ERR(D) for larger doses, and sub-multiplicity of smoking and radiation risks.
In a recent descriptive analysis of the Mayak data, Gilbert et al. found a linear dose response, modified significantly only by a drop-off with attained age [37], see also Appendix. The value at age 60, ERR(t = 60a) * 7D/Gy, is somewhat higher than for the cohort average of the mechanistic models presented here. By contrast, the multi-stage models exhibit a strongly nonlinear dose dependence especially for higher doses. Furthermore, they typically display a decrease with attained age only for large enough ages, most pronounced after the end of exposure. Initially, an increase with age is seen due to exponential clonal growth, at least for high enough doses.
In contrast to descriptive models, where an exposure modifier may not be significant when parametrized explicitly, mechanistic models implicitly make predictions for the risk dependence on any exposure scenario. This is exemplified by the age-at-exposure dependence or the inverse dose-rate effect shown by several mechanistic models (Fig 5). Furthermore, a non-significant trend in Ref. [37] indicated a sub-multiplicative interaction between Plutonium dose and smoking. This is in agreement with results from the 2-stage model, which further offers a mechanistic interpretation in terms of exponential cell growth, combined with earlier malignancies for smokers (see Results).

Implications for lung carcinogenesis
We have shown that several three-stage models give an improved description of the data compared to one involving two stages. Evidently, this stochastic inference is based solely on the lung-cancer endpoint and cannot replace experimental insight into the dynamics of intermediate stages. Even so, it is worth emphasizing that the Mayak data do not fully allow to distinguish general mechanisms for carcinogenesis. Rather, the results here mostly rely on the radiation-associated risk. In fact, the deviances of the various mechanistic and descriptive baseline models (risk factors age and smoking) do not differ noticeably-in contrast to the models including radiation (Table 1).
That statement seemingly contradicts the fact that there is only a fraction * 25% of excess cases relative to the baseline (as can be seen from ERR $ 5 D=Gy, the whole-cohort average dose being D $ 0:05Gy). This translates to 60-70 excess cases (Table 2), compared to about 320 baseline cases. However, truly spontaneous baseline cases (30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40) are outnumbered by smoking-related ones by a factor of * 10. Moreover, the smoking variable is only binary (and noisy). This makes it hard to resolve the actual baseline risk accurately, especially since the multi-stage models do not differ in their qualitative behavior except for young ages. By contrast, the models do differ markedly for different irradiation scenarios as found in the Mayak data.

Conclusion
We have shown that carcinogenesis models extended to three stages offer an improved description of the Mayak lung-cancer data, as compared to the two-stage model. All favored 3-and 2-stage models indicate a radiation-enhanced proliferation rate of premalignant cells, suggesting that this is a robust finding-at least in the framework of multi-stage models with clonal expansion-and not limited to the case of two stages. Despite that structural similarity, the models make qualitatively different predictions for the risk following certain exposure scenarios. For instance, those models whose radiation impact is on an earlier stage exhibit a strongly suppressed risk for older ages at exposure. As opposed to the two-stage case, all three-stage models reveal a critical dose(rate) above which the excess risk increases sharply. Moreover, while an inverse dose-rate effect is predicted by all models, only those with three stages also display a pronounced direct effect for longer exposure durations.
One aim of this study has been to elucidate which aspects of carcinogenesis models are persistent themes or rather model-specific features. Such a better understanding should facilitate the development of improved carcinogenesis models, both mechanistic but also descriptive ones. There is still some way to go toward a more realistic description. Ultimately, this would involve biological input on rates or premalignant stages so as to cut the number of undetermined parameters. Closer at hand, a natural next step may be validating the current models on other data sets. An extension to more than three stages conceivably provides a further improvement. Moreover, to get a more accurate description of the biological mechanisms, it is desirable to develop models specifically for the different histological cancer subtypes. Table 2. Observed numbers of lung-cancer cases by dose categories, compared with those predicted by the descriptive and multi-stage models (in brackets: excess cases). As a reference, we also give the person years (py), i.e., the number of years spent in the respective (sub)cohort summed over all persons.

Appendix: Parameter estimates
We now present some details on the model results. The reference descriptive model closely follows Ref. [37]. The main risk factor is smoking, which is included simply as a factor in Eq (1), e c smk ¼ e 2:3AE0:3 % 10. The alcohol status further increases the baseline risk by e c alc ¼ e 0:6AE0:1 % 1:8. In addition, the birth year was found to elevate the risk between years 1915 and 1935 by about e 0.4±0.1 % 1.5. With 3 extra parameters, this birth-cohort effect is moderately significant (ΔD = −11). It might be interpreted as a smoking modifier (as which it yields only a slightly higher deviance), possibly related to the changed smoking levels between the wars. The effect is irrelevant for the radiation risk, and generally no birth-cohort distinction is made in the scenarios in Figs 3-6. We mention that a calendar-year dependence has been tested but discarded: It oscillates strongly but shows no conclusive trend, nor does it influence the dose response.
The dose dependence is modeled as a linear function, ERR = cD, c = (4.7 ± 0.9)Gy −1 , with no evidence for threshold or quadratic terms. Note that the lag time entering the dose, D D(t −t lag ), does not alter the deviance significantly between 0-10a. We fix it at t lag = 5a, also for all multi-stage models. No significant time-dependent dose modifiers are found. However, a trend suggests a decrease of ERR/D with attained age, modeled as (t/60a) −2.4±2.5 , and a marginal decrease with median age at exposure. Even though the reference model (Eq 1) implies multiplicative risks of smoking and radiation, a non-significant trend indicates sub-multiplicity, with an ERR for smokers reduced by a factor e −0.9(±1.5) * 0.4.
The maximum-likelihood parameters of the TSCE model are shown in Table 3. We have omitted two covariables not relevant for the radiation risk: alcohol and birth year. Both are included as addends to the growth rate γ, consistent with an interpretation as smoking surrogates. We have also tested an explicit age dependence of the rate parameters, which indicated a reduction of μ 1 (t) or, equivalently, γ(t) above age * 50a. However, these effects were numerically unstable and have thus been discarded. To convey an impression of the fit quality, Table 2 shows a juxtaposition of observed cancer cases and those predicted by the various models across the dose range.
A comment is in order on the parameter estimates of the three-stage models ( Table 4). Some of the (structurally) identifiable parameter combinations may be practically indeterminable insofar as they leave the minimum deviance D virtually unchanged. In model A, for instance, D is independent of μ 1 μ 2 ! 0, which has been fixed at an arbitrarily small value. Furthermore, the best estimate of g ð0Þ 1 ¼ a 1 À b 1 % À0:1a À1 is not significant, and g ð0Þ 1 is set to zero. Worse yet, for model B (C), the background rate g ð0Þ 2 (g ð0Þ 1 ) is unstable, tending to Beyond Two-Stage Models for Lung Carcinogenesis unrealistically large negative values. Since only the smokers' value is stable, g ð0Þ 2 þ Dg 2 ðs ¼ 1Þ $ 0:1a À1 , we have set g ð0Þ 2 0, even at the price of a higher deviance. It is also for these intricacies that we have opted to present several of the best three-stage models, rather than relying strictly on the weight suggested by Akaike's index. This way, a more plausible impression is given of the uncertainty of model predictions. Table 4. Same as Table 3, but for the parameters of the best three-stage models.