Many tumors are hierarchically organized and driven by a subpopulation of tumor-initiating cells (TIC), or cancer stem cells. TICs are uniquely capable of recapitulating the tumor and are thought to be highly resistant to radio- and chemotherapy. Macroscopic patterns of tumor expansion before treatment and tumor regression during treatment are tied to the dynamics of TICs. Until now, the quantitative information about the fraction of TICs from macroscopic tumor burden trajectories could not be inferred. In this study, we generated a quantitative method based on a mathematical model that describes hierarchically organized tumor dynamics and patient-derived tumor burden information. The method identifies two characteristic equilibrium TIC regimes during expansion and regression. We show that tumor expansion and regression curves can be leveraged to infer estimates of the TIC fraction in individual patients at detection and after continued therapy. Furthermore, our method is parameter-free; it solely requires the knowledge of a patient's tumor burden over multiple time points to reveal microscopic properties of the malignancy. We demonstrate proof of concept in the case of chronic myeloid leukemia (CML), wherein our model recapitulated the clinical history of the disease in two independent patient cohorts. On the basis of patient-specific treatment responses in CML, we predict that after one year of targeted treatment, the fraction of TICs increases 100-fold and continues to increase up to 1,000-fold after 5 years of treatment. Our novel framework may significantly influence the implementation of personalized treatment strategies and has the potential for rapid translation into the clinic. Cancer Res; 76(7); 1705–13. ©2016 AACR.

Major Findings

Under the cancer stem cell hypothesis, a tumor population is driven by a fraction of cells with the capacity to self-renew. Absolute and relative size of this population in human cancers remains unknown at any stage of the disease. We develop a parameter-free mathematical model that allows estimation of cancer stem cell fractions from longitudinal measurements of tumor burden, which are often available from imaging or liquid biopsies. Our novel framework provides critical information for personalized treatment strategies and only requires routinely available clinical information.

Quick Guide to Equations and Assumptions

The continuous proliferation and production of cells in a stem cell–driven tissue with multiple cell differentiation stages is well described by a hierarchical compartment model (24–26). Each compartment represents a certain stage of cell differentiation. Tumor-initiating cells (TIC), or cancer stem cells (CSC), are at the root of the hierarchy (Fig. 1). TICs divide at a rate S and either self-renew (producing two initiating cells) with probability p or differentiate (producing two differentiated cells) with probability 1−p. We neglect the possibility of asymmetric cell divisions, which would not contribute to the growth or shrinkage of the initiating cell population and only influence the time scale of the dynamics but not the underlying growth patterns. In addition, initiating cells might die at a rate d per cell division. Differentiated cells proliferate at a rate D, also die at a rate d per cell division and undergo a maximum of m cell doublings before they enter cell senescence, resembling a Hayflick limit (27, 28). Under these assumptions, the system takes the form of a hierarchically ordered, coupled set of ordinary differential equations that count the influx and outflux of cells of each compartment

formula

Here, N0(t) counts the number of TICs at time t, N1(t) corresponds to the number of differentiated cells that still have m cell replications left before they undergo senescence, and Nm(t) is the number of differentiated cells that have exhausted their replication capacity and are removed from the system (reached the Hayflick limit).

TIC fraction during tumor growth

The system of differential equations (1) can be solved recursively for general initial conditions. If we set Ni(0) as the initial number of cells at time t = 0 in compartment i, we find N0(t) = N0(0)eαt, where α = (1 + d – 2p)S is the net growth of the TIC population. Higher compartments (i > 0) of differentiated cells evolve according to

formula

The outflow of each differentiated compartment is β=(1+d)D, where γ=βα corresponds to the differential outflow of stem and nonstem compartments.

The tumor population's transition from a fast into a slower growth regime is determined by the signs of α and β. As D and d are strictly positive, all terms in Eq. B that contain eβt vanish in the long run. If ρ > (1 + d)/2, α is positive and determines long-term tumor growth. Therefore, if we start from a single TIC, Ni(0) = 1 only, the total number of all differentiated cells grows by

formula

and follows an exponential growth curve with an offset ag that involves all microscopic parameters. The fraction r(t) of TICs is then formally given by

formula

but in the second growth phase described by Eq. C, this expression simplifies considerably:

formula

The fraction of TICs becomes time independent and stays constant during further tumor growth. This does not imply that the number of TICs is constant. Both the total tumor mass and the TIC population continuously expand with time. However, after a sufficiently long time, the growth characteristics of the total tumor population reaches a steady state governed by self-renewal.

TIC fraction during treatment

During successful treatment, tumor reduction in equilibrium can be written as |$\bar N(t) = {a_t}\,\,{\rm{e}^{- {\alpha _g}t}}$|⁠. In the reduction phase, the TIC fraction contains an additional term due to a change in initial conditions. Instead of a single TIC, we have |${\rm{e}^{- {\alpha _g}T}}$| TICs at time of diagnosis T. Consequently, the fraction of TICs is calculated as |r(t) = {\hbox{$1$} \!\mathord {/{{{({1 + {a_t}\,\,{\rm{e}^{- {\alpha _g}T}}})}}}\kern-\nulldelimiterspace} $}}|⁠. By resubstituting the age of the tumor, we find an expression for the fraction of TICs under treatment

formula

Again, the TIC fraction transitions into a dynamic equilibrium during continued treatment.

