## Abstract

Chronic myeloid leukemia (CML) is characterized by a specific chromosome translocation, and its pathobiology is considered comparatively well understood. Thus, quantitative analysis of CML and its progression to blast crisis may help elucidate general mechanisms of carcinogenesis and cancer progression. Hitherto, it has been widely postulated that CML blast crisis originates mainly via cell-autonomous mechanisms such as secondary mutations or genomic instability. However, recent results suggest that carcinogenic transformation may be an inherently multicellular event, in departure from the classic unicellular paradigm. We investigate this possibility in the case of blast crisis origination in CML. A quantitative, mechanistic cell population dynamics model was employed. This model used recent data on imatinib-treated CML; it also used earlier clinical data, not previously incorporated into current mathematical CML/imatinib models. With the pre-imatinib data, which include results on many more blast crises, we obtained evidence that the driving mechanism for blast crisis origination is a cooperation between specific cell types. Assuming leukemic–normal interactions resulted in a statistically significant improvement over assuming either cell-autonomous mechanisms or interactions between leukemic cells. This conclusion was robust with regard to changes in the model's adjustable parameters. Application of the results to patients treated with imatinib suggests that imatinib may act not only on malignant blast precursors, but also, to a limited degree, on the malignant blasts themselves. *Cancer Res; 71(8); 2838–47. ©2011 AACR*.

A comprehensive mechanistic model gives evidence that the main driving mechanisms for CML blast crisis origination are interactions between leukemic and normal cells.

## Introduction

Chronic myeloid leukemia (CML) is regarded as one of the best understood cancers, being characterized by Ph^{+} cells, that is, cells having a Philadelphia (*BCR-ABL)* chromosome translocation (1). In the 1990s, the tyrosine kinase inhibitor imatinib mesylate (imatinib), which suppresses *bcr-abl* oncoprotein action, was found to improve prognosis dramatically. However, imatinib does not definitively cure the disease. Many patients may need to continue taking the drug indefinitely. In some cases, the treatment fails and the disease progresses (2). The use of other tyrosine kinase inhibitors has mitigated but not fully solved this problem (3).

CML has a simpler etiology than most cancers (4); moreover, its time course is comparatively easy to monitor in the clinic (5–7). Consequently, despite being considerably less prevalent than major solid tumors, CML has often been regarded as a kind of “model organism” for quantitative modeling of carcinogenesis (8); as such it is potentially important for many areas, including radiation risk estimation. Relevant to our quantitative approach here is a recent cell population dynamical model of CML and its response to imatinib treatment (9). This model was grounded in the biology of hematopoietic differentiation [reviewed in (10, 11)], extended by parameter calibration using data from CML clinical trials with imatinib (12), and used to analyze treatment protocols (13). Figure 1 summarizes some of its methods. Many other mathematical/computational CML/imatinib models have been suggested in the past decade [reviewed, e.g., in (5, 8, 14)]. In the Discussion section below, we will compare our present approach to some of these models.

Our mathematical model consists of the following:

Coupled, nonlinear, ordinary differential equations for time-dependent cell numbers in a cell population consisting of various normal (Ph

^{−}) and leukemic (Ph^{+}) hematopoietic cell subpopulations;Parameter values and initial conditions for the differential equations;

A stochastic formalism for determining blast crisis incidence from solutions of the differential equations.

The differential equations for Ph^{−} and Ph^{+} cells are respectively:

In Equations (A) and (B), and in Fig. 1:

*t*is time, with the origin of time taken as the time of clinical CML presentation (Fig. 2A).Variables

*x*with subscripts denote numbers of normal (i.e., Ph^{−}) cells (Fig. 1). For example,*x*_{q}(*t*) is the number of normal quiescent stem cells at time*t.*Leukemic (Ph

^{+}) cell numbers are denoted similarly, by*y*.*n*_{0}≡*x*_{0}+*y*_{0}is the total number of cycling stem cells at time*t*, leukemic plus normal, used to model a nonlinear coupling between Equations (A) and (B).The remaining letters are adjustable parameters which characterize cellular proliferation, differentiation, density effects, apoptosis, and transitions between compartments (Fig. 1). See Materials and Methods ("Parameter interpretation and values") for details.

