## Abstract

The genetic and phenotypic heterogeneity of cancers can contribute to tumor aggressiveness, invasion, and resistance to therapy. Fluorescence imaging occupies a unique niche to investigate tumor heterogeneity due to its high resolution and molecular specificity. Here, heterogeneous populations are identified and quantified by combined optical metabolic imaging and subpopulation analysis (OMI-SPA). OMI probes the fluorescence intensities and lifetimes of metabolic enzymes in cells to provide images of cellular metabolism, and SPA models cell populations as mixed Gaussian distributions to identify cell subpopulations. In this study, OMI-SPA is characterized by simulation experiments and validated with cell experiments. To generate heterogeneous populations, two breast cancer cell lines, SKBr3 and MDA-MB-231, were co-cultured at varying proportions. OMI-SPA correctly identifies two populations with minimal mean and proportion error using the optical redox ratio (fluorescence intensity of NAD(P)H divided by the intensity of FAD), mean NAD(P)H fluorescence lifetime, and OMI index. Simulation experiments characterized the relationships between sample size, data standard deviation, and subpopulation mean separation distance required for OMI-SPA to identify subpopulations.

© 2015 Optical Society of America

## 1. Introduction

Solid tumors are highly heterogeneous, both across patients and within individual tumors. Tumor heterogeneity may contribute to tumor aggression, invasion, metastases, and therapy resistance [1–3]. Studies of tumor heterogeneity are challenging because traditional cell and tissue analyses require pooling of proteins, RNA, or DNA from hundreds to thousands of cells, which provides information on the protein and genetic expression of the majority of cells but may mask unique expression profiles and phenotypes of heterogeneous populations. Therefore, a technology capable of resolving cancer cell behaviors at single cell resolution is imperative for studies of tumor heterogeneity.

Breast cancers are often classified into three subtypes based on the expression or lack of expression of oncogenic proteins within the malignant cells. In particular, estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2), are used to classify breast cancers and inform clinical therapy selection. Breast cancer clinical subtypes include: triple negative (not expressing any of the three receptors), ER positive (expressing ER but lacking HER2 overexpression), and HER2 overexpressing (overexpressing HER2 receptors, ER and PR may be expressed or not). Primary tumors are heterogeneous, containing subpopulations of cells with different genetic profiles and protein expression [1, 4, 5]. Therapy selection is often based on the expression or lack of expression of ER and HER2 in a small sampling of cells, a biopsy. However, due to heterogeneity across large tumors and within the biopsy itself, populations of cells with different receptor expression and/or adaptations that make them resistant to therapy, may escape clinical detection [1]. This suggests that different portions of a tumor may respond differently to treatment, and while therapies are selected to target the majority of cells, there may be subpopulations of malignant cells not adequately targeted. Subpopulations of cells with stem cell like behaviors may contribute to drug resistance and recurrent tumor growth [6]. Cancer stem cells comprise between 3 and 35% of the number of cells within a solid breast tumor [6].

Fluorescence microscopy is uniquely suited to study tumor heterogeneity in cells and tissues due to the high resolution capabilities of microscopy and molecular specificity attained by probing fluorophores. Dynamic profiles of cellular metabolism can be obtained with optical metabolic imaging (OMI). OMI probes the endogenous fluorescence properties of two coenzymes involved in metabolism, reduced nicotinamide adenine dinucleotide (NAD(P)H) and flavin adenine dinucleotide (FAD). NAD(P)H and FAD are used in multiple metabolism processes including glycolysis and oxidative phosphorylation. The endpoints of OMI include the redox ratio, NAD(P)H fluorescence lifetime, FAD fluorescence lifetime, and a combination variable, the OMI index. The redox ratio is the intensity of NAD(P)H fluorescence relative to the intensity of FAD fluorescence and provides relative information on the global metabolism of the cell [7, 8]. The redox ratio is sensitive to shifts in metabolic pathways [7, 9]. The fluorescence lifetimes report changes in the microenvironment of NAD(P)H and FAD and are especially sensitive to the binding state of the fluorophore, as well as local temperature, pH, and proximity to quenchers such as molecular oxygen [10]. Both NAD(P)H and FAD fluorescence lifetimes can be either short or long, depending on the binding state of NAD(P)H and FAD (free or bound to an enzyme complex) [11, 12]. Previous studies have shown that OMI endpoints are sensitive to metabolism differences between cancer subtypes [9, 13, 14]. Additionally, the OMI endpoints provide dynamic readouts of cellular metabolism and detect pre-malignant transformations within tissues [15, 16], classify subtypes of breast cancer cells [9, 13], and detect response to anti-cancer drugs [14].