Cancer comprises a group of diseases that involve abnormal and uncontrolled proliferation of cells. These aberrant properties are induced by alterations in genes that control cell-regulatory mechanisms, microenvironmental response, and cell–cell signaling (1). Large-scale genomic studies have revealed the spectrum of genomic profiles in many cancers (2) but accumulating evidence shows that cancers are characterized by extensive inter- and intratumor heterogeneity (3, 4). Tumor heterogeneity is a consequence of tumor evolution (5), bridges multiple scales (6, 7), and may imply complex interactions (8–10).

Tumor cells are often found in distinct stages of differentiation (11). This phenotypic diversity likely is a remnant of the hierarchical organization of the tissue of origin. In most healthy tissues, stem cells maintain tissue homeostasis. A certain number of differentiation compartments give rise to the production of distinct mature cell types (12, 13). The finding that tissue organization can be maintained in tumors has led to the postulation of the existence of tumor-initiating cells (TIC; cancer stem cells). Under this hypothesis, a fraction of tumor cells are uniquely able to seed, maintain, and reseed tumors (11). First identified in leukemia (12), TICs have since been shown to drive a number of solid tumors, including colon (14–16), brain (17), breast (18), head and neck (19), and lung (20), among others.

The implicit self-renewing feature of TICs has led to the conclusion that the effective treatment of hierarchically organized cancers can only be achieved if all TICs are eradicated (21–23). However, the fraction of TICs is unknown at diagnosis. Currently, the only method that can infer this information for specific patients is from direct biopsies before and after treatment. This method has major limitations due to the poor reliability of stem cell markers and important sampling biases. Ideally, we require a continuous measure of the TIC fraction that can be obtained by minimally invasive procedures.

Here we combined mathematical modeling and clinical data to make predictions about the fraction of TICs. We used a compartment model approach that describes cell population dynamics in hierarchically organized tumors (26, 29). It has been observed that during treatment, hierarchical tumors tend to transition from a fast into a relatively slow response phase, particularly in the targeted treatment of leukemias (1, 27, 30). These transitions are tied to differential dynamics between stem-like and differentiated compartments. Our results allowed us to exploit a universal property of hierarchical tumor organization. One can estimate the fraction of tumor-driving TICs from purely macroscopic observables, for example, tumor burden, routinely gleaned from medical imaging or liquid biopsies.

Our estimates of the tumor initiating cell fraction in individual patients are parameter-free; they do not depend on microscopic parameters of cell-cycle dynamics and cellular interactions. Knowledge of cell proliferation or death rates is not required. We used our framework to infer patient-specific TIC fractions from two independent chronic myeloid leukemia (CML) datasets. Each dataset recorded patient treatment responses over time during continued treatment (24, 30). Our analysis allowed estimation of the actual number of TICs before and after treatment. This information can be invaluable to gauge treatment response and inform future personalized treatment strategies.

We used a system of inhomogeneous linear ordinary differential equations to describe the dynamics of tumor cells. We thereby described two distinct tumor cell populations. The first cell population consists of stem-like TICs. These cells can potentially divide infinitely many times. At each point in time, a TIC can either self-renew into two TICs, differentiate into two differentiated cells, or die. The second cell population consists of differentiated cells that can only divide a finite number of times, or die. If differentiated cells can divide m number of times, the system consists of m+1 ordinary differential equations. We solved this system by variation of parameters and wrote down an exact mathematical expression for the fraction of TICs. We observed that the exact analytic form of this expression followed a simple mathematical form. This form does not require knowledge of the proliferation and death rates of stem-like and differentiated cells, but only the macroscopic behavior of tumor burden (Fig. 2). Hence, observing a patient's tumor trajectory right before and after sufficiently long treatment, we could leverage the result of our mathematical approach to estimate the fraction of TICs. We used two previously published datasets by Michor and colleagues (1) and Roeder and colleagues (30), which followed CML patients' tumor burden over time during continued treatment. The workflow of data analysis is given in of the Supplementary Fig. S1 in Supplementary Information. Supplementary Figures S2 and S3 present individual patient trajectories used for our numerical estimates. From this analysis, we calculated distributions of the TIC fraction before and after treatment, which differed significantly (Fig. 3); we used Wilcoxon signed rank test as implemented in the software Wolfram Mathematica 10.

Figure 1.

Model schematic, showing key parameters governing the mathematical model based on the assumption that there exists a population of cells that have the capacity to self-renew (tumor-initiating cells or cancer stem cells). The tumor-initiating cell population of size N0 is exclusively able to maintain the tumor due to occasional self renewal. Transit amplifying cells (N1,…, Nm) undergo m cell division before they enter cell senescence. Tumor initiating cells (TIC) proliferate with a rate S, self-renew with probability p and die at a rate S d. Transit-amplifying cells proliferate with rate D and die at a rate D d. The tumor-initiating cell pool expands via symmetric stem cell divisions and the hierarchy is filled with transit amplifying cells via cell differentiation. Typically, the tumor is initiated with a single TIC in compartment 0.

Figure 1.

Model schematic, showing key parameters governing the mathematical model based on the assumption that there exists a population of cells that have the capacity to self-renew (tumor-initiating cells or cancer stem cells). The tumor-initiating cell population of size N0 is exclusively able to maintain the tumor due to occasional self renewal. Transit amplifying cells (N1,…, Nm) undergo m cell division before they enter cell senescence. Tumor initiating cells (TIC) proliferate with a rate S, self-renew with probability p and die at a rate S d. Transit-amplifying cells proliferate with rate D and die at a rate D d. The tumor-initiating cell pool expands via symmetric stem cell divisions and the hierarchy is filled with transit amplifying cells via cell differentiation. Typically, the tumor is initiated with a single TIC in compartment 0.

Close modal
Figure 2.

