In cancer, a gain of stemness can have profound implications on tumor aggressiveness, drug response and clinical outcome (Figure 1A). Here, we provide a systems biology overview of how the immune cell niche, cellular contexts, and molecular or genetic perturbations contribute to stem cell-like properties of malignant cells. We start by describing the sources of tumor heterogeneity and defining the cellular niche as a dynamic spatial domain harboring cancer stem cells (CSCs). We then elaborate on how this cellular niche is critical for communication between cell populations and move through the biological pathways endowing these stemness properties. We outline how systems biology approaches present an important strategy for identifying the crosstalk between immune cells, the bulk cancer, and CSCs. Finally, we discuss how these systems biology approaches open new avenues to dissect and manipulate cellular states and niches to impact cancer progression.

Tumor complexity not only arises from the intrinsic diversity of the tumor cells themselves, but also from their interactions with other non-tumor cell types. The intratumoral heterogeneity of the cancer cells arises through mechanisms such as genomic instability, clonal evolution, cellular differentiation, and reversible plasticity. These contributions have been well studied, but most analyses fail to systematically characterize the contributions of other cells to cancer cell heterogeneity. Immune cells, which are also part of the tumor niche, present over 350 CD (cluster of differentiation) antigens, secrete close to 100 cytokines and chemokines, and belong to subtypes that express thousands of unique gene signatures (Davis et al., 2017). These immune cells, as well as other non-malignant stromal cells, co-exist with cancer cells, thus creating complex biological entities. Moreover, these stromal cells interact with and shape the evolution of the cancer cell populations (Figure 1B).

From a systems biology perspective, the contribution of each source of heterogeneity can be understood by dissecting cell populations, cellular niches, niche dynamics and the crosstalk between them. This will ultimately help devise disease spatial resolution maps that shed light on mechanisms of cancer progression and drug resistance. Such complex systems exhibit emergent properties that are not just the aggregation of unconnected individual behaviors, but are instead unique phenomena that reflect synergy in interactions between cells. The scientific discipline of systems biology aims to provide a quantitative and dynamic framework that leverages mathematical and computational techniques to unravel the complexity of cellular processes. This approach is applied with the view that a holistic view of systems will not only inform how a gene connects to a protein and an activity, but also help elucidate how a cellular phenotype is contingent on the biological context (where and when) (Kirschner, 2005). As such, systems biology approaches have the power to dissect the interplay between tumor and immune cells to characterize higher order properties and systems interconnectivity. For example, interrelationships between ligands and receptors have been an active area of research where a variety of system methodologies have been applied to define ligand/receptor co-expression, formation of multimeric complexes and cellular activity (see section entitled Probing the cellular niche: Technological advances). The large breadth of ‘omics data now available is enabling systems approaches to give new insights. Here, we will highlight the potential of systems biology to dissect the determinants of cancer stemness present in the cellular niche.

Stemness, defined as the potential for self-renewal and differentiation from the cell of origin, was originally attributed to normal embryonic stem cells that give rise to all cells in adult organisms. However, it has long been known that cancer cells show markers and properties of embryonic stem cells (Friedmann-Morvinski and Verma, 2014). Cancer progression often involves the gradual loss of a differentiated phenotype and acquisition of progenitor or stem cell-like features. Undifferentiated primary tumors are more likely to result in metastasis to distant organs, causing disease progression and poor prognosis, as well as resistance to available therapies. Within a tumor, phenotypic diversity and spatial cellular variation may also impact expression of key embryonic stem cell regulators, resulting in distinct paths to cancer stemness (Figure 1C).

A challenge with CSCs is that only a small number of stem cell markers have been identified. This limited set of markers reflects the small proportion of CSCs generally present in tumors and the poor conservation of CSC surface markers across cancers. Therefore, from a systems perspective, an ab initio approach is best for the identification of CSCs.

A diverse array of high-throughput technologies has been used to zoom in on CSCs and their niches within cancer patients. The oldest and most common modality is RNAseq. A limitation of traditional bulk RNAseq, however, is that it averages expression across thousands of cells within a sample. Because of this limitation, understanding and constructing regulatory networks for distinct types such as CSCs has been difficult with bulk RNAseq alone.