This study describes and validates OMI subpopulation analysis (OMI-SPA) for identification and quantification of cellular heterogeneity. OMI can be performed at high resolution to allow reporting of OMI endpoints at the single-cell level. Heterogeneity of the cell population can then be assessed by subpopulation analysis (SPA). SPA uses mixed distribution Gaussian models with multiple components to fit the population density representation of the data. In this study, we used both computational simulation and co-cultured experiments to validate OMI-SPA for identification and quantification of cellular subpopulation heterogeneity. Simulation experiments demonstrate the population characteristics required for robust SPA identification of two subpopulations. These experiments were validated by OMI-SPA of heterogeneous samples created by co-culturing two different breast cancer cell lines, a triple negative breast cancer (MDA-MB-231) and a HER2 + cell line (SKBr3) at varying proportions. These two cell lines were chosen to represent cell populations responsive (SKBr3) and resistant (MDA-MB-231) to the anti-HER2 antibody, trastuzumab. The OMI endpoints and morphology of these two cells are sufficiently different to allow computational separation and heterogeneity analysis.

## 2. Methods

#### 2.1 SKBr3 and MDA-MB-231 Specific Simulations

The average and standard deviation for the experimentally measured redox ratio, NAD(P)H mean lifetime, and FAD mean lifetime for the SKBr3 cells and MDA-MB-231 cells cultured independently were determined and used to generate model data sets. A random number generator (Matlab) was used to generate a data set of repeated observations within the distributions of the OMI endpoints for a population of *N* cells, containing *a* percent SKBr3 cells and *b* percent MDA-MB-231 cells. *N*, *a*, and *b* were varied to investigate the ability of OMI-SPA to identify two populations under a range of conditions. *N* ranged from 25 to 1000 cells and *a* varied from 0 to 1 (*b* = 1-*a*, and varied from 1 to 0). These simulate data sets were then evaluated with subpopulation analysis (Methods 2.2).

#### 2.2 Subpopulation Analysis (SPA)

Each cell population is modeled as a Gaussian mixture distribution model [9, 14, 17], $f\left(y;{\text{\Phi}}_{g}\right)=\underset{}{\overset{}{{\displaystyle {\displaystyle {\sum}_{i=1}^{g}{\pi}_{i}\varphi (y;{\mu}_{i},{V}_{i})}}}}$,where g is the number of components, $\varphi (y;{\mu}_{i},{V}_{i})$ represents a normal probability density function with mean μ_{i} and variance V_{i}, and π_{i} is the mixing proportion. Φ_{g} represents the unknown parameters, (π_{i}, μ_{i}, V_{i}): i = 1…g in a g-component model. The mixture model is fitted by maximum likelihood using the expectation maximization algorithm to determine the optimum parameters, (π, μ, V), (Matlab). In the simulation and co-culture experiments, data is modeled three times as Gaussian mixture distribution models with 1-3 components (g = 1, 2, or 3). Fit parameters for each model including the Akaiki information criteria (AIC), population means, population standard deviations, and proportions, were recorded. The AIC information criteria is a measure of model goodness of fit and is minimized in the optimal model [18]. The most representative model of the data was selected as the model with the lowest AIC.

#### 2.3 Cell culture

MDA-MB-231 and SKBr3 cells were grown separately in DMEM with 10% fetal bovine serum and 1% penicillin: streptomycin. Cells were plated at a density of 10^{6} cells per 35-mm glass bottom petri dish (MatTek Corp.), 48 hours before imaging, in the following proportions based on cell count (Table 1): 100% (10^{6}) MDA-MB-231 cells; 70% (7 x 10^{5}) MDA-MB-231 cells and 30% (3 x 10^{5}) SKBr3 cells; 50% (5 x 10^{5}) MDA-MB-231 cells and 50% (5 x 10^{5}) SKBr3 cells; 30% (3 x 10^{5}) MDA-MB-231 cells and 70% (7 x 10^{5}) SKBr3 cells; and 100% (10^{6}) SKBr3 cells. Two 35-mm glass bottom petri dishes were plated per group.

#### 2.4 Fluorescence lifetime instrumentation

Fluorescence lifetime imaging was performed on a custom built multi-photon microscope (Bruker), which has previously been described [9, 14, 19]. Excitation and emission light are coupled through a 40X oil immersion (1.3 NA) objective of an inverted microscope (Nikon, TiE). For NAD(P)H excitation, a titanium:sapphire laser (Coherent Inc.) was tuned to 750 nm (average power 7.5-7.9 mW). For FAD excitation, the laser was tuned to 890 nm (average power 8.4-6 mW). Bandpass filters, with a 440-480 nm passband for NAD(P)H and a 500-600 nm passband for FAD, isolated emission light. A pixel dwell time of 4.8 μs was used to acquire 256x256 pixel images. Each fluorescence lifetime image was collected using time correlated single photon counting electronics (SPC-150, Becker and Hickl) and a GaAsP PMT (H7422P-40, Hamamatsu). Photon count rates were maintained above 5x10^{5} for the entire 60 second image acquisition time, ensuring no photobleaching occurred and adequate signal for fluorescence lifetime decay fits.

