Systems biology approaches to measure and model phenotypic heterogeneity in cancer

The recent wide-spread adoption of single cell profiling technologies has revealed that individual cancers are not homogenous collections of deregulated cells, but instead are comprised of multiple genetically and phenotypically distinct cell subpopulations that exhibit a wide range of responses to extracellular signals and therapeutic insult. Such observations point to the urgent need to understand cancer as a complex, adaptive system. Cancer systems biology studies seek to develop the experimental and theoretical methods required to understand how biological components work together to determine how cancer cells function. Ultimately, such approaches will lead to improvements in how cancer is managed and treated. In this review, we discuss recent advances in cancer systems biology approaches to quantify, model, and elucidate mechanisms of heterogeneity.


Introduction
Cancers are marked by substantial genetic, epigenetic, and phenotypic heterogeneity. With the advent of single cell technologies, our ability to quantify cell-to-cell variability is advancing rapidly; however, the biological interpretation of these complex data is only at a nascent stage. A systems biology approach, in which quantitative technologies, comprehensive experimental measurements, and computational analyses are carefully integrated, is required to understand intratumoral heterogeneity and extract biological meaning. Moreover, computational frameworks will aid in the elucidation of the underlying biological mechanisms of heterogeneity. Such integrated approaches will ultimately lead to a deeper understanding of cancer and how to manage it to improve patient outcomes. Here, we cover recent advancements in cancer systems biology that have improved our understanding of heterogeneity in genetic, phenotypic, signaling, and pharmacokinetic realms ( Figure 1).

Genetic heterogeneity
It is well-established that the initiation and progression of cancers depends on the acquisition of multiple driver mutations that activate oncogenic pathways and lead to the induction of cancer hallmarks 1,2 . Moreover, the growth of tumors likely reflects a complex process of clonal evolution that ultimately results in considerable genetic heterogeneity within tumors 3 . A recent pan-cancer analysis of primary tumors revealed that the majority of cancers harbor multiple genetically distinct clones, and tumors with more clones harbored more driver mutations and showed greater nuclear heterogeneity 4 . Moreover, increased genetic heterogeneity is associated with poor outcome, indicating that measures of intratumoral heterogeneity may serve as useful biomarkers 4 . With the widespread adoption of single-cell sequencing technologies, we are poised to develop deeper and more precise views of intratumoral heterogeneity 5 .
Recent years have seen a substantial advancement in our understanding of genetic heterogeneity, particularly as it relates to therapeutic response. Using cell line model systems, Ramirez and colleagues studied drug-tolerant persister cells and found that they may provide a latent reservoir for the emergence of drug-resistance. Specifically, long-term drug exposure leads to survival and expansion of cells through a drug-tolerant state that can mediate genetically-driven resistance mechanisms 6 . Consistent with this idea, single-cell DNA and RNA sequencing of triple negative breast cancer patient samples profiled before and after treatment demonstrated that resistant genotypes were pre-existing and adaptively selected after treatment, while new transcriptional profiles were acquired by reprogramming in response to therapy 7 . Taken together, these studies further motivate the importance of understanding and managing intratumoral genetic heterogeneity. Laajala, et al (this issue) reviews systems biology approaches to study cancer genetic heterogeneity in greater depth.