The underlying assumptions Equations (A) and (B) implement are stated and motivated in (9) and its references. In brief, with appropriate parameter values the solutions of the differential equations can explain clinically observed behavior of imatinib-treated CML in terms of different proliferation properties for different cell types in the hematopoietic hierarchy.

Initial conditions for the differential equations are obtained as follows (see Fig. 2A). In the absence of Ph^{+} cells there is an equilibrium solution to Equation (A) that is approached after a long time and corresponds to dynamical homeostasis of a normal mature blood system. At negative time *t* = −*L*** _{CML}**, equilibrium with no Ph

^{+}cells present is assumed; then 1 Ph

^{+}cycling stem cell is introduced. Here

*L*

**is the “CML latency time,” calculated as follows: the system is tracked using Equations (A) and ( B) and**

_{CML}*L*

**is the interval until there are 10**

_{CML}^{12}mature Ph+ cells (“manifest CML” in Fig. 2A). That manifest CML time is identified with the time of clinical presentation, and, for treated patients, with the start of imatinib treatment; it is chosen as the time origin

*t*= 0.

The blast crisis model compares mechanisms involving pairwise cell–cell cooperation for blast crisis origination (e.g., Fig. 2B2) with standard 1-cell mechanisms (e.g., Fig. 2B1), using stochastic irreversible mass-action kinetics. Specifically, our basic assumption is:

Blast crisis origination is a time-inhomogeneous Poisson process with rate λ*Λ*(*t*−*L*_{BC}), where

Here the following interpretations and comments hold.

“Origination” of blast crisis is defined as the appearance of a blast clone sufficiently large that its probability of stochastic extinction is negligible. Blast crisis origination (labeled by

*t = U*in Fig. 2A) could occur probabilistically at any time*t*> −*L*._{CML}The interplay between stochastic versus deterministic calculations is based, as often (5), on small versus large subpopulation numbers; it is summarized in the Supplementary Materials section S2.

Λ(

*t*) is the mass-action term for the incoming cell types involved. For example, if blast cells are produced by a 1-cell mechanism, as in Fig. 2B1, then*z*_{b}= 1 and Λ(*t*) is simply the appropriate cell number, for example, Λ =*y*_{0}(*t*) for the mechanism of Fig. 2B1. If the mechanism for producing blast cells is a cell–cell cooperation event, then Λ(*t*) in Equation (C) is bilinear or quadratic. For example, for the mechanism in Fig. 2B5,*y*_{a}=*y*_{0}and*z*_{b}=*x*_{1}in Equation (C).*λ*is an adjustable parameter for each mechanism, referred to henceforth as the “intensity parameter.” The dimensions of*λ*are inverse time.After blast crisis origination, there is an unknown (and with present techniques not directly observable) “blast crisis latency time”

*L*until “manifest” blast crisis, that is, until accelerated phase or blast crisis can be diagnosed by standard criteria. We treat_{BC}*L*as a parameter to be adjusted from the data. Manifest CML and manifest blast crisis are different events, and blast crisis latency time_{BC}*L*is less than CML latency time_{BC}*L*(Fig. 2A)._{CML}

In Fig. 2B6, we consider nonmyeloid and nonleukemic cells, assuming their number is time-independent; this number can be absorbed into the intensity parameter *λ* and then the mathematical treatment of the cases in Fig. 2B1 and B6 becomes identical, despite the biological differences, so we shall not separately consider the case of Fig. 2B6 or similar cases in the further calculations. As will be seen below, this procedure leaves 45 different blast crisis mechanisms to consider.

