Missense, nonsense, and frameshift mutations, and dysregulated alternative splicing in GCSFR have been isolated in patients with MDS/AML. In the study of Beekman et al. (2012), nonsense and missense mutations in GCSFR arose during the course of the disease. In the model we envision, population genetics, and population dynamics of proliferating bone marrow cells are closely intertwined.

Proliferating healthy cells in the bone marrow mutate at random times, possibly influenced by super-pharmacological doses of GCSF. A summary of possible mutations and their consequences for proliferation dynamics of granulocyte precursors is depicted in Figure 1. GCSF signaling occurs through its cognate receptor, GCSFR. It involves both proximal signaling networks consisting of signaling molecules such as Lyn, Jak, STAT, Akt, and ERK, and distal gene regulatory networks consisting of transcription factors. Together, these signaling networks promote proliferation, survival, and differentiation. In patients with SCN, the earliest known mutation to contribute to transformation to secondary MDS or AML is a nonsense mutation in the GCSFR gene. This mutation leads to a truncated receptor, GCSR delta 715 (Glaubach and Corey, 2012, and reference therein).

It follows from a simple calculus of mutation events that as long as the cell population size is kept in check, the rate at which new mutant clones appear in the population is rather low. When the population expands, new mutant clones arise faster (see further on).

In our model we take the view that carcinogenesis is driven by a succession of small-scale (e.g., point) mutations in specific loci. Other viewpoints (epigenetic effects, karyotypic alterations, intercellular interactions, etc.) have been suggested. In treatment-related MDS some drugs (e.g., many alkylating agents) induce t-MDS primarily via large scale alterations that lead to karyotypic instability (Bhatia, 2011).

Limited mutation load at the SCN phase causes neutropenia and fluctuations of cell population size. With time, accumulation of driver mutations causes expansion of mutant clones, which however are not yet expanding at a dramatic rate. At some point in time, mutations accumulate sufficiently to cause a major change in the proliferation law and the now malignant cell population starts rapidly expanding.

Our model is based on the following hypotheses (Figure 2):

At the time of diagnosis of SCN, GCSF therapy is initiated, which induces an initial series of X driver mutations, occurring at random times.

In the model, the increasing proliferation rate of successive mutant clones causes acceleration of growth of the malignant bone marrow stem cell population, which shortens the time interval to appearance of new clones, which in turns increases proliferation rate, and so forth; this results in a positive feedback. As we will see, the stochastic nature of the process (the times to appearance of each next mutant are random) causes a spread of the timing of the subsequence mutations, particularly the first X mutations during the SCN phase. This may result in the transition to MDS not manifesting itself for a very long time in a fraction of cases.

We explain some other intuitions underlying the model. For a new subclone, stochastic theory is used to estimate extinction probability, with extinction after more than a few cell generations being negligible in view of the growth advantage of the new clone. However the time at which the next mutation occurs in a cell clone is also stochastic and it is as a rule more dispersed for the slower-growing clones. Therefore the time to reaching the threshold number of bone marrow stem cells (which in our model defines the time at sAML diagnosis), is a random variable. One of the questions we ask is if dispersion of this time matches the wide distribution of the times at diagnosis (Rosenberg et al., 2010).

The population-genetic effect of population size-dependent accumulation of mutations occurs as a natural consequence of the proliferation law in the form of a multitype Galton–Watson branching process:

Consecutively arising surviving mutant clones are numbered with the index k, ranging from 1 to K; time interval between the appearance of the k-th and k + 1-st surviving mutant clones is denoted by τ_k. k-th mutant cells have accumulated k driver mutations (assuming the clone in SCN bone marrow at diagnosis has a single cell with one driver mutation, which seems a defendable idealization).

These three rules allow one derive the probability distributions of time intervals τ_k, probabilities of survival of each clone, and expected growth laws of each clone. Mathematical details follow from the theory of Galton–Watson branching process; see for example the monograph by Kimmel and Axelrod (2002).

We assume that cell division is effective with probability b, i.e., the probability generating function (pgf) of the number of progeny cells per parent cell has the form f(s) = bs² + (1 − b). The extinction probability q is the smaller solution of the equation q = f(q), which is less than 1 if the process is supercritical. In our case, (1)q=bq2+1-b⇒q=b-1-1;b∈(0.5,1].

Similarly, the expected number of progeny of a cell is equal to f ′(1−) = 2b, hence the expected number of cells at time t is equal to N(t) = (2b)^(t/T), which yields the value of λ (2)expλt=2bt∕T⇒λ=ln2b∕T.

We will use “continuous” time t for notational convenience. However, we consider generations of cells dividing at discrete times t_i = iT, where T is the average cell cycle time. As it is known, the expected (mean) growth law in the Galton–Watson process has the form (3)E# cells, at time t, in a clone started at time t0=defNt=expλt-t0, as t→∞.

To determine the distribution of time to a mutation creating a new non-extinct clone, we consider a newborn cell. In this cell, mutation may occur with probability u, and if the extinction probability of the mutant clone is q′, then the probability that the cell does not produce a new mutant clone is equal to 1 − u(1 − q′). Until time t_i = iT, approximately ∑j=0iNtj=Nt0expλTi+1-1expλT-1-1 new cells are born, and assuming independence, we obtain (4)Pr[no mutant initiating nonextinct clone appears untiltime ti=iT]=Pr[τ=:time to nearest nonextinctmutant clone > ti=iT]=1-u1-q′∑j=0iNtj=1-u1-q′Nt0expλTi+1-1expλT-1=adexpct-1. where, for the k-th mutant population a=1-1-qk+1uk=1-uk2bk+1-1∕bk+1,c=λk=ln2bk∕T,d=expλkT-1-1=2bk-1-1.