Phenotypic heterogeneity
Cancers display substantial phenotypic heterogeneity and rapid adaptation in response to stimuli that occur on timescales that cannot be explained by genetic evolution or clonal selection, indicating a substantial contribution of the epigenome to cell-to-cell variation. A landmark study by Sharma and colleagues revealed that small sub-populations of cancer cells could mediate therapeutic resistance on time-scales too rapid to be explained by genetic evolution and selection 8 . Moreover, the acquisition of this drug-tolerant phenotype was found to be transient and result from chromatin changes, indicating a dynamic regulation of phenotypic heterogeneity which stands in contrast to long-term "hard-wired" genetic changes. In support of this hypothesis, a recent study of therapeutic response revealed substantial cell-to-cell transcriptional variability that involves infrequent, semi-coordinated transcription of a number of resistance markers in a small percentage of "jackpot" cells 9 . Such findings indicate that heterogeneity is mediated at multiple molecular levels and therefore requires integrated experimental approaches that can link across modalities at the single cell level.
Several groups have developed computational frameworks to shed light on the molecular basis and dynamic nature of phenotypic heterogeneity. Gupta and colleagues found that subpopulations of sorted phenotypically-distinct triple-negative breast cancer cells could return to equilibrium proportions over time 10 . This observation could be explained by a Markov model in which cells transition stochastically between phenotypically distinct states, indicating phenotypic plasticity even in the absence of an external pressure such as therapy. Building on this idea, Risom and colleagues used an ordinary differential equation (ODE) model to show that drug tolerant states can arise through differentiation state transitions rather than Darwinian selection of preexisting subpopulations 11,12 . Such switching between differentiation states may enable cells to rapidly evade therapy. These studies highlight the importance of coupling experimental observation with computational frameworks to gain insight into the underlying mechanisms of cell-to-cell variation and plasticity.
The wide-spread adoption of single cell sequencing methods provides insights into the phenotypic state of hundreds to thousands of individual cells and has spurred novel analytics to aid in interpretation of these complex data [13][14][15] . In particular, some of the challenges of single cell sequencing include both technical artifacts related to differences in read-depth, dropouts (i.e., stochastic false-negative gene expression measurements), and variation in total transcript abundance 16,17 . However, many of these technical challenges are being overcome with novel analytical approaches. For example, Azizi and colleagues developed the BISCUIT algorithm, a Bayesian approach for simultaneous clustering and imputing single cell RNAseq data 13 . Seurat is a suite of tools for analysis and visualization of single cell RNAseq data sets based on common sources of variation, which enables identification of shared populations across data sets; its major strengths are ease-of-use and extensive data visualization tools 14 . As new technologies develop, so too does the interest in integrating data from across various modalities to develop a complete picture of the biological mechanisms associated with heterogeneity 18 .
Ultimately, effective cancer therapy regimens will need to consider heterogeneity, and several methods to design therapies combatting cell-to-cell variation have emerged in recent years [19][20][21][22] . For example, DRUG-NEM is an integrated experimental-computational approach to rationally identify drug combinations that account for intratumoral heterogeneity 19 . This approach applies cyTOF (Mass Cytometry Time-of-Flight), a variant of flow cytometry that replaces fluorescence detection with metal ion mass spectrometry for greatly improved multiplexing 23 . cyTOF enables the simultaneous measurement of many intracellular and surface markers for hundreds of thousands of single cells in a sample 23 . The authors assess single-cell proteomic responses to individual drugs and then use a nested effects model to prioritize drug combinations that produce the maximal desired intracellular effects. Such nested effects models can identify common and unique signaling changes in response to drug treatments and this information can be used to optimize the minimum number of drugs that induce desired signaling changes 24 . Functional approaches such as these are quite advanced for blood cancers, where primary cancer cells are relatively accessible 25 .

Heterogeneity in cell signaling and pathway activity
While statistical models like those described above have refined our view of how cell-to-cell variability drives resistance, mechanistic information about relevant resistance pathways enables in silico experimentation of how individual molecular events influence overall behavior. Cells respond dynamically to treatment, and so genetically-encoded single cell reporters and in vitro live-cell time-lapse imaging, paired with dynamical systems models, have become a mainstay of studying this variability [26][27][28][29][30] . The continued development and use of these experimental and theoretical tools together demonstrate the power of this approach. From foundational studies of variability in TRAIL-induced apoptosis, the key elements of studying cell heterogeneity in specific pathways have remained largely unchanged [31][32][33][34] . Common to all these studies has been that cell-to-cell variability is explained by variation in the abundance of a few key proteins, rather than variation across many pathway components or intrinsic noise. Critically, computational models identified key components that both predict variation in phenotypic response across cells and pinpoint how to control this variability.
Imaging paired with dynamical systems analysis will continue to be used in studies of pathway activity variability and drug resistance. However, new experimental and computational tools will certainly expand the impact this pairing provides. Studies of pathway-specific cell-to-cell heterogeneity will always depend upon detailed knowledge of pathways and their relevant regulatory processes. Therefore, general improvements in our biological knowledge will expand the range of pathways and processes that can accurately be modeled. While imaging provides uniquely rich, dynamic data on single cell responses, it remains limited in the number, types, and throughput of molecular measurements one can perform. Methods to interrogate molecular state, such as endogenous protein tagging, will widen the accessible range of pathways and resistance mechanisms 35 . The number of measurements that can be multiplexed in single living cells is unlikely to considerably improve, and cell-to-cell variability dictates non-destructive techniques for sampling multiple measurements in the same cell. However, dynamic measurements can be paired with destructive measurements as an endpoint, and data can then be integrated to develop an understanding of the biological mechanisms that may be driving differences in dynamic responses 36 . Modeling developments to parameterize dynamical models on mixtures of population and individual cell measurements will be critical to enhance the value of these mixed measurements 37 . Finally, individual pathways operate within a larger context of broader molecular changes. Methods to integrate "phenotypic" and mechanistic cell variability will ultimately be necessary to recognize the simultaneous contributions of pathway-specific and global changes. These burgeoning improvements in measurement and modeling show that dynamical models are poised for expanded use in studies of cell signaling.