Untreated CML usually progresses to the accelerated phase and blast crisis, for which the prognosis is poor. Malignant blasts arise that apparently differ irreversibly from any of the cell types found in the chronic stage (15–17). The malignant blasts often have additional, heterogeneous karyotype alterations (7, 8, 18), and exhibit genomic instability (19, 20). Morphologically and immunophenotypically, they generally show characteristics of both leukemic stem cells and some other cell type: in most cases, granulocyte-monocyte progenitor cells (17); in a substantial minority of cases, immature members in the lymphoid, rather than myeloid, branch of the hematopoietic developmental tree (17, 21); and in a few cases cell types in the myeloid branch different from granulocyte-monocyte progenitors [e.g., (22–24)]. The mechanism of malignant blast origination remains unclear. Often, an additional point mutation or other genetic/epigenetic alteration in a Ph^{+} cell is postulated (15). Such alterations have been modeled quantitatively as autonomous “one-cell” processes [e.g., (16)].

However, there is very substantial evidence that intercellular interactions play an important role in tumor development for both leukemias and solid tumors. Among relevant growth-modulating intercellular interactions reported are the following: cell–cell signaling [e.g., (25–28)], which has been suggested as underlying a feedback mechanism responsible for CML blast crisis (29); tumor interactions with the immune system [e.g., (30, 21)]; and niche-epithelial cell contact effects [e.g., (32)], in particular for the hematopoietic system (33, 34). Mounting evidence (35) suggests that multiple cell types may play a direct role in tumor clone generation, a departure from the classic concept that tumorigenesis is an autonomous one-cell event. Here, we investigate the possibilities that a blast clone may originate from cooperating Ph^{+} cells, or even from cooperating Ph^{+} and normal cells. For concreteness we consider only pairwise cell–cell cooperation (Fig. 2); our results do not exclude 3-way or higher-order intercellular interactions as blast originating mechanisms.

Current computational CML/imatinib models deemphasize clinical CML data from the pre-imatinib era. However, earlier data remain important. They not only contribute to understanding untreated CML as a cancer paradigm, they can also, as we shall show, help calibrate the parameters of recent models. A 1987 data compilation (36), not previously incorporated into the recent models, includes many more observed blast crises than do the imatinib era data sets, and thus gives better statistics. We first use this pre-imatinib data to help establish the most likely mechanisms of blast crisis origination, then briefly compare resulting predictions to data on blast crises in imatinib-treated patients.

## Methods

### Parameter interpretation and values

In Equations (A) and (B) and Fig. 1, the following intuitive remarks apply.

The parameter

*p*determines how much the cycling stem cell proliferation rate is decreased by crowding, with leukemic cells tolerating overcrowding better than normal ones (Table 1).The parameter ratio α/β determines the ratio of quiescent to cycling stem cells at equilibrium. Biologically, α/β > 1 (i.e., more quiescent than cycling stem cells) is preferred (37).

With no leukemic cells and the parameter values used in our main calculation (Table 1), dynamical equilibrium for Equation (A) occurs when the number of quiescent stem, cycling stem, progenitor, differentiated, and mature cells is, respectively, 9 × 10

^{6}, 10^{8}, 10^{10}, and 10^{12}.The parameter

*r*is a maximum rate for cycling stem cell proliferation by symmetric division, estimated to be larger for leukemic than for normal cells in the absence of treatment (Table 1). This parameter is closely related to the calculated CML mean latency time*L*(Supplementary Materials S1 gives some details)._{CML}The

*e*_{j}parameters are compartment outflows due to both apoptosis and (except in the case of*e*_{3}) cell differentiation.Parameters

*a*and*c*are the effective input rates for formation of more differentiated by less differentiated cells, taking into account multiplication within downstream compartments. For example, the effect of divisions within the compartment of myeloid progenitors, involving an increase in numbers coupled with transition to more differentiated subcompartments, is included in the parameter*a*. In the model, Ph^{+}cells have somewhat larger values for*a*and*b*than do normal cells. However, for Ph^{+}cells, imatinib treatment decreases*a*(Table 1) putatively due to increased apoptosis of Ph^{+}cycling stem cells and/or less proliferation within the Ph^{+}progenitor cell compartment. Similarly imatinib decreases*b*for Ph^{+}cells (Table 1).Overall, imatinib treatment is modeled as follows (Table 1): a drastic decrease in some leukemic cell compartment input parameters reflecting both apoptosis and decreased downstream proliferation; and a moderate decrease in the maximum symmetric proliferation rate for cycling leukemic stem cells.

