Kinetic Modeling of ABCG2 Transporter Heterogeneity: A Quantitative, Single-Cell Analysis of the Side Population Assay

The side population (SP) assay, a technique used in cancer and stem cell research, assesses the activity of ABC transporters on Hoechst staining in the presence and absence of transporter inhibition, identifying SP and non-SP cell (NSP) subpopulations by differential staining intensity. The interpretation of the assay is complicated because the transporter-mediated mechanisms fail to account for cell-to-cell variability within a population or adequately control the direct role of transporter activity on staining intensity. We hypothesized that differences in dye kinetics at the single-cell level, such as ABCG2 transporter-mediated efflux and DNA binding, are responsible for the differential cell staining that demarcates SP/NSP identity. We report changes in A549 phenotype during time in culture and with TGFβ treatment that correlate with SP size. Clonal expansion of individually sorted cells re-established both SP and NSPs, indicating that SP membership is dynamic. To assess the validity of a purely kinetics-based interpretation of SP/NSP identity, we developed a computational approach that simulated cell staining within a heterogeneous cell population; this exercise allowed for the direct inference of the role of transporter activity and inhibition on cell staining. Our simulated SP assay yielded appropriate SP responses for kinetic scenarios in which high transporter activity existed in a portion of the cells and little differential staining occurred in the majority of the population. With our approach for single-cell analysis, we observed SP and NSP cells at both ends of a transporter activity continuum, demonstrating that features of transporter activity as well as DNA content are determinants of SP/NSP identity.


Introduction
The side population (SP) assay is used to identify stem cells by flow cytometry through the characteristic of enhanced dye efflux mediated via ATP-binding cassette (ABC) transporters [1]. The SP was first identified by Goodell et al. as hematopoietic stem cells in samples of murine bone marrow aspirate [2]. The role of the SP has since expanded to serve as a means to identify stem cell populations based, primarily, on ABCG2 activity [3], though additional ABC transporters, such as P-glycoprotein/ABCB1, can also mediate formation of a SP [4]. ABCG2, also known as breast cancer resistance protein (BRCP), can mediate multidrug resistance (MDR) in breast [5][6][7][8][9] and other cell lines [10][11][12][13][14]. The SP has been implicated in numerous cancers as a harbinger of MDR-mediated chemoresistance [15][16][17][18] and cancer stem cells (CSCs) [19][20][21][22] in in vitro cancer cell lines; thus the presence of a SP is understood as an undesirable indicator.
SPs are identified by splitting samples into conditions with and without an ABC transporter inhibitor followed by Hoechst staining, which enables population-level comparison of differences in cell staining due to ABC transporter activity between the two conditions ( Fig 1A). Blocking of transporter mediated Hoechst efflux by the inhibitor serves as a basis for comparison of cell staining in the condition without the transporter inhibitor ( Fig 1B). When comparing the two conditions, SP cells are observed as a population with decreased staining in the lower left of the Hoechst Red and Blue staining plot ( Fig 1A).
The basis of differential staining is thought to be driven by impaired dye efflux in the presence of ABC transporter inhibitor, with SP cells exhibiting high-ABC transporter activity and decreased staining compared to low-ABC transporter-activity NSP cells with uninhibited transporter activity [3,23,24]. Although the kinetic, ABC transporter-mediated mechanism is universally accepted as the basis for differential staining of cell populations in the +inhibitor and -inhibitor conditions of the SP assay, differences in cell staining due to transporter activity have not been demonstrated from a kinetic perspective nor at the single-cell level. This gap in knowledge persists due to a technical limitation that precludes an individual cell from Hoechst staining in both +inhibitor and -inhibitor conditions, thereby preventing any measurement of shifts of individual cells within the population distribution. Therefore, it is unclear how the heterogeneity of ABC transporter activity within a population influences staining characteristics in the +inhibitor and -inhibitor conditions of the SP assay and unclear how the heterogeneity in transporter activity is reflected in individual cells of the SP and NSP. We hypothesize that We developed an algorithmic approach for objective quantification of SP size using a Hoechst staining threshold to define the gate that delineates SP and NSP regions (S3 Fig). Although suggestions have been made for standardizing the reporting of SP assay results [1], gating protocols vary considerable in the literature. Here, we defined Hoechst staining intensity as the x-y projection of z-score transformed Hoechst Red and Blue raw signals with a threshold set at the 1 st percentile in the +FTC condition. This projection gating approach demonstrated a strong correlation (R 2 = 0.98, S4B Fig) between manual and automated algorithmic measurement of percentage SP.
To further elucidate the role of TGFβ on SP size, we made serial measurements of SP size using our projection gating method. A549 cells were expanded in culture for 4 days without TGFβ treatment and then passaged. On the subsequent day (Day 0), treatment with 0, 1, 10, and 100 pM TGFβ was started and SP size was measured at 2 day intervals ( S4A Fig). On the day of passage, Day 4, the SP constituted 20% of the cell population; however, following cell , and APC/ABCG2 (B) antibodies on live A549 cells following 4-day treatment with 0, 1, 10, and 100 pM TGFβ. Surface marker staining data were obtained from 3-color staining, with compensation, from 3 biological replicates. Values are plotted as the geometric mean ± 95% confidence interval of the sample geometric mean fluorescence intensities (GMFI). A significant difference from the 0 pM TGFβ condition was determined with a two-way ANOVA, p<0.05, and are indicated by asterisks (*). (C) Mean %SP versus ABCG2 expression for the Day 4 condition from the SP time course experiment, along with a best-fit line from a linear regression, R 2 shown, and 95% confidence interval. (D) SP size measured by projection gating with increasing time in culture (up to 6 days) and TGFβ. Performed in triplicate and plotted as mean ± SEM. passage, the SP consisted of 2% of the population (Fig 2D). In the following days, SP size increased and plateaued near 20% in the control condition. With increasing amounts of TGFβ exposure, the repopulation of the SP attenuated in a dose-dependent manner (Fig 2D).
We derived clonal populations of A549 cells by expanding individual low-and high-ABCG2 expressing cells to determine whether subclones within the A549 cell line define either the SP or NSP subpopulations. Cells were labeled with anti-ABCG2 antibody and 96 high- Collectively, these results indicate that the SP within the A549 cell line is dynamic with respect to cell phenotype (Fig 2A), cyclically variable in culture (Fig 2D), and is not genetically distinct from the NSP (S5 Fig).