There is, however, precedent for studying distinct tumor components using bulk RNAseq. CIBERSORT (Newman et al., 2015) has been used extensively to computationally dissect immune cell populations from these bulk samples; and several methods exist for constructing tissue level networks. CIBERSORT uses a knowledge-based signature matrix that encapsulates major 22 functionally defined human hematopoietic subsets to deconvolute and infer cell composition from gene expression profile data, utilizing support vector machine (SVM) regression methods. This approach has played a key role in decoding immune cell population, particularly in bulk solid tumors where tissue dissociation protocols and cellular enrichment techniques limit the efficacy of single-cell methods and provide only a partial view of the wider cell heterogeneity. Extending this approach in a new direction, a groundbreaking study by Thorsson et al. used TCGA RNAseq data across 33 cancer types and more then 10,000 patients to identify signatures for six TME subtypes and demonstrate that these TME subtypes can be associated with prognosis as well as genetic and immune modulatory alterations (Thorsson et al., 2018). Reflecting the importance of the TME to tumor behavior, the TME subtypes identified by Thorsson et al. were reportedly able to predict disease outcomes and help guide novel treatments. More recently, Bagaev et al. derived 29 expression signatures to establish four TME subtypes (immune enriched fibrotic, immune enriched non-fibrotic, fibrotic and depleted) that correlate with immunotherapy efficacy in melanoma (Bagaev et al., 2021). As demonstrated by these studies, it is possible to use RNAseq to study a subpopulation of cells (e.g., immune cells or stromal cells) if the transcriptional profile is detectable and distinctive. However, given the rarity of CSCs, it has been difficult to study CSCs themselves using bulk RNAseq. Moreover, while publicly available multi-omics signatures from non-neoplastic stem cells were previously used to derive a stemness score in tumors (Malta et al., 2018) that can then be applied using deconvolutional methodologies to identify CSCs in bulk RNAseq, a key caveat is that CSCs do not necessarily share the same transcriptional programs as normal stem cells counterparts. Taking a different approach, Aran et al. probed the normal tissue adjacent to the tumor (NAT), which is morphologically similar but phenotypically different from the tumor and focused on stromal pathways to query the interaction between adjacent tissues (Aran et al., 2017). These authors uncovered NAT-specific characteristics, namely activation of pro-inflammatory immediate-early response genes concordant with endothelial cell stimulation. Furthermore, previous studies on breast NAT suggested that the microenvironment surrounding the tumor, not the epithelial cells, is essential for understanding disease recurrence and developing surgical strategies (Graham et al., 2011). These studies have shown that by zooming into the tumor niche and its adjacent boundaries we can acquire deeper understanding of cell-cell interplay and dissect the phenotypic states that sustain stemness.

With the advent of single-cell RNAseq (scRNA-seq), new avenues to investigate and quantify molecular features at single-cell resolution have emerged. Single-cell technologies have the ability to directly evaluate cell states, heterogeneity, and lineages (Saelens et al., 2019). However, more effort is required to provide understanding of cell-cell interactions as well as their co-evolution and adaptation to perturbations in the cellular milieu, which is needed for the construction of a comprehensive cellular interaction map (Elmentaite et al., 2019). Current efforts to build a human body atlas (Regev et al., 2017) will provide important information on the baseline stromal, tumor and immune tissue heterogeneity. Tissue specificity is particularly important in immune cells because they both circulate and remain tissue resident; profiling of both populations can enhance understanding of which cells are recruited to become part of the TME. Despite the promise of scRNA-seq, it also has a key limitation: Enrichment and dissociation strategies that are used to separate single cells for profiling and capture cells of interest disrupt tissue organization and result in loss of spatial information (Elmentaite et al., 2019).