Parameters (days^{−1})
. | . | Normal cells . | Leukemic cells^{a}
. | . | |
---|---|---|---|---|---|

. | . | . | Imatinib^{−}
. | Imatinib^{+}
. | |

Cycling/quiescent stem cell transitions | α | 9.0E-4 | |||

β | 1.0E-4 | ||||

Cell number outflow | e_{0} | 3.0E-3 | |||

e_{1} | 8.0E-3 | ||||

e_{2} | 5.0E-2 | ||||

e_{3} | 1.0 | ||||

Cell number inflow and proliferation | a | 0.8 | 1.6 | 1.6E-2 | |

b | 5.0 | 10 | 1/75 | ||

c | 1.0E2 | ||||

Stem cell proliferation | r | 5.0E-3 | 8E-3 | 2.0E-3 | |

Crowding^{b} | p | 6.7E-7 | 1.7E-7 |

Parameters (days^{−1})
. | . | Normal cells . | Leukemic cells^{a}
. | . | |
---|---|---|---|---|---|

. | . | . | Imatinib^{−}
. | Imatinib^{+}
. | |

Cycling/quiescent stem cell transitions | α | 9.0E-4 | |||

β | 1.0E-4 | ||||

Cell number outflow | e_{0} | 3.0E-3 | |||

e_{1} | 8.0E-3 | ||||

e_{2} | 5.0E-2 | ||||

e_{3} | 1.0 | ||||

Cell number inflow and proliferation | a | 0.8 | 1.6 | 1.6E-2 | |

b | 5.0 | 10 | 1/75 | ||

c | 1.0E2 | ||||

Stem cell proliferation | r | 5.0E-3 | 8E-3 | 2.0E-3 | |

Crowding^{b} | p | 6.7E-7 | 1.7E-7 |

NOTE: Values shown were taken from reference (9). The blast crisis latency time was *L*** _{BC}** = 0. The CML latency time (

*L*

**= 11.83) years and the blast crisis intensity parameters (λ, listed in Supplementary Materials S1) were calculated from the model.**

_{CML}^{a}Parameter symbols have primes (e.g., α′). Empty cells have same values as the cells to their left.

^{b}p is dimensionless.

### The Sokal data set

In a data set presented by Sokal in 1987 (36), the number of blast crisis cases—1,635—was markedly larger than in any recent data set. Moreover, serendipitously or by foresight, Sokal's approach was very well suited for our present purposes. He emphasized cell population dynamics, discussed why he thought the clinical treatments involved did not substantially affect blast crisis onset kinetics, estimated the time between manifest blast crisis and death, and discussed issues related to heterogeneity in time of actual clinical diagnosis compared with the time when the disease could in principle have been diagnosed. All these are important to our use of the data.

Unfortunately, the database Sokal used to prepare the paper was lost several years ago during a computer upgrade. We thus resorted to automated analysis of Sokal's Figure 1 in his 1987 paper, a semilogarithmic actuarial survival curve for patients diagnosed with CML but not blast crisis or accelerated phase. Sokal interpreted subsequent deaths in his data set as due to blast crises and, in this aggregate analysis, took the time between manifest blast crisis and death to be 6 months, the same value as suggested by a more recent estimate (38).

Using CurveUNscan curve recognition software (www.curveunscan.com), we scanned Sokal's figure to estimate empirical survival fractions. The output consisted of “survival” (defined here as being alive and not having manifest blast crisis or accelerated phase) for 1,635 equispaced time intervals from year 0 (clinical presentation) to 15 years. Data cleansing involved the following: changing a few initial survival fractions so that survival never increased as time progressed; translating the curve by 6 months corresponding to Sokal's estimate of the time interval between manifest blast crisis and death; and determining an overall scale from the condition that survival (as defined above) was 100% for this cohort at the time of clinical CML presentation.

### Blast crisis mechanisms

