Novel Hybrid Phenotype Revealed in Small Cell Lung Cancer by a Transcription Factor Network Model That Can Explain Tumor Heterogeneity Free

Small cell lung cancer (SCLC), accounting for ∼13% of lung cancers (1), is exceptionally aggressive. Patients with extensive disease die ∼1 year from diagnosis, and patients with limited disease experience a dismal 20% cure rate (2–4). Standard of care (2), confined to chemotherapy and radiotherapy for half a century, is largely ineffective as SCLC patients exhibit high initial response rates rapidly followed by treatment-refractory relapse.

Expression-based subtyping, impactful in other cancers (5), may be effective in SCLC because of phenotypic variability (2, 6) with respect to neuroendocrine features of its cell of origin (7, 8). A recent study (9) identified two transcriptional SCLC subtypes distinguishable by Notch pathway activity and aggressiveness, but without mutational differences. In genetic mouse models of SCLC, Calbo and colleagues showed that spontaneously occurring neuroendocrine and non-neuroendocrine cell phenotypes coexist and cooperate to promote metastasis (8).

These reports indicate that a deeper understanding of cellular phenotypes could produce insights into biology and evolution of SCLC. A limitation of previous studies (8, 9) is that analyses were based on population averages, whereas variability in tumors should be considered at the single-cell level (5, 10, 11). It also remains unclear why this heterogeneity emerges.

To fill these knowledge gaps, we investigate SCLC phenotypic heterogeneity at the single-cell level using an integrative computational and experimental approach. Consistent with previous reports, we found two transcriptional subtypes at the population level in cell lines and patient specimens, characterized by gene coexpression modules enriched in neuroendocrine/epithelial (NE) and mesenchymal-like (ML) features. To understand how these phenotypes may arise in the absence of driving mutations (9), we hypothesized that they are attractors of a regulatory transcription factor (TF) network. This approach is grounded in the mathematical interpretation of Waddington's epigenetic landscape (12, 13), whereby attractors correspond to biological differentiation states or stable phenotypes. Based on this view, it has previously been proposed that malignant phenotypes in cancer correspond to attractors (14, 15), and some have suggested “differentiation therapy” from malignant to benign attractors as a possible treatment strategy (14–17).

To this end, we construct an SCLC master regulatory network of TFs from NE and ML gene-expression signatures. We then adopt a discrete Boolean modeling approach to simulate the behavior of this TF network and evaluate its ability to dynamically control NE and ML phenotypes. Discrete models are well suited to provide insight into complex TF networks by identifying steady state TF patterns of expression, termed attractors. Mathematically, attractors represent the stable configurations of the dynamic TF network. Biologically, attractors correspond to transcriptional steady states of an epigenetic landscape formed by active (ON) and silent (OFF) TFs regulating each other. While discrete models are coarse approximations, they are nevertheless informative and circumvent the obstacle of unfeasible parameter acquisition (18, 19).

Simulations of our SCLC TF network predict attractors corresponding to the NE and ML SCLC subtypes. Furthermore, by distilling the NE and ML states to their core driving TFs, the model highlighted a shortcoming of the two-subtype classification, as several cell lines and patient samples did not match any attractors. Western blots revealed that these samples expressed similar levels of both NE and ML markers. Flow cytometry revealed that this “double-positive” phenotype reflected the character of individual single cells, confirming the existence of a previously unreported “hybrid” single-cell phenotype in SCLC. Exposure to cytotoxic and epigenetic drugs caused NE and ML cells to transition toward the hybrid state, implicating it as a refuge for survival of treated SCLC tumors.

The Cancer Cell Line Encyclopedia (CCLE) dataset (20) was downloaded from the Broad Institute as CEL files. Data were normalized and median centered using quantile RMA normalization using Affy Bioconductor package (21) in R. Probe-level data for all the datasets were converted to gene-level data by probe merging using collapseRows (22). Probes with no known gene symbols were removed.