#### 2.5 Cell imaging

Cells were imaged directly through the bottom of 35mm glass-bottom petri dishes (MatTek Corp). For each dish, six representative fields of view were imaged, for a total number of ~200 cells per group. First, an NAD(P)H lifetime image was acquired, and then an FAD lifetime image was acquired from the same field of view. Sequential fields of view were separated by at least 1 field of view, 270 μm.

#### 2.6 Generation of redox ratio images

A fluorescence intensity image was generated by integrating the fluorescence lifetime decay over time for each pixel in the lifetime image. The total number of NAD(P)H photons per pixel was then divided by the total number of FAD photons per pixel and used to create a redox ratio image for each field of view (Matlab, Mathworks). NAD(P)H and FAD fluorescence specific to cellular metabolism is localized in the cytoplasm and mitochondria. Therefore, the redox ratio image was thresholded to remove background and nuclear fluorescence, and the average redox ratio for each remaining cell cytoplasm was computed (ImageJ). Cells were manually segmented in ImageJ.

#### 2.7 Generation of NAD(P)H and FAD lifetime images

For each pixel, the photon counts for the 9 surrounding pixels were binned (SPCImage). Fluorescence lifetime components for each pixel were computed by de-convolving a measured system response curve from urea crystals and fitting the fluorescent decay to a two component model, $I\left(t\right)={\alpha}_{1}\mathrm{exp}\left(-\frac{t}{{\tau}_{1}}\right)+{\alpha}_{2}\mathrm{exp}\left(-\frac{t}{{\tau}_{2}}\right)+C$. In this model, I(t) is the fluorescence intensity at time t after the laser excitation pulse, α_{1} and α_{2} are the fractional contributions of the free and bound molecules (i.e. α_{1} + α_{2} = 1), τ_{1} and τ_{1} are the fluorescence lifetimes of the short and long lifetime components, and C is a constant that accounts for background light. Matrices of the lifetime components were exported as ascii files for further processing in ImageJ. The mean lifetime, τ_{m} = α_{1} τ_{1} + α_{2} τ_{2}, was computed for each pixel to create a mean lifetime image for each field of view. The NAD(P)H τ_{m} and FAD τ_{m} for each cell cytoplasm was computed and recorded. The OMI index [14] is a linear combination of mean centered redox ratio, NAD(P)H τ_{m}, and FAD τ_{m} data and can be computed per cell as follows: $OMIIndex=\frac{R{R}_{i}}{\u3008RR\u3009}+\frac{NAD\left(P\right)H{\tau}_{mi}}{\u3008NAD\left(P\right)H{\tau}_{m}\u3009}-\frac{FAD{\tau}_{mi}}{\u3008FAD{\tau}_{m}\u3009}$.

#### 2.8 OMI-SPA behavior simulations

Simulation experiments were performed to model the behavior of OMI-SPA for generalized data. In the first simulation, populations of cells were generated with a normalized mean of 1 and equal proportions of the two populations. The standard deviation and distance between the means was varied from 0 to 1 and 0-2, respectively. In the last simulations, the mean distance was varied from 0 to 2, the population proportions varied from 0 to 1, and the standard deviation of the populations was 0.05, 0.1, 0.25, or 0.5. These normalized standard deviation values span the range typically observed in OMI data [9, 14]. Simulated cell populations were generated for each mean distance, population proportion, and standard deviation at increasing sample sizes (up to 10,000) to determine the minimum sample size necessary to resolve two subpopulations (AIC_{2} < AIC_{1} −20).

## 3. Results

#### 3.1 MDA-MB-231 and SKBr3 simulations

First, SPA was performed on simulated cell populations using experimentally determined mean and standard deviation values of the redox ratio of SKBr3 and MDA-MB-231 cells, 1.92 +/− 0.39 and 1.06 +/− 0.22, respectively. The AIC represents a model goodness of fit and is minimized in the optimal model [18]. The AIC includes penalties for an increased number of model components, to account for increased fitness of models with increased components [18]. The AIC for two components is less than the AIC for the single component model (AIC_{1}-AIC_{2} > 0) at large sample sizes and SKBr3 proportions less than 0.8 (Fig. 1a). The error of the SKBr3 redox ratio mean is greatest (>10%) at high and low proportions and low sample sizes (Fig. 1b). SKBr3 redox ratio mean error is minimized (<5%) at high sample sizes and proportions near 50% (Fig. 1b). Likewise, MDA-MB-231 mean error is maximal (>10%) at small sample sizes and low proportions, and minimized (<5%) with high sample sizes and proportions around 50% (Fig. 1c). The large errors of MDA-MB-231 mean values at SKBr3 population proportions greater than 0.8, correspond to data points where the two-component model fit the data less well than the single component model. The population proportion error, or the error of how many cells are attributed to each cell subpopulation, is minimal (<5%) at high sample sizes and at SBRr3 proportions less than 80% (Fig. 1d).