SP & NSP Are Not Distinct Subpopulations
The study of side populations has been limited by the approach to define the SP, which relied on manually defined gates to define a boundary between SP and NSP regions within the Hoechst Red and Hoechst Blue plot. In the development of our automated projection gating approach to measure SP size, we closely examined numerous aspects associated with SP data.
A key observation was that the SP and NSP are not clearly separated within the-FTC condition. We observed a continuous staining distribution in both Hoechst Red and Blue channels of the-FTC condition, which is redistributed towards lower signal intensities compared to the +FTC condition (S6A and S6B Fig). The decrease of mean Hoechst signal intensity in the-FTC condition, compared to the +FTC condition, correlated with increased SP size (S6C Fig).
The intensity of both Hoechst Red and Blue staining is greatly reduced in the conventional binary assignment of cells into SP or NSP subpopulations. We sought to develop analytical methods to enable characterization of staining data from the SP assay with greater preservation of staining information. We processed the Hoechst Red and Blue flow cytometry signals into 2-dimensional population density functions (PDF) for both the +FTC (PDF +FTC ) and-FTC (PDF -FTC ) conditions, which converted raw data into normalized data sets (S3 Fig). This allowed us to directly compare spatial differences in staining intensity in the +FTC and-FTC conditions by subtracting the PDF +FTC from the PDF -FTC to define the ΔFTC density for a given sample (Fig 3A). Conversion of the raw data into PDFs was critical because it allowed for direct comparison of flow cytometry samples of un-equal event counts. An increased cell population density in the-FTC condition compared to the +FTC condition is visualized as a red signal whereas blue indicates decreased population density. By converting raw data into PDFs and thereby normalizing, we able to make quantitative comparisons in spatial staining intensities between +FTC and-FTC conditions of a particular sample. Furthermore, we are also able to compare the differences in staining redistribution (ΔFTC) between different samples by taking the difference between respective ΔFTC densities to compute the ΔSP density ( Fig 3B). These methods enable visualization of the influence of ABCG2 inhibition on staining while preserving the quantitative aspects of staining intensity and density.
We analyzed the influence of tert-butylhydroquinone (tBHQ) on ABCG2 activity and SP size using an imaging cytometer using these data analysis steps. Following 48 hours of 50 μM tBHQ treatment, A549 cells (4 days post-passage) were found to have and increased SP size from 6% to 9% (S7 Fig). The red signal within the ΔFTC plots reflected the presence of the SP in each sample; similarly, the red signal within the ΔSP plot reflected the increased SP size in the tBHQ-treated sample compared to the control (Fig 3C). Comparing differences in the TGFβ-treated samples shown in Fig 2 (S8 Fig), we observe decreased staining present in the-FTC condition without arbitrary segregation SP and NSP subpopulations from a continuous staining distribution.

Modeling Transporter Heterogeneity in an In Silico Population
We developed an approach for modeling of the SP assay to investigate the role of transporter heterogeneity within the cell population on the generation of SP and NSP responses of individual cells. At the cellular scale, we considered the influence of transporter activity in a massaction kinetic model to simulate Hoechst transport dynamics and cell staining. At the population scale, we considered the role of heterogeneity in the cell population, defined by variation in transporter properties and cell morphology. We defined an ensemble as the pairing of a particular set of kinetic constraints with all of the staining simulations for an in silico population. Having sampled M = 10,000 sets of kinetic rate constants (K), we have then had 10,000 ensembles with in silico staining results that we then compared to our experimental data to assess the ensemble for a SP response at the population level. We then analyzed the influence of transporter properties on staining results at the single-cell level for ensembles demonstrating a SP response.
Each of the cells in the N = 1,000 cell in silico population was uniquely defined by parameters for cell volume, cell membrane surface area, nuclear volume, nuclear membrane surface area, DNA content, and transporter heterogeneity (Fig 4A). Morphological parameters for each cell were sampled from corresponding experimental distributions using Latin hypercube sampling (LHCS), which ensured that the resulting in silico population faithfully represented the experimental distributions ( S9 Fig and S10 Fig). Staining of the in silico population was carried out for each individual cell under 4 different levels of transporter expression, T 1 -T 4 , corresponding to distributions derived from 0, 1, 10, and 100 pM TGFβ-treated cells (S10 Fig), under +FTC and-FTC conditions. Single-cell staining simulations consisted of a compartmental mass-action kinetic model of intracellular dye transport processes, governed by kinetic rate First, a heterogeneous in silico cell population was generated from experimental distributions, which was then paired with each of the kinetic parameter sets that were derived from Latin hypercube sampling of parameter space. The cell population was implemented in 4 different transporter-variant versions, where relative transporter expression was derived from experimental data. An ensemble is defined as a pairing of a particular parameter set with the variety of populations for which it is simulated. Three different models of transporter activity heterogeneity were implemented. (B) For each ensemble, each of the 4 transporter-variant populations were underwent simulated Hoechst staining at the single-cell level both with (+FTC) and without (-FTC) transporter inhibition. Using in silico flow cytometry, the Hoechst concentrations following kinetic staining simulation were converted into Hoechst Red and Blue signals constants for the given ensemble and the morphologic parameters associated with the cell ( Fig  4A and Fig 1C). Following simulation of each of the cells in a population, the final dye concentrations of the samples were processed via in silico flow cytometry to determine Hoechst Red and Blue signals (Fig 4B and S11B, S11C and S11D Fig). Simulated +FTC and-FTC conditions were used to measure SP size, PDF +FTC , PDF -FTC , ΔFTC, and ΔSP distributions for each of the 4 transporter expression levels (T 1 -T 4 ; Fig 4C), which were then compared to corresponding values from the experimental data ( Fig 4D). Ensembles were determined to have an SP response if each of the following criteria were met: 1) the change in mean Hoechst Red Signal (ΔHRS mean ) and Hoechst Blue Signal (ΔHBS mean ) were negative when comparing the-FTC to the +FTC condition (indicating decreased staining in the-FTC condition); 2) PDF +FTC , PDF -FTC distributions overlapped (indicating the presence of a NSP); 3) the cross-correlation between experimental and in silico ΔFTC and ΔSP distributions was positive (indicating similarity in spatial response distributions); and 4) in silico SP sizes correlated with experimental SP sizes for the 4 transporter expression levels (indicating a consistent decrease in SP size with decreasing transporter expression levels). Ensembles producing a SP response were further analyzed at the single-cell level to investigate the influence of transporter properties on SP/NSP status and differential staining the +FTC and-FTC conditions. We independently implemented 3 different modes of transporter heterogeneity within the in silico cell population to assess the role of transporter heterogeneity in the formation of SP responses in the SP assay. In the first mode of transporter heterogeneity (I), the concentrations of transporter across the population is uniform, which due to the heterogeneity in cell size, leads to a distribution of transporter numbers within the population. In the second mode of heterogeneity (II), the number of transporters per cell across the population is uniform, generating a distribution of transporter densities within the population. For the first two modes, uniform number and densities of transporters for the T 1 -T 4 expression levels were determined by mean values of ABCG2 expression from flow cytometry data ( Fig 2B). We explicitly designed the third mode of heterogeneity (III) to be more expansive than the first two modes. In this mode, transporter expression within the populations were sampled from experimental distributions of ABCG2 expression (S10A Fig) as part of the LHCS approach in generating the in silico population. Further, we enabled a non-linear relationship between transporter expression levels and transporter activity levels at the single cell-level. This mimics cooperativity in which higher transporter levels led to greater transporter activity than would be accounted for by a simple 1-to-1 correspondence between expression and activity. We implemented this modeling approach to assess the relationship between transporter heterogeneity and SP responses for the 3 modes of transporter heterogeneity.