Single-cell computational frameworks have thus far focused in dissecting the ligand-receptor interactome based on expression of ligands and receptors across cell populations. The various computational frameworks take different approaches (Zhou et al., 2017; Cohen et al., 2018; Kumar et al., 2018). CellPhoneDB (Efremova et al., 2020) uses scRNA-seq to decode intercellular communication networks by inferring multimeric ligand-receptor interactions between cell states based on expression of a receptor by one cell state and a ligand by another cell state. By assessing which structural complexes are ubiquitously expressed in a tissue and not varied between cellular states, ligand-receptor cell specificity is assured. Because multi-subunit heteromeric complexes often swap subunits to ensure ligand specificity, expression of all subunits of a functional complex is required for downstream signaling. Vento-Tormo first applied CellPhoneDB to the study of early human pregnancy and interaction between fetal and maternal cells with a focus on trophoblast-decidua interactions that underlie common diseases of pregnancy such as pre-eclampsia and still births (Vento-Tormo et al., 2018). The authors uncovered receptor-ligand complexes that can modulate trophoblast differentiation. Moreover, they identified an association between three major decidual natural killer (NK) cell populations, their cell markers, and their expression of cytokines, chemokines and receptors, with different reproductive outcomes, providing a cell atlas to understand normal and pathological pregnancies (Vento-Tormo et al., 2018). In the context of CSCs, such a strategy could potentially be applied to dissect receptor-ligand interactions that help sustain stemness if achievable scRNA-seq resolution allows for the reliable identification of CSCs.

Using RNA velocity measurements that predict the fate of cell populations (splicing and non-splicing ratio) and provide a measure of stemness, we can potentially identify CSCs in a non-biased manner. Such strategy is particularly attractive when we consider the potential heterogeneity of the CSC compartment. This strategy was employed by Gautam et al. (2021) to generate a multispecies cell roadmap for human and porcupine ocular compartments. The ability to use organoids and enrich for stem cell like progenitors offers new possibilities to further decode CSCs.

Boisset et al. have manually microdissected 727 mouse interacting cellular structures and used scRNA-seq to create an interacting cellular network (Boisset et al., 2018). Using permutations, they have created a null distribution of random interactions between cell populations that permits identification of enriched or depleted interactions as compared to the background model. More recently, the authors of NicheNet have assembled prior knowledge on ligand-receptors and their signaling pathways using public resources and used expression datasets to generate weighted networks with regulatory scores that infer activity (Browaeys et al., 2020). Overall, these studies provide immense insight into the cell-cell interactome, but these studies all make several key assumptions: 1) That transcripts encoding receptors are translated into proteins that are translocated to the membrane, 2) that ligands are successfully transported out of the cells, and 3) that cells are in proximity with interacting partners in tissue space (Elmentaite et al., 2019). These approaches are now starting to be applied to extremely rare populations such as CSCs and yield great promise to depict their secretome.

Mass cytometry (cytometry by time-of-flight, CyTOF) has revolutionized human immune cell profiling by allowing the simultaneous measurement of >40 surface markers on a single-cell basis (Bracci et al., 2020). CyTOF is particularly valuable when analyzing samples with a limited number of cells such as tumor biopsies, longitudinal studies or response to therapy, e.g., in clinical trials (Gadalla et al., 2019). The high degree of multiplexing in the CyTOF is possible because antibodies are coupled to rare metal isotopes that provide a unique mass tag to each marker to overcome spectral overlap commonly observed in multiparameter flow cytometry. Recently, Gallad et al. validated CyTOF panel data against flow cytometry and demonstrated equivalent ability to identify cell subsets or perform phenotyping, but with smaller cell numbers (Gadalla et al., 2019). CyTOF barcoding strategies enable the marking cells from separate samples, which allows for the simultaneous measurement of multiple samples, thus eliminating multiple sources of sample-to-sample variability (Behbehani et al., 2014). CyTOF has been applied to characterize CSCs specifically in leukemias and lymphomas; however, this requires a priori knowledge of the proteins at the cell surface and available antibodies amenable to conjugation. More recent studies have used CyTOF in combination with CITE-Seq (cellular indexing of transcriptomes and epitopes by sequencing), thus creating a multimodal approach with simultaneous quantification of single-cell transcriptomes and surface proteins. For example, Yao et al. have demonstrated that B cell, monocyte/macrophage, and plasmacytoid dendritic cell abundance across three methods (scRNA-seq, CyTOF and CITE-Seq) is consistent, but T cell measurements have greater variability. Additionally, this trimodal analysis indicates that there is a good correlation across scRNA-seq and protein expression for highly expressed cell type markers (Yao et al., 2020). Because these three single-cell approaches enable the simultaneous identification of cell types, cell states and characterization of cellular heterogeneity at transcriptomic and/or protein levels, understanding the concordance of the measurements among these three modalities is of great interest. Collectively, these observations suggest that the of use multimodal omic datasets can open new avenues to dissect CSC-immune cell interactions.