With *y* and *z* as in Equation (C), there are 5 × 11 = 55 products and 10 redundancies (e.g., *y*_{0}*y*_{1} ≡ *y*_{1}*y*_{0}) so the number of distinct autonomous 1-cell plus pairwise interaction mechanisms is 45 (examples are shown in Fig. 2B). Many of these mechanisms give quite similar results computationally, reflecting close coupling among some of the cell subpopulation numbers.

Our main assumption implies by Poisson statistics that the survival *S*(*t*) (i.e., the probability that manifest blast crisis has not occurred by time *t*) involves “blast integrals” |$\mu (t) = \int_{ - L_{\rm BC} }^{t - L_{\rm BC} } {\Lambda (\tau )d\tau }$| (with integrand *Λ* given by Equation (C)), as the following calculation shows

For a given *L*** _{BC}**,

**μ**was calculated directly from the model characterized earlier [in Equations (A) and (B) and Table 1]. Then

*λ*was estimated by fitting to Sokal's data, using a nonlinear least squares method. Details on the statistical techniques used are provided in Supplementary Materials S3.

## Results

We first performed the main calculation as follows: (i) we calculated the time-dependent number of cells of various types; (ii) given the cell numbers, we calculated the timing of blast crisis origination for 45 different candidate blast crisis mechanisms (a few of which are shown in Fig. 2B), using comparatively extensive data from the pre-imatinib era to determine intensity parameters; and (iii) we identified preferred mechanisms. After carrying out the main calculation, we checked robustness of the results by a number of auxiliary calculations. Finally our results were applied to estimate the (small) chance of blast crisis for imatinib-treated CML, and briefly compared with current data.

### Main calculation

Our main calculation used the parameter values specified in Table 1, for example, assumed that in equilibrium prior to initiation of CML there were more quiescent than cycling normal stem cells. The calculation used blast crisis latency time *L*** _{BC}** = 0, corresponding to very rapid blast expansion after blast crisis origination (Fig. 2A). Equations (A) and (B) gave the cell numbers shown in Fig. 3 and a CML latency time of

*L*

**= 11.83 years.**

_{CML}Using the cell numbers in Fig. 3, we compared 45 possible blast crisis mechanisms (some of which are shown schematically in Fig. 2B) to the data of reference (36) by using Equations (C) and (D). The results are shown graphically in Fig. 4.

Using nonlinear regression and a log-rank test to estimate goodness of fit, as described in Supplementary Materials S3, we assigned, to each of the 45 mechanisms, *P* values for the null hypothesis that the predicted survival and Sokal data do not differ. The *P* values are organized in Table 2 according to the cell–cell associations; for example, the (*y*_{0}*,x*_{1}) box in the table refers to associations between cycling leukemic stem cells (*y*_{0}) and normal myeloid progenitor cells (*x*_{1}) in Equation (C). It was apparent throughout that in such tables the entries involving *y*_{2}, *y*_{3}, *x*_{2}, and/or *x*_{3} are always quite similar to entries involving progenitor cells instead. The intuitive reason is that changes in the differentiated cell numbers mirror changes in progenitor cell numbers with only a comparatively small time delay. We omitted the columns and rows for differentiated cells in Table 2 to focus attention on essential differences, and to take into account the fact that stem and progenitor cells, unlike more differentiated cells, are found predominantly in the bone marrow, where there are niches and more opportunities for orchestrated cell–cell associations. Tables without the omissions, from which Table 2 was condensed, are given in the Supplementary Materials S1 and S4 for completeness.

. | y_{q}tblfn3
. | y_{0}
. | y_{1}
. | 1 . | x_{q}
. | x_{0}
. | x_{1}
. |
---|---|---|---|---|---|---|---|

y_{q} | <E-22 | <E-15 | <E-16 | <E-11 | <E-10 | <E-5 | <E-5 |

y_{0} | ≤E-7 | <E-8 | 3E-3 | 0.01 | 0.76 | 0.89 | |

y_{1} | <E-9 | <E-3 | 2E-3 | 0.40 | 0.54 |

. | y_{q}tblfn3
. | y_{0}
. | y_{1}
. | 1 . | x_{q}
. | x_{0}
. | x_{1}
. |
---|---|---|---|---|---|---|---|

