The Cancer Stem Cell Fraction in Hierarchically Organized Tumors Can Be Estimated Using Mathematical Modeling and Patient-Specific Treatment Trajectories

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.

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.

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.

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 r_g 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

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.

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

See Fig. 2A.

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

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

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:

which corresponds to logistic growth. Here, T is the tumor age at diagnosis and s₂ 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, s₁ and s₂. First, we directly measured the initial fast reduction under treatment σ₁, and the rate of tumor reduction in dynamic equilibrium, σ₂, 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 s₁ long before detection. Assuming similar effects of treatment on all cancer cells, we consider s₁∝ σ₁, s₂∝ σ₂, and |${s_2} \approx {s_1}\ {\rm{}}{\sigma _{\rm{2}}}/{\sigma _{\rm{1}}}$|⁠. These approximations likely overestimate the initial slope s₁ and lead to a slight overestimate of TIC fractions. All individual regressions are shown in Supplementary Figs. S2 and S3, where σ₁, σ₂ correspond to the decline rates of tumor burden under treatment, and s₁ 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⁻¹⁰ TICs at diagnosis. Tumor burden typically is in the order of approximately 10¹² 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⁻⁸ 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 10⁹ 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 p_res = (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⁻⁹, 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 p_res = (1−10⁻⁵)¹⁰⁰⁰⁻¹ = 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⁻⁶ to 10⁻⁴ 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.