## 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 well-established mouse model of Kras-induced lung cancer (*Kras ^{G12D/+}*). 5-Bromo-2-deoxyuridine (BrdUrd) was administered to tumor-bearing mice, and samples were collected at defined times during pulse and chase phases. Mathematical and statistical modeling of the label-retention data during the chase phase supported the existence of a slowly cycling label-retaining 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 cell-cycle kinetics is difficult and may help evaluate tumor cell subpopulations with distinct cell-cycling rates.

*Cancer Res; 73(12); 3525–33. ©2013 AACR*.

Mathematical modeling of label-retention data from a mouse model of lung cancer provides quantitative evidence for the presence of two populations of tumor cells with distinct cell-cycling kinetics.

**Equations**

The two mathematical models compared in the analysis of the chase data are the following:

Where

*W*is the natural logarithm of the percentage of BrdUrd^{+}cells within a tumor, and .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 non-negative 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 slow-cycling population in tumors.

At given time

*t*during the chase phase, gives the observed percentage of slow-cycling cells in BrdUrd^{+}cells.

Major Assumptions of the Model

Proliferation rates are proportional to the number of cells present;

Each tumor cell has a cross-sectional 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 caspase-3 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 log-transformed 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, “label-retention” studies have been used to identify and characterize tissue-specific stem cells for decades, following the pioneering work by Potten and colleagues in the intestine (3). The existence of label-retaining cells has been proposed to be important for radiation response (3, 4). Label-retention 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 “tumor-initiating cells (TIC) share with normal stem cells the ability to self-renew 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 cell-cycle 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 cell-cycle 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*model. Previous work using this model identified a population of CD45

^{G12D/+}^{−}PECAM

^{−}Sca1

^{+}cells enriched for double-positive 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*model (27). The fact that only a small percentage of cells in the

^{G12D/+}*Kras*model can transplant tumors suggests that there is a subpopulation of cells with CSC characteristics.

^{G12D/+}“Label retention” of slowly cycling cells provides an alternative method for fluorescence-activated 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 label-retention studies is that they allow for rigorous testing of the hypothesis that a slowly cycling population exists without the need for complex crosses with lineage-tracing reporters and other genetic tools (18, 19, 28). *A priori*, it cannot be assumed that CSCs are label retaining. However, if a label-retaining population can be identified by label-retention 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 label-retaining population in breast tissue supports this approach (29).