Similar simulations demonstrate the ability of OMI-SPA to identify SKBr3 and MDA-MB-231 subpopulations using NAD(P)H τ_{m} values (Fig. 2). The SKBr3 cells have a mean NAD(P)H τ_{m} of 1.31 +/− 0.13 ns and the MDA-MB-231 cells have a mean NAD(P)H τ_{m} of 0.85 +/− 0.11 ns, determined from experiments. In the simulation, OMI-SPA identifies two populations in all simulated data sets with a sample size greater than 75 (Fig. 2a) and the error of the SKBr3 and MDA-MB-231 subpopulation means identified by the models is within 1% at sample sizes greater than 100 and proportions between 20% and 80% (Fig. 2b-c). Likewise, the proportion error is minimized (<5%) at sample sizes > 200 and SKBr3 proportions between 20 and 80%.

Finally, the simulations were repeated using the mean and standard deviations of the SKBr3 and MDA-MB-231 mean FAD lifetimes, which were experimentally measured as 1.09 +/− 0.10 ns and 1.12 +/− 0.13 ns, respectively. For most of the simulated cell populations, the two component model does not have an AIC value less than the single component fit (Fig. 3a). The error of the FAD lifetime values for the modeled SKBr3 and MDA-MB-231 subpopulations was within 10% of the true FAD lifetime values except at sample sizes < 100 (Fig. 3b-c); however, the error of the proportion of cells assigned to each subpopulation was greater than 10% at almost all sample sizes and proportions (Fig. 3d).

#### 3.2 MDA-MB-231 and SKBr3 co-culture experiments

High-resolution images (Fig. 4) demonstrate the variability of redox ratios, NAD(P)H mean lifetimes, and FAD mean lifetimes between the two cell lines. MDA-MB-231 cells have a lower redox ratio, a shorter mean NAD(P)H lifetime, and a longer mean FAD lifetime, than SKBr3 cells. The images from the co-culture experiments demonstrate heterogeneity in cellular morphology and fluorescence imaging endpoints (Fig. 4). The images from the co-culture experiments were examined manually and each cell designated as SKBr3 or MDA-MB-231 based on morphology. Any cells that could not be identified were excluded from the analysis. While the cells were plated at the following proportions, 30/70%, 50/50%, and 70/30%, the manual analysis verified the heterogeneity of the co-cultured experiments and reports the actual proportions of the imaged cells for each co-culture group (Table 1).

OMI-SPA of the co-culture experiment reveals that the single cell data sets are modeled best by a one component Gaussian distribution model and data from the co-culture groups are best modeled by two component Gaussian distribution models. OMI-SPA of the redox ratio correctly identifies two populations in all three co-culture experiments (Fig. 5). The mean redox ratio values of the SKBr3 modeled populations are within 6% of the true mean redox ratio (defined as the mean redox ratio of manually classified SKBr3 cells) (Table 2). Likewise, the mean error of the MDA-MB-231 cells is 5% or less for all co-culture experiments (Table 2). The error of the model-estimated proportion of each subpopulation of cells is 5% or less for all three co-culture groups (Table 2). The error of the standard deviations for the identified SKBr3 populations ranges from 5% for the 100% SKBr3 culture to 25% for the 30% SKBr3 culture, while the standard deviation errors for the MDA-MB-231 model redox ratios is 10% or less for all experiments (Table 2).

The mean redox ratio for each population of cells, SKBr3 and MDA-MB-231, was compared across all experimental groups to evaluate the stability of the redox ratio when cells are grown in different environments. The MDA-MB-231 cells grown independently of the SKBr3 cells had a mean redox ratio of 1 and a standard deviation of 0.23. The mean MDA-MB-231 redox ratio was not significantly different for any of the co-culture experiments: 0.97+/− 0.22 for 70% MDA-MB-231 cells, 0.99+/−0.2 for 50%, and 0.82 +/− 0.25 for 30%. Likewise, the redox ratio for the SKBr3 cells was similar when the cells were cultured without and with the MDA-MB-231 cells. The only statistically significant (p<0.05) difference in redox ratio was observed for the SKBr3 cells grown in the 30% SKBr3/70% MDA-MB-231 culture versus the 100% SkBr3 culture, 1.2+/−0.5 vs. 1.8 +/− 0.39, respectively. These results suggest that the redox ratio remains relatively stable irrespective of whether cells are grown independently or in co-culture.