Heterogeneity of Kinetic Transport Processes Generate SP Responses
In our modeling approach, the primary aim was to identify the ensembles in which an SP response resulted. Again, an ensemble consisted of the simulated staining of the in silico cell population using a specific set of governing kinetic rate constants (K j ) for a particular mode of using a linear transformation that incorporates spectral excitation and emission properties of simulated Hoechst dye and flow cytometers. (C) In silico flow cytometry signals were converted into Hoechst Scores and Hoechst Score PDFs. Projection gating objectively assess SP size. Hoechst Score PDFs were used to measure ΔFTC and ΔSP distributions. (D) Hoechst staining metrics of Scores, PDFs, ΔFTC, and ΔSP were used to identify models exhibiting SP responses by comparison to corresponding experimental metrics. Ensembles meeting qualitative selection criteria are then scored according to their similarity in %SP to experimental data. Ensembles demonstrating SP responses were then analyzed on a single-cell level and the distribution of single-cell SP responses analyzed.
transporter heterogeneity in both +FTC and-FTC conditions for each of the T 1 -T 4 transporter expression levels. For each mode of heterogeneity, we assessed M = 10,000 kinetic rate constant sets (K) for an SP response with less than 5% of the simulation ensembles aborted due to timeout for long simulation time and with SP responses only identified in a small subset,~5%, of successfully completed ensemble simulations (S15 Fig). In addition, the range of SP responses for each of the modes exhibited a similar quality of fit to experimental SP size (S16 Fig) and a wide range of kinetic rate constants were permissive of SP responses (S17 Fig). An increased frequency of ensembles with SP responses was observed in the concentration (I) and number (II) modes of heterogeneity with larger transporter expression slopes (k 8 ; S17A and S17B Fig), in which a larger k 8 value corresponds to greater discrepancy in expression between T 1 -T 4 conditions. Similarly, in the experimental distribution mode of heterogeneity (III), an increased frequency of SP responses resulted in ensembles with greater non-linearity (k 8 , k 9 ; S17C Fig).
Having met the selection criteria, ensembles exhibiting an SP response were analyzed at the single-cell level for features.
Experimentally, SP responses are apparent when a subpopulation of cells stains less intensely in the-FTC condition compared to the +FTC condition. Using our ensemble modeling approach, we were able to simulate SP responses in heterogeneous in silico cell populations. The simulated flow cytometry data generated during the simulations was analyzed using the same approaches for experimental SP data, yielding outputs of SP size, PDF +FTC , PDF -FTC , ΔFTC, and ΔSP distribution data for the T 1 -T 4 conditions ( Fig 5A). We observed subtle differences in the staining patterns of SP cells, in which some responses exhibit relatively fewer SP cells with a more significant decrease in cell staining (subsequently defined as a Subpopulation Response). In contrast, other responses had more SP cells with a less significant decrease in cell staining (subsequently defined as a Full Response) between +FTC and-FTC conditions. Experimentally, such differences could only be investigated at the population level due to technical limitations; however, we designed our modeling approach to circumvent these limitations.
A key advantage to our modeling approach was in the consistency of the cell population for staining simulations, in which an identical cell population could be seeded for both +FTC and-FTC conditions. Furthermore, the same in silico population was consistent within an ensemble, with the exception of transporter expression by T 1 -T 4 conditions, and from ensemble to ensemble, differing only by the respective set of governing kinetic rate constant sets. This permitted a level of comparison that would be impossible in an experimental setting. For example, we were able to reduce Hoechst Red and Blue scores onto an x-y projection, as is done in the projection gating method, to obtain a single staining metric for each cell. Because the same in silico cell population was used for both +FTC and-FTC conditions, we were then able to compute the difference in Hoechst staining projection scores (-ΔH proj ) on a cell-by-cell basis. A larger -ΔH proj corresponds to a greater decrease in staining in the-FTC condition compared to the +FTC condition, thus indicating a greater influence of transporter activity on staining intensity.
Next, we investigated the relationship between the distribution of single-cell -ΔH proj responses with a population and the staining pattern of SP cells for a given ensemble. In plotting histograms of -ΔH proj values we were able to appreciate patterns in responses ( Fig 5B). In one pattern, we observed a -ΔH proj distribution with the greatest number of -ΔH proj values near zero and a tail of increasing -ΔH proj values (Fig 5B top). This indicates that transporter inhibition had little to no effect on staining in a majority of the cells but greatly influenced staining in a portion of the cells. We described such a scenario as a Subpopulation Response. Contrastingly, in the Full Response, the entire population exhibits non-zero -ΔH proj values and a more normal-like distribution ( Fig  5B bottom), indicating that staining of each cell in the population was affected by inhibition of transporter activity as well as consistency in the magnitude of this effect.
The nature of our modeling approach enabled us to characterize a great number of relationships not possible experimentally. We compared the influence of transporter inhibition on staining intensity (-ΔH proj ) against staining intensity in the-FTC condition (H proj -FTC; Fig  5C). In the Subpopulation Response, we observed SP cells with little influence of transporter inhibition on staining as well as NSP cells with significant influence of transporter inhibition (Fig 5C top). In the Full Response, we observed NSP cells with greater influence of transporter inhibition on staining (Fig 5C bottom). Similar comparisons were made to +FTC conditions and by transporter numbers (S19 Fig).
While both were capable of generating SP sizes consistent with experimental data, the Subpopulation Responses (S13 Fig For most ensembles, the -ΔH proj response distributions were not clearly Subpopulation or Full Response in nature, rather most exhibited aspects of both. To more objectively characterize the -ΔH proj responses, the distributions were evaluated by standardized skewness and bimodality coefficient (Fig 6). Standard skewness is a metric that quantifies asymmetry of the distribution while the bimodality coefficient quantifies the similarity between the distribution and a purely bimodal distribution (S20 Fig). Both transporter number (I) and concentration (II) modes of heterogeneity had a greater tendency to generate -ΔH proj distributions more closely resembling the Full Response with greater distribution symmetry and bimodality coefficients similar to a normal distribution (blue and red; Fig 6). Ensembles with SP responses generated in the experimental distribution of transporter heterogeneity (III) were more varied, ranging from Full Responses to Subpopulation Responses (green; Fig 6). The overall spread of distribution mappings in this plot demonstrates a consistency of SP responses arising from -ΔH proj distributions that are roughly normal and have a right-sided tail of variable magnitude.

Kinetic Determination of SP Response
SP cells arise through ABC transporter activity; however, the technical limitations that prevent staining of specific cells in both inhibitor-free and inhibitor-containing conditions have prevented direct measurement of the influence of transporter activity on Hoechst staining a the single-cell level. Our novel computational approach to in silico staining of a cell population in numerous conditions, differing only by transporter activity, revealed that specific distributions of kinetic transport heterogeneity within a population can alone account for the formation of a SP, independent of genetic or other phenotypic traits, such as CSC. This biological insight was obtained by the pairing of flow cytometry distribution data with a novel computational approach for modeling dye kinetics at the single-cell level. This builds on similar approaches by You et al that demonstrated cell-to-cell variability and heterogeneity in gene expression can govern biphasic responses to extracellular cues [27] and that variability of bacterial uptake by hosts are governed in a probabilistic manner determined by host receptor expression levels [28]. It is increasingly understood, even within clonal/isogenic populations, that noisy or variable gene expression leads to heterogeneity within a population and can contribute to differences in phenotype [29][30][31][32]. A role for noisy gene expression in generating diverse phenotypic responses in clonal population has been supported both experimentally and computationally; further, it has been postulated as a mechanism of "bet hedging", thereby conferring survival advantages in stressful or alternate environments [33][34][35][36][37].
The kinetic aspects involved in the SP assay caused us to more carefully consider the significance of SP/NSP discrimination. SP/NSP phenotype are not strictly defined by inheritance or genetic factors, as is indicated by the fact that SP and NSP cells emerge in clonal populations derived from single-cell sorting (S5 Fig) and by reports that describe the interconversion between SP and NSP phenotypes from isolated cell populations [20,38,39]. Further, in our modeling, we observed a range of transporter activities in SP and NSP cells (Fig 5C), suggesting that not all SP cells have high transporter activities and not all NSP cells have low transporter activities. This can be attributed to the interplay between heterogeneity in transporter activity within the population as well as relative differences in DNA content across the population, which alters that staining potential in the SP assay [40]. Thus, the SP phenotype of a specific cell is not discretely defined by membership in a homogenous subpopulation; rather it is a kinetically defined property, inextricably defined by the specific experimental conditions.

Variability in SP Size
We observed SP sizes ranging from 2% to 20% in the A549 cell line (Fig 2D), which is consistent with reported SP sizes of 2.5% to 30% for A549 cells in the literature [21,38,[41][42][43][44][45]. We observed consistent measures of SP size with reproducible trends with increasing time/passage in culture and decreasing exposure to TGFβ ( Fig 2D). It should be noted that SP size is determined by the specific experimental conditions in the SP assay, the analytic technique employed by the investigators [1,[46][47][48], and the conditions in which samples are maintained in culture [39,45]. Additional variability may arise from discrepancies in gating SP/NSP regions in cytometry data [21,38,[41][42][43][44][45], which is highly subjective and inconsistent and complicated by the fact that SP and NSP staining regions are continuous with no self-evident border of separation (S6A and S6B Fig).
Proper gating to exclude erroneous events is key to accurate SP measurement as cellular debris or dead cells may be mistaken as SP cells. Apoptotic cells, for instance, lie near the SP staining region with greater decreased Hoechst red signal than Hoechst blue signal [1]. Failure to properly exclude these events could result in misclassification of this population as SP when using manual gating methods thereby falsely elevating SP size. The automated SP measurement we described provides an objective measure of the SP without user bias; however, the reliability of the output, %SP, is highly dependent on the input, Hoechst red and Hoechst blue intensities of cytometry events. To ensure the quality of our input we sequentially excluded cellular debris, multiple events, and dead cells in the gating tree to define the final Hoechst Red and Hoechst Blue stained cells as recommended by Golebiewska et al. [1]. Further method development may be necessary to consider scenarios with large apoptotic populations in order to properly exclude such events from the final analysis.
Encouragingly, the automated methods we have described may have potential to measure SP size despite the presence of apoptotic cells as these cells would be expected to stain similarly in -inhibitor and +inhibitor conditions. While an apoptotic population may underestimate SP size by falsely lowering the 1 st percentile level in the projection gating method, the staining differences between +inhibitor and -inhibitor conditions would be negated in the ΔFTC distribution. Therefore, the Hoechst staining density distributions may be further developed for SP measurements. While automated methods hold promise for objective measurement of SP size, we emphasize that additional studies will be necessary to establish validity of such measurements in the presence of errant populations.
Variability in biological, experimental, and analytical conditions likely account for the much of the variability in literature-reported SP sizes for the A549 cell line, which highlights the importance of reporting detailed experimental conditions and calls to question the utility of comparing SP size in different experimental and analytical conditions. For example, SP size in untreated A549 cells 4 days after passage had SP sizes of 20% ( Fig 2D) and 6% (S7C Fig) when measured using a BD LSR II flow cytometer or Amnis FlowSight imaging cytometer, respectively, despite having identical Hoechst staining conditions and SP data analysis. Instrument settings such as excitation, emission, Hoechst Red and Blue channel windows, detector sensitivity, and detector scaling all greatly influence the recorded Hoechst signals.

Caveats to the Kinetic Interpretation of SP/NSP Identity
The kinetic interpretation of the SP may be applied to some immortalized cell lines, such as the one used in this study; however, SP may not arise strictly from kinetic properties in all situations, especially in samples of multicellular composition in situ. For example, in the original description of the SP, the SP corresponded to hematopoietic stem cells within bone marrow aspirate [2] and the SP in glioblastoma-derived samples were tumor stromal cells while the glioblastoma cells, including glioblastoma CSCs, did not contribute to the SP [49]. Our investigation adds to this body of evidence as it validates an ABCG2-dependent kinetic basis for the formation of the SP and provides a novel perspective on the distribution of single-cell ABCG2 activity across a population as well as within the SP itself.
Independent of cell line, SP size is also a function of passage frequency, Hoechst staining conditions, fluorescence excitation/emission settings, and SP/NSP gating strategies. Therefore, to focus our investigation on the kinetic aspects of the SP assay, we performed the SP assay with attention to consistency of the culture conditions, Hoechst staining conditions, cytometry settings, and unbiased SP analysis. While this approach enabled investigation of assay kinetics based on differences in ABCG2 expression, other kinetic factors, such as Hoechst staining concentration and duration [46,48], influence SP measurement. A more refined and comprehensive modeling approach would be necessary to determine precise kinetic rate parameters to more definitively characterize the transport phenomena generating the SP.

The SP Phenotype: Heterogeneity in Therapeutic Targeting
The SP has been the subject of many investigations, many with conflicting results. Some studies claim the SP to be a population of CSCs [20,50], or that ABCG2 is necessary to maintain a stem cell phenotype [51]. In contrast, other studies have identified stem cells lacking a SP phenotype or ABCG2 expression [49,[52][53][54] or identified SP lacking stem cell phenotypes [55][56][57][58].
In A549 cells, TGFβ simultaneously leads to enhanced expression of CSC properties [38,43,59] and decreased SP size (Fig 2B) [38]. These findings suggest that SP and CSC phenotypes are nonequivalent; however, the phenotypes appear to be frequently co-expressed.
The co-expression of SP and CSC phenotypes may have a significant role in tumorigenesis as ABCG2 expression can confer a MDR phenotype to CSCs [60]. Additional cellular processes, related to ABCG2 transporter activity, may contribute to resistance in the MDR phenotype. For instance, activity of Nrf2, a master regulator of the cellular redox environment is also a key regulator of ABCG2 expression [61][62][63][64]. ABCG2 has been implicated in antioxidant processes [65][66][67][68][69][70][71] and TGFβ signaling, down-regulates ABCG2 expression (Fig 2B) in addition to down-regulating multiple antioxidants [38,[72][73][74][75][76]. Furthermore, population-level studies have demonstrated increased antioxidant expression in SP cells compared to NSP cells. [77] Therefore, the intracellular redox processes may act in concert with ABCG2-mediated drug efflux to promote the MDR phenotype. Targeting functional properties specific to high ABCG2-expressing cells may be a promising approach to overcome a MDR phenotype. For instance, addition of the tyrosine kinase inhibitor axitinib to topotecan enhanced topotecan-mediated apoptosis in A549 cells through inhibition of ABCG2-mediated transport, independent of ABCG2 expression. [78] Despite the fact that numerous chemotherapeutics are substrates of ABC transporters, evaluation of ABC transporter inhibitors in clinical trials has failed to demonstrate added benefit due to drug toxicities and the inability to achieve sufficient concentrations to effectively inhibit transporters [79]. Further investigation of the kinetic processes involving ABCG2 will refine our understanding of the functional consequences of ABCG2 activity in cancer cells and potentially inspire novel chemotherapeutic approaches.

Conclusion
This investigation is the first to demonstrate the kinetic mechanisms that form the basis of the SP assay. We validated the ABCG2 transporter-mediated differential staining between +FTC and-FTC conditions across multiple transporter concentrations in a heterogeneous cell population, accounting for differences in responses for SP and NSP cells. Our modeling approach leveraged these experimentally determined dynamic distributions of ABCG2 expression levels as a function of TGF-β treatment and culture time. The computational model tested influences of transporter properties on cell staining across a heterogeneous population, which would otherwise be impossible to achieve due to the technical limitations of the SP assay. Our results suggest that in particular distributions of transporter kinetics within a population, a subset of cells within the population exhibit marked enhancement of transporter activity compared to the main cell population. Analysis of thousands of single cell simulations provided unique insight that NSP and SP cells both lie along a spectrum of ABCG2 activity. Collectively, these results support our hypothesis that specific single-cell distributions of ABC transporter activity yield differential staining and serve as a kinetic basis in forming side populations.

Side Population Assay
Hoechst staining and flow cytometry to measure the SP in A549 cells closely followed the protocol described in [23,48]. Cells were trypsinized and resuspended in CO2 conditioned DMEM +, consisting of high-glucose DMEM without phenol red, 10 mM HEPES, 2% FBS, and 2 mM EDTA, at a concentration of 1x106 cells/ml. Samples were split into +FTC and -FTC conditions, supplemented with DMSO-mobilized FTC (EMD Millipore) at a final concentration of 10 μM or DMSO alone, respectively. The solutions were incubated in a 37°C water bath for 30 minutes, after which they were supplemented with Hoechst 33342 (Life Technologies) at a final concentration of 5 μM for 90 minutes, with mixing at 30-minute intervals. The staining solutions were then centrifuged at 1,000 RCF for 10 minutes and resuspended in HBSS+, consisting of HBSS without phenol red, 10 mM HEPES, 2% FBS, and 10 mM EDTA. Cells were incubated with the viability stain SYTOX Blue (Life Technologies, 1:1000) for 5 minutes prior to fluorescence measurement via flow cytometry. Positive controls for dead cell staining were obtained by incubating cells at 56°C for 45 minutes followed by SYTOX Blue staining. The BD LSR II flow cytometer was monitored using fluorescent beads to ensure optimization of system optics and detectors for quality control in polychromatic settings. Samples were excited with a 355 nm UV laser and the Hoechst Red signal was measured in a λem = 675/50 nm channel with linear scaling and the Hoechst Blue signal measured in a λem = 450/50 nm channel with linear scaling, collecting 100,000 events per flow sample. The Hoechst stained -FTC condition was initially used to tune Hoechst Red and Hoechst Blue photomultiplier tube settings, which were held constant for all subsequent studies. Additionally, samples were excited with a 445 nm violet laser with SYTOX blue emission measured in a λem = 473/10 nm channel with logarithmic scaling. Sample gating proceeded as follows: 1) Debris exclusion with FSC-area/SSCarea (λex = 488 nm) gating; 2) Single-cell selection with FSC-height/FSC-area; 3) Live cell selection with the violet-λem = 473/10 nm channel. Events retained through all 3 gates were used for subsequent SP analysis in the Hoechst Red and Blue channels. Manual selection of SP gates was determined using the +FTC conditions where a quadrant gate was placed as tight as possible such that greater than 99% of the cells in the +FTC condition were located in the upper right quadrant. The same gates were then applied to the -FTC condition where the two left gates were considered to be SP gates and the right two gates considered to be NSP gates. The measured %SP in the manual gating approach is the sum of the percent of parent population in the SP gates. For a given sample, the %SP was determined using the specific +FTC and -FTC conditions for that sample.

Hoechst Staining Signal Processing
Hoechst Score Transformations. Side population flow cytometry data for a particular sample consists of individual events with associated Hoechst Red and Blue signals for both +FTC and -FTC conditions. The Hoechst Red and Blue signals are expressed in independent arbitrary units. To generalize the interpretation of side populations from Hoechst signals, independent of raw signal units, Hoechst signals were converted into Hoechst Scores. Hoechst Scores are based upon the standard score, or z-score, in which a distribution is mean-centered and normalized to the standard deviation. To compute the Hoechst Score for the Hoechst Red channel, the mean and standard deviation of the Hoechst Red signal from the +FTC condition were calculated. Next, the +FTC Hoechst Red mean was subtracted from the Hoechst Red signals from each of the events in the +FTC and -FTC conditions. Similarly, each event in the +FTC and -FTC conditions were divided by the standard deviation from the +FTC condition. The resulting event data, +FTC mean centered and +FTC standard deviation normalized, constituted the Hoechst Red Scores for the two conditions. Hoechst Blue Scores were derived in an analogous fashion. Hoechst Scores were calculated at a per sample basis between each pairing of +FTC and -FTC condition data.
Projection Gating for %SP Measurement. Hoechst Scores Projections were derived from Hoechst Red and Blue Scores. Projections were derived from Score data and not signal data due to the arbitrariness of signal magnitude in the signal data. Score data from Hoechst Red (HRS) and Blue (HBS) channels have common units and relative magnitudes. Projection values were derived for each event within a sample, in which the Hoechst Score Projection (H proj ) was defined: Hoechst Projections were used to set a threshold, or gate, intensity at the lower limit of the NSP. The 1 st percentile mark of the Hoechst Projection data from the +FTC condition was used define the threshold and applied to the -FTC condition. The percent of events falling below the threshold in the -FTC condition was set as the %SP in this projection gating approach.
Hoechst Scores PDF Distributions. For a given Hoechst condition (+FTC or -FTC), Hoechst Red and Blue scores were provided transformed into a 2D probability density function (PDF) on the Hoechst Red and Hoechst Blue plane with the frequency, or cell density, defined at each paired Hoechst Red and Blue coordinate. Smoothed surfaces over the Hoechst Red and Blue plane were derived by the method for smoothing scatter plot data described by Eilers and Goeman [80]. Next, the area under the surface was calculated and normalized to one for each condition.
In its most basic form, flow cytometry data is a set of coordinate data/scatter, with each event represented by an intensity value along each of the measured dimensions. Comparison of data between two sets relies on the comparison of some sort of statistical transformation of the data (i.e. mean or median). However, such transformations result in loss of spatial information. Histograms allow for comparison of spatial information, but differences in event number complicate interpretations. In order to permit more rigorous comparisons, we converted Hoechst staining coordinate data into probability density functions along both Hoechst Red and Blue Score dimensions. This was similar, in effect, to constructing a 2D histogram with smoothing and normalization to account for differences in event number between data sets. In this format, spatial differences in staining distributions were easily computed by taking the difference between PDFs.
ΔFTC Distributions. For a given sample, a ΔFTC distribution was calculated by taking the point-wise difference between the PDF of the -FTC condition (PDF -FTC ) and PDF +FTC across the Hoechst Red and Blue plane to calculate the difference in normalized cell density.
ΔSP Distributions. To compare the difference between ΔFTC distributions between two samples we calculated a ΔSP distribution. To compare a test sample to a control sample the ΔSP was derived as ΔFTC test -ΔFTC control , where the point-wise difference in ΔFTC intensities was calculated at each Hoechst Red and Blue pair.

FlowSight Imaging Cytometer
We imaged samples at 20X magnification using a FlowSight imaging cytometer with the Quantitative Imaging Upgrade (Amnis, Seattle, WA). Single color controls were used to set-up compensation matrices. For the SP assay, the compensation matrix was manually edited to allow collection of the Hoescht Blue (470/35 nm) and Red (694/51 nm) signals using the 405 nm laser. Images were analyzed with IDEAS analysis software (Amnis). Using the gradient root mean square feature for the brightfield channel, "Focused cells" were selected according to the manufacturer's recommendation. Debris was eliminated by gating single cells using the area and aspect ratio features for the brightfield channel. Live cells were gated using the intensity feature in the green channel (533/27 nm) for SYTOX Green staining. For ABCG2 analysis, the intensity feature of the APC channel (694/51 nm; 642 nm excitation laser) was used to quantify the expression of ABCG2. For the side population assay, double positive cells from were selected by gating in the Hoechst Red and Hoechst Blue channels.

Ensemble Modeling of Side Population Responses
General Overview of Approach. Heterogeneous cell populations were simulated in an array of experimental conditions across a wide range of kinetic conditions to investigate the influence of transporter function on Hoechst staining kinetics that give rise to SP phenotypes. Three models were generated describing numerical (Model 1), concentration (Model 2), and experimentally derived (Model 3) transporter distributions. Individual cells (P i ) within a population (P) of size N = 1,000 are described by a set of morphological parameters (volumes, surface areas, & DNA content, Table 1). The parameter values for the population were assigned via LHCS of PDFs derived from experimental distributions of cell radii, nuclear radii, and DNA content. Kinetic parameter sets (M = 10,000) were obtained from LHCS of uniform distributions in log space (Fig 4A). For each parameter set, Hoechst staining is simulated in the cell population across multiple transporter conditions and with and without transporter inhibition, simulating +FTC and -FTC conditions. Each population simulation consists of N singlecell mass-action ODE simulations of Hoechst staining (Fig 4B). Following the kinetic simulations, Hoechst concentrations within individual cells were converted to Hoechst Red and Blue signals via linear transformation with a signal matrix accounting for the spectral excitation and emission properties of free and DNA-bound Hoechst dyes as well as an in silico flow cytometer (Fig 4B). Simulated flow cytometry signals, in arbitrary units, were then converted to Hoechst Score PDFs where projection gating was applied to measure the %SP (Fig 4C). Hoechst scores PDFs were used to calculate ΔFTC and ΔSP distributions, allowing for visualization of differences in population simulations with and without transporter inhibition as well as across different transporter conditions (Fig 4C). Hoechst score metrics from in silico flow cytometry populations were compared to metrics derived from experimental data to gauge the extent of SP response in the populations. Parameter sets with in silico populations meeting the selection criteria for identifying a SP were then accepted and ranked according to similarity to %SP measured experimentally according to the normalized root mean-squared error (RMSE, Fig 4D). Accepted sets were then analyzed at the single-cell SP response. The distribution of responses to inhibition were used to classify the homo/heterogeneity of response magnitudes and the uniformity/bimodality of response frequency (Fig 4D).
Whole Cell and Nuclear Radii Probability Density Functions. FlowSight imaging cytometry was used to measure whole cell and nuclear radii of Hoechst stained cells in the presence of FTC as per the imaging cytometer instructions. Cells and nuclei had aspect ratios > 0.9 and were assumed to be spherical for the purposes of simplification. Whole cell and nuclear areas, reported in micrometers, were then used to estimate cell and nuclear radii for the population of A549 cells (S9A Fig). The size of the nucleus was partially correlated to the whole cell size (S9B Fig); therefore, instead of treating each as in independent distribution, a 2D PDF was constructed in a manner analogous to the used to calculate PDFs for Hoechst Red and Blue Scores.

Different Transporter Expression & Distributions.
Transporter Expression: Different models of transporter distribution across a population were implemented; however, at the single-cell level, the Hoechst staining kinetic model was identical. Therefore, what differed between the models was how different cells within a given population were assigned transporter expression. In each of the models, experimentally-derived expression levels were used to inform model expression. In Models 1 & 2, the relative geometric mean expression of ABCG2 from TGFβ -treated A549 cells (Fig 2B) were used while in Model 3, the flow cytometry staining distribution served as a probability density function from which transporter expression frequency was sampled (S10A Fig). Model 1. Number Distribution/Equal Concentration: Each population was associated with the geometric mean as the relative transporter level for the population. This factor served as the scaling factor relative to the maximum geometric mean intensity from the untreated control condition. The magnitude of the relative differences between transporter levels of the samples was set by k 8 . The relative differences of transporter levels were reconfigured for each kinetic parameter set. Once the relative transporter level was determined, it was set as the maximum transporter concentration in the cytosol or nucleus where the relative intensities of cytosolic to nuclear transporter activities were expressed as: Each cell within the population was assigned total cytosolic and nuclear transporter activities of T CTA & T NTA , respectiviely.
Model 2. Concentration Distribution/Equal Number: Transporter activity for Model 2 proceeded in the same manner as Model 1; however, after calculation of T CTA & T NTA , the average number of molar equivalents in cells of the distribution was calculated. Each of the cells in the population were then assigned the same number of molar activity equivalents, which was then converted to concentration activity equivalents on a cell-by-cell basis using the cytosolic and nuclear volumes.
Model 3. Experimental Distribution: Flow cytometry ABCG2 surface marker staining data from TGFβ -treated A549 cells were exported from FlowJo as compensated fluorescence intensities. The distributions were then loaded into MATLAB where they were converted to PDFs for each individual replicated. PDFs were generated with the ksdensity function. For a particular sample, a final PDF was taken as the unit normalized average PDF of three experimental replicates. The distributions were then scaled to fall between 0 and 1 (S10A Fig). Therefore, values within the distribution reflect relative expression within the distribution.
Relative transporter levels within a distribution (T Total ) were randomly selected values (P Ti ) from ABCG2 expression PDF distributions where i = [1,4], corresponding to ABCG2 distributions from 0, 1, 10, & 100 pM TGFβ treatments, respectively. T Total levels are expressed in units of T Level . Next, the relative transporter levels in cytosolic (T CT ) and nuclear (T NT ) compartments were calculated under the assumption that the total transporter level is split into cytosolic and nuclear compartments at a fixed ratio, which remains constant during the simulation.
Total Transporter Level: Cytosolic/Nuclear Transporter Levels: Total transporter activity levels were then calculated from the Hill equation (θ c ), reflecting the cooperative interactions of transporters within each compartment: Total Cytosolic Transporter Cooperatively: Total Nuclear Transporter Cooperatively: Transporter Cooperatively: Transporter Activity Level Inhibition. Prior to kinetic modeling of Hoechst staining, each cell was assigned a total cytosolic and nuclear transporter activities, T CTA & T NTA . The absolute transporter activity in the kinetic simulation was then determined by the absolute scaling factor k 6 as well as the degree of transporter inhibition (i T ), where i T = 0.99 in the inhibited condition (+FTC) while i T = 0 in the uninhibited condition (-FTC).
Total Cytosolic Transporter Activity: Total Nuclear Transporter Activity: Hoechst/DNA-Binding Site Expression Probability Density Function. DNA content distributions within cell populations can be measured with Hoechst staining. [81] Therefore we took the +FTC Hoechst stained samples to represent the relative distribution of DNA content within a population. Hoechst Blue signals for each of the +FTC conditions in the SP time course study were loaded into MATLAB and converted into a PDF using the ksdensity function. An overall distribution was constructed from the average of the 39 individual PDFs (S10B Fig). The distribution was then normalized to the mode so that the distribution represented a distribution relative to the mode. The mode was assumed to represent a cell in the G 0 /G 1 phase and have a relative DNA content of the size of 1 genome for an aneuploid A549 cell. We used this distribution to estimate the number of Hoechst-binding sites in DNA (S10B Fig). For a given cell, the relative DNA level sampled as DNA L with units of DNA Level . Within the population, each of the DNA intensities was identically scaled, though maintaining their relative distribution, to determine absolute binding sites and converted to molar binding sites. Finally, binding site number was factored by nuclear volume to derive a molar concentration of binding site number.
DNA Total ¼ DNA L Á ðGenome SizeÞ Á ðBase pairs per siteÞ Á k 2 =AvgN Sampling PDFs to Construct In silico Cell Populations. To sample PDFs and produce in silico populations, PDFs were converted, approximately, to cumulative distribution functions (CDF) by taking the cumulative sum of a PDF. The CDFs were then normalized to a range of 0 to 1. To generate a population of N cells from the PDF, N random numbers were drawn from the interval of 0 to 1 and mapped to the CDF to find the corresponding expression value. In our observation, random sampling over the entire interval produced highly variable effects. To circumvent this issue, we implemented Latin hypercube sampling (LHCS) of the CDF, which more uniformly sampled the distribution. To sample the radii, Nx2 random numbers between 0 and 1 were generated. The whole cell PDF was sampled as an independent PDF. Next, the nuclear PDF for the given whole cell radii was sampled to sample the conditional PDF for nuclear radii size. Reconstructed cell populations of various sizes are shown in S9C Fig. After the radii were sampled for a given cell, the radii were used to calculate cell and nuclear volumes and surface areas, assuming spherical morphology. Cytosolic volumes were taken as the difference between whole cell and nuclear volumes. Notably, the same LHCS vector was used to sample the different CDFs of transporter expression. Therefore, the sampling of transporters between populations is consistent.
A549 Genome Size = 7.3x10 8 Base Pairs [82] Base Pairs Per Hoechst Binding Site = 80 [83] Binding  Table 2) were assigned via LHCS, which segments a dimension of parameter space into uniform segments. Within each segment, a parameter value is selected from a uniform random distribution. Thus, LHCS generates a collection of randomly chosen parameter choices with nearly uniform sampling of the parameter space. Within MATLAB, the LHCSdesign function was used to generate an MxQ LHCS matrix of M samples within the interval (0,1) for each of the L parameters (Models 1 & 2, Q = 8; Model 3, Q = 9). The criterion correlation and maxmin were enabled and 50 iterations were permitted to reduce correlation and maximize point-to-point distance within the LHCS matrix.
To convert the LHCS of the parameter space ranges in Log 10 space for a given parameter k q , the q th column of the LHCS matrix was scaled by the Log 10 of the range size and increased by Log 10 of the lower limit of the range. Finally, the parameter value k mq was obtained by taking the antilog of the m,q th entry of the transformed LHCS matrix. In Model 3, the Hill Half-Maximal Level, k 9 , was sampled uniformly from 0 to 1.
Initial modeling included k off within the parameter search space; however, early interrogation of the system demonstrated insensitivity to variation in k off . Therefore we maintain a fixed k off relative to k 1 in all ensembles based on the reported K D of 10 −7 . [83] Kinetic Model of Hoechst Staining: Simulation of Hoechst staining took place at the single-cell level. Within each population, cells were assigned variable cell and nuclear sizes, DNA content, and, in Model 3, relative transporter levels. Across the set of kinetic parameters were common across the entire population. Hoechst staining within a single cell was modeled using mass-action kinetics to describe the rates of reaction the transport across plasma and nuclear membranes (S11A Fig). Each single cell system was modeled with three spatial compartments and was simulated with 90 minutes of staining.
In each simulation, cells were initialized with no Hoechst species within the cell. The extracellular compartment was assumed to be so large so as to not experience changes in Hoechst concentration throughout the simulation. We assumed total DNA binding sites, cytosolic transporter, and nuclear transporter levels were conserved during the time course of the staining and, using conservation of mass, we algebraically reduced the order of the system, setting Table 2. Single-Cell ODE Model Kinetic Parameters. Sets of random kinetic parameters were assigned using LHCS to uniformly sample the parameter range in Log 10 -spaced intervals. The set of kinetic parameters is a uniformly applied across a population.

Parameter Symbol Range Units
Hoechst-DNA Association Rate the order of the system at 5 differential variables ( Table 3). The set of kinetic reactions were used to compose the set of differential equations, which governed the dynamics of Hoechstassociated species within the kinetic model. For each single-cell simulation, variables were assigned initial conditions and submitted with system reaction equations to the ode15s solver in MATLAB. Constant variables remain unchanged during the course of the simulation. Algebraic variables were derived from constant and differential variables using algebraic conservation equations at each time point in the solver. Differential variables were solved at each time point in the solver according to the set of differential equations.
Mass Conservation Equations: During a simulation of Hoechst staining in a single cell, the amount of transporter and Hoechst/DNA-binding sites are assumed conserved and un-changed in total quantity. Algebraic terms accounting for this conservation are substituted into the model for simplification and to reduce the order of the model.
Hoechst/DNA-Binding Sites: Cytosolic Transporter: Nuclear Transporter: Kinetic Reaction Equations Plasma Membrane Diffusion: Hoechst-DNA Association: Hoechst-DNA Dissociation: Cytosolic Hoechst-Transporter Association: Cytosolic Hoechst-Transporter Dissociation & Efflux: Nuclear Membrane Diffusion: Nuclear Hoechst-Transporter Association: Nuclear Hoechst-Transporter Dissociation & Efflux: Differential Equations: Cytosolic Hoechst: DNA-Bound Hoechst: Cytosolic Transporter-Bound Hoechst: Nuclear Hoechst: Nuclear Transporter-Bound Hoechst: In silico Flow Cytometry Simulation. Conversion of simulated Hoechst staining into Hoechst Red and Blue signals was mediated by a linear transformation of DNA-bound Hoechst and non-DNA-bound (free) Hoechst species within each cell in a process we refer to as in silico flow cytometry. Following the kinetic simulation of Hoechst staining, the molar quantity of total free Hoechst and DNA-bound Hoechst are determined for each cell. The quantities of these dyes were used to calculate a corresponding Hoechst Red and Hoechst Blue signal. Hoechst Red and Blue signals result from the combination of Hoechst Red and Blue emission from both DNA-bound and free Hoechst species (S11B Fig). DNA-bound and free Hoechst dyes possess different spectral properties, including quantum yield, excitation maxima, and emission maxima. [84] These differences manifest as differences in relative excitation efficiency and emission strength in the Hoechst Red and Blue emission channels (S11C Fig). Accounting for these factors, we are able to formulate a signal transformation matrix with which we can transform quantities of DNA-bound and free Hoechst into relative Hoechst Red and Blue signals.
Most flow cytometry techniques aim to isolate the signal from an individual fluorophore to a single detection channel. The SP assay, however, relies on spectral spillover of the Hoechst emission into two detectors, Hoechst Red and Hoechst Blue. Inherent spectral differences in DNAbound (H b ) and non-DNA-bound/free (H f ) Hoechst dyes would indicate that the two forms of Hoechst can influence the detected signal in each of these channels (S11C Fig). For example, the quantum yield of DNA-bound Hoechst (0.38) is roughly 10-fold larger than that of free Hoechst (0.034) [84], and more readily induced to emit fluorescent light upon excitation. Further, the excitation/emission maxima of Hoechst in the DNA-bound form differs from the free form. [84] The differences in emission spectra result in differential emission contributions to each of the detection channels (S11D Fig). In this schema, a number of factors influence the magnitude of the Hoechst Red and Blue emission signals. Nonetheless, it is a somewhat constrained system in that the Hoechst Red signal is composed of emission from both DNA-bound and free Hoechst and the Hoechst Blue signal is composed of emission from both.
Hoechst signals are calculated based upon the quantities of DNA-bound and free Hoechst species within the cell, the spectral properties of the Hoechst species, and the spectral properties of the simulated flow cytometer used to excite and measure Hoechst fluorescence. Hoechst Red signal (HR sig ) is the sum of the emission from DNA-bound (H b ) Hoechst in the Hoechst Red channel and from free Hoechst (H f ) in the Red Channel. The emission from H b in the Hoechst Red channel is proportional to its excitability (quantum yield, Q b ), relative excitation efficiency (E b ), the area of spectral emission overlap with the Hoechst Red channel (R b ), and the amount of H b . Likewise, the emission from H f in the Hoechst Red channel is proportional to its corresponding Q f , E f , R f , and H f ( Table 4). Signal for the Hoechst Blue channel can be similarly constructed.
We are then able to modify the representation of to obtain the linear transformation matrix S.
The Hoechst Red and Blue signals resulting from linear transformation with the signal matrix are arbitrary in that the units do not have a specific meaning. Nonetheless, within a range of Hoechst signals produced under the same circumstances, differences in magnitude reflect differences in the quantities of Hoechst dyes used to generate them. Therefore, the Hoechst signals can be compared for staining within populations. Similarly, with identical conditions, comparisons can be made across populations. Conversion of Hoechst signals into Scores is an approach to make measured changes in Hoechst staining more applicable in a broader, less experimentally specific sense.
Hoechst Score & PDF Conversion of Flow Cytometry Signals. Following the calculation of in silico flow cytometry values of Hoechst Red and Blue signals for each cell in a population, Hoechst Scores data was derived from Hoechst signal data. Processing of data from the in silico flow data was identical to that of data processed in real flow cytometry data. For a given in silico sample (inhibition and no inhibition pairing), PDF +FTC , PDF -FTC , and ΔFTC distributions were calculated. The SP size was measured using the projection gating approach. Finally, across the four conditions, the ΔFTC distributions were compared to the ΔFTC distribution from the highest transporter sample, corresponding to untreated control, to calculate ΔSP distributions.
Data-Driven Qualitative Selection of SP Responses. Simulation of each of the ensembles produced the following data: differences in Hoechst Scores staining metrics, Hoechst Score PDF+FTC, Hoechst Score PDF-FTC, ΔFTC, and ΔSP relative to the control, and %SP for each in silico sample. In the analogous experimental conditions, Day 4 of the SP time course with 4 differing TGFβ sample conditions, we possess equivalent experimental data (Fig 2D, S3 Fig and  S8C Fig). We used the experimental data to guide selection of models in terms of SP response. A series of selection check points were setup, which each ensemble was required to meet all selection criterion in order to be accepted as exhibiting a SP response.
First, the two highest transporter conditions were required to have negative ΔHRS mean and ΔHBS mean values, reflecting an overall decrease in Hoechst staining. Next, the PDF +FTC and PDF -FTC were not allowed to share any less than 25% overlap, indicating that the entire range of population did not shift without inhibition. The responses across all of the transporter conditions were required to reflect that of the experimental data. Ensembles were required to demonstrate a positive correlation for both the ΔHRS mean with experimentally observed ΔHRS mean values and ΔHBS mean with experimentally observed ΔHBS mean values. Experimental PDF +FTC , PDF -FTC , ΔFTC, and ΔSP distributions for all conditions were exported from the SP time course study. The normalized, aligned 2D cross correlation was calculated for each simulation/ experimental pairing. Within each of the categories, PDF +FTC , PDF -FTC , ΔFTC, and ΔSP, the average cross correlation of all of the conditions was required to be positive. Finally, ensembles were required possess a %SP of at least 5% for the control sample and a differential %SP of 2.5% between the control and lowest transporter expressing sample. Ensembles meeting the selection criteria were then scored according to the normalized root mean-square error (RMSE) of the %SP and differences in %SP with the experimental data using the MATLAB function goodnessOfFit with a normalized root mean square error (NRMSE) cost function. Fit values using this approach range from 1, with a perfect fit of simulated data to experimental data, to negative infinity with increasingly poor fit to experimental SP data.

Analysis of Single-Cell Side Population Response Distributions.
Within an ensemble, each of the cell populations is identical to one another, except for the relative amount of transporter activity. Unlike experimental assays, cells within the in silico assay are indexed and differences between samples perfectly controlled for. Therefore, we can examine how the exact same cell will stain differently under very tightly controlled alternate scenarios. Because of this feature, we can tabulate the difference in Hoechst Score projection in the inhibited and uninhibited Hoechst staining simulations. Thus, for each cell in a sample, we determine the difference between these two conditions as -ΔH proj , in which a larger value corresponds to a larger single-cell SP response. For each ensemble passing the qualitative selection process, the distributions of -ΔH proj values for each of the samples was further analyzed to interrogate the "shape" of the distribution. For each sample distribution, the 3 rd standardized moment (skewness) and 4 th standardized moment (kurtosis) was calculated. The skewness and kurtosis were used to calculate the bimodality coefficient: Model Implementation. At the start of each simulation, an in silico population of N cells was generated and M kinetic parameter sets constructed using LHCS. In parallel, kinetic parameter sets were submitted with cell populations to conduct the simulations within an ensemble. Many parameter sets were stiff to numerical solving. To prevent stalled simulation of the overall model, populations or single-cells that failed to solve within and allotted time window were aborted. Upon completion each ensemble was checked for a SP response. Upon completion of all of the models, the single-cell SP response distributions were analyzed for each of the passing ensembles in a model and aggregated for comparison.

Software
Flow cytometry and Flow Sight imaging cytometry data were processed and analyzed using FlowJo for Mac OS X version 10.0.7, Tree Star, Inc. Statistical analyses of experimental data were performed within Graphpad Prism for Mac OS X version 6.0e. Imaging Cytometry images were segmented using ImageJ and the SCIPY platform in python. Cytometry distribution analyses were performed using MATLAB version 2014a (64-bit), MathWorks Inc, in 64-bit Windows 8.1. Side population simulations were implemented in MATLAB version 2014a for Linux and run in parallel on the PACE cluster at Georgia Tech, which consisted of 64 single core 3.8 GHZ AMD processors with over 240 GB of total RAM available (10 GB per node). The following MATLAB File Exchange entries (accessed on 11/5/14) were implemented in MATLAB to analyze or display cytometry or simulation data: smoothhist2d (13352) [85], tight_subplot (27991) [86], suplabel (7772) [87], redblue (25536) [88], progress monitor (32101) [89], and distributionPlot (23661) [90]. In silico flow cytometry SP plots from an ensemble from Model 2, in which cells within the population express variable transporter concentrations. A) Schematic representation of the distribution of single-cell SP responses in a Full Type response. B) Hoechst Scores PDF +FTC and PDF -FTC plots for different transporter samples whose means were drawn from 0, 1, 10, and 100 pM TGFβ experimental conditions. Projection gating (gray line) was used to measure the %SP. C) %SP from the ensemble is compared to the means of experimental conditions. D) Differences in PDF -FTC and PDF +FTC distributions from A are shown as ΔFTC for each transporter sample. ΔSP distributions are differences in ΔFTC distributions compared to the 0 pM sample. (EPS)