Consensus clustering was performed using ConsensusClusterPlus v1.24.0 package in R v3.2.3 (23) on both the 53 SCLC cell line dataset from CCLE (20) and 28 SCLC patients from the Clinical Lung Cancer Genome Project (CLCGP; ref. 24), with 80% subsampling of both genes and samples, 1,000 repetitions, 1 - Pearson correlation, and k-means. Both the CCLE and CLCGP datasets were subsetted to only include genes measured in the GSE6044 (25) dataset, in order to maximize overlap with our previous work (26).

Coexpression network analysis was performed in R v3.2.3 using the weighted gene coexpression network analysis (WGCNA) package v1.49 (27). As with the consensus clustering analysis, the CCLE SCLC dataset was subsetted to only include genes measured in GSE6044 (25). We used 1 - Pearson correlation to build a coexpression based dissimilarity matrix. Modules were generated using unsupervised average-linked hierarchical clustering with a static height of 0.95. We required that each module contain at least 100 genes.

Comparative Pathway enrichment analysis of the Blue and Turquoise modules was performed using BINGO and Enrichment map (28) and visualized in Cytoscape. Pathway enrichment for a cell line was computed by transforming expression values to Z-scores and computing the average expression of all genes within a given pathway.

To generate an SCLC-specific transcriptional network, we applied the bootstrap version of ARACNE (29) on gene expression profiles from the 53 CCLE SCLC cell lines using the following parameters: P = 10⁻⁷, dpi = 0 and 100 bootstraps, resulting in 27,224 interactions among 8,706 nodes. To evaluate if genes in Blue and Turquoise module are enriched for targets of any specific TF, we used Fisher exact test (30) to compute TF enrichment with the Blue and Turquoise module genes. We selected all TFs to be candidate master regulators if the Fisher exact test P value was ≤0.05, leaving 96 and 207 TFs for Blue and Turquoise modules, respectively. Of these, 23 TFs were common to both modules (Supplementary Fig. S7D).

Modeling gene regulatory networks that were inferred entirely from the available data can suffer from circular reasoning: a dataset generates a network, which then predicts the dataset. To avoid this fallacy, we built the network topology strictly using sources external to the datasets of interest. Thus, we first filtered ARACNE TF predictions based on gold-standard TF-target binding site databases CHEA, ENCODE, TRANSFAC, JASPAR using EnrichR (31–35), and literature databases such as Pubmed and Glad4U (36). These filtration steps produced a list of 76 likely TF regulators of NE and/or ML differentiation. We took only heterogeneously expressed TFs (median absolute deviation above the 50th percentile) yielding a list of 38 TFs that we used to build a Boolean network for SCLC. Next, we extracted directed interactions between these TFs using only information from the above gold-standard references and literature. Where possible, interactions were classified as activating or inhibiting by manually searching the literature. Interactions that we could not find in the literature were classified as activating if the TFs are positively correlated across the CCLE dataset, or inhibitory given negative correlation (Supplementary Fig. S7B). Five TFs did not have outgoing edges and were eliminated, leaving 33 TFs.

The TF network was simulated as a Boolean network where each node was either ON (active) or OFF (silent). Nodes were updated using the random order asynchronous method (18). We considered two distinct approaches updating TFs: (i) threshold updating, in which the total number of ON inputs are compared such that the target node is switched ON when there are more activators, and OFF when there are more inhibitors, and (ii) inhibitory dominant, in which having any inhibitor ON is sufficient to turn the target node OFF. The threshold update rule is referred to as a Hopfield neural network in some literature. Because the network's state space is so large, we only simulated a random subsample of the states. For both update rules, a state transition network was seeded with 8,000 randomly generated initial states, in which each TF had a 50% chance of being ON or OFF. When TFs are updated in asynchronous random order, it is possible that the state at time t may have several possible outgoing trajectories to time t+1. To account for this, each state that we observed was initialized and updated 30 times to sample distinct trajectories to the next state that may be influenced by the update order. Each newly observed state was queued to be updated in this fashion, until there were no new states identified. Attractors were identified by applying the attracting components algorithm from NetworkX (https://networkx.github.io/) to the state transition graph for each update rule. Using threshold updates, we found 57 fixed-point attractors, and no oscillating ones. We observed the same set of attractors using synchronous updates with 2¹⁷ and 2¹⁸ initial states. Using the inhibitory dominant update rule, we found 6 fixed-point attractors, and 5 two-state oscillating attractors, finding the same set of attractors using synchronous updates from 2¹² to 2¹⁸ initial states.