A commonly used method to identify label-retaining cells is to mark proliferating cells using the thymidine analogue 5-bromo-2-deoxyuridine (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 “label-retaining” cells. Mathematical modeling of BrdUrd label incorporation and retention can quantify the behavior of the labeled cells and enable identification of populations with different cell-cycle kinetics. Rigorous mathematical approaches have been developed that use BrdUrd label retention to determine cell-cycle 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 (Sigma-Aldrich) 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 Z-fix (Anatech Inc), and embedded in paraffin. Immunohistochemical or immunofluorescence staining was conducted on 4-μm sections. For primary antibodies, anti-BrdUrd (BD Pharmigen clone 3d4, 1:200 dilution), anti-TTF1 (Abcam, 1:200 dilution), anti-E-cadherin (BD Transduction Laboratories, clone 36, 1:200 dilution), and anti-cleaved caspase-3 (Cell Signaling, 1:300 dilution) were used. Detection was conducted with goat anti-rat Alexa Fluor 594 (Invitrogen), donkey anti-mouse Alexa Fluor 488 (Invitrogen), and donkey anti-rabbit Alexa Fluor 488 (Invitrogen) for immunofluorescence, or horseradish peroxidase-conjugated goat anti-mouse and ABC goat-anti-rabbit kit for immunohistochemistry.

### Quantification of immunohistochemistry

The area of each tumor cross-section was measured using Bioquant software, and the number of BrdUrd-positive (BrdUrd^{+}) cells within each tumor was assessed by manual counting. The number of BrdUrd^{+} cells per tumor cross-sectional 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 tumor-bearing 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 label-retaining 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 goodness-of-fit of the models. The first was an exponential decay model, which assumes that there is a single cell-cycling rate. The second was a biexponential decay model, which assumes that there are 2 distinct cell-cycling 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 2-sided 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.

### Goodness-of-fit 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. Goodness-of-fit 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 anti-BrdUrd antibody in paraffin-fixed 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 pulse-chase 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 label-retaining 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 anti-BrdUrd antibody and an antibody against lung epithelial marker E-cadherin or TTF-1 (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, slow-cycling 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 label-retaining 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 goodness-of-fit of the models (Fig. 2C and Table 1). The first was an exponential decay model, which assumes that there is a single cell-cycling rate (Eq. A). The second was a biexponential decay model, which assumes that there are 2 distinct cell-cycling 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 2-sided 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.

. | A (%BrdUrd^{+} cells)
. | α (1/week) . | B (%BrdUrd^{+} cells)
. | β (1/week) . | MSE . | AIC . |
---|---|---|---|---|---|---|

Exponential model | 23.4 (7.2%) | 0.269 (4.4%) | — | — | 0.476 | 496.46 |

Biexponential model | 20.3 (43%) | 0.549 (49%) | 9.14 (110%) | 0.176 (50%) | 0.456 | 490.25 |

. | A (%BrdUrd^{+} cells)
. | α (1/week) . | B (%BrdUrd^{+} cells)
. | β (1/week) . | MSE . | AIC . |
---|---|---|---|---|---|---|

Exponential model | 23.4 (7.2%) | 0.269 (4.4%) | — | — | 0.476 | 496.46 |

Biexponential model | 20.3 (43%) | 0.549 (49%) | 9.14 (110%) | 0.176 (50%) | 0.456 | 490.25 |

NOTE: Percent relative standard error (%RSE) for each parameter estimate is included in parentheses on the left. Measures of “goodness-of-fit” on the right are derived from either MSE or AIC method.

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 S-phase 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, slower-cycling cells may require fewer divisions for the label to become undetectable. This could lead to an undercounting of slower-cycling 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 area-based 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*model are highly proliferative and yet show very few apoptotic cells. Cell death rates as determined by caspase-3 staining are indeed quite low (<1%) in this setting (1). Therefore cell death was not included in the mathematical models.^{G12D/+}We also fit exponential and biexponential models to the pulse data and compared the goodness-of-fit 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 goodness-of-fit, 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).

### Cross-validation supports that the 2-compartment model is a superior model

To further verify that the biexponential model was indeed a better fit to the data, a number of cross-validation tests were conducted. The cross-validations give a measure of the sensitivity of the models to the data by testing whether the models overfit the data. For each cross-validation, 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 cross-validations is described below.

#### Cross-validation 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.

#### Cross-validation 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.

#### Cross-validation 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.

#### Cross-validation 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.

#### Cross-validation 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 cross-validations 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 cross-validations.

Percentage of training data sets for which the biexponential model fit was better . | ||
---|---|---|

Cross-validation . | According to MSE . | According to AIC . |

1 | 100% | 100% |

2 | 100% | 98.2% |

3 | 100% | 97.9% |

4 | 100% | 75.3% |

5 | 100% | 100% |

Percentage of validation data sets for which the biexponential model fit was better | ||

Cross-validation | According to MSE | |

1 | 53% | |

2 | 79.2% | |

3 | 80.9% | |

4 | 86.5% | |

5 | 100% |

Percentage of training data sets for which the biexponential model fit was better . | ||
---|---|---|

Cross-validation . | According to MSE . | According to AIC . |

1 | 100% | 100% |

2 | 100% | 98.2% |

3 | 100% | 97.9% |

4 | 100% | 75.3% |

5 | 100% | 100% |

Percentage of validation data sets for which the biexponential model fit was better | ||

Cross-validation | According to MSE | |

1 | 53% | |

2 | 79.2% | |

3 | 80.9% | |

4 | 86.5% | |

5 | 100% |