Spatial heterogeneity
Properties of the in vivo tumor environment introduce sources of phenotypic and molecular heterogeneity beyond those outlined above. In contrast to most in vitro models of cancer, tissues present variation in the extracellular matrix, resident host cells, and barriers to drug access 38 . Each of these can promote resistance; indeed, pharmacokinetic/pharmacodynamics (PK/PD) limitations have long been known to drive resistance across disease indications, particularly for tissues such as the central nervous system 39 . Spatial variation in these properties may co-exist in the same tumor and hinder drug function even in normally accessible tissues 40 .
In vivo and fixed tissue imaging modalities are uniquely powerful for capturing spatial variation in tumor properties. Mass spectrometry coupled to localized tissue sampling methods has recently enabled drug distribution imaging, alongside traditional immunohistochemistry techniques 41 . These methods also extend whole organism imaging methods in resolution and sensitivity to identify drug biodistribution within individual tissues. While limited to specific tissue sites, fluorescence intravital imaging has also provided a window into the spatial variation of tumors 42 . For example, even within a single patient, individual tumor nodes can have striking variation in immune infiltrates 43 . A continuing challenge will be the many length scales of spatial variation, from poor drug distribution across tissues to drug sequestration by neighboring cells 44 .
Pharmacokinetic models are widespread and accepted for modeling the diffusion and transport of molecules within cells and tissues. PK/PD models provide precise predictions linking molecular mechanism to cell response but are only accurate if built with a detailed understanding of the relevant properties to model. Randall et al recently integrated in vivo imaging of drug distribution and phosphoproteomic evaluation of signaling effect to reveal that intracranial EGFR inhibitor availability cooperates with other mechanisms of resistance in glioblastoma cells 45 . This integration is especially powerful for distinguishing between resistance driven by drug availability and other mechanisms.
Despite spatial effects being long-recognized as a source of intratumoral heterogeneity, more robust experimental tools and models are needed to help address how this variation influences cancer biology and treatment. In particular, the range of perturbations possible and the scale on which they can be performed is quite limited in the in vivo environment. In vitro methods to recapitulate critical features of the in vivo environment will be critical to enable mechanistic studies. For example, Wu et al. have applied a microfluidic device to study the interaction between breast cancer cell drug resistance and migration 46 . Others have demonstrated that hydrogel biomaterials can mimic the tumor extracellular matrix and enable detailed study of individual components 47,48 . A key requirement of such mimetic systems is the inclusion of salient variables and features, which are not always well understood.

Future prospects for understanding and managing intratumoral heterogeneity
In reality, these distinct forms of heterogeneity likely operate to varying extents in parallel within a single tumor. Therefore, applying new findings from model systems will require methods to evaluate which forms of heterogeneity are most relevant to therapeutic response. Integrating methods to study these forms of heterogeneity will also reveal how they can impact each other. For example, Keren et al. recently showed that imaging cyTOF can provide both high dimensional and spatial information in breast tumors, and that the spatial organization of tumors can predict prognosis 49 . Forms of single-cell gene expression analysis 50 that preserve spatial information are now also enabling similar analysis, along with highly-multiplexed immunofluorescence 51,52 . These approaches afford assessment of multiple biomarkers in a single cell, all while maintaining spatial information about cell-tocell variability and organizational structure. Of course, the studies described above have universally relied on a theoretical model tailored to the experimental system to understand the heterogeneity involved. As new experimental systems are developed, new computational models will need to be developed to leverage these data. Finally, a deep understanding of heterogeneity will ultimately require the integration of data across multiple physical and time scales, which can be achieved with multiscale modeling 53 .

Summary and conclusions
Genetic heterogeneity perhaps offers a glimpse into the future of how we will understand and manage other forms of intratumoral heterogeneity. Observation and characterization have given way to broad behavioral principles, such as that of clonal selection and expansion. Ultimately, all these forms of heterogeneity dynamically respond to treatment and therefore clinical management will require strategies to monitor patient's tumors over time and in response to therapy 7,54 . An integrated cancer systems biology approach will provide the framework for understanding the molecular basis of heterogeneity and how to manage it clinically.

•
Cell-cell variation manifests at multiple levels and scales in cancer  Classes of therapeutically-relevant heterogeneity. Variation may arise from cell-to-cell differences in mutational spectra, phenotype, cell signaling pathway activity, or pharmacokinetics. Measurement technologies and modeling approaches to study these various forms of heterogeneity are indicated. Meyer Curr Opin Syst Biol. Author manuscript; available in PMC 2020 September 03.