To score the correlation between samples and attractors, CCLE and CLCGP expression data were independently scaled between 0 and 1, and Pearson r was calculated pairwise between attractors and cell lines/patients. Statistical significance was determined by considering the highest correlation between each attractor and cell line, and comparing against a null distribution, obtained by generating 10,000 random TF vectors. The Mann–Whitney U test was used to compare the score distributions. The model attractors had higher correlation with cell lines than the random attractors (threshold: P = 9.5e−34, inhibitory-dominant: P = 7.6e−9), and higher correlation with patients than with the random attractors (threshold: P = 1.7e−23, inhibitory-dominant: P = 1.7e−6). Individual attractors were assigned a P value for their highest correlation with each sample by direct comparison with the null distribution (Supplementary Table S5).

Cells need to be able to robustly guide their differentiation choices depending on driving signals, and therefore we would expect trajectories toward cell attractors to be robust. Structural coherence (37) is a topological metric that measures how reliably an initial condition evolves toward its appropriate attractor given a perturbation. This metric requires an estimate of the total size of the basin of attraction, so we were only able to apply it to attractors for which we could reliably estimate this size (a few basins were too small to reliably estimate their size). Values of random, maximum, observed, and structural coherence are reported in Supplementary Table S4. We also calculated Derrida curves (Supplementary Fig. S9; ref. 38), as the average growth or decay of an initial perturbation after one step, for perturbation sizes ranging from 1 to 32 TFs.

Antibodies used for Western blotting include EphA2, LEF1, c-MYC, NFkB1, ZEB1, MITF, CD56 and CD44 (Cellsignal), OVOL2, SOX2, ASCL1, GAPDH, and β-actin (Sigma), FOXA2 and SMAD3 (BD Biosciences). Fluorophore-conjugated primary antibodies were used for flow cytometry—CD56 BV605, CD151 PE, CD24 BUV395 (BD Biosciences), CD44 Pacific blue, CADM1 A647 (MBL), EPHA2 A488 (R&D Systems), CD133 PE-Cy7 (Biolegend).

Authenticated cell lines were obtained from ATCC from 2012 to 2015, authenticated by DNA STR profiling, morphology, and mycoplasma detection. Cell lines were grown in company recommended media. New cell lines were obtained as needed every 1 to 2 years. Cell lines were expanded in culture for less than 2 months, then frozen in aliquots for subsequent use. Cell lines are passaged no more than 30 times before being discarded (approximately 4 months). Any contaminated cell lines were discarded, and new aliquots were thawed and cultured. Mycoplasma test was performed on all cell lines in culture every 2 weeks.

One to two million cells were plated in T75 or T150 flasks the previous day and collected for flow experiment the next day as described below. For drug treatment experiments, cells were plated same as above the previous day, followed by drug addition the next day. Cells were incubated with drugs for 3 days at 37 degrees then collected for flow experiments.