Since the distribution tail of random variable τ_k has the form Prτk>τ=adexpct-1=exp-lna-1dexpcτ-1. it can be algorithmically generated using the inverse tail method (5)τk=c-1lnlnr∕dlna+1, where r is a pseudo-random number uniformly distributed from 0 to 1. In this framework, a sample path of the number of cells in the k-th mutant clone (which contains cells with k mutations accumulated) is equal to (6)Nkt=0t≤∑j=1k-1τjexpλkt-∑j=1k-1τjotherwise

The derivations presented are quite similar to those of Bozic et al. (2010), except that in that paper, expected times E(τ_k) to the next mutation have been used. Here, we are interested in exposing stochastic variability in the time course of the SNC → sMDS → sAML transition. Another refinement would be to use distributions of cell counts instead of expected values N_k(t). This would result in serious computational problems, arguably without much impact on the results.

Equations 1 and 2 allow generating realizations of times to successive driver mutations under different values of mutation rates and proliferative characteristics of the mutant clones. We make the following assumptions:

Transition to sMDS requires one or two somatic driver mutations, whereas the transition to sAML requires at least three somatic driver mutations (cf. Table 1).

We carried out computations based on two literature sources and then used rounding to the nearest order of magnitude to obtain a working threshold number of the leukemic initiating cells (LIC) (Bonnet and Dick, 1997). In both cases we assume that the volume of human bone marrow is equal to V = 1700 ml and that LIC cells constitute a fraction ψ = 10⁻⁶ of leukemic bone marrow mono-nucleated cells (BMMNC). We also assume that in sAML, fraction ρ = 0.8 of BMNNC is constituted by leukemic cells.

Under given values of parameters κ and ε as well as mutation rates u_k, the time at diagnosis of sAML, defined as the time T from initiation of GCSF treatment such that ∑kNkT=L depends on parameter A according to an approximate power law T=fA=expαAβ, where β < 0. This dependence, which was obtained via simulation studies (not shown), allows finding the distribution of A that leads to a clinically observed distribution of the time of sAML diagnosis according to the following expression for distribution tails F̄Aa=1-F̄Tfa, where F̄Tt=PrT>t is the tail of the distribution of time T. This in turn allows generating pseudo-random realizations of A according to the expression (8)A=f-1F̄T-1R=exp-αF̄T-1R1∕β, where R is a pseudo-random number from the uniform distribution on the (0, 1) interval.

We need to approximate the tail of the distribution of the time at diagnosis of sAML. A recent source is the paper by Rosenberg et al. (2010). These authors reported results of a prospective study of 374 SCN patients, and included estimates of hazard rates and cumulative probability of sMDS/sAML as a function of time after GCSF treatment. Hazard rate grows for the first 5 years and then plateaus. To simplify computations we adopted a piecewise constant estimate of the hazard rate h_T(t) hTt=0.01t∈[0,3)0.02t∈[3,∞) with time in years. Comparing with Figure 1A in Rosenberg et al. (2010) we see that h_T(t) remains within the confidence band computed based on the prospective study. Using the expression F̄Tt=exp-∫0thTτdτ and inverting the tail function F̄Tt we complete the derivation of expression Eq. 8 (elementary details not shown).

The form of expression Eq. 8 and plausible estimates of parameters κ and ε as well as of mutation rates u_k, are difficult to be uniquely determined with the data available at the present time. We used the following heuristic procedure:

Driver mutation rates increase from the reference value by a factor of 5, starting mutation 3, so that u_1,2 = u but u_3,4,5 = 5u. The increase is needed for the later mutations to occur in quick succession, so that mutation 3 occurs before ∑kNkT>L, with L = 10⁴ being a relatively low value.

Figure 3 depicts the impact of successive driver mutations on the natural course of the SCN → sMDS → sAML transition. Figure 3A depicts counts N_i(t) of cells in successive mutant clones as a function of time, under model as in Eq. 7 with A = 0.005, ε = 0.02, and κ = 2. Straight lines with increasing slopes are counts of cells in successive mutant clones. We observe that the time intervals separating the origins of successive clones are decreasing with each mutation event. Thick dashed line represents the total mutant cell count. It is also interesting to observe that clones with increasing numbers of mutations dominate transiently, until they are replaced by other clones with higher proliferative capacity (selective value). Figure 3B depicts relative proportions n_i(t) = N_i(t)/Σ_jN_j(t) of cells belonging to successive mutant clones.

It is somewhat surprising that under any combination of coefficients A and k, the range of simulated times at sAML diagnosis is rather narrow. Figure 4B depicts ranked simulated times at sAML diagnosis under model as in Eq. 7 with A = 0.005, ε = 0.02, and κ = 2. Spread of these values is narrow, with interquartile range between 15 and 21. Systematic simulation experiments demonstrate that this is the case for a wide range of A and κ parameter values. This outcome is in contrast to the wide spread of times at diagnosis summarized in Table 1 and that based on Rosenberg et al. (2010).

Simulation-estimation experiment outlined in the Methods demonstrates that distribution of simulated times (counting form initiation of CGSF treatment) at sAML (Rosenberg et al., 2010) is reproduced by our model. Figure 4A cumulative distribution of the times at sAML diagnosis under model as in Eq. 7 with κ = 2, ε = 0.02, and A generated from the distribution in Eq. 8 with α = − 0.655 and β = −0.912.