The mean NAD(P)H τ_{m} values of the two cell lines grown independently are 1.3 ns for SKBr3 cells and 0.85 ns for the MDA-MB231 cells (Fig. 6a-b). SPA of the NAD(P)H τ_{m} data from the co-culture experiments reveals two populations for all three co-culture groups (Fig. 6c-e). The errors of the means of the computationally determined SKBr3 and MDA-MB-231 subpopulations are within 7% of the respective, true NAD(P)H τ_{m} means (Table 2). Likewise, the proportion of cells attributed to each population as determined computationally by the SPA models, is within 10% of the true populations and the standard deviation errors within 28% (Table 2). To demonstrate the cellular heterogeneity visually, the average NAD(P)H τ_{m} value of each cell was used to classify cells as SKBr3 or MDA-MB-231, and cells are color coded red if SKBr3 and blue if MDA-MB-231 (Fig. 7). An NAD(P)H τ_{m} cut-off value of 1.06 ns was chosen because it is equal-distance, 1.9 standard deviations, from each mean and minimized classification error. This NAD(P)H τ_{m} threshold value yielded an overall accuracy of 95.2% to correctly classify all cells, with 94.2% of SKBr3 cells correctly classified and 95.5% of MDA-MB-231 cells correctly classified.

SPA of the heterogeneous co-cultures is less successful at separating the two populations using the FAD τ_{m} (Fig. 8, Table 2). While two populations are identified in the 30/70% experiments, the model fails to identify two populations in the 50/50% group (Fig. 8d).

The OMI index is a linear combination of mean centered redox ratio, NAD(P)H τ_{m}, and FAD τ_{m} values computed per cell [14]. The mean OMI index values for SKBr3 and MDA-MB-231 cells are 1.67 +/− 0.31 and 0.54 +/− 0.22, respectively. OMI-SPA of the OMI index for the co-culture experiments resolves a single population for the 100% SKBr3 and 100% MDA-MB-231 groups and two populations for each of the co-culture groups (Fig. 9). The error of the SKBr3 subpopulation OMI index means is <6% for all co-culture experiments (Table 2). The mean error for the OMI index for the modeled populations of MDA-MB-231 cells are all less than 4% (Table 2), and the proportions for the SKBr3 and MDA-MB-231 populations within 1% (Table 2). The standard deviation error for the modeled SKBr3 populations is 20% or less for all experiments and 11% or less for the MDA-MB-231 populations (Table 2).

#### 3.3 Behavior of AIC

Simulations were performed to determine the characteristics of data sets from which two distinct subpopulations can be resolved. In the first experiments, populations of size N = 300 were generated with a random number generator and each population proportion was 0.5. All populations had a normalized mean of 1, and the distance between the means was varied from 0 to 2. The standard deviation of the populations varied from 0 to 1. For these simulations, the two component model is most representative of the data (AIC_{2}<AIC_{1}) at low standard deviations and at greater distances between the means (Fig. 10a). If a cutoff of AIC_{2}<AIC_{1} is used to select the most representative model, a two component model is most representative when the standard deviation < 0.46*(distance between means) + 0.0821. However, if the difference AIC_{1}-AIC_{2} must be greater than 5% of AIC_{1}, then standard deviations < 0.2563*(Distance between means) + 0.0832 describes the populations with two components that can be identified by OMI-SPA. However, even though the AIC selects a two component model for these populations, the error of the computed population mean values, standard deviations, and proportions increases with increased subpopulation standard deviation and decreased distance between the means (Fig. 10b-d).

Similar simulations were performed to determine the sample size necessary to resolve two populations of known, normalized mean distances, standard deviations, and proportions (Fig. 11). In Fig. 11, the simulations, in which the normalized distance varied from 0 to 2 and the proportion of each population varied from 0 to 1, are repeated at increasing standard deviation values (a) 0.05, (b) 0.1, (c) 0.25, and (d) 0.5. These simulations demonstrate that SPA requires a larger sample size to identify two subpopulations at smaller distances between the means (Fig. 11). As the standard deviation increases, an increased sample size is required to identify subpopulations. SPA fails to identify two subpopulations within a sample size of 10,000, at low mean distances: mean distances < 0.3 for a standard deviation of 0.05, mean distances < 0.4 for a standard deviation of 0.1, mean distances < 0.6 for a standard deviation of 0.25, and mean distances < 0.8 for a standard deviation of 0.5. Additionally, SPA fails to identify two populations in data sets up to 10,000 cells, at percentages of 0 and 1, as expected for uni-modal populations. Increased sample sizes are required to identify populations of very small proportions, in the 0.02-0.12 range (Fig. 11c). The dashed red circle (Fig. 11b) and red dashed line (Fig. 11d) encompass measured subpopulations [9, 14] in breast cancer (see Discussion). Altogether, these simulations on normalized data sets demonstrate the relationships between mean distance, standard deviation, population proportion, and sample size necessary to resolve two subpopulations.