Cells were dissociated using TryplE (GIBCO) for 10 to 15 minutes followed by staining with Alexa 700 dye (Molecular Probes) for 5 minutes at 37 degrees. Cells were then washed and fixed with 2% paraformaldehyde (10 minutes at room temperature), followed by surface marker staining or permeabilization with ice-cold 100% methanol at −20 for 30 minutes. Cells were then stained with fluorescent conjugated antibodies for 30 minutes in dark at room temperature. Samples were washed with PBS and run on BD 5-laser instrument at the Vanderbilt Flow cytometry core. Fluorescent channels were compensated using anti-mouse IgK beads (BD Biosciences) that were tagged with fluorescent antibodies. First intact cells were gated on the Forward Scatter Area (FSC-A) and Side-Scatter Area (SSC-A) plot, where debris has a low FSC-SSC ratio. This was followed by gating for Alexa700-negative viable cells as described previously (Box 1 in ref. 39), where A700-positive populations are dead/dying cells. Finally, A700-negative viable cells were further gated to include only singlet populations. At least 20,000 to 30,000 viable singlet cells were collected per sample. All subsequent gates on fluorescent markers are made on viable singlet cells. The raw cytometer intensity readouts for fluorescent channels were first converted to log scale by using the asinh() function with a cofactor of 150. Gating was conducted in Cytobank (40).

Cell lines were plated for 2 days in complete medium to equilibrate. Lysates were prepared by spinning cells at 4°C, aspirating the media, and adding M-PER lysis buffer (Pierce) containing 1X phosphatase inhibitors 2 and 3 and protease inhibitor (Sigma-Aldrich). Lysates were incubated for five minutes at room temperature, vortexed for 30 seconds and centrifuged at 15,000 rpm for 15 minutes (at 4°C). The protein concentration was quantified using BCA assay (Pierce). Lysates were boiled for 10 minutes at 100 degrees with 1X NuPage sample buffer (Molecular Probes) and run on 8% or 4% to 12% Tris-glycine gels (Molecular Probes). Semi-dry transfer was followed by blocking with 1X Casein-TBS. Blots were imaged by chemiluminescence.

We applied consensus clustering (23) to the 53 SCLC cell lines in the CCLE database (20), and to 28 tumor samples from the CLCGP (24), with the expectation to detect previously observed neuroendocrine and nonneuroendocrine phenotypes (8, 9). Consensus clustering analysis, which provides a rationale for determining the number of robustly separated subtypes, indicated that the CCLE cell lines and the CLCGP patient specimens are most consistently separated into two distinct clusters (termed clusters 1 and 2; Fig. 1 and Supplementary Fig. S1, Supplementary Table S1).

To characterize these clusters, we applied WGCNA (27) to the 53 CCLE cell lines. WGCNA identified 13 gene coexpression modules (Supplementary Fig. S2A and Supplementary Table S2). To summarize the overall expression of each module, eigengenes (first principal component of all genes within a module) were computed (Fig. 2). The Blue module eigengene was upregulated (Bonferroni adjusted P < 0.001) within cluster 2, while the Turquoise module eigengene was upregulated (Bonferroni adjusted P < 0.001) within cluster 1 (Fig. 2 and Supplementary Fig. S2B). No other module eigengenes showed statistically significant differences between the consensus clusters (Supplementary Fig. S2B).

WGCNA modules are derived agnostically and reveal coexpressed genes participating in similar biology. Blue and Turquoise module genes therefore provide information about the overall character of the consensus clusters. Visual inspection reveals that hubs (highly connected genes) in the Blue and Turquoise modules are well-known NE or mesenchymal biomarkers, respectively (Supplementary Figs. S3 and S4). Statistical analyses by Gene Ontology and Enrichment Map (28) confirm enrichment of the Blue module for neuroendocrine and epithelial differentiation processes, while the Turquoise module is enriched for pathways involved in mesenchymal phenotype and epithelial-to-mesenchymal transition (EMT; Supplementary Figs. S5 and S6A). Additionally, we performed gene-set enrichment analysis (GSEA; ref. 41) with published signatures of proneural, mesenchymal, and proliferative glioblastoma subtypes (30) on the 12 cell lines with highest Blue module eigengene expression, and 9 cell lines with highest Turquoise module eigengene expression. The proneural (Supplementary Fig. S6B, bottom) and mesenchymal (Supplementary Fig. S6B, top) signatures were enriched in the cell lines with high Blue module and Turquoise module expression, respectively. As a control, the proliferative signature was significantly enriched in neither (P > 0.1, data not shown).

