Abstract
Slowly cycling tumor cells that may be present in human tumors may evade cytotoxic therapies, which tend to be more efficient at destroying cells with faster growth rates. However, the proportion and growth rate of slowly cycling tumor cells is often unknown in preclinical model systems used for drug discovery. Here, we report a quantitative approach to quantitate slowly cycling malignant cells in solid tumors, using a wellestablished mouse model of Krasinduced lung cancer (Kras^{G12D/+}). 5Bromo2deoxyuridine (BrdUrd) was administered to tumorbearing mice, and samples were collected at defined times during pulse and chase phases. Mathematical and statistical modeling of the labelretention data during the chase phase supported the existence of a slowly cycling labelretaining population in this tumor model and permitted the estimation of its proportion and proliferation rate within a tumor. The doubling time of the slowly cycling population was estimated at approximately 5.7 weeks, and this population represented approximately 31% of the total tumor cells in this model system. The mathematical modeling techniques implemented here may be useful in other tumor models where direct observation of cellcycle kinetics is difficult and may help evaluate tumor cell subpopulations with distinct cellcycling rates. Cancer Res; 73(12); 3525–33. ©2013 AACR.
Major Findings
Mathematical modeling of labelretention data from a mouse model of lung cancer provides quantitative evidence for the presence of two populations of tumor cells with distinct cellcycling kinetics.
Quick Guide to Equations and Assumptions
Equations
The two mathematical models compared in the analysis of the chase data are the following:
Exponential model (1compartment model):
Where W is the natural logarithm of the percentage of BrdUrd^{+} cells within a tumor, and .
Biexponential model (2compartment model)
Where W is the natural logarithm of the percentage of BrdUrd^{+} cells within a tumor and .
The first equation uses an exponential decay function to model the percentage of tumor cells that are BrdUrd^{+}. The parameters A and α are non–negative constants, and ε is an error term from a normal distribution with mean zero and SD σ. The second equation models the same data as a biexponential decay function. The parameters A, α, B, and β are nonnegative constants, and ε is an error term from a normal distribution with mean zero and SD σ.
Parameter descriptions and interpretation

The parameter t denotes the time elapsed from time zero (after all BrdUrd administration ended) during chase phase.

The parameters α and β model how fast the cycling cells lose BrdUrd labeling.

The parameters A and B give the initial percentage of BrdUrd^{+} cells for each population.

The parameter ratio A/B gives the relative ratio of fast and slowcycling population in tumors.