## 4. Discussion

Recent evidence suggests that tumor heterogeneity is a major source of drug resistance [1]. Cancer stem cell populations may contribute to drug resistance and tumor recurrence [6]. Within breast cancers, stem cells can compromise 3-35% of the tumor cells [6]. Therefore, non-invasive single cell analysis technologies that can monitor changes in cellular subpopulations over time are critically needed to understand and combat drug resistance in tumors. Optical imaging of NAD(P)H and FAD is highly sensitive to metabolic differences among breast cancer subtypes [9, 14], and can be performed at high resolution for single-cell analysis. Here, we demonstrate the capabilities and limitations of OMI-SPA to quantify subpopulations of cells by co-culture and simulation experiments.

Simulations of the OMI-SPA behavior for co-cultured SKBr3 and MDA-MB-231 cells (Fig. 1-3) demonstrate that for these two cell lines, OMI-SPA will be able to identify the two populations, with minimal mean and proportion error, using the NAD(P)H τ_{m} values at population sizes greater than 50 cells and at proportions between 0.1 and 0.9 (Fig. 2). The normalized distance between mean redox ratios for SKBr3 and MDA-MB-231 is greater than that of NAD(P)H τ_{m} (0.57 vs. 0.42), but the normalized average standard deviation is twice as great for the redox ratios than NAD(P)H τ_{m} (0.22 vs. 0.12). Therefore, the redox ratio simulations demonstrate that OMI-SPA can resolve two populations with minimal mean and proportion error at larger sample sizes (>100 total cells) and proportions between 20 and 80% SKBr3 (Fig. 1). While the normalized average standard deviation of the SKBr3 and MDA-MB-231 FAD τ_{m} values is 0.10, the smallest of these three OMI endpoints, the normalized mean FAD τ_{m} distance for SKBr3 and MDA-MB-231 cells is only 0.06, and OMI-SPA is unable to identify two populations in the FAD τ_{m} simulations due to this small difference in means of the two populations (Fig. 3). This highlights the fact that OMI-SPA can identify two populations only if the two populations have sufficiently large mean separation and small standard deviations. Otherwise, OMI-SPA cannot identify both subpopulations.

Co-culture experiments of SKBr3 and MDA-MB-231 cells were performed to test OMI-SPA on experimental data. Due to the large differences in mean redox ratio and NAD(P)H τ_{m}, the redox ratio and NAD(P)H τ_{m} images of co-cultured SKBr3 and MDA-MB-231 cells appear heterogeneous (Fig. 4). This highlights the advantage of microscopy techniques that provide images for qualitative and quantitative assessment. Furthermore, false-coloring of the images based on the cellular NAD(P)H τ_{m} values visually demonstrates the presence of SKBr3 cells and MDA-MB-231 cells in the images of co-cultured cells (Fig. 7). OMI-SPA of the redox ratio and NAD τ_{m} values of the co-cultured SKBr3 and MDA-MB-231 cells identified two subpopulations with minimal error (Fig. 5, 6; Table 2). OMI-SPA of the FAD τ_{m} data failed to resolve two populations in the 50%/50% co-culture experiment, which was expected from the simulation data that suggested the FAD τ_{m} means were too similar for OMI-SPA to distinguish (Fig. 8).

The OMI index is a linear combination of the redox ratio, NAD(P)H τ_{m}, and FAD τ_{m}, and provides a single optical variable that is highly sensitive to cellular metabolic changes with receptor expression and drug response in breast cancer [14]. Subpopulation analysis of the OMI index of co-cultured SKBr3 and MDA-MB-231 cells identifies two subpopulations when grown at 30/70%, 50/50%, and 70/30% with minimal error in population means, standard deviations, and subpopulation proportions (Fig. 9, Table 2). In comparison with the three individual OMI index constituents, the OMI index has the greatest normalized mean distance between SKBr3 and MDA-MB-231 cells, 1.02, and a normalized standard deviation, 0.24, indicating robust SPA for these two cell lines.

Generalized simulations were performed to assess the characteristics of data sets that can be accurately modeled as two subpopulations by OMI-SPA. These generalized simulations can be used to predict SPA performance for any two subpopulations by computing a subpopulation mean separation distance and standard deviation normalized to the mean of the data set. The results of these simulations (Fig. 10-11) demonstrate the relationships between sample size, population standard deviation, and distance between means necessary to resolve two populations. A larger standard deviation is tolerated if the distance between means is greater (Fig. 10) or the sample size is larger (Fig. 11). OMI-SPA subpopulation mean, standard deviation, and proportion errors increase with increased population standard deviation and decreased distances between means (Fig. 10). Cancer stem cells compromise 3-35% of the cells within a tumor [6], and the generalized OMI-SPA results (Fig. 11) demonstrate the capabilities of OMI-SPA to resolve these small populations provided that the stem cells have different OMI endpoints from the non-stem cell population.