Thus, the Blue module genes are a signature of the canonical neuroendocrine state of SCLC, while the Turquoise module genes are a signature of nonneuroendocrine/mesenchymal cells. We therefore refer to cluster 2 (in which Blue module genes are high) as the SCLC NE subtype, and to the cluster 1 (in which Turquoise module genes are high) as the SCLC ML subtype (Fig. 1B and C).

Because NE and ML subtypes are not associated with driver mutations (9), we hypothesized they may be regulated epigenetically, as during normal cell differentiation. Cell identity in differentiation is largely controlled by regulatory networks of TFs that coordinate expression of each other and of target genes (42, 43). To understand specification of NE or ML cell identity, we derived a network of TFs that regulate expression of genes within the Blue and Turquoise modules (see Materials and Methods; Supplementary Fig. S7). By pruning nodes that had no outgoing edges, we reduced this network to a core set of 33 TFs and 361 interactions (Fig. 3A and Supplementary Table S3). We used random asynchronous order Boolean simulation with both threshold and inhibitory dominant update rules to identify attractors. Using the threshold update rule, we found 57 stable fixed-point attractors (Fig. 3B) and no oscillating attractors. With the inhibitory dominant rule, we found only 6 fixed-point attractors, but 5 two-state oscillating attractors as well. Both the threshold and inhibitory dominant attractors revealed similar features in all subsequent analyses, suggesting that these results are robust to the precise nature of the regulatory interactions (Supplementary Fig. S8). We also inferred network robustness to perturbations using structural coherence (Supplementary Table S4; ref. 37) and Derrida curves (Supplementary Fig. S9; ref. 38). These results indicate that the basins of attraction are more robust than those of a random network, and that the network dynamics are ordered, rather than chaotic. Together, these observations suggest that the internal structure of interactions imposed to the TF network (Fig. 3A) is non-random and naturally leads to well-regulated stationary states, and that our coarse modeling approach is acceptable for this system.

Each attractor is a 33-dimensional vector of TFs, differing by the overall TF ON–OFF expression pattern (columns in Fig. 3B). Many of the attractors differ only by one or a few TFs. Hierarchical clustering segregated the threshold attractors into four distinct clusters (Fig. 3B). Results from the inhibitory dominant network are qualitatively similar and are shown in Supplementary Fig. S8. Based on the TF ON–OFF patterning in each attractor within the clusters, attractors 49 to 25 (along the dendrogram, Fig. 3B) have active TFs known to be either neuroendocrine (INSM1, POU3F2, SOX2, SOX11) or epithelial (FOXA2, OVOL2). In contrast, attractors 35 to 34 (along the dendrogram, Fig. 3B) contain active TFs (MYC, NFKB1, SMAD3) involved in mesenchymal differentiation or EMT.

We computed a correlation score between attractors and both cell lines and patients by scaling the CCLE (Fig. 3C) and CLCGP (not shown) expression data from 0 to 1 and calculating Pearson correlation coefficient pairwise between each attractor and sample. Several attractors exhibited high correlation with cell lines within the NE consensus cluster, while others exhibited high correlation with cell lines within the ML consensus cluster (Fig. 3D). Similar results were observed with patients (Fig. 3E), and in both cases the model's attractors were found to be significantly more correlated to the samples than random (Supplementary Fig. S10 and Supplementary Table S5).

These results confirm that the simulated dynamics of the 33 TF network agree well with the NE and ML nature of the cell lines, suggesting that the gene expression signatures of these phenotypes is driven by these underlying NE and ML TF attractors. Nevertheless, a few attractors showed no significant correlation with any cell line or tumor. These attractors were not further pursued, as we do not observe them biologically. More significantly, however, several cell lines and patient samples did not have significant correlation with any attractor. These samples may represent mixed populations of NE and ML cells, leading to poor correlation with either subtype at the population level, or alternatively may represent populations of “hybrid” cell types.