S15 Fig. Simulation Outcomes by Transporter Condition.
For each of the transporter conditions (A-transporter number distribution (mode i); B-transporter concentration distribution (mode ii); C-experimental transporter distribution (mode iii)), a total of M = 10,000 kinetic parameter sets were generated by LHCS and cell staining simulated for the corresponding ensembles. For each of the 10,000 parameter sets, started for each condition, >95% of the sets successfully simulated the corresponding ensemble whereas <5% of the sets had to be aborted while simulating the corresponding ensembles due to failure to meet simulation time-out restrictions (left). Of the successfully simulated parameter sets, >95% failed to produce a SP response as determined by passing the qualitative selection process (right).  (Fig 2). For each of the parameter sets having generated ensembles with SP response, the %SP fit RMSE was included in the histogram. RMSE was calculated as a normalized root mean square error using the MATLAB function goodnessOfFit, which range from 1, perfect fit, to negative infinity with progressively worse fit to experimental data. Distributions of SP response (-ΔH proj ) in individual cells within the populations of ensembles meeting SP selection criteria were characterized by standardized skewness and bimodality. Simply put, skewness measures the degree of asymmetry of a distribution around the mean with a value of 0 corresponding to symmetry and positive values corresponding to distributions with larger ranges in the distribution above the mean. Likewise, negative skew values correspond to distributions with a larger range in the distribution below the mean than above it. The bimodality coefficient is calculated from the standardized skewness and standardized kurtosis. It has a range from 0 to 1, in which a value of 0 reflects a distribution with a single value while 1 corresponds to a distribution with exactly two values. Distributions with two fairly distinct modes score closer to 1 while distributions with a singular mode with a higher frequency score closer to 0. The mappings of a wide variety of example distributions are depicted along with representations of the range single-cell SP responses in an example cell population. Lower case letters correspond to positioning on the Response Distribution Map.