The utility of OMI-SPA is to identify subpopulations of cells within unknown, heterogeneous populations. OMI-SPA has been used in xenograft and patient derived- primary tumor organoids to identify subpopulations with varied response to anti-cancer drug treatment [14]. The normalized mean separation for the OMI index of trastuzumab responsive (BT474, ER+/HER2 + ) cells and resistant cells (HR6, ER+/HER2 + ) is 1.16 and the normalized OMI index standard deviation is 0.42 for responsive cells and 0.28 for resistant cells [9]. The normalized simulation experiments (red dashed line; Fig. 11d) suggest that OMI-SPA will be able to resolve heterogeneous populations of responsive and non-responsive cells, if there are at least 500 cells in the population. Analysis of primary-tumor derived organoids with heterogeneous drug-response, identified subpopulations of cells with varied OMI index values [14]. The normalized mean separation for OMI-SPA subpopulations within patient primary tumor derived organoids ranged from 0.4 to 1.4, population percentages ranged from 0.15 to 0.85, and the normalized standard deviations ranged from 0.02 to 0.36. As indicated by the red dashed circle (Fig. 11b), OMI-SPA can resolve these subpopulations with high mean separations with as few as 100 cells, and can resolve the majority of these subpopulations with less than 500 cells.

Heterogeneous tumors may contain more than two subpopulations and the characteristics of these subpopulations may be unknown in advance. Furthermore, the metabolism of these subpopulations may change dynamically over time and in response to stimuli. SPA, as described here, can be used to evaluate the presence of multiple subpopulations. Additionally, the simulation experiments (Fig. 10 and 11) generalize the OMI-SPA technique and show this is a robust method for populations of cells with varying OMI endpoint means and standard deviations and not limited to the two specific cell lines used in the co-culture experiments. Furthermore, one advantage of OMI is the non-invasive capabilities which enable imaging of *in vitro* samples over time to study subpopulation and individual cell dynamics. Altogether, OMI-SPA is a robust, attractive platform for evaluating and monitoring subpopulation heterogeneity.

OMI-SPA assumes that the OMI endpoints distributions for each population exhibit normality. A Wilk-Shapiro test for normality revealed that the NAD(P)H and FAD mean lifetimes of both cell lines exhibited normality but the OMI index did not for either cell line and neither did the redox ratio for the SKBr3 cells. The lack of normality of these data sets may introduce error when modeled as Gaussian curves. Even with this know error, however, the errors for the mean OMI index modeled for both cell lines and the errors for the SKBr3 redox ratio in the co-culture experiments was less than 6%, suggesting OMI-SPA performs well with these distributions represented as Gaussian curves. SPA can be improved to account for different data distributions by using additional distributions that better represent the homogenous population.

The results of this experiment demonstrate that OMI-SPA can be used to identify cancer cell subpopulations based on the OMI endpoints: redox ratio, NAD(P)H mean lifetime, FAD mean lifetime, and OMI index. Furthermore, these results characterize the relationships between sample size, standard deviation, and mean distance required for OMI-SPA to accurately describe the two populations. Our previously published analyses of cellular subpopulations within *in vivo* tumors [9] and patient-derived tumor organoids [14] also indicate that OMI-SPA can accurately identify cell subpopulations with a sample of as few as 300 cells.

## Acknowledgements

We acknowledge Miranda Kunz for assistance with cell culture and imaging. Funding sources include Vanderbilt Provost Fellowship (AJW), the DOD BCRP (DOD-BC121998), the NSF (AJW; DGE-0909667), the NIH/NCI (NIH R00-CA142888, R01-CA185747), the Mary Kay Foundation (067-14), and the NCI Breast Cancer SPORE (P50-CA098131).

## References and links

**1. **R. Fisher, L. Pusztai, and C. Swanton, “Cancer heterogeneity: implications for targeted therapeutics,” Br. J. Cancer **108**(3), 479–485 (2013). [CrossRef] [PubMed]

**2. **D. J. Kiviet, P. Nghe, N. Walker, S. Boulineau, V. Sunderlikova, and S. J. Tans, “Stochasticity of metabolism and growth at the single-cell level,” Nature **514**(7522), 376–379 (2014). [CrossRef] [PubMed]

**3. **K. J. Cheung, E. Gabrielson, Z. Werb, and A. J. Ewald, “Collective invasion in breast cancer requires a conserved basal epithelial program,” Cell **155**(7), 1639–1651 (2013). [CrossRef] [PubMed]