NOTE: The top part of the table shows the percentage of cases in which the biexponential model was superior to the exponential model, when both were fit to the training data sets. The bottom part of the table shows the percentage of cases in which the biexponential model was superior to the exponential model, as determined by lower value of MSE on the testing data set for each model determined by fitting to the validation data sets.

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 cross-validation 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 best-fitting 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 2-sided 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).

Time (wks) . | 0 . | 1 . | 2 . | 3 . | 4 . | 6 . | 8 . | 10 . | 12 . |
---|---|---|---|---|---|---|---|---|---|

Predictive check 1 . | |||||||||

Data vs. log exponential | 32.3 | 4.9 | 2.9 | 2.7 | 8.4 | 7.1 | 6.9 | 0.2 | 42 |

Data vs. log biexponential | 5.9 | 2.9 | 3.3 | 2.6 | 5.0 | 0.8 | 14.4 | 0.8 | 12.7 |

Log exp vs. log biexponential | 11.8 | 4.7 | 5.6 | 6.2 | 8.7 | 10.1 | 5.3 | 3.6 | 11.9 |

Predictive check 2 | |||||||||

Data vs. log exponential | 98.8 | 0.1 | 0 | 0 | 0.6 | 0.8 | 0 | 0 | 100 |

Data vs. log biexponential | 0 | 0 | 0 | 0 | 0.4 | 0 | 2.3 | 0 | 15.4 |

Log exp vs. log biexponential | 99.5 | 36.3 | 10.8 | 61.9 | 89.2 | 92.2 | 40.3 | 14.2 | 97.5 |

Time (wks) . | 0 . | 1 . | 2 . | 3 . | 4 . | 6 . | 8 . | 10 . | 12 . |
---|---|---|---|---|---|---|---|---|---|

Predictive check 1 . | |||||||||

Data vs. log exponential | 32.3 | 4.9 | 2.9 | 2.7 | 8.4 | 7.1 | 6.9 | 0.2 | 42 |

Data vs. log biexponential | 5.9 | 2.9 | 3.3 | 2.6 | 5.0 | 0.8 | 14.4 | 0.8 | 12.7 |

Log exp vs. log biexponential | 11.8 | 4.7 | 5.6 | 6.2 | 8.7 | 10.1 | 5.3 | 3.6 | 11.9 |

Predictive check 2 | |||||||||

Data vs. log exponential | 98.8 | 0.1 | 0 | 0 | 0.6 | 0.8 | 0 | 0 | 100 |

Data vs. log biexponential | 0 | 0 | 0 | 0 | 0.4 | 0 | 2.3 | 0 | 15.4 |

Log exp vs. log biexponential | 99.5 | 36.3 | 10.8 | 61.9 | 89.2 | 92.2 | 40.3 | 14.2 | 97.5 |

NOTE: Predictive check 1 uses a sample size equal to the number of samples in the original data set at each time point. A sample size of 500 at each time point was used for predictive check 2.

### Predicted values

The parameter estimates for the 2-compartment 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 label-retention 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 cross-validation 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 label-retaining 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 label-retaining cells are indeed tumor-initiating CSCs.

## Conclusions

Mathematical modeling was applied to BrdUrd-labeled tumor cells in a mouse model of lung cancer. Measures of goodness-of-fit 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 “label-retaining” tumor cells. The mathematical modeling also allowed us to estimate the proliferation rate of this label-retaining population as well as to estimate the proliferation rate of the faster-cycling, non–label-retaining 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 2-compartment 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. Sweet-Cordero

**Development of methodology:** Y. Zheng, H. Moore, E.A. Sweet-Cordero

**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. Sweet-Cordero

**Writing, review, and/or revision of the manuscript:** Y. Zheng, H. Moore, E.A. Sweet-Cordero

**Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases):** Y. Zheng, E.A. Sweet-Cordero

**Study supervision:** Y. Zheng, E.A. Sweet-Cordero

## Grant Support

E.A. Sweet-Cordero was funded by a Research Scholar Grant from the American Cancer Society and by NCI 5 R01 CA157510-03. Y. Zheng was supported by a Senior Research Fellowship from the American Lung Association (RT-82480-N)

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.