We probed the expression of cell lines for 10 of the 33 TFs. In general, NE TFs were expressed at a higher level in NE than ML cell lines, and vice versa (Fig. 4A). Consistent with the finding that some cell lines did not correlate with either the NE or ML attractors (Fig. 3D and Supplementary Fig. S8B), we observed cell lines that simultaneously expressed both NE and ML TFs. Similarly, CD56 (NE marker) and CD44 (ML marker) were found coexpressed in 3 out of 10 patient samples (pt1112-1, pt216-1, and pt460-1, Fig. 4B). The other tumor samples had mutually exclusive expression of either CD56 or CD44.

These double-positive cell lines may be either composed of mixed populations of NE and ML cells, or hybrid single cells simultaneously coexpressing features of both phenotypes. We investigated this using single-cell flow cytometry with well-established NE and mesenchymal differentiation biomarkers (Fig. 5 legend). To aid visualization, we defined overall NE and ML scores as the unweighted sum of NE and ML biomarkers, respectively.

In biaxial plots of these scores, NE cell lines consist primarily of NE⁺ML⁻ single cells, while ML cell line single cells were NE⁻ML⁺ (Fig. 5A and C). However, in several cell lines, single cells simultaneously expressed similar levels of both NE and ML biomarkers (Fig. 5B). These results confirm the existence of a novel SCLC NE⁺ML⁺ hybrid phenotype comprised of both NE and ML characteristics.

Other reports of hybrid phenotypes in cancer have associated the hybrid cells with a more stem-like phenotype (44, 45). We measured single-cell expression of CD133, a cancer stem-cell marker, and found no significant difference between the NE, ML, and hybrid cell lines (Supplementary Fig. S11).

Next, SCLC cell lines were treated with etoposide and cisplatin, and with epigenetic modulators valproic acid (HDAC inhibitor) and 5-azacytidine (DNA methylation inhibitor). Flow cytometry measurement of NE and ML markers (see Fig. 5 legend) was used to characterize how treatment shifted the phenotypic identity of SCLC cells. Biaxial plots of NE vs. ML scores showed phenotypic shifts at the single-cell level for all perturbations (Fig. 6), converging toward hybrid populations. The hybrid phenotype was thus an end state for SCLC cells subjected to stress.

In this work, we identified signatures of phenotypic heterogeneity in SCLC, and a set of TFs regulating expression of these genes. Through discrete Boolean model simulations, we showed that a master SCLC TF network naturally settles into states that were identified as NE or ML. This suggests that the NE and ML cell identities can naturally emerge from regulatory dynamics of TF networks, rather than being driven by genomic mutations. Most significantly, by distilling the NE and ML phenotypes to their essential TF drivers, the model also exposed some cell lines that had been improperly classified as NE or ML through transcriptome-based clustering. In vitro validation confirmed existence of SCLC cell lines and patient samples that simultaneously express both NE and ML markers and TFs. We verified this hybrid phenotype at the single-cell level in cell lines and showed that cell lines transition toward this hybrid state for survival upon drug treatment.

As resolving cancer heterogeneity can have profound impact on patient care and outcomes (46), a deeper understanding of SCLC heterogeneity, i.e., NE, ML, hybrid and beyond, may translate to benefits in the clinic. Historically, attempts to separate SCLC into subtypes in the clinic were abandoned due to poor reproducibility among pathologists and unclear clinical relevance (47). However, our findings indicate that SCLC heterogeneity is dynamic, because a core TF network specifies both the NE and ML phenotypes, which may confound static associations. The hybrid phenotype may be an additional confounder, exacerbated by the tendency of both NE and ML cells to move toward the hybrid phenotype under treatment. We argue that SCLC heterogeneity is best interpreted as states within a phenotypic landscape, and understanding phenotypic mobility within this landscape could provide ties to clinical relevance.