y_{q} | <E-22 | <E-15 | <E-16 | <E-11 | <E-10 | <E-5 | <E-5 |

y_{0} | ≤E-7 | <E-8 | 3E-3 | 0.01 | 0.76 | 0.89 | |

y_{1} | <E-9 | <E-3 | 2E-3 | 0.40 | 0.54 |

^{a}Columns *y*_{q} to *y*_{1} specify pairwise leukemic–leukemic cell–cell association mechanisms, the column labeled “1” specifies 1-cell mechanisms, columns *x*_{q} to *x*_{1} specify pairwise leukemic–normal cell–cell association mechanisms. Underlined values are the preferred ones for these 3 cases. Here E refers to powers of 10, for example E-3 = 1/1000.

Results in Table 2 show statistically significant differences between observations and predictions (*P* < 0.05) for all but the 4 leukemic–normal cell interactions in the lower right. Of these 4 mechanisms, it is seen that particularly high *P* values are obtained for the 2 mechanisms involving cycling leukemic stem (i.e., *y*_{0}) cells. For leukemic–leukemic or 1-cell mechanisms, the mechanism with the highest *P* value also includes *y*_{0} cells, but the *P* values are markedly smaller. Underlined entries in Table 2 were chosen as preferred representatives of the 3 different colors in Fig. 4. The difference between colors is striking even allowing for Bonferroni-like corrections for multiple comparison bias (39). We estimate that at most 9 different possibilities, corresponding to *y* = (*y*_{q}, *y*_{0}) and *z* = (*y*_{q}, *y*_{0}, 1, *x*_{q}, *x*_{0}) in Equation (C), are effectively independent for our purposes.

### Effect of changes in the blast crisis latency time

We next considered results obtained by replacing blast crisis latency time *L*** _{BC}** = 0 with various time lags. We conjectured that large times

*L*

**> 3 years would lead to poor fits, the rationale being the following.**

_{BC}*L*

**> 3 years would imply that more than half of the patients who present with CML and no manifest blast crisis already have cryptic malignant blast clones, since absent treatment most manifest blast crises in such patients occur within 3 years of manifest CML (36). Given this pattern, imatinib treatment would not be able to stave off blast crisis for long because the drug is not fully effective against malignant blasts (3). But the imatinib era data show that, to the contrary, imatinib is very effective in preventing or long-postponing blast crises in a large majority of cases. This contradiction indicates that**

_{BC}*L*

**< 3 years. For small values of**

_{BC}*L*

**, the contradiction does not arise since the action of imatinib on blast cell precursors could explain its effectiveness.**

_{BC}Table 3 shows results for 3 non-zero values of *L*** _{BC}**, using the same format as Table 2. Table 3 shows that leukemic–normal pairwise associations are robustly preferred as blast crisis origination mechanisms even for

*L*

**≠ 0. It is also seen that, as conjectured, large blast crisis lag times are inconsistent with data on untreated CML. This confirmation of our conjecture highlights how a quantitative model can integrate disparate data sets (here pre- and post-imatinib era data) in unexpected ways to elucidate the observations.**

_{BC}L_{BC}
. | . | y_{q}
. | y_{0}
. | y_{1}
. | 1 . | x_{q}
. | x_{0}
. | x_{1}
. |
---|---|---|---|---|---|---|---|---|

4 Months | y_{q} | <E-24 | <E-17 | <E-18 | <E-12 | <E-11 | <E-6 | <E-6 |

y_{0} | <E-9 | ≤E-10 | ≤E-3 | 2E-3 | 0.51 | 0.66 | ||

y_{1} | <E-11 | <E-4 | .02 | 0.20 | 0.30 | |||

11 Months | y_{q} | <E-28 | <E-22 | <E-23 | <E-15 | <E-14 | <E-8 | <E-7 |

y_{0} | <E-14 | ≤E-15 | ≤E-5 | <E-4 | 0.16 | 0.23 | ||

y_{1} | <E-17 | <E-7 | <E-6 | 0.03 | 0.05 | |||