Inferring the fraction of TICs from tumor growth curves and treatment response. A and B, the solid line shows one realization of the model for tumor growth (from a single cell), followed by an example of treatment response. The dashed lines correspond to linear regression of the log-linear tumor size over time. We specifically marked the offsets, which allow calculation of the TIC fraction without the need for detailed knowledge of the cell properties. If the offsets ηg and ηt can be estimated from the regression during growth (g) and treatment (t), respectively, and the tumor size at the beginning of treatment is |$\bar N$|⁠, then one can infer the equilibrium fraction of TICs during tumor growth [phase (ii)] and treatment response [phase (iv)] from purely macroscopic observables. C, four CML patients from Michor and colleagues (1) and a pooled group of 68 patients from Roeder and colleagues (D; ref. 30), all treated with imatinib. In both cases, the treatment response transitions from a faster into a slower decline. This is an innate property of hierarchically organized tumors.

Figure 2.

Inferring the fraction of TICs from tumor growth curves and treatment response. A and B, the solid line shows one realization of the model for tumor growth (from a single cell), followed by an example of treatment response. The dashed lines correspond to linear regression of the log-linear tumor size over time. We specifically marked the offsets, which allow calculation of the TIC fraction without the need for detailed knowledge of the cell properties. If the offsets ηg and ηt can be estimated from the regression during growth (g) and treatment (t), respectively, and the tumor size at the beginning of treatment is |$\bar N$|⁠, then one can infer the equilibrium fraction of TICs during tumor growth [phase (ii)] and treatment response [phase (iv)] from purely macroscopic observables. C, four CML patients from Michor and colleagues (1) and a pooled group of 68 patients from Roeder and colleagues (D; ref. 30), all treated with imatinib. In both cases, the treatment response transitions from a faster into a slower decline. This is an innate property of hierarchically organized tumors.

Close modal
Figure 3.

Estimated TIC fraction of chronic myeloid leukemia at diagnosis and treatment. Individual fractions of TICs right before treatment were estimated for each patient of the Michor and colleagues dataset (ref. 1; full disks). We followed the linear regression method described in the main text (also see Quick Guide to Equations and Assumptions). The median TIC fraction was 10−10 (gray line), but variation between individuals was large. Treatment typically selects for TICs: TIC fractions after treatment (circles) led to a median of 10−8 (P = 0.02, Wilcoxon signed-rank test between the two TIC fraction distributions). Interestingly, the variation between individuals was reduced under treatment, potentially due to selection against BCR-ABL cells. Combined with the average decrease of tumor burden, this suggests a slowly shrinking TIC population for continued treatment of several years. Nine of the 46 individual cases show a decrease in TIC fraction. This could be explained either by our conservative fitting procedure or by the fact that the tumor under treatment has not yet reached its equilibrium phase.

Figure 3.

Estimated TIC fraction of chronic myeloid leukemia at diagnosis and treatment. Individual fractions of TICs right before treatment were estimated for each patient of the Michor and colleagues dataset (ref. 1; full disks). We followed the linear regression method described in the main text (also see Quick Guide to Equations and Assumptions). The median TIC fraction was 10−10 (gray line), but variation between individuals was large. Treatment typically selects for TICs: TIC fractions after treatment (circles) led to a median of 10−8 (P = 0.02, Wilcoxon signed-rank test between the two TIC fraction distributions). Interestingly, the variation between individuals was reduced under treatment, potentially due to selection against BCR-ABL cells. Combined with the average decrease of tumor burden, this suggests a slowly shrinking TIC population for continued treatment of several years. Nine of the 46 individual cases show a decrease in TIC fraction. This could be explained either by our conservative fitting procedure or by the fact that the tumor under treatment has not yet reached its equilibrium phase.

Close modal

We model hierarchical tumor organization using a compartment approach (Fig. 1). Each compartment represents cells at a certain differentiation or proliferation stage (2, 3). One can observe that in such systems the tumor growth curve decomposes into two regimes. In the first regime, differentiated compartments are filled by TICs. In the second regime, a dynamic equilibrium is reached: the drive from stem cells is balanced by loss of differentiated cells due to senescence or other factors. In addition, one can infer the TIC fraction during the second regime of tumor expansion, or during the second regime of tumor reduction in response to treatment.

Tumor growth

To model tumor growth, we initiated the system with a single TIC. This is in line with the cancer stem cell hypothesis (31) and does not necessarily imply that the cell of origin was a tissue-specific stem cell. Stem-like properties can be acquired in transit amplifying stages, which is, for example, common in different acute leukemias (17, 28). The proliferation parameters of TICs determine the long-term behavior of the tumor. A tumor grows continuously (and potentially becomes a detectable cancer), if the self-renewal probability p fulfills p > (1 + d)/2. The probability of TIC self-renewal needs to be sufficiently large to compensate loss of TICs due to cell death. The tumor would vanish for insufficient self-renewal rates. We assume that self-renewal is sufficiently common to allow for a growing tumor.

In the scenario of a growing tumor, the TIC population expands. In addition, TICs differentiate into transit-amplifying cells, which form the bulk of the tumor. After sufficient time, one observes a balance between loss of differentiated cells due to cell death or cell senescence, and gain of differentiated cells due to doublings of transit-amplifying cells. In this balanced regime, tumor growth is limited by the expansion rate of TICs, see Fig. 2.

TIC fraction