While SCLC has a high mutation rate, no correlations between mutations and distinct subtypes have been observed (9). In this work, we show a non-mutational mechanism for regulation of distinct SCLC cell identities, as attractors of a TF regulatory network. This framework of equating cell types to attractors was first advanced by Waddington in his eponymous epigenetic landscape (12), and has because been expounded upon by many mathematical biologists (13). It has further been argued that malignant cancer states and cancer heterogeneity are best described by the concept of attractors (14–16). Significantly, this framework has been shown to provide possible therapeutic strategies by identifying targets in silico, which most significantly perturb the attractor states (17). We do not exclude the possibility that mutations can drive SCLC heterogeneous phenotypes, but our results indicate that epigenetic causes should also be considered.

The Boolean modeling approach used here is an established way of obtaining a coarse-grained picture of large network behavior, and is well suited to identifying the stable states of TF networks (18, 19). Our approach was novel in constructing the TF network in a blind fashion from SCLC gene expression datasets, TF databases and literature, and was validated by its ability to reproduce correlations with known phenotypes. Significantly, our model revealed the existence of a previously unrecognized hybrid phenotype, by showing that some cells do not correlate to either the NE or ML attractors. However, we did not identify attractors correlated with the hybrid phenotype. Hybrid EMT phenotypes have been previously reported in NSCLC (48) and lung adenocarcinoma (49), and other groups have recently reported computational modeling of hybrid EMT phenotypes by driving EMT networks with external stimuli (44, 50). Additionally, the Boolean modeling approach cannot capture intermediate levels of expression, and therefore attractors corresponding to the hybrid state may not be identifiable using this method. Continuum modeling approaches may be needed to better understand the hybrid state. Identification of gene coexpression modules enriched in the hybrid phenotype may also reveal additional relevant TFs. Considering all of these possibilities, current ongoing work in our laboratory is directed at identifying hybrid attractors.

Our study establishes the hybrid phenotype as a refuge for drug-treated SCLC cell lines. We anticipate that this phenotype may play a significant role in the evolution of SCLC tumors under treatment, and possibly in relapse. SCLC TF networks may serve as a guide for interventions aimed at preventing phenotypic transitions into resistant states. While we focused on SCLC here, such transcriptional regulation may play a similar role in maintaining non-mutational heterogeneity in other cancer types, and our approach should be generally useful in uncovering underlying mechanisms.

A. Califano is cofounder and chief scientific advisor at DarwinHealth, Inc. and has ownership interest (including patents) in the same. J.M. Irish is co-founder and board member of Cytobank, Inc. and received research support from Incyte Corp. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.R. Udyavar, D.J. Wooten, L. Estrada, P.P. Massion, V. Quaranta

Development of methodology: A.R. Udyavar, D.J. Wooten, P.P. Massion

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A.R. Udyavar, M. Hoeksema

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A.R. Udyavar, D.J. Wooten, M. Bansal, S. Schnell, J.M. Irish, P.P. Massion

Writing, review, and/or revision of the manuscript: A.R. Udyavar, D.J. Wooten, A. Califano, S. Schnell, J.M. Irish, P.P. Massion, V. Quaranta

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A.R. Udyavar, P.P. Massion

Study supervision: A.R. Udyavar, A. Califano, P.P. Massion, V. Quaranta

We acknowledge the Vanderbilt Flow Cytometry shared resource lab for the help and assistance in training for the use of flow cytometry instruments. We thank Dr. Charles Rudin (Memorial Sloan Kettering Center, New York, NY) and Dr. Julien Sage (Stanford University, Stanford, CA) for providing us with patient-derived xenograft tumors.

This work was supported by NIH/NCIU54 CA113007-09, Integrative Cancer Biology Program (ICBP, V. Quaranta), Grant 1I01CX000242 from the Department of Veterans Affairs, CA90949 from the NCI SPORE program (P.P. Massion), NIH/NCIR00 CA143231 (J.M. Irish), NIHT32 CA009592 (J. Chen), National Center for Advancing Translational Sciences UL1-TR000445-06 G. Bernard, and the Vanderbilt-Ingram Cancer Center (P30 CA68485 to J. Pietenpol).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.