19 Months | y_{q} | <E-32 | <E-27 | <E-29 | <E-18 | <E-18 | <E-10 | <E-10 |

y_{0} | ≤E-21 | <E-23 | ≤E-9 | <E-9 | 7E-3 | 0.01 | ||

y_{1} | <E-25 | <E-12 | <E-11 | <E-3 | <E-3 |

L_{BC}
. | . | y_{q}
. | y_{0}
. | y_{1}
. | 1 . | x_{q}
. | x_{0}
. | x_{1}
. |
---|---|---|---|---|---|---|---|---|

4 Months | y_{q} | <E-24 | <E-17 | <E-18 | <E-12 | <E-11 | <E-6 | <E-6 |

y_{0} | <E-9 | ≤E-10 | ≤E-3 | 2E-3 | 0.51 | 0.66 | ||

y_{1} | <E-11 | <E-4 | .02 | 0.20 | 0.30 | |||

11 Months | y_{q} | <E-28 | <E-22 | <E-23 | <E-15 | <E-14 | <E-8 | <E-7 |

y_{0} | <E-14 | ≤E-15 | ≤E-5 | <E-4 | 0.16 | 0.23 | ||

y_{1} | <E-17 | <E-7 | <E-6 | 0.03 | 0.05 | |||

19 Months | y_{q} | <E-32 | <E-27 | <E-29 | <E-18 | <E-18 | <E-10 | <E-10 |

y_{0} | ≤E-21 | <E-23 | ≤E-9 | <E-9 | 7E-3 | 0.01 | ||

y_{1} | <E-25 | <E-12 | <E-11 | <E-3 | <E-3 |

### Other robustness calculations

We performed further robustness checks by varying parameters (in a “1 at a time” approach), as described in Supplementary Materials. For example, analyses varying the leukemic stem cell proliferation parameter *r*′ showed that a leukemic–normal cell association remains the preferred blast crisis origination mechanism in all cases (Supplementary Materials S1). This result was particularly important because the value used for *r*′ in our main calculation is somewhat problematic. Originally, *r*′ was estimated from results on the Japanese atomic bomb survivors using older data (12), which have meanwhile been revised (40–42); moreover the estimate of *r*′ neglected some complications due to possible stochastic extinction of small Ph^{+} clones (see Supplementary Materials S2). Robustness of our main result under variations of *r*′ shows that these problems are technical rather than crucial.

We also examined how changing the statistical methods used affected our conclusions (see Supplementary Materials S3). We found our main results and conclusions were robust with respect to alterations of the statistical approach.

### Blast crises in imatinib-treated CML

Starting from a formalism originally designed for analyzing CML response to imatinib treatment, we found that data on blast crises in the pre-imatinib era can be used to characterize and calibrate a model that analyzes the cause and estimates the time of blast crisis origination. We therefore used the calibrated results for estimates of blast crisis incidence in imatinib-treated patients, assuming imatinib prevents or delays blast crisis by decreasing the number of leukemic cells that can act as precursors of the blast cells, rather than by acting directly on blast cells. The estimates were compared with recent data (43), with no further adjustable parameters introduced. It was found that, despite the absence of adjustable parameters, the predictions and data have the same order of magnitude and similar shapes. Because the shapes correspond, introducing 1 extra adjustable parameter for each curve, that is, allowing λ to be modified by imatinib treatment (corresponding to some imatinib action directly on the malignant blasts, not just on their precursors), can give a close match between predictions and observations (results not shown).

## Discussion

### Recapitulation

To assess the role of cell–cell cooperation in CML blast crisis formation, we analyzed alternative mechanisms for blast crisis origination. Throughout, we emphasized data for the actual end points of interest—human hematopoiesis and CML—rather than surrogate (i.e., *in silico, in vitro*, animal, or loosely related human) end points. We used data, not previously incorporated explicitly into current quantitative CML models, from the pre-imatinib era, when many more blast crises were observed. Assuming a pairwise cooperativity mechanism between leukemic and normal cells for blast crisis origination robustly gave a statistically significant improvement over assuming either a pairwise cooperation between leukemic cells or an autonomous 1-cell action. Applying the model with pairwise cooperation between normal progenitor and leukemic cycling stem cells to blast crises for imatinib-treated patients gave predictions qualitatively consistent with the data. Thus, within the framework of our model, we found that cooperation between leukemic and normal cells is the indicated mechanism for origination of blast crisis. We suggested this result could provide fundamental insights into the development of other cancers.

