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{Φ}}_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⁶ 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⁶) MDA-MB-231 cells; 70% (7 x 10⁵) MDA-MB-231 cells and 30% (3 x 10⁵) SKBr3 cells; 50% (5 x 10⁵) MDA-MB-231 cells and 50% (5 x 10⁵) SKBr3 cells; 30% (3 x 10⁵) MDA-MB-231 cells and 70% (7 x 10⁵) SKBr3 cells; and 100% (10⁶) SKBr3 cells. Two 35-mm glass bottom petri dishes were plated per group.

Table 1. Experimental groups for the SKBr3 and MDA-MB-231 co-culture experiments.

View Table | View all tables in this article

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⁵ 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\exp \left(-\frac{t}{{\tau }_1}\right)+{\alpha }_2\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, α₁ and α₂ are the fractional contributions of the free and bound molecules (i.e. α₁ + α₂ = 1), τ₁ and τ₁ 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 = α₁ τ₁ + α₂ τ₂, 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: $OMI Index= \frac{R{R}_i}{〈RR〉}+\frac{NAD\left(P\right)H {\tau }_{mi}}{〈NAD\left(P\right)H {\tau }_m〉}-\frac{FAD {\tau }_{mi}}{〈FAD {\tau }_m〉}$.

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₂ < AIC₁ −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₁-AIC₂ > 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).

 

Fig. 1 Three dimensional graphs of the behavior of OMI-SPA to model simulated data using the redox ratio mean and standard deviation of SKBr3 and MDA-MB-231 cells. (a) Difference in the AIC values of 1 and 2 component Gaussian models of the simulated cell populations as a function of varying sample sizes and varying proportions (0% to 100%) of SKBr3 cells. AIC₁-AIC₂ > 0 indicates two component model is a better fit than the one-component model. (b) Error of the mean redox ratio computed for the SKBr3 subpopulation, and (c) MDA-MB-231 population. (d) Error of the estimated proportion of the SKBr3 subpopulation.

Download Full Size | PDF

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%.

 

Fig. 2 The behavior of OMI-SPA to model data simulated from the mean and standard deviations of NAD(P)H τ_m of SKBr3 and MDA-MB-231 cells. (a) Difference in the AIC values of 1 and 2 component Gaussian models of the simulated cell populations as a function of varying sample sizes and varying proportion of SKBr3 cells. (b) Error of the mean NAD(P)H τ_m computed for the SKBr3 subpopulation, and (c) MDA-MB-231 population. (d) Error of the estimated proportion of the SKBr3 subpopulation.

Download Full Size | PDF

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).

 

Fig. 3 The behavior of OMI-SPA to model data simulated from the mean and standard deviations of the FAD τ_m of SKBr3 and MDA-MB-231 cells. (a) Difference in the AIC values of 1 and 2 component Gaussian models to the simulated cell populations as a function of varying sample sizes (N) and varying proportion of SKBr3 cells. (b) Error of the mean FAD τ_m computed for the SKBr3 subpopulation and (c) MDA-MB-231 population. (d) Error of the estimated proportion of the SKBr3 subpopulation.

Download Full Size | PDF

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).

 

Fig. 4 Representative images of the redox ratio (NAD(P)H/FAD), NAD(P)H mean lifetime, and FAD mean lifetime of MDA-MB-231 cells, SKBr3 cells, and mixed populations.

Download Full Size | PDF

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).

 

Fig. 5 Histograms of redox ratios, quantified per cell, of a population of approximately 200 cells of varying percentages of MDA-MB-231 cells and SKBr3 cells measured experimentally. The solid blue line represents best mixed-model Gaussian distribution fit. The red dashed curves represent the two component contributions, if a two component model is optimal. Histograms are normalized to have a total area of 1.

Download Full Size | PDF

Table 2. Mean, proportion (P), and standard deviation and % errors computed from the optimal fitting Gaussian distribution model of the co-culture experimental data.

View Table | View all tables in this article

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.

 

Fig. 6 Histograms of NAD(P)H τ_m from populations of approximately 200 cells of varying percentages of MDA-MB-231 cells and SKBr3 cells. The solid blue line represents the best mixed-model Gaussian distribution fit. The red dashed curves represent the two component contributions, if a two component model is optimal. Histograms are normalized to have a total area of 1.

Download Full Size | PDF

 

Fig. 7 Representative images with cells color coded red if the cell mean NAD(P)H τ_m value is greater than 1.06 ns and blue if the NAD(P)H τ_m value is less than 1.06 ns.

Download Full Size | PDF

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).

 

Fig. 8 Histograms of FAD τ_m from populations of approximately 200 cells of varying percentages of MDA-MB-231 cells and SKBr3 cells. The solid blue line represents the best mixed-model Gaussian distribution fit. The red dashed curves represent the two component contributions, if a two component model is optimal. Histograms are normalized to have a total area of 1.

Download Full Size | PDF

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).

 

Fig. 9 Histograms of OMI Index from populations of approximately 200 cells of varying percentages of MDA-MB-231 cells and SKBr3 cells. The solid blue line represents the best mixed-model Gaussian distribution fit. The red dashed curves represent the two component contributions, if a two component model is optimal. Histograms are normalized to have a total area of 1.

Download Full Size | PDF

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₂<AIC₁) at low standard deviations and at greater distances between the means (Fig. 10a). If a cutoff of AIC₂<AIC₁ 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₁-AIC₂ must be greater than 5% of AIC₁, 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).

 

Fig. 10 Behavior of OMI-SPA for generalized data. Simulation initial conditions include a population of size N = 300 cells, each population proportion was 0.5, and all populations had a normalized mean of 1. The distance between the means was varied from 0 to 2. The standard deviation of the populations varied from 0 to 1. (a) AIC difference for the 1 and 2 component models to fit the simulated data. (b) Error of the mean, (c) error of the variance, and (d) error of the subpopulation proportions computed for the two-component models of the simulated data.

Download Full Size | PDF

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.

 

Fig. 11 Minimum sample size (N), if less than 10,000 cells, required to resolve two populations with minimal population mean and proportion error (AIC₂<AIC₁-20) at varied mean distances (0-2) and population proportions (0-1). Standard deviation of the populations was varied: (a) 0.05, (b) 0.1, (c) 0.25, and (d) 0.5. The red, dashed circle in (b) encompasses the majority of normalized OMI endpoint mean distances and proportions of subpopulations observed in patient-derived organoids. The red, dashed line (d) represents the normalized distance (1.16) and standard deviation between two HER2 positive cell lines, one responsive to trastuzumab and one resistant.

Download Full Size | PDF

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.