A major limitation of these single-cell techniques is the lack of spatial transcriptomic information that is key for understanding the CSC niche. 3D intact-tissue RNA sequencing methods exemplified by STARmap allow for mapping of more than 1,000 genes in sections of mouse brain to define cell types and establish a roadmap for cell organization principles (Wang X et al., 2018). Alternatively, imaging cytometry (IMC) is emerging as a transformative technique that allows for multiparametric analysis of >40 protein markers in frozen and FFPE (formalin-fixed, paraffin-embedded) tissue samples (Chang et al., 2017; Bouzekri et al., 2019). IMC provides a unique window into structural features of the tissue under investigation, validates spatial ligand-receptor interactions and identifies the distribution of cell types within the tissue to build a spatial map (Bodenmiller, 2016). In breast cancer, for example, clear differences have been noted in markers like cytokeratin, HER2, E-cadherin and c-MYC within single tumors, highlighting the intra- and interpatient heterogeneity that can be captured with IMC (Giesen et al., 2014; Wagner et al., 2019). As a result, IMC has the ability identify each cell type in its environmental cellular context, thereby providing a unique spatially resolved view of the tissue and allowing more accurate inference of cellular functional states. The downside is that this approach is not yet suitable for high-throughput screening of ligand-receptor interactions due the limited scanned tissue area (1 μm), high sample cost, and paucity of computational platforms available for data analyses. Bodenmiller et al. have developed multivariate computational tools to visualize and analyze multiplexed images of human tissue sections generated by IMC (Schapiro et al., 2017). Another platform recently proposed by Greenald et al. utilizes a library containing a large collection of manually curated cells, TissueNet (one million cells from six organs and different imaging platforms) and a deep learning model to achieve human-like cell and nuclear segmentation (Greenwald et al., 2021).

Future IMC studies focusing on a well-defined tissue sections that can be layered with ‘omic information will be of great interest for dissecting the interplay of between cell populations (Bodenmiller, 2016). Whether IMC layered with ‘omic data can also inform studies of rare cells such as CSCs remains to be seen.

With the sea of available ‘omic data, there is an urgent need for multi-omic data integration to provide systems-level information. Clustering of multi-omics has been implemented at various steps in the data analysis pipeline (e.g., early vs. late integration) [see reviews (Ramazzotti et al., 2018; Rappoport and Shamir, 2018)]. A drawback to all these strategies is that all data are treated equally, which may not reflect the biology of a disease (Ramazzotti et al., 2018). Recently, Stuart el al. suggested harmonizing single cell data across distinct modalities by selecting anchors (a common set of features) between datasets to recover matching cell states (Stuart et al., 2019; Stuart and Satija, 2019). A second key challenge is that some methods implicitly assume that data heterogeneity arises mostly from technical variation and is of no biological importance. To overcome this limitation, LIGER (linked inference of genomic experimental relationships) uses integrative non-negative matrix factorization (iNMF) for data reduction combined with joint clustering to define a cell label and assemble a neighborhood graph (Welch et al., 2019). Gao et al. have subsequently expanded LIGER capabilities to iteratively incorporate multiple data modalities and massive datasets (Gao et al., 2021). Overall, these strategies recover cell states but do not zoom in deep enough to capture the interplay between cells or levels of ‘omic data.

To dissect the cellular niches, future strategies will need to integrate ‘omic data in a way that achieves a better understanding of cellular states and their contributions to the diversity of CSC phenotypes. Such innovations will be key if we are to design novel therapeutic strategies that can disrupt key cell-cell communications and lead to the demise of CSC populations.