### Other computational models of CML

One limitation of the present treatment is a need for more systematic calibration of the parameters in the underlying CML/imatinib model summarized above in Equations (A) and (B); at present we have not found a convincing way to assign error bars to the values in Table 1 and thus have not developed an estimate of whether the model is overparameterized. However, a more general limitation is the fact that robustness of our basic conclusions under variations in the underlying CML/imatinib model has not been demonstrated.

In other quantitative CML/imatinib models that have been suggested in the last decade (5, 8), the mathematics involved is usually quite straightforward, but there are marked differences in the underlying assumptions and in the parameter estimates. A series of papers [e.g., (14, 44, 45)] emphasize additional modern clinical data, cell-cycle effects, and the possibility of stochastic extinction of small Ph^{+} clones. Another series of papers incorporate additional stochastic effects and more specific biological information about hematopoiesis [reviewed in (46)]. Spatial effects, potentially important to any model that includes cell–cell association effects, have also been considered recently [e.g., (47)], as have immune system responses to leukemic cells [e.g., (30, 31)]. A model which considers feedback from differentiated to pluripotent cells (29) emphasizes intercellular interaction blast crisis mechanisms. This model could not be compared in detail to our approach because its parameters have not been calibrated with clinical data.

These and other recent CML/imatinib models disagree substantially on the relevant cell population dynamics—for example, on how many Ph^{+} cycling stem cells there are; whether these have a growth advantage over their Ph^{−} counterparts; whether imatinib affects them; and how rapidly alterations in progenitor cell numbers are reflected in the number of differentiated cells. Blast crisis origination models also differ on key points. For example, several [e.g., (29, 48)] regard blast crisis as a different dynamic state of the same cell types present in the chronic phase of CML rather than assuming the alternative (15–17), that malignant blasts differ irreversibly. Here we assumed the latter and suggest the blast phenotypic diversity reviewed in the Introduction perhaps indicates origination by interactions of Ph^{+} cells with diverse cell types.

## Prospects

In view of the wide range of suggested formalisms, and disagreements among them, deeper quantitative analysis of blast crisis now likely depends on first choosing and calibrating a consensus mathematical/computational model for CML under imatinib treatment. Developing a consensus model will require use of large, up-to-date clinical databases and careful selection of which biological effects to deemphasize to avoid over-parameterization. Critically comparing various current models that emphasize cell population dynamics and intercellular interactions while using subcellular data as auxiliary rather than primary is perhaps the most promising starting point.

It is already clear from the results of the present paper that a major effort to record detailed information on the proportions of various cell subpopulations at the time of CML presentation and shortly thereafter would be an important aid to such modeling. Almost all current models consider the time-course of untreated CML. With rare exceptions, corresponding longitudinal studies are no longer carried out in the clinic. However, due to variability in presentation times, the cross-sectional information available from blood and marrow samples at the time of presentation, together with modeling, could help fill the resulting gap.

The primary reason for constructing a credible, generally accepted quantitative model of the action of imatinib and tyrosine kinase inhibitors on CML would of course be finding insights that improve CML treatment. The results of the present paper suggest strongly that such a model would also improve our understanding of how blast crisis originates—whether it is basically a unicellular or, as our findings suggest, a multicellular process. This improvement in turn should help clarify analyses of how other blood and solid cancers, more complex than CML, evolve.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Acknowledgments

## Grant Support

This study was supported by NASA NSCOR NNJ06HA28G (L. Hlatky), NIH ICBP 1U54CA149233 (L. Hlatky), and DOE grant DE-SC0001434 (P. Hahnfeldt).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked *advertisement* in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.