The fraction of TICs during tumor growth is |{r_g}(t) = {\hbox{${{N_0}(t)}$} {/{{{{\bar{N}}(t)}}|⁠, where |${N_0}(t)$| corresponds to the number of cells with self-renewal potential and |$\bar N(t)$| is the total size of the tumor cell population. The TIC fraction decreases monotonically in the first phase of tumor growth. Then, it evolves towards an equilibrium state for any possible combination of microscopic cell proliferation parameters (Fig. 2). The numerical value of rg in dynamic equilibrium depends on the sign of the differential flow between the initiating and the noninitiating compartments. This flow describes the difference between the net growth of the tumor-initiating cell compartment and the net loss in any given differentiated (noninitiating) compartment due to cell differentiation or death.

If the differential flow is negative, more TICs are lost by differentiation than gained by self-renewal. In this case, the initiating cell fraction tends to zero. If the differential flow is exactly zero, TICs furnish half of the tumor population. For positive values of the differential flow, we have a surplus in the production of differentiated tumor cells as compared with their losses during further differentiation. In this third case, the time-dependent components converge to a constant value and the fraction of TICs takes a nontrivial constant value between 0 and 1. Thus, the relative composition of the growing tumor remains constant in the second phase of tumor expansion and the fraction of TICs is given by

formula

We analytically calculated the value of the constant in the denominator and showed that it involves all model parameters. Thus, from the mathematical model's perspective, a detailed microscopic knowledge of a tumor's properties seems to be a prerequisite to estimate TIC fractions. For example, all proliferation parameters of the different cancer cell types should play a role. Yet, this detailed knowledge is difficult or even impossible to be obtained in a clinical setting. In the following, we propose an alternative and parameter-free method that determines the constant in the denominator of Eq. G. This approach therefore allows the estimation of the fraction of TICs in hierarchically organized tumors.

Estimating the TIC fraction during growth

In dynamic equilibrium, tumor growth is exponential, governed by microscopic parameters of the dynamics. We find that the coefficient α of the exponential growth coincides with the constant in the exact expression of the fraction of TICs. We can thus write |${1}/{{(1 + {\rm{\ralpha})}}}$| for the equilibrium TIC fraction, which allows inference of TIC fractions from tumor growth curves directly. Instead of calculating α analytically, one can fit an exponential growth law to an equilibrium tumor growth curve. This fit gives a slope s and an offset ηg, and requires no microscopic knowledge about the tumor. Moreover, the offset ηg of the regression corresponds to α and leads to

formula

See Fig. 2A.

Estimating the TIC fraction under treatment

We can model treatment strategies by altering the death rates (or other parameters) of cancer cells. A hierarchically organized tumor shrinks continuously under treatment, if the death rate of cancer cells exceeds the self-renewal capability of stem cells: d > 2p − 1. Under continuous treatment, similar to unperturbed growth, we observe that tumor burden reduces fast initially and transitions into slower response after a characteristic time. The initial reduction phase is dominated by the death events in differentiated cells. In this phase, treatment selects for tumor-initiating cells (Fig. 2B). The tumor then reaches a dynamic equilibrium characterized by a balance between cell renewal and loss. The TIC fraction remains constant, despite a continuous decrease in tumor size and the initial treatment effect diminishes.

The active stem cell fraction during treatment response can be estimated by linear regression to the tumor kinetics after transformation, |${\rm {log}} \ \bar N({{\it t}}) = {\rm{log}}\ {\eta _{\rm{t}}} + \sigma t$|⁠, but the initial condition (at time of detection) is not a single seeding TIC. We have to take into account the entire tumor cell population. The offsets ηt (treatment) and ηg (growth) and the total tumor size at treatment initiation time |$T, \bar N({{T}})$|⁠, suffice to accurately estimate the fraction of TICs under treatment

formula

This is exact and parameter-free but requires an estimate of the age of the tumor.

The TIC fraction in chronic myeloid leukemia

We applied our framework to estimate TIC fractions of chronic myeloid leukemia (CML) during tumor growth and treatment response. CML is induced by translocations of chromosome 9 and 22 that result in the expression of the BCR-ABL fusion gene in hematopoietic stem cells (32). CML patients often present with abnormal accumulation of undifferentiated granulocyte progenitor cells in the bone marrow and blood. Targeted anticancer treatment strategies were first introduced in CML. For about 15 years, patients have stayed in remission under continued administration of tyrosine kinase inhibitors, such as imatinib, nilotinib, or bosutinib (33, 34).

We analyzed two independent and previously published datasets (1, 30) of CML treatment response to Imatinib. Relative tumor burden (BCR-ABL/ABL) was measured using reverse transcriptase PCR at diagnosis and was followed at fixed time intervals. The dataset of Michor and colleagues (1) contained 103 patients followed every 3 months for 2 years. In Fig. 2C we show four representative patients, Supplementary Figs. S2 and S3 show the selected patients with clear biphasic response. Figure 2D shows the median treatment response of these patients over time. The median tumor burden shows a fast initial decline followed by slower decline after approximately 6 to 12 months of continued therapy. The Michor and colleagues data permitted individualized estimates of stem cell fractions (Fig. 3). We disregarded individuals with incomplete follow-up, as well as 7 individuals with relapse due to emergence of a resistant subclone. This left us with 46 individual treatment response curves that meet our criterion of biphasic tumor burden reduction. The Roeder and colelagues median data were used independently.

Our framework requires estimates of the offset of the tumor growth curve during treatment, ηt, and during initial growth, ηg. The former can be inferred from exponential fits on each patient-specific dataset directly. The latter requires further calculations. Typical datasets lack explicit information about the initial tumor growth. However, we can use |$\bar N(t) = {\eta _g}{\rm{exp\{\it T}}{{{s}}_{\rm{2}}}{\rm{\}}}$| and substitute for ηg in Eq. I:

formula

which corresponds to logistic growth. Here, T is the tumor age at diagnosis and s2 is the growth rate of the tumor in equilibrium before detection. For CML, it was estimated to take 5 to 7 years from the initiating mutational hit to diagnosis (1), supported by a peak in CML incidence approximately 6 years after of the atomic bomb in Hiroshima (35). However, we account for the uncertainty of tumor age and show that our results are robust for a wide range of tumor ages (see Supplementary Fig. S5). Furthermore, we needed to estimate the growth rates of the tumor during untreated growth, s1 and s2. First, we directly measured the initial fast reduction under treatment σ1, and the rate of tumor reduction in dynamic equilibrium, σ2, from treatment response curves. Second, the Michor and colleagues' data (1) contained 7 cases of relapse due to a resistant subclone. We used the slopes of these 7 relapses to estimate the initial fast growth s1 long before detection. Assuming similar effects of treatment on all cancer cells, we consider s1∝ σ1, s2∝ σ2, and |${s_2} \approx {s_1}\ {\rm{}}{\sigma _{\rm{2}}}/{\sigma _{\rm{1}}}$|⁠. These approximations likely overestimate the initial slope s1 and lead to a slight overestimate of TIC fractions. All individual regressions are shown in Supplementary Figs. S2 and S3, where σ1, σ2 correspond to the decline rates of tumor burden under treatment, and s1 is the growth rate of relapse. The respective distributions of these regression parameters are shown in Supplementary Fig. S4.

In our entire patient cohort, we found a median fraction of approximately 10−10 TICs at diagnosis. Tumor burden typically is in the order of approximately 1012 cells (1). This suggests a population of a few hundred to up to a few thousand TICs in CML at diagnosis. The same order of magnitude was also found by extensive sorting of patient-derived blood samples for quiescent (G0) BCR-ABL–positive cells, combined with in vitro and in vivo cell culture expansions (36).

Note that our method itself is independent of microscopic details of the hierarchy, for example, the proliferation parameters of cancer cells at different stages of the hierarchy. However, our estimate requires the tumor to be in a dynamic equilibrium. The transition time to this dynamic equilibrium depends on the microscopic properties of each individual cancer. It is therefore possible that we overestimated the fraction of TICs if the tumor trajectories used had not yet reached equilibrium at diagnosis. In addition, we idealized the tumor into two major compartments. A detailed description would require many more compartments (27). However, during late tumor growth and after sufficiently long continued therapy, the growth of hierarchically organized tumors is always dominated by the slowest, stem-like compartment.

We found variation within the patient population (Fig. 3). The treatment phase changed the composition of the tumor and we found a median TIC fraction of approximately 10−8 in equilibrium phase of continued treatment. Thus, treatment of one year increased the relative fraction of TICs by a factor of 100. This finding is robust to a wide variation of fitting parameters (see Supplementary Fig. S5). The tumor burden after one year of treatment is between 1% and 0.1% of the burden at diagnosis, suggesting a standing total number of a few hundred TICs. Slow decline of TICs during treatment is likely a consequence of the slow turnover rate inherited from hematopoietic stem cells (37), as opposed to an insensitivity of TICs to imatinib. The longer follow up of Roeder and colleagues showed a continued decrease of tumor burden to a total tumor population of 109 cells after 6 years. There could have been 10 to 100 remaining TICs in those patients, but the elimination of the last TICs is highly stochastic. Therefore, the extinction time is expected to have a wide distribution, causing significant differences between patients (27).

Our results explain why targeted treatment strategies are successful in CML and why patients stay in remission without emerging resistant subclones. This may be in contrast to targeted therapies of many solid tumors, such as in colon, lung, or breast, where resistance evolution or preexisting resistance is more common (38, 39).

The probability that no resistant mutant is present in a cell population of size N is pres = (1−μ)N−1, where μ is an effective rate that accounts for the emergence of resistant mutations (40). The cell population at risk (TICs) is small in CML. If we use our estimate of N ≈ 1,000 TICs, a point mutation rate of 10−9, 100 different mutations that can confer resistance, and a ratio of one symmetric tumor-initiating self-renewal in a hundred asymmetric stem cell divisions, we estimate that the probability of no resistance is pres = (1−10−5)1000−1 = 0.99. This suggests that only 1% of CML patients would have a preexisting resistant subclone within the population of tumor-initiating cells. Even if we assume 10,000 TICs at diagnosis, 90% of patients would not present with resistance in TICs, which is supported by clinical observations (33). This is opposed to solid tumors, where stem cell fractions might be in the order of 10−6 to 10−4 and preexisting resistance is highly likely, explaining the relapse of most CML patients after a long time of discontinued therapy: a considerable number of TICs survive therapy. However, most patients remain sensitive to targeted treatment and remission can often be reintroduced, because the TIC number is low overall.

The cancer stem cell hypothesis has attracted much attention, but also many critics and much skepticism since its formulation (41). While the existence of TICs in some tumors is well established, the situation in other cancers remains somewhat unclear (18, 31). However, knowledge about the size of the population of tumor initiating (or reinitiating) cells undoubtedly impacts cancer treatment approaches.

Numerous mathematical models have shed light on clinical phenomena from a theoretical perspective, and helped to explain patterns of treatment response and evolution of resistance (1, 28, 38, 42). Unfortunately, many models require involved parameterization, difficult to obtain in clinical settings. In this work, we presented a simple but general method to estimate the fraction of tumor-driving initiating cells with stem-like features. This estimate can be made exclusively from the shape of a tumor's growth or reduction curve. It only consists of a single exponential fit (linear regression) in the case of tumor growth, and two such regressions of longitudinal tumor size data during treatment. It is important to note that we do not provide a method to estimate model parameters by linear regression. Rather, we point to a direct functional link between two tumor properties, namely tumor expansion/reduction and the fraction of TICs. It is also important to point out that our method does not necessarily estimate fractions of tissue-specific stem cells. Instead, it estimates the cancer-driving reservoir of self-renewing cells, which can potentially be induced by oncogenic transformation at any stage of the hierarchy of the healthy tissue of origin (27, 28).

Our method of data analysis is parameter-free as it requires no knowledge about microscopic properties of the tumor. It only requires longitudinal measurements from techniques that are already used routinely in clinical care. For example, one can use high-resolution imaging or liquid biopsies. Thus, our method could readily complement current treatment protocols and inform clinicians about the relative size of the active pool of TICs. Patients with only few remaining TICs might benefit from treatment discontinuation, while patients with a remaining high number of TICs would benefit from changing the treatment strategy.

We found that the pool of self-renewing cells in CML is in the order of a few hundred cells, and thus relatively small. Furthermore, targeted treatment for 1 year increases the fraction of self-renewing cells by a factor of 100, although the total number of TICs likely declines together with the tumor bulk. This effect has also been observed in solid tumors experimentally, for example, in mammary tumors in mice, where platinum treatment increased the fraction of self-renewing cells 3-fold, compared with treatment-naïve cases (43).

The chance of a preexisting treatment–resistant tumor-initiating cell or the emergence of such resistance during therapy depends critically on the size of the TIC pool (39, 40) and thus can be assessed by our method. Our estimate of a relatively low number of a few hundred TICs in CML suggests a low probability of preexisting-resistant subclones. This is in agreement with clinical observations. Patients stay in remission on targeted treatment for years, in contrast to many solid tumors, where resistant subclones expand rapidly after treatment initiation. Possible late resistance cannot entirely be excluded, but its risk continuously decreases (44, 45).

Our model neglects a spatial component of tumor growth (46). This assumption leads to exponential growth in equilibrium, a situation well met in most leukemias (28, 47, 48). In some cases, the spatial component might be of importance and tumor growth becomes polynomial, rather than exponential (49). Another mechanism that can lead to subexponential growth of the tumor are positive or negative feedback loops due to signals secreted from differentiated cells that inhibit stem cell division and stem cell self-renewal (28, 50, 51). In such cases, our method provides only an approximation of the TIC fraction. Yet, the divergence might be small compared with the unavoidable measurement errors.

We provided a testable hypothesis concerning the fraction of TICs in both liquid and solid tumors. The ability to infer individual TIC fractions at diagnosis and during treatment can influence treatment strategies, as it informs about the tumor's self-renewal capability. Aggressiveness and duration of treatment might critically depend on the amount of TICs. Furthermore, the composition of drug cocktails, as well as treatment timing and scheduling, could be adjusted according to the knowledge gained about the TIC population.

No potential conflicts of interest were disclosed.

Conception and design: B. Werner, J.G. Scott, A. Sottoriva, P.M. Altrock

Development of methodology: B. Werner, J.G. Scott, P.M. Altrock

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): B. Werner, A. Traulsen, P.M. Altrock

Writing, review, and/or revision of the manuscript: B. Werner, J.G. Scott, A. Sottoriva, A.R.A. Anderson, A. Traulsen, P.M. Altrock

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P.M. Altrock

Study supervision: B. Werner, A. Sottoriva, A.R.A. Anderson, P.M. Altrock

This paper was conceived during the third annual Integrated Mathematical Oncology workshop on Personalized Medicine at the Moffitt Cancer Center, Tampa, FL.

B. Werner is supported by the Geoffrey W. Lewis Post-Doctoral Training Fellowship. J.G. Scott thanks the NIH for its generous loan repayment grant. A. Sottoriva is supported by The Chris Rokos Fellowship in Evolution and Cancer. This work was also supported by the Wellcome Trust (105104/Z/14/Z). P.M. Altrock acknowledges generous financial support from Deutsche Akademie der Naturforscher Leopoldina, grant no. LPDS 2012-12.

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.

1.
Hanahan
D
,
Weinberg
RA
. 
Hallmarks of cancer: the next generation
.
Cell
2011
;
144
:
646
74
.
2.
Consortium T1GP
. 
An integrated map of genetic variation from 1,092 human genomes
.
Nature
2012
;
490
:
56
65
.
3.
Vogelstein
B
,
Papadopoulus
N
,
Velculescu
VE
,
Zhou
S
,
Diaz
LA
 Jr
. 
Cancer genome landscapes
.
Science
2013
;
339
:
1546
58
.
4.
Burrell
RA
,
McGranahan
N
,
Bartek
J
,
Swanton
C
. 
The causes and consequences of genetic heterogeneity in cancer evolution
.
Nature
2013
;
501
:
338
45
.
5.
Greaves
M
,
Maley
CC
. 
Clonal evolution in cancer
.
Nature
2012
;
481
:
306
13
.
6.
Merlo
LMF
,
Pepper
JW
,
Reid
BJ
,
Maley
CC
. 
Cancer as an evolutionary and ecological process
.
Nat Rev Cancer
2006
;
6
:
924
35
.
7.
Byrne
HM
. 
Dissecting cancer through mathematics: from the cell to the animal model
.
Nat Rev Cancer
2010
;
10
:
221
30
.
8.
Marusyk
A
,
Tabassum
DP
,
Altrock
PM
,
Almendro
V
,
Michor
F
,
Polyak
K
. 
Non-cell-autonomous driving of tumour growth supports sub-clonal heterogeneity
.
Nature
2014
;
514
:
54
8
.
9.
Archetti
M
,
Ferraro
DA
,
Christofori
G
. 
Heterogeneity for IGF-II production maintained by public goods dynamics in neuroendocrine pancreatic cancer
.
Proc Natl Acad Sci
2015
;
112
:
1833
8
.
10.
Gerlee
P
,
Altrock
PM
. 
Complexity and stability in growing cancer cell populations
.
Proc Natl Acad Sci
2015
;
112
:
E2742
3
.
11.
Meacham
CE
,
Morrison
SJ
. 
Tumour heterogeneity and cancer cell plasticity
.
Nature
2013
;
501
:
328
37
.
12.
Bonnet
D
,
Dick
JE
. 
Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell
.
Nat Med
1997
;
3
:
730
7
.
13.
Reya
T
,
Morrison
SJ
,
Clarke
MF
,
Weissman
IL
. 
Stem cells, cancer, and cancer stem cells
.
Nature
2001
;
414
:
105
11
.
14.
O'Brien
CA
,
Pollett
A
,
Gallinger
S
,
Dick
JE
. 
A human colon cancer cell capable of initiating tumour growth in immunodeficient mice
.
Nature
2006
;
445
:
106
10
.
15.
Ricci-Vitiani
L
,
Pallini
R
,
Biffoni
M
,
Todaro
M
,
Invernici
G
,
Cenci
T
, et al
Tumour vascularization via endothelial differentiation of glioblastoma stem-like cells
.
Nature
2010
;
468
:
824
8
.
16.
Yeung
TM
,
Gandhi
SC
,
Wilding
JL
,
Muschel
R
,
Bodmer
WF
. 
Cancer stem cells from colorectal cancer-derived cell lines
.
Proc Natl Acad Sci
2010
;
107
:
3722
7
.
17.
Singh
S
,
Hawkins
C
,
Clarke
I
,
Squire
J
. 
Identification of human brain tumor initiating cells
.
Nature
2004
;
432
:
396
401
.
18.
Al-Hajj
M
,
Wicha
MS
,
Benito-Hernandez
A
,
Morrison
SJ
,
Clarke
MF
. 
Prospective identification of tumorigenic breast cancer cells
.
Proc Natl Acad Sci
2003
;
100
:
3983
8
.
19.
Prince
ME
,
Sivanandan
R
,
Wolf
GT
,
Kaplan
MJ
,
Dalerba
P
,
Weissman
IL
, et al
Identification of a subpopulation of cells with cancer stem cell properties in head and neck squamous cell carcinoma
.
Proc Natl Acad Sci
2007
;
104
:
973
8
.
20.
Eramo
A
,
Lotti
F
,
Sette
G
,
Pilozzi
E
,
Biffoni
M
,
Di Virgilio
A
, et al
Identification and expansion of the tumorigenic lung cancer stem cell population
.
Cell Death Differ
2007
;
15
:
504
14
.
21.
Dingli
D
,
Michor
F
. 
Successful therapy must eradicate cancer stem cells
.
Stem Cells
2006
;
24
:
2603
10
.
22.
Lenaerts
T
,
Pacheco
JM
,
Traulsen
A
,
Dingli
D
. 
Tyrosine kinase inhibitor therapy can cure chronic myeloid leukemia without hitting leukemic stem cells
.
Haematologica
2010
;
95
:
900
7
.
23.
Werner
B
,
Gallagher
RE
,
Paietta
EM
,
Litzow
MR
,
Tallman
MS
,
Wiernik
PH
, et al
Dynamics of leukemia stem-like cell extinction in acute promyelocytic leukemia
.
Cancer Res
2014
;
74
:
5386
96
.
24.
Michor
F
,
Hughes
TP
,
Iwasa
Y
,
Branford
S
,
Shah
NP
,
Sawyers
CL
, et al
Dynamics of chronic myeloid leukaemia
.
Nature
2005
;
435
:
1267
70
.
25.
Dingli
D
,
Traulsen
A
,
Pacheco
JM
. 
Compartmental architecture and dynamics of hematopoiesis
.
PLoS ONE
2007
;
4
:
e345
.
26.
Werner
B
,
Dingli
D
,
Lenaerts
T
,
Pacheco
JM
,
Traulsen
A
. 
Dynamics of mutant cells in hierarchical organized tissues
.
PLoS Comput Biol
2011
;
7
:
e1002290
.
27.
Rodriguez-Brenes
IA
,
Wodarz
D
,
Komarova
NL
. 
Minimizing the risk of cancer: tissue architecture and cellular replication limits
.
J R Soc Interface
2013
;
10
:
20130410
0
.
28.
Werner
B
,
Beier
F
,
Hummel
S
,
Balabanov
S
,
Lassay
L
,
Orlikowsky
T
, et al
Reconstructing the in vivo dynamics of hematopoietic stem cells from telomere length distributions
.
eLife
2015
;
10.7554
:
08687
.
29.
Werner
B
,
Dingli
D
,
Traulsen
A
. 
A deterministic model for the occurrence and dynamics of multiple mutations in hierarchically organized tissues
.
J R Soc Interface
2013
;
10
:
20130349
.
30.
Roeder
I
,
Horn
M
,
Glauche
I
,
Hochhaus
A
,
Mueller
MC
,
Loeffler
M
. 
Dynamic modeling of imatinib-treated chronic myeloid leukemia: functional insights and clinical implications
.
Nat Med
2006
;
12
:
1181
4
.
31.
Visvader
JE
,
Lindeman
GJ
. 
Cancer stem cells in solid tumours: accumulating evidence and unresolved questions
.
Nat Rev Cancer
2008
;
8
:
755
68
.
32.
Daley
G
,
Van Etten
R
,
Baltimore
D
. 
Induction of chronic myelogenous leukemia in mice by the P210bcr/abl gene of the Philadelphia chromosome
.
Science
1990
;
247
:
824
30
.
33.
Druker
BJ
,
Guilhot
F
,
O'Brien
SG
. 
Five-year follow-up of patients receiving imatinib for chronic myeloid leukemia
.
N Engl J Med
2006
;
355
:
2408
17
.
34.
Redaelli
S
,
Piazza
R
,
Rostagno
R
,
Magistroni
V
,
Perini
P
,
Marega
M
, et al
Activity of bosutinib, dasatinib, and nilotinib against 18 imatinib-resistant BCR/ABL mutants
.
J Clin Oncol
2008
;
27
:
468
9
.
35.
Bizzozero
JJO
,
Johnson
KG
,
Ciocco
A
,
Hoshino
T
,
Itoga
T
,
Toyoda
S
, et al
Radiation-related leukemia in hiroshima and nagasaki 1946-1964: distribution, incidence and appearance time
.
N Engl J Med
1966
;
274
:
1095
101
.
36.
Holyoake
TL
,
Jiang
X
,
Drummond
MW
,
Eaves
AC
,
Eaves
CJ
. 
Elucidating critical mechanisms of deregulated stem cell turnover in the chronic phase of chroni myeloid leukemia
.
Leukemia
2002
;
16
:
549
58
.
37.
Busch
K
,
Klapproth
K
,
Barile
M
,
Flossdorf
M
,
Holland-Letz
T
,
Schlenner
SM
, et al
Fundamental properties of unperturbed haematopoiesis from stem cells in vivo
.
Nature
2015
;
518
:
542
6
.
38.
Diaz
LA
,
Williams
RT
,
Wu
J
,
Kinde
I
,
Hecht
JR
,
Berlin
J
, et al
The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers
.
Nature
2012
;
486
:
537
40
.
39.
Bozic
I
,
Reiter
JG
,
Allen
B
,
Antal
T
,
Chatterjee
K
,
Shah
P
, et al
Evolutionary dynamics of cancer in response to targeted combination therapy
.
eLife
2013
;
2
:
e00747
.
40.
Kessler
DA
,
Levine
H
. 
Large population solution of the stochastic Luria-Delbruck evolution model
.
Proc Natl Acad Sci
2013
;
110
:
11682
7
.
41.
Wicha
MS
,
Liu
S
,
Dontu
G
. 
Cancer stem cells: an old idea–a paradigm shift
.
Cancer Res
2006
;
66
:
1883
90
.
42.
Komarova
NL
,
Wodarz
D
. 
Drug resistance in cancer: principles of emergence and prevention
.
Proc Natl Acad Sci
2005
;
102
:
9714
9
.
43.
Shafee
N
,
Smith
CR
,
Wei
S
,
Kim
Y
,
Mills
GB
,
Hortobagyi
GN
, et al
Cancer stem cells contribute to cisplatin resistance in brca1/p53-mediated mouse mammary tumors
.
Cancer Res
2008
;
68
:
3243
50
.
44.
Gugliotta
G
,
Castagnetti
F
,
Breccia
M
,
Levato
L
,
D'Adda
M
,
Stagno
F
, et al
Long-term outcome of a phase 2 trial with nilotinib 400 mg twice daily in first-line treatment of chronic myeloid leukemia
.
Haematologica
2015
;
100
:
1146
50
.
45.
Hochhaus
A
,
Brien
SGOA
,
Guilhot
F
,
Druker
BJ
,
Branford
S
,
Foroni
L
, et al
Six-year follow-up of patients receiving imatinib for the first-line treatment of chronic myeloid leukemia
.
Leukemia
2009
;
23
:
1054
61
.
46.
Talkington
A
,
Durett
R
. 
Estimating tumor growth rates in vivo
.
Bull Math Biol
2015
;
77
:
1934
54
.
47.
Lapidot
T
,
Sirard
C
,
Vormoor
J
,
Murdoch
B
,
Hoang
T
,
Caceres-Cortes
J
, et al
A cell initiating human acute myeloid leukamia after transplantation int SCID mice
.
Nature
1994
;
367
:
645
8
.
48.
Guibal
FC
,
Alberich-Jorda
M
,
Hirai
H
,
Ebralidze
A
,
Levantini
E
,
Di Ruscio
A
, et al
Identification of a myeloid committed progenitor as the cancer-initiating cell in acute promyelocytic leukemia
.
Blood
2009
;
114
:
5415
25
.
49.
Gerlee
P
. 
The model muddle: in search of tumor growth laws
.
Cancer Res
2013
;
73
:
2407
11
.
50.
Rodriguez-Brenes
IA
,
Komarova
NL
,
Wodarz
D
. 
Evolutionary dynamics of feedback escape and the development of stem-cell–driven cancers
.
Proc Natl Acad Sci
2013
;
108
:
18983
8
.
51.
Stiehl
T
,
Baran
N
,
Ho
AD
,
Marciniak-Czochra
A
. 
Cell division patterns in acute myeloid leukemia stem-like cells determine clinical course: a model to predict patient survival
.
Cancer Res
2015
;
75
:
940
9
.