Single cell experiments allow the study of tumors and their heterogenous cell populations that often coexist in an evolutionary continuum. Pseudotime trajectory methods (also known as trajectory inference methods) (Trapnell et al., 2014; Haghverdi et al., 2016; Setty et al., 2016; Saelens et al., 2019; Wolf et al., 2019) allow for the placement of cell populations within a trajectory based on their expression profiles. Because most single cell datasets are static and represent a unique snapshot in time, these cellular trajectories represent more a cellular states axis rather than real time, while encapsulating information from past and present. A challenge with the current pseudotime methodologies is that multiple dynamics can lead to the same cellular “state.” For example, cellular fate decisions such as cell death, proliferation, or differentiation can result from multiple causes. In order to understand the origin of these fates, we need to know the directionality of the biological process. Our ability to merge several biological snapshots in time, for example, from time series data, can help to decode the range of “real” existing states. Otherwise, the exploration of cellular states is currently being addressed in the field using dynamic inference methods. For example, Weinreb et al. used a dynamic inference framework named population balance analysis (Weinreb et al., 2018) that aims to limit the state search space by introducing biophysical constrains (cell density) to define cellular states. Conversely, Fischer et al. (2019) used RNA velocity information (spliced to unspliced ratio) as a mean to provide directionality to cellular trajectories. RNA velocity provides a measure of a cell’s internal compass (La Manno et al., 2018; Bergen et al., 2021) by quantifying nascent and mature messenger RNA. Here, we can also use matched multi-omic data integration (scRNA-seq and scATAC-seq) to identify key biological entities in each cell population and constrain the range of truly available cellular states to infer cell population dynamics. Future integration of multi-modal omics with spatial information derived from IMC will help to increase our spaciotemporal resolution, infer more accurate cell relatedness measures, and tackle key questions such as the extent of diversity of the CSC niche and which candidate targets are promising to reduce or re-wire CSC populations.

Figure 2 summarizes representative systems biology approaches that can decipher stemness of CSCs. Current and new technological advances are essential to probe different layers of CSC regulation. Moreover, there is an urgent need for deep cell phenotyping using insights gained from integration of available systems approaches. Cell subsets can be investigated for in silico isolation of CSCs and used to explore the degrees of stemness defined by the combined actions of signaling pathway activities, cell-cell interactions, and TF regulatory networks to formulate phenotypic recipes that drive cellular states. To mirror the fact that CSCs are exposed to a myriad of cellular environmental stimuli, perturbation experiments open the possibly to study the impact of genetic changes, zoom into cellular rewiring, and identify cellular vulnerabilities. A combination of these approaches may overcome current limitations in targeting CSCs and allow the identification of novel targetable genes.

In summary, CSCs reside in tumor ecosystems composed of a plethora of cell types communicating in ways that drive cellular phenotypes. Therapeutics that target cell-cell interactions have become a useful tool in clinical practice and can be considered for CSC targeting (He and Xu, 2020). Key examples include ipilimumab, which targets the CD28/cytotoxic T lymphocyte antigen 4 (CTLA4) interaction, as well as pembrolizumab and nivolumab, which both target the programmed cell death 1 (PD1)/PD1 ligand 1 (PD-L1) interaction (Topalian et al., 2012). In addition to these FDA-approved immune checkpoint inhibitors, other inhibitory immunoreceptors that are being evaluated as potential targets for clinical intervention include LAG3, TIM3, TIGIT, B7H3, CD93, CD73, adenosine A2A receptor and BTLA (He and Xu, 2020). Despite the clear success of immune checkpoint inhibitors in several tumor types, response rates are limited (10–30%), especially in solid tumors (Hamanishi et al., 2015; Sarkar et al., 2016; Dempke et al., 2017; Disis et al., 2019; Palaia et al., 2020). The limited response rates seen with these agents are likely due to the complex network of interactions involving multiple cell types present in the TME, including CSCs, which leads to phenotypic heterogeneity we do not yet adequately understand (Sarkar et al., 2016). CSCs themselves exhibit immunosuppressive properties such as expression of inhibitory checkpoint ligands, low expression of MHC-I molecules and natural killer cell (NK cell) receptors, and low or absent expression of innate immune receptors, which renders the CSCs resistant to killing by a variety of immune cells (Bayik and Lathia, 2021). Inhibitory activity of immune checkpoints are determined by the cell surface levels of inhibitory receptors, ligand interactions with those receptors, turnover processes, and posttranslational modifications that regulate signal transduction (He and Xu, 2020). Therefore, most recently, CSC targeting has been carried out with combinations of dendritic cell-based vaccines, oncolytic viruses, and immune checkpoint blockades (for review see (Badrinath and Yoo, 2019)). Better understanding of potentially targetable ligand-receptor interactions involving CSCs has the potential to tailor these strategies. Furthermore, such understanding has the potential to expand the targetable interactions beyond tumor-immune cell interactions to include tumor-nonimmune stroma interactions as well.