At given time t during the chase phase, gives the observed percentage of slowcycling cells in BrdUrd^{+} cells.
Major Assumptions of the Model
Proliferation rates are proportional to the number of cells present;
Each tumor cell has a crosssectional area of approximately 121 μm^{2} (measurements of cell counts and areas of a random sample of 10 tumors gave a median value of 121 μm^{2}, with a SD of 25 μm^{2}; Supplementary Table S1);
Cell death rates are negligible for the time period of the experiment (we previously confirmed that cell death rates as determined by caspase3 staining are indeed quite low (<1%) in this setting) (ref. 1; also Supplementary Fig. S1, Supplementary Table S2);
The variance of the residuals of the logtransformed data is approximately constant;
Proportional error structures best characterize the residual errors for both models (based on tests of proportional and additive error structures);
During BrdUrd pulse, every cell that undergoes cell division takes in enough BrdUrd to be detectable after only one cell division;
The incorporation of BrdUrd has negligible effect on the survival of cells and their rate of cycling;
All cells require the same number of cell divisions to reach undetectable label levels (which can actually undercount cells that slowly cycle);
The rate of loss of BrdUrd label is proportional to the rate of proliferation of labeled cells;
All cells in the region labeled as tumor are actual tumor cells;
Tumors at the same time point but from different mice are comparable in age and size distributions;
Tumor measurements at the same time point are independent from one sample to another whether from the same mouse or from different mice;
Selection of tumor slices to sample is random; and
Tumor growth favors cells that cycle quickly over cells that cycle more slowly, in that any cells that cycle more slowly will become a smaller proportion of the total number of cells over time (for this analysis, the implication is that any population of slowly cycling cells that is detected at later time points was likely a larger proportion of all cells at earlier time points).
Introduction
Normal adult stem cells are thought to be relatively quiescent, a property that protects them from proliferative exhaustion (2). Because of this property, “labelretention” studies have been used to identify and characterize tissuespecific stem cells for decades, following the pioneering work by Potten and colleagues in the intestine (3). The existence of labelretaining cells has been proposed to be important for radiation response (3, 4). Labelretention approaches have also been used to identify stem cells in the interfollicular epidermis (5–9) and the hematopoietic system (10–12). In the hematopoietic stem cell (HSC) compartment, some studies have suggested the existence of a slowly cycling stem cell population (9, 10, 12), whereas other investigators have not found label retention in this compartment (11). In cancer research, increasing attention has focused on the heterogeneity of tumor cells present within the tumor mass (distinct from the heterogeneity of nontumor cells due to the presence of a tumor microenvironment) due to the hypothesis that certain subpopulations of tumor cells have increased capacity to propagate. These “cancer stem cells” (CSC) or “tumorinitiating cells (TIC) share with normal stem cells the ability to selfrenew and “differentiate” into committed cell types with more limited proliferative capacity (13). Cancer stem cell populations have been identified and extensively characterized in several solid tumors (14–17). However, it is not clear whether these CSC populations share the property of “label retention” that has been ascribed to some normal stem cells. Recent data in glioblastoma and in skin cancer suggest that indeed such slowly cycling CSCs may exist (18, 19). However, whether this is true for other tumor types remains to be determined. A limitation to carrying out a rigorous analysis of the cellcycle kinetics of solid tumor populations is the lack of robust and rigorous quantitative methods.
CSCs in solid tumors are traditionally identified by differential cell sorting with specific cell surface markers (20). The operational “stemness” property is defined as an increased ability to transplant tumors into an immunocompromised host. However, this definition of CSCs remains controversial, as transplantation success rates can vary widely depending on a wide range of technical factors (21). In human lung cancer, several reports have suggested that a CSC population can be enriched using cell surface markers, yet the cellcycle kinetics of these subpopulations (or indeed, whether slowly cycling cells are synonymous with cancer stem cells defined by transplantation assays) has not been explored (22, 23).
The Kras^{G12D/+} mouse model has been well validated for the study of lung cancer initiation and progression in vivo (24). Previous work has shown that this model closely recapitulates the molecular changes seen in human non–small cell lung cancer (NSCLC; ref. 25). In this model, an oncogenic Kras allele is expressed only after delivery of Cre recombinase. Cre recombinase is delivered via nasal instillation using an adenoviral vector (AdCre). Thus, timing of tumor initiation can be well controlled. After Cre delivery, the lung epithelium begins to proliferate and over time the lung parenchyma is occupied by multiple adenomas, some of which then progress to adenocarcinoma. In the studies described here, BrdUrd incorporation was measured in tumors within the lungs in this Kras^{G12D/+} model. Previous work using this model identified a population of CD45^{−}PECAM^{−}Sca1^{+} cells enriched for doublepositive SPC/CC10 cells [bronchioalveolar stem cells (BASC)] that are activated after lung injury and thus thought to be CSCs (26). However, this population was not enriched for transplantation ability into immunocompromised mice in the Kras^{G12D/+} model (27). The fact that only a small percentage of cells in the Kras^{G12D/+} model can transplant tumors suggests that there is a subpopulation of cells with CSC characteristics.
“Label retention” of slowly cycling cells provides an alternative method for fluorescenceactivated cell sorting of surface markers to identify populations of cells within tumors with important biological properties that may bypass the limitations of the transplantation approach. An advantage of labelretention studies is that they allow for rigorous testing of the hypothesis that a slowly cycling population exists without the need for complex crosses with lineagetracing reporters and other genetic tools (18, 19, 28). A priori, it cannot be assumed that CSCs are label retaining. However, if a labelretaining population can be identified by labelretention studies, this could justify further experimentation to allow for isolation of these cells. In addition, it may be important in some settings to determine whether the CSC population is synonymous with the slowly cycling population or whether these are overlapping but distinct entities. Recent data isolating and characterizing a labelretaining population in breast tissue supports this approach (29).
A commonly used method to identify labelretaining cells is to mark proliferating cells using the thymidine analogue 5bromo2deoxyuridine (BrdUrd). After administration of BrdUrd ends, cycling cells will lose this DNA label as it will be diluted by half of each complete cell cycle. Therefore, both BrdUrd uptake (“pulse”) and BrdUrd dilution (“chase”) data are useful in quantifying cell turnover rates. A key element of such an approach is to develop appropriate statistical and mathematical modeling techniques to show the presence of “labelretaining” cells. Mathematical modeling of BrdUrd label incorporation and retention can quantify the behavior of the labeled cells and enable identification of populations with different cellcycle kinetics. Rigorous mathematical approaches have been developed that use BrdUrd label retention to determine cellcycle kinetics in other settings. For example, Bonhoeffer and colleagues used a mathematical model of BrdUrd incorporation to describe the kinetics of CD4^{+} and CD8^{+} lymphocytes circulating in uninfected and simian immunodeficiency virus–infected macaques (30). The model was used to estimate proliferation and death rates of these specific lymphocyte subsets. Similarly, Kiel and colleagues modeled BrdUrd incorporation and retention data from mice to test the “immortal strand” hypothesis of asymmetric chromosome segregation during division of HSCs and concluded that asymmetric segregation of chromosomes does not occur in the division of HSCs (11). Here, we use mathematical and statistical analysis of in vivo BrdUrd labeling data to determine whether tumors in a mouse model of lung cancer contain a population of slowly cycling tumor cells. Our results strongly suggest that a population of slowly cycling tumor cells does exist in the Kras^{G12D/+} model. These studies establish the necessary justification for further analysis of this slowly cycling population in this model. More generally, we propose that the mathematical analysis of BrdUrd labeling in tumor models applied here may be broadly useful for other similar tumor models to elucidate the heterogeneity of tumor cell kinetics.
Materials and Methods
Mice
For all studies, the Kras^{G12D/+} mouse model was used. Mice were genotyped as previously described (24). AdCre was delivered by intranasal instillation between 4 and 6 weeks of age. Mice were then allowed to age for 12 weeks before BrdUrd administration. All animal experiments were approved by the Stanford University School of Medicine Administrative Panel on Laboratory Animal Care (APLAC; Stanford University, Stanford, CA).
BrdUrd labeling
BrdUrd (SigmaAldrich) was administered to the mice for up to 15 days for the pulse experiment and 7 days for the chase experiment. Administration was done simultaneously both through drinking water (1 mg/mL) and daily intraperitoneal injections (100 mg per kilogram body weight) to maximize delivery.
Immunohistochemistry and immunofluorescence
After sacrifice, mouse lungs were removed, fixed with Zfix (Anatech Inc), and embedded in paraffin. Immunohistochemical or immunofluorescence staining was conducted on 4μm sections. For primary antibodies, antiBrdUrd (BD Pharmigen clone 3d4, 1:200 dilution), antiTTF1 (Abcam, 1:200 dilution), antiEcadherin (BD Transduction Laboratories, clone 36, 1:200 dilution), and anticleaved caspase3 (Cell Signaling, 1:300 dilution) were used. Detection was conducted with goat antirat Alexa Fluor 594 (Invitrogen), donkey antimouse Alexa Fluor 488 (Invitrogen), and donkey antirabbit Alexa Fluor 488 (Invitrogen) for immunofluorescence, or horseradish peroxidaseconjugated goat antimouse and ABC goatantirabbit kit for immunohistochemistry.
Quantification of immunohistochemistry
The area of each tumor crosssection was measured using Bioquant software, and the number of BrdUrdpositive (BrdUrd^{+}) cells within each tumor was assessed by manual counting. The number of BrdUrd^{+} cells per tumor crosssectional area was recorded. A manual count of 10 randomly selected tumors gave an estimate of 121 μm^{2} area per tumor cell (SD 24.9 μm^{2}). This estimate was used to convert BrdUrd^{+} cells per area into units of BrdUrd^{+} cells per total tumor cells. Two mice at each time point were used to generate the pulse data. At each time point in the chase period, approximately 5 or 6 mice were selected at random (one time point had 7 mice, and another had 3). For each selected mouse, multiple slices of the tumorbearing lung were obtained and placed at random positions on slides. Approximately 5 tumors were selected randomly from the available slices for each mouse for analysis (1 mouse with 1 tumor, 1 mouse with 2 tumors, 9 mice with 4 tumors, and 39 mice with 5 tumors).
Mathematical modeling
At the end of the 12 weeks of the chase, labeled cells were still detected. This suggested the possibility of a labelretaining subpopulation with a slower cycling time than other cells. To test this hypothesis, we fit 2 mathematical models to the chase data and compared the goodnessoffit of the models. The first was an exponential decay model, which assumes that there is a single cellcycling rate. The second was a biexponential decay model, which assumes that there are 2 distinct cellcycling rates. Each model included an error structure that was selected to give the best fit for that particular model to account for residual variability. Different error structures were tested on each model, and a 2sided Kolmogorov–Smirnov test was conducted to confirm that a proportional error structure resulted in residuals that were indistinguishable from a normal distribution for each of the models.
Because of the random selection of the slices and tumors for analysis (described above), we assume that the %BrdUrd values from a single mouse form a random sample. This is also supported by Supplementary Fig. S2, which shows chase data plotted by mouse. Supplementary Fig. S2 also provides graphical evidence that at each time point, the distributions of the %BrdUrd^{+} values are similar between mice. In addition, the Kruskal–Wallis test, an analogue of the ANOVA test for non–normally distributed data, was conducted to compare the median %BrdUrd^{+} values for all of the mice at a given time point. With a significance level of 0.05, 5 of the 9 time points showed no significant difference between those median values. For the remaining 4 time points, there was at least one mouse whose median differed from those of the other mice (although only one of those time points would have shown a difference at a significance level of 0.01). To further investigate possible differences between mice, we used the Wilcoxon rank sum test, a nonparametric analogue of the t test, to compare the median %BrdUrd^{+} values for each pair of mice at a time point. A Bonferroni correction was used to account for the multiple pairwise tests conducted at each time point. For 118 of the 119 possible comparisons, the median values for mice at a given time point were not pairwise different at a significance level of 0.05 (for the other comparison, the P value was 0.045). On the basis of all of these results, we assume that all %BrdUrd^{+} values at a given time point form a random sample.
Goodnessoffit testing
Both mathematical models were fit to the data using nonlinear regression. In both cases, a natural logarithm transformation of the data was used to account for the heteroscedasticity (nonuniform variance) of the data. Additive and proportional error structures were then tested on the models. For both models, proportional error structures gave the best fit and were used throughout the analysis. Goodnessoffit comparisons of the 2 models were made using mean square error (MSE) and Akaike information criterion (AIC) values.
Results
Determination of minimal sufficient length of labeling phase
To determine how long the BrdUrd labeling period should be for the chase experiment, an extended pulse experiment was conducted. Kras^{G12D/+} mice were treated with AdCre to initiate expression of oncogenic Kras. Twelve weeks after introduction of AdCre, when tumors were well established, mice were given BrdUrd to label proliferating tumor cells. To determine the minimal sufficient length of the labeling phase, an initial “pulse” of BrdUrd was delivered to mice for up to 15 days (see Materials and Methods). Mice were sacrificed on days 2, 3, 5, 7, 9, 11, 13, and 15 during the labeling period and BrdUrd label incorporation was assessed by immunohistochemistry using an antiBrdUrd antibody in paraffinfixed sections of lung. As shown in Fig. 1, the cell population was effectively fully labeled by day 7 (chase time “0 week”), with a median day 7 level of 29.4%. Therefore, for the chase data collection, a pulse phase of 7 days was used for the labeling of the mice.
Detection of BrdUrd^{+} cells through a pulsechase experiment
For the chase experiment, BrdUrd was administered over 7 days to a second cohort of Kras^{G12D/+} mice 12 weeks after tumor initiation (n = 50). After 7 days of BrdUrd “pulse,” mice were sacrificed at weeks 0, 1, 2, 3, 4, 6, 8, 10, and 12 (Fig. 2). Four to 6 mice from this cohort were sacrificed at each specified time, and the remaining number of BrdUrd^{+} cells was determined again by immunohistochemistry (Fig. 2A). To confirm that BrdUrd labelretaining cells, and particularly those present at the end of the chase experiment, were indeed lung tumor cells and not cells from the surrounding stroma, we used double labeling with immunofluorescence with an antiBrdUrd antibody and an antibody against lung epithelial marker Ecadherin or TTF1 (Fig. 2B). The results show that the vast majority of BrdUrd^{+} cells are indeed of epithelial origin and therefore are likely to be lung tumor cells.
Mathematical modeling supports the existence of a distinct, slowcycling population
At the end of the 12 weeks of the chase, labeled tumor cells were still detected (Fig. 2A, 12 weeks; and Fig. 2B). This suggested the possibility of a labelretaining tumor cell subpopulation with a slower cycling time than that of other tumor cells. To test this hypothesis, we fit 2 mathematical models to the chase data and compared the goodnessoffit of the models (Fig. 2C and Table 1). The first was an exponential decay model, which assumes that there is a single cellcycling rate (Eq. A). The second was a biexponential decay model, which assumes that there are 2 distinct cellcycling rates (Eq. B). Each model included an error structure (as described in Materials and Methods) that was selected to give the best fit for that particular model to account for residual variability. Different error structures were tested on each model, and a 2sided Kolmogorov–Smirnov test was conducted to confirm that a proportional error structure resulted in residuals that were indistinguishable from a normal distribution for each of the models.
We considered a few refinements for these 2 models but finally used the simple models described above:
One refinement considered was allowing cells to take differing numbers of cell divisions to reach an undetectable level of label. Kiel and colleagues (11) estimated that approximately 3 rounds of cell division are needed to decrease the level of BrdUrd present in a cell below the threshold of detection after full labeling. Bonhoeffer and colleagues (30) estimated that the number of cell divisions needed to have the BrdUrd level drop below the lower limit of detection to be 5 to 6 divisions. We assumed that all cells require the same number of divisions during the chase phase to reach an undetectable level of label, which allowed us to use the simple models. This is a reasonable assumption if all cells begin the chase phase with equally high levels of label. However, if there are slower cycling cells, they could begin the chase phase with lower levels of label than other cells owing to fewer cell divisions during the labeling period. Because BrdUrd is taken up into DNA during the Sphase of cell division, the number of times a cell divides during the BrdUrd pulse phase affects the amount of BrdUrd incorporated at the end of the pulse phase. Cells that divide more rapidly can take up more label per time period than cells that divide more slowly (31). If the label period is not sufficiently long, slowercycling cells may require fewer divisions for the label to become undetectable. This could lead to an undercounting of slowercycling cells at later times during the chase phase, if such cells exist.
In contrast with the situation in normal tissues, the size of the compartment for solid tumors changes over time. To account for this, we normalized counts of BrdUrd^{+} cells by an areabased estimate of the total number of cells in each tumor sample. We did not include cell death in our mathematical models because tumors in the Kras^{G12D/+} model are highly proliferative and yet show very few apoptotic cells. Cell death rates as determined by caspase3 staining are indeed quite low (<1%) in this setting (1). Therefore cell death was not included in the mathematical models.
We also fit exponential and biexponential models to the pulse data and compared the goodnessoffit of the models. However, the data were insufficient to be able to distinguish one model as significantly better than the other. Therefore we did not include a model for the pulse data in our results.
MSE and AIC were used to compare the models. Both MSE and AIC are measures of goodnessoffit, with a smaller value indicating a better fit. However, AIC penalizes for extra parameters, taking into account that a model with more parameters has more degrees of freedom and will give a better fit to data for this reason.
As shown in Table 1, the MSE values were 0.476 and 0.456 for the exponential and biexponential models, respectively. The AIC values were 496.46 and 490.25 for the exponential and biexponential models, respectively. Therefore the mathematical model that assumed 2 distinct cycling periods gave a measurably superior fit to the data than the model with a single cycling period using either criterion. Both models incorporated variability, as described in Materials and Methods, so the difference was not simply due to random variability favoring a model with multiple cycling periods. In addition, a Lowess fit was made for the residuals with a proportional error structure for each model, and visual comparison confirmed that the residuals for the biexponential model were more uniformly distributed than the residuals for the exponential model (Fig. 3A and B).
Crossvalidation supports that the 2compartment model is a superior model
To further verify that the biexponential model was indeed a better fit to the data, a number of crossvalidation tests were conducted. The crossvalidations give a measure of the sensitivity of the models to the data by testing whether the models overfit the data. For each crossvalidation, a subset of data (the validation data) was removed from the full data set. The remaining data (the training data) were fit with the 2 models described in the Quick Guide to Equations and Assumptions. The structure of the training data sets for the various crossvalidations is described below.
Crossvalidation 1 (omit 1 mouse).
A training data set was formed by removing all data taken from a single mouse. This procedure was repeated so that each mouse was left out once, for a total of 50 cases. The fitting algorithm for the biexponential model did not converge in one of these cases, so there were only 49 cases in which the 2 models could be compared.
Crossvalidation 2 (omit 1 mouse per time point).
A training set of data was formed by removing all data from one mouse at each time point. This was done 1,000 times, with a random selection of the mice.
Crossvalidation 3 (omit 1 tumor per mouse).
Each training set of data was formed by randomly removing all of the data from one tumor from each mouse. This was done 1,000 times.
Crossvalidation 4 (omit 2 tumors per mouse).
Training sets were formed by removing 2 samples at random from the data from each mouse. This was done 1,000 times.
Crossvalidation 5 (omit 20% of data per time point).
Training sets were formed by removing approximately 20% of the samples at each time point at random. This was done 1,000 times.
The MSE and the AIC were computed for each of the models when they were fit to a training data set. The percentages of times the biexponential model gave a superior fit (as measured by lower MSE or AIC) are recorded in Table 2. For all 5 types of crossvalidations conducted, the biexponential model gave a better fit than the exponential model as determined by both MSE and AIC. The MSE was lower for the biexponential model in 100% of cases. The AIC was lower for the biexponential model in more than 75% of cases for each of the 5 crossvalidations.
For each of the fitted models derived from the training data, the MSE was computed for the validation data. The percentages of cases in which the MSE value was better (lower) for the biexponential function than the exponential function are recorded in Table 2. In each of the crossvalidation methods conducted, the biexponential model had lower MSE values than the exponential model for the majority of cases (more than 50%), and was therefore concluded to be the superior model.
Predictive checks supported selection of the biexponential model
As further model validation, 2 predictive checks were conducted for the models. The bestfitting exponential and biexponential models to the full data set were used to produce simulated data.
Predictive check 1.
The same number of data points at each time point as contained in the original data set were simulated for each time point, and this was done 1,000 times to give 1,000 full datasets of simulated data for each of the 2 models.
Predictive check 2.
A total of 500 data points were simulated at each time point. This was repeated 1,000 times to give 1,000 simulated data sets (each with 500 data points at each time) for each of the 2 models.
For each predictive check, the distributions of the simulated and real data at each time point were compared. The 2sided Kolmogorov–Smirnov test was applied to determine statistical significance for differences between the distributions, with the null hypothesis being that any 2 distributions being compared cannot be distinguished (with a significance level of 0.05). Table 3 summarizes the results of the tests, with the entries giving the percentages of cases (out of 1,000 simulations) in which the null hypothesis was rejected and the 2 distributions being compared were distinguishable according to the Kolmogorov–Smirnov test. In addition to the comparisons of simulated and real data, comparisons were made between simulated data sets for the 2 different models. Those results are included in the tables as well. In both predictive checks, for the majority of the points, the biexponential model was indistinguishable from the data more often than the exponential model was. This was true for 6 of 9 time points in the first test (Table 3, top) and 8 of 9 time points in the second test (Table 3, bottom).
Predicted values
The parameter estimates for the 2compartment model in Table 1 indicate that approximately 31% of the cells were initially in the slow population and that the mean turnover rates for the fast and slow populations were 1.8 and 5.7 weeks per cell, respectively. Figure 4 shows the biexponential function prediction of the slowly cycling population over the course of the chase experiment. It was created by first bootstrapping the data sets to create 1,000 new datasets and fitting each data set to get a new set of parameter estimates. Those 1,000 sets of parameter estimates were then used to calculate the fraction of slowly cycling cells at each experimental time point. The 2.5, 50, and 97.5 percentiles of those values were plotted.
Discussion
We conducted mathematical modeling and statistical analyses of BrdUrd labelretention data in tumors using a mouse model of human lung cancer. In every test conducted, a model with 2 distinct cycling rates was found to fit the data better than a model with a single cycling rate. Multiple methods of crossvalidation and predictive checks confirmed this conclusion. Residuals for both models were fit by Lowess curves (Fig. 3A and B), visually showing a more uniform distribution of residuals for the biexponential model. The modeling and analysis provide quantitative evidence for the existence of a labelretaining population of cells and thus indicate that there are likely 2 distinct proliferation rates in this solid tumor model. Further experimentation will be required to determine whether these labelretaining cells are indeed tumorinitiating CSCs.
Conclusions
Mathematical modeling was applied to BrdUrdlabeled tumor cells in a mouse model of lung cancer. Measures of goodnessoffit for the competing mathematical models showed that 2 distinct proliferation rates are more likely than a single proliferation rate, providing quantitative evidence of a population of “labelretaining” tumor cells. The mathematical modeling also allowed us to estimate the proliferation rate of this labelretaining population as well as to estimate the proliferation rate of the fastercycling, non–labelretaining population. In addition, these estimated rates allow for the determination of the fraction of tumor cells that are label retaining at various times during the chase phase. This provides an upper bound on the fraction of tumor cells in the samples that may have stem cell–like properties and an estimate of their proliferation rate. The parameter estimates for the 2compartment model in Table 1 indicate that approximately 31% of the cells were initially in the slow population and that the mean turnover rates for the fast and slow populations were 1.8 and 5.7 weeks per cell, respectively.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: Y. Zheng, H. Moore, E.A. SweetCordero
Development of methodology: Y. Zheng, H. Moore, E.A. SweetCordero
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Zheng, T. Solis
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Zheng, H. Moore, A. Piryatinska, E.A. SweetCordero
Writing, review, and/or revision of the manuscript: Y. Zheng, H. Moore, E.A. SweetCordero
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Zheng, E.A. SweetCordero
Study supervision: Y. Zheng, E.A. SweetCordero
Grant Support
E.A. SweetCordero was funded by a Research Scholar Grant from the American Cancer Society and by NCI 5 R01 CA15751003. Y. Zheng was supported by a Senior Research Fellowship from the American Lung Association (RT82480N)
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.
Footnotes
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
 Received November 24, 2012.
 Revision received February 17, 2013.
 Accepted March 15, 2013.
 ©2013 American Association for Cancer Research.