**4. **V. Almendro, Y. K. Cheng, A. Randles, S. Itzkovitz, A. Marusyk, E. Ametller, X. Gonzalez-Farre, M. Muñoz, H. G. Russnes, A. Helland, I. H. Rye, A. L. Borresen-Dale, R. Maruyama, A. van Oudenaarden, M. Dowsett, R. L. Jones, J. Reis-Filho, P. Gascon, M. Gönen, F. Michor, and K. Polyak, “Inference of tumor evolution during chemotherapy by computational modeling and in situ analysis of genetic and phenotypic cellular diversity,” Cell Rep **6**(3), 514–527 (2014). [CrossRef] [PubMed]

**5. **K. Polyak, “Tumor Heterogeneity Confounds and Illuminates: A case for Darwinian tumor evolution,” Nat. Med. **20**(4), 344–346 (2014). [CrossRef] [PubMed]

**6. **J. E. Visvader and G. J. Lindeman, “Cancer stem cells in solid tumours: accumulating evidence and unresolved questions,” Nat. Rev. Cancer **8**(10), 755–768 (2008). [CrossRef] [PubMed]

**7. **I. Georgakoudi and K. P. Quinn, “Optical imaging using endogenous contrast to assess metabolic state,” Annu. Rev. Biomed. Eng. **14**(1), 351–367 (2012). [CrossRef] [PubMed]

**8. **B. Chance, B. Schoener, R. Oshino, F. Itshak, and Y. Nakase, “Oxidation-reduction ratio studies of mitochondria in freeze-trapped samples. NADH and flavoprotein fluorescence signals,” J. Biol. Chem. **254**(11), 4764–4771 (1979). [PubMed]

**9. **A. J. Walsh, R. S. Cook, H. C. Manning, D. J. Hicks, A. Lafontant, C. L. Arteaga, and M. C. Skala, “Optical metabolic imaging identifies glycolytic levels, subtypes, and early-treatment response in breast cancer,” Cancer Res. **73**(20), 6164–6174 (2013). [CrossRef] [PubMed]

**10. **J. Lakowicz, *Principles of fluorescence spectroscopy* (Plenum Publishers, New York, 1999).

**11. **J. R. Lakowicz, H. Szmacinski, K. Nowaczyk, and M. L. Johnson, “Fluorescence Lifetime Imaging of Free and Protein-Bound NADH,” Proc. Natl. Acad. Sci. U.S.A. **89**(4), 1271–1275 (1992). [CrossRef] [PubMed]

**12. **F. Tanaka, N. Tamai, and I. Yamazaki, “Picosecond-resolved fluorescence spectra of D-amino-acid oxidase. A new fluorescent species of the coenzyme,” Biochemistry **28**(10), 4259–4262 (1989). [CrossRef] [PubMed]

**13. **A. Walsh, R. S. Cook, B. Rexer, C. L. Arteaga, and M. C. Skala, “Optical imaging of metabolism in HER2 overexpressing breast cancer cells,” Biomed. Opt. Express **3**(1), 75–85 (2012). [CrossRef] [PubMed]

**14. **A. J. Walsh, R. S. Cook, M. E. Sanders, L. Aurisicchio, G. Ciliberto, C. L. Arteaga, and M. C. Skala, “Quantitative optical imaging of primary tumor organoid metabolism predicts drug response in breast cancer,” Cancer Res. **74**(18), 5184–5194 (2014). [CrossRef] [PubMed]

**15. **M. C. Skala, K. M. Riching, D. K. Bird, A. Gendron-Fitzpatrick, J. Eickhoff, K. W. Eliceiri, P. J. Keely, and N. Ramanujam, “In vivo multiphoton fluorescence lifetime imaging of protein-bound and free nicotinamide adenine dinucleotide in normal and precancerous epithelia,” J. Biomed. Opt. **12**(2), 024014 (2007). [CrossRef] [PubMed]

**16. **M. C. Skala, K. M. Riching, A. Gendron-Fitzpatrick, J. Eickhoff, K. W. Eliceiri, J. G. White, and N. Ramanujam, “In vivo multiphoton microscopy of NADH and FAD redox states, fluorescence lifetimes, and cellular morphology in precancerous epithelia,” Proc. Natl. Acad. Sci. U.S.A. **104**(49), 19494–19499 (2007). [CrossRef] [PubMed]

**17. **W. Pan, J. Lin, and C. T. Le, “Model-based cluster analysis of microarray gene-expression data,” Genome Biol. **3**(2), H0009 (2002). [CrossRef] [PubMed]

**18. **H. Akaike, “A new look at the statistical model identification,” IEEE Trans. Automatic Control **19**(6), 716–723 (1974). [CrossRef]

**19. **A. J. Walsh, K. M. Poole, C. L. Duvall, and M. C. Skala, “Ex vivo optical metabolic measurements from cultured tissue reflect in vivo tissue status,” J. Biomed. Opt. **17**(11), 116015 (2012). [CrossRef] [PubMed]