Successful application of singe cell technologies and the integration of different modalities are paramount for better understanding of the factors that determine whether CSCs will respond to various therapies or not. Improved ability to quantify ligand-receptor interactions, dissect their spatial relationship (co-localization) and assess their association with outcome appear to be needed if we are to devise individualized CSC targeting. Indeed, current analyses suggest that for tumor-immune targeting there may not be a single predictor of clinical response. Instead, it is possible that the strengths of multiple ligand-receptor interactions will need to be assessed in addition to immune checkpoint components in order to predict response to immune checkpoint blockade (Ribas and Wolchok, 2018). Systems biology techniques are capable of extracting this biological information from ‘omic data and helping derive a systems view.

In contrast to CSC cell surface markers that are expressed on stem cells, intracellular stem cell markers such as aldehyde dehydrogenase (ALDH) enzymes are an intracellular target amenable of intervention that can also enable some level of treatment individualization. ALDH activity identifies CSCs in numerous cancers (Deng et al., 2010; Landen et al., 2010; Silva et al., 2011). Recently, Raghavan et al. demonstrated that ovarian cancer spheroids derived from cells that survived chemotherapy (cisplatin) displayed lower ALDH expression, complete loss of CD133 expression, and resistance to cisplatin/ALDH inhibitor combination treatment while spheroids that were resistant to cisplatin/JAK2 inhibitor combinations contained an increased number of ALDH⁺ cells (Raghavan et al., 2017; Chefetz et al., 2019). Thus, stratification of patients according to their CSC type (i.e., CSC markers and expression levels) pre- and post-therapy can potentially lead to more personalized treatment approaches (Patsalias and Kozovska, 2021). Ultimately this approach will not only allow us to better stratify patients for targeted and combinatorial therapies, but also potentially identify novel targets in CSCs.

Targeting CSCs, as well as the phenotypic and functional heterogeneity of the bulk tumor cells derived from them, remains an unsolved clinical challenge. Phenotypic plasticity of CSCs fuels adaptive phenotypes that contribute to tumor chemoresistance and poor clinical outcomes. A second layer of complexity arises from spatially distinct CSC niches. In this review we focused on interrogating the tumor niche to zoom in on the resident CSCs and provide insights on potential pharmacological approaches. To date no scholarly review covers this specific topic under a unified systems biology framework and new insights on targeting CSCs and strategies to engineer TME are still limited. The use of ‘omics technologies as well as systems biology and computational methodologies has the potential to revolutionize the process of identifying and studying the CSC niche. While a wide array of methods have been developed to interrogate tumor heterogeneity, available methods still do not fully characterize the phenotypic states of CSCs or asses their plasticity and its impact on drug sensitivity. Therefore, there is an urgent need to develop tools that better integrate multimodal datasets to capture these cellular complexities. A crucial next step in the study of CSCs is to characterize their contributions to intratumoral heterogeneity, at various regulatory levels, from genotype to phenotype. If we are successful, the development of models to spatially localize CSC populations within tumors and dissect their contributions to tumor initiation, progression, and drug resistance will ultimately allow us to devise innovative strategies that target CSCs and improve cancer therapy.