Modeling heterogeneity in the pluripotent state: A promising strategy for improving the efficiency and fidelity of stem cell differentiation.

Pluripotency can be considered a functional characteristic of pluripotent stem cells (PSCs) populations and their niches, rather than a property of individual cells. In this view, individual cells within the population independently adopt a variety of different expression states, maintained by different signaling, transcriptional, and epigenetics regulatory networks. In this review, we propose that generation of integrative network models from single cell data will be essential for getting a better understanding of the regulation of self-renewal and differentiation. In particular, we suggest that the identification of network stability determinants in these integrative models will provide important insights into the mechanisms mediating the transduction of signals from the niche, and how these signals can trigger differentiation. In this regard, the differential use of these stability determinants in subpopulation-specific regulatory networks would mediate differentiation into different cell fates. We suggest that this approach could offer a promising avenue for the development of novel strategies for increasing the efficiency and fidelity of differentiation, which could have a strong impact on regenerative medicine.

(human/mouse) embryonic stem cell

hematopoietic stem cell

(pluripotency) gene regulatory network

(induced) pluripotent stem cell

transcription factor

Introduction

Experimental studies at the single cell level have revealed that embryonic stem cells (ESCs), and more generally pluripotent stem cells (PSCs), exhibit significant gene expression heterogeneity 1, 2, 3, 4, 5. The heterogeneity of gene expression in the pluripotent state has been studied in different model systems including mouse, rat, primates, and human ESCs 2, 6, 7, 8, 9. Although there is a lot of debate on whether this heterogeneity is an inherent feature of PSCs 2, the pluripotent space is characterized by a continuum of cellular subpopulations in vitro and in vivo, in which different subpopulations express different levels of pluripotency regulators, such as Nanog, Myc, Dppa3, and Rex1 2, 3, 5, 10. Single cell studies have allowed gene clustering, depending on the levels of heterogeneity in the gene expression landscape in PSCs. In mouse ESCs genes exist that are uniformly expressed in most cells and exhibiting a unimodal distribution (Oct4, Rest, Tcf3, Sal4); other genes exhibit bimodal expression and are expressed in some populations but not in others (Nanog, Rex1, Tet1, Esrrb), and yet another group of genes display sporadic expression (Neurod1, Klf4, Otx2, Pax6) and are undetected in most cells but highly expressed in some specific subpopulations 3, 5. Similar results were obtained studying heterogeneity in human ESCs and induced PSCs, in which pluripotency regulators such as Nanog and Oct4, exhibit similar expression patterns across different subpopulations to the expression patterns of their orthologous counterparts in mouse 2, 11. Gene expression heterogeneity has also been studied in hematopoietic stem cells (HSCs), in which gene expression uni‐ and bimodality have been observed for different genes in different hematopoietic progenitor cells 1, 12. Interestingly, it has also been demonstrated that genes exhibiting bimodal gene expression tend to be co‐expressed, and regulators such as Rex1, Nanog, and Esrrb display strong correlations 3, 5, which is related to the dynamics in the interplay among TFs for regulating pluripotency genes in ESCs 13. On the other hand, genes displaying sporadic expression exhibit more complex correlation patterns 5, suggesting that there exist no clear functional relationships among genes expressed in burst in some specific subpopulations.

The differential expression patterns of some pluripotency regulators in different subpopulations have been linked to fine tuning the balance between self‐renewal and cell fate commitment 7, 14, 15, 16, 17. The mechanisms underlying gene expression heterogeneity in the pluripotent state are still not well understood, and evidence suggests that gene fluctuations could arise from stochastic noise in gene expression regulatory mechanisms 18, 19, 20, 21. These fluctuations in gene expression are possible due to a permissive chromatin configuration across the genome of PSCs 22, allowing stochastic expression of lineage‐specific and differentiation‐promoting genes even before cell‐fate commitment 3. A central challenge nowadays is to characterize the transcriptional states of different phenotypic stem cell subpopulations, and the mechanisms that control their stability and promote transitions among them 5. Molecular studies of the regulatory mechanisms controlling pluripotency have shed some light on the pluripotency gene regulatory networks (PGRNs) 4, 23 formed by canonical pluripotency regulators such as Nanog, Oct4, Sox2 among others. Since the characterization of the initial PGRNs, more complex regulatory networks have been described including many other transcription factors (TFs) (Sall4, Zscan10, and Dax1) as key regulators of self‐renewal and differentiation of ESCs into different cell‐fates 4, 24. The chromatin landscape and signaling pathways also play a key role in regulating pluripotency and differentiation, although the cross talk between these regulatory levels during lineage specification has been much less studied 5, 25, 26.

Based on these considerations, it is feasible to view pluripotency as a functional property of cell populations rather than a well‐defined characteristic of single‐cells 27. Recent reports highlight the need for obtaining a deeper understanding of the mechanisms through which the microenvironment influences self‐renewal or priming of different subpopulations to differentiate into specific cell types 2, 3, 4. This is necessary in order to realize the biomedical potential of stem cells in regenerative medicine 4. In this regard, it will be essential to develop computational formalisms for performing systems level studies that go beyond the analysis of independent genes or groups of genes, but take advantage of the wealth of information generated by high throughput technologies. Despite the novelty of this systems‐level description of pluripotency, it is important to further explore how this view could help to address key questions in the field of cellular differentiation. Namely, which are the gene regulatory and signaling networks, and the epigenetic mechanisms controlling self‐renewal and cell‐fate commitment; and which are the mechanisms mediating the interaction of the niche with stem cells. Here, we argue that the implementation of integrative computational models of pluripotency is essential to get a better understanding of the regulation of self‐renewal and differentiation, especially taking advantage of recent developments in single‐cell technologies. We propose that the integration of signaling, transcriptional, and epigenetics levels into network models will be key to obtaining an accurate identification of cell fate determinants and cellular mechanisms regulating priming and self‐renewal. Finally, we discuss the importance of these systems‐level studies, which in combination with experimental approaches will be indispensable for designing novel strategies for increasing fidelity and efficiency of cell fate determination. These studies will have a strong impact on regenerative medicine.

Single‐cell profiling is key for studying pluripotent state gene expression heterogeneity

The groundbreaking developments of recent years for profiling gene expression and epigenetics in individual cells have significantly revolutionized molecular biology, making it possible to elucidate the mechanisms of cell‐to‐cell variability, and their implication in complex biological processes, such as differentiation. This field has been propelled by remarkable technological advances for efficiently isolating single cells 28, 29, and the reduction of the detection levels of protocols for profiling gene expression 30, 31, histone marks 25, 26, or chromatin accessibility 32, 33 in individual cells. These techniques allow researchers to overcome the limitations of population studies, in which averaging the measurements over a heterogeneous population of cells masks the variability at the epigenetics and transcriptional levels among cells within a culture or tissue 26, 34, 35. Moreover, the analysis of large numbers of single cells provides statistical power for predicting functional correlations among genes 36 in different cells within a heterogeneous population, hence facilitating the reconstruction of subpopulation‐specific gene regulatory networks that determine priming for commitment to different cell fates.

Single‐cell profiling has been extensively used for studying multiple biological systems to try to understand the underlying mechanisms regulating gene expression heterogeneity in the pluripotent state. Recently, single‐cell studies have helped unravel lineage specifiers triggering differentiation, pinpointing the genes taking part in lineage commitment during hematopoiesis 1. Other reports underscore the role of cell‐to‐cell gene expression variability in cell commitment in different hematopoietic progenitor cells, resulting in the independent activation of regulator genes in the absence of a coordinated lineage program 36, 37, which suggests that cell fate commitment can occur through multiple alternative pathways. Similar results have been obtained in other systems, such as the study of murine lung development, in which single cell transcriptomics data revealed cell‐type specific transcriptional regulators that discriminate between different populations that define the cellular hierarchy of the distal mouse lung epithelium 38. Single‐cell studies of the sub‐regions of the embryo, have also provided key insights of the initial phases of multicellular organisms development, allowing the identification of regulators triggering segregation between cell populations in early mouse embryos 39, and the delineation of gene regulatory mechanisms underlying progressive development of early mammalian embryos 40, 41. Profiling of the epigenetics landscape of individual cells has shown that lineage‐specific master regulators are associated with single‐cell epigenomic variability across several cell types, suggesting that control of single‐cell variance is a fundamental characteristic of different biological states 32. Moreover, another study compiled a comprehensive catalog of enhancer histone marks at the single‐cell level, and unveiled that there exists significant de novo establishment of lineage‐specific enhancers during hematopoiesis 25. Similarly, based on single‐cell profiling of histone modifications dynamics in different tissue‐resident macrophages populations, the same group showed that a combination of tissue‐ and lineage‐specific transcription factors form the regulatory networks controlling chromatin specification in tissue‐resident macrophages 26.

Pluripotent state gene expression heterogeneity is tightly regulated at different regulatory levels

A deeper analysis of the functional categories of genes in different expression groups in PSCs subpopulations shows that genes expressed in a unimodal and bimodal fashion tend to be enriched in housekeeping and metabolic functions 3. On the other hand, genes expressed in bursts are highly enriched in signaling proteins 2 and the Polycomb family of epigenetic regulators, and some of these genes are expressed at levels as high as known pluripotency regulators 3. These findings are supported by the results from another single‐cell study based on a recently developed experimental technique (inDrop) 42. This study shows that among the variable genes in mESCs are included pluripotency factors previously reported to fluctuate in pluripotent cells (Nanog, Rex1/Zfp42, Dppa5a, Sox2, Esrrb), and more strikingly, the most highly variable genes included known markers of Primitive Endoderm fate (Col4a1/2, Lama1/b1, Sox17, Sparc), markers of Epiblast fate (Krt8, Krt18, S100a6), and key epigenetic regulators of the ESC state (Dnmt3b) 42. The high variability in the epigenetic landscape has recently been studied at the single‐cell level in HSCs 25, in which significant chromatin reorganization in different subpopulations plays a key role during cell‐fate commitment. Similar observations have been made in human ESCs, where the dynamics of histone chromatin marks and DNA methylation is strongly associated to the binding of specific TFs, such as Sox17, Otx2, and Gata6, which defines and stabilizes the phenotypes corresponding to different germ layers 13. Moreover, the acquisition of a poising state – i.e. H3K4 nucleosome monomethylation – at enhancers in specific genes is essential for the ability of endodermal intermediates to respond to inductive signals during pancreatic and hepatic differentiation of human ESCs 43. In order to obtain a deeper understanding of the mechanisms regulating heterogeneity in PSCs, different studies have been performed based on the imposition of different growth conditions – e.g. inhibitors of signaling proteins, epigenetic regulators or gene knockout variants – to heterogeneous PSCs populations 3, 4, 5, 44, 45. Culturing PSCs in 2i medium – i.e. with inhibitors of Erk and Gsk3 signaling – significantly reduces expression variability in many genes, either by eliminating bimodality or by increasing their burst frequency 5, an effect also observed upon impairment of microRNA production 3. On the other hand, impairment of PRC2 function through the loss of the Polycom‐group protein Eed results in greater population heterogeneity across most genes. Furthermore, Dgcr8 and Dicer knockout mouse ESCs, in which microRNA maturation is impaired, are assigned to the ground or self‐renewing state 3. Among genes whose expression levels change more significantly in Dgcr8 and Dicer knockout mESCs cultured in serum + LIF (leukemia inhibitory factor), and wild‐type mESCs in 2i + LIF in comparison to mESCs in other conditions, are Myc, and methyltransferases Dnmt3b and Dnmt3l 3. Their reduced expression may be linked to the reduction of DNA methylation in PSCs 5. ESCs cultured in these conditions also exhibit significantly lower levels of H3K27me3 at promoters 3. Thus, these results underscore the crosstalk between the transcriptional, microRNA, and epigenetic levels in the regulation of heterogeneity in PSCs.

The interaction between signaling, transcriptional, and epigenetics levels plays a key role in regulating self‐renewal and cell‐fate commitment. Hence, it is essential to get a better understanding of the molecular mechanisms mediating the integration of signals from the niche, and how these signals modify chromatin and gene expression landscapes. The cross‐talk between signaling and transcriptional levels has been widely studied in vitro in PSCs cultured in different conditions. LIF cause the activation of the LIF/JAK/Stat3 signaling pathway, through the activation of some Janus tyrosine kinases (JAK), triggering the phosphorylation of Stat3 and the activation of the expression of several genes required for mESCs and hESCs self‐renewal 46, 47, 48, 49. The combination of inhibitors in 2i and 3i culture mediums, which act on some specific genes from the Wnt/β‐catenin signaling pathway, promote pluripotency in mESCs through the stabilization of β‐catenin that override the inhibitory effect of Tcf3 on some important pluripotency regulators, such as Oct4 and Nanog 50, 51. The TGF‐β/SMAD signaling pathway maintains self‐renewal in mESC through the BMP/SMAD signalling activation of Id family genes 52, while in hESCs is the Activin/Nodal/SMAD2/3 cascade the one responsible for promoting pluripotency 53. Interestingly, several studies have demonstrated that pluripotency regulators form intricate circuits at the transcriptional level 13, 36, 53, 54, and many of the TFs (Essrb, Klf4, Stat3, Tcf3) in these regulatory motifs are downstream effectors of the signaling pathways regulating self‐renewal and differentiation. Although the complete spectrum of signaling pathways regulating pluripotency has not been fully described 55, these results demonstrate the confluence of different environmental signals from the microenvironment for the regulation of pluripotency. It has recently been shown that within the heterogeneous population of Neural Stem Cells, there exist dormant subpopulations of cells that are able to enter a primed‐quiescent state before activation, which is accompanied by down‐regulation of glycolytic metabolism, Notch and BMP signaling, and a concomitant up‐regulation of lineage‐specific transcription factors and protein synthesis 56. There is also evidence of the cross‐talk between signaling and epigenetics layers mediated by post‐translational modifications. Some members of the JAK tyrosine kinase family, activated under LIF conditions, control histone H3Y41 phosphorylation at Nanog promoter, and also have global effects on heterochromatin in ESCs 57. The activation of Erk signaling pathway triggers a cascade of post‐translational modifications that potentiate changes in Polycomb‐2 complex occupancy and poising of RNA‐polymerase specifically on pluripotency regulators in ESCs 58. Recently, researchers have described the role of Activin–SMAD2/3 signaling pathway that cooperates with Nanog to recruit histone modifiers onto key developmental genes, which promote histone H3K4 trimethylation in pluripotency genes (Nanog, Oct4, Wnt3). This process is key for regulating self‐renewal and differentiation in human ESCs 59. In summary, these observations highlight the repertoire of mechanisms by which signals from the microenvironment are transmitted through the signaling layer to coordinate changes in the transcriptional and chromatin landscape in PSCs.

Subpopulation‐specific gene regulatory networks can be inferred from differential network analysis

One of the key challenges in stem cell research is to try to understand how different environmental cues are integrated at the transcriptional level for controlling self‐renewal and differentiation in PSCs. It has been suggested that different topologies of pluripotent gene regulatory networks maintain the heterogeneous molecular states characteristic of individual PSCs 60. Recently, different studies have addressed this important question by integrating single cell expression and perturbation data for identifying the underlying PGRNs 4, 16, 22, 61, which are involved in the stabilization of the pluripotent state and lineage commitment 62, 63. In most of these studies, topological analyses of PGRNs have identified feed forward loops (FFLs) linked to regulation of the balance between self‐renewal and differentiation. These network motifs, which combine positive and negative regulatory interactions – i.e. activation and inhibition, respectively – have great information processing potential, allowing the generation of alternative output solutions. The feedback structure of the PGRNs have been reported to be partially destroyed, with main variations in the loop between Nanog and other pluripotency factors, leading to a transient activation of different motifs that trigger differentiation 16. Incoherent FFLs have been inferred in mESCs subpopulations linking Oct4 and Nanog with their gene targets, suggesting an antagonistic interaction between these TFs for the regulation of pluripotency and self‐renewal 4. Based on quantitative modeling of these network motifs in the context of the PGRNs, the authors propose that any stalling or block during differentiation, cause a reversion of the primed cell back to the ground pluripotent state. In this model, the authors also found other important genes associated to pluripotency (Tbx3, Klf4, Esrrb, and Sal4), which constitute checkpoints in the PGRN, demarcating the transition toward differentiation 4. In another study 22, the authors analyze how PSCs filter out stochastic gene expression fluctuations, arising in the heterogeneous pluripotent state that include expression of lineage‐specific genes alongside pluripotency regulators 3. A detailed quantitative analysis of the response of the regulatory circuit formed between Nanog, Oct4, and Sox2 to the fluctuations of Brn2, which determines Nanog expression for triggering neural cell commitment 22, shows that Nanog responds to the fluctuations in Brn2 expression in different subpopulations of ESCs like a two‐state switch. In pluripotent cells expressing high levels of Nanog, fluctuations below a 100‐fold increase of Brn2 are considered as noise and do not cause a response of the cells. However, Brn2 induction above this threshold disrupts the Oct4/Nanog/Sox2 pluripotency complex, which generates a sharp response causing a significant silencing of Nanog, priming cells toward differentiation 22. It was also found that this circuit integrates Brn2 fluctuations based upon both magnitude and duration of the input, and the rapid response is highly dependent on Nanog transcript lifetime 22. An in‐depth analysis of the PGRNs derived from independently generated stem cell data sets (iPSCs and hESCs) has shown that gene expression variability is highest in network regions with fewer connections, and conversely, highly connected network regions exhibit the most stable, least variable expression pattern 64. These observations suggest that transition from self‐renewal to lineage commitment is accompanied by changes in the underlying network structure, such that genes become increasingly co‐regulated as the population become more sensitive to differentiation signals 64. Together, these results suggest that heterogeneity in PSCs plays a key role promoting transitions between metastable self‐renewing and lineage‐primed subpopulations, and in the integration of signals of the microenvironment or gene expression fluctuations to respond accordingly to environmental cues during differentiation.

Classical approaches to infer cell type specific gene regulatory networks rely on cell population data – i.e. gene expression and ChiP‐seq data – 65, 66, 67, and therefore, they give little insights into differences in regulatory interactions that underlie cellular heterogeneity in the pluripotent state. Indeed, the application of differential network analysis to study cellular differentiation is fundamental, since regulatory interactions inferred using population studies might not actually occur in individual cells, in which cell‐fate decisions take place 36, 68. In contrast, single‐cell based differential network analysis should give insights into cellular heterogeneity and possible regulatory interactions likely to be crucial for cell fate determination of individual progenitor cells. Nevertheless, as previously discussed, cellular states are characterized by different combinations of gene expression profiles of individual cells present in the population 2, 3, 4, 5. Indeed, in recent years there has been an increasing interest in dynamical differential network analysis approaches that are starting to replace static descriptions of biological networks 67, 69, 70, 71, 72, 73, 74. These approaches have been mainly used to infer and compare cell‐type/condition specific networks; and they rely on multiple considerations, such as combining literature‐based information with cell‐type/condition specific data – i.e. gene or protein expression, protein‐protein/DNA interactions, epigenetic states – 67, 70, 73. In particular, methods for differential gene regulatory network analysis have focused not only on detecting differences in individual regulatory interactions, but also in network motifs and regulatory modules across different cell types or cellular conditions 67, 75. More recently, differential network analysis methods, which rely on single‐cell expression data, have also been implemented. These methods have aimed at studying stem/progenitor cell populations during differentiation 12, 76, and at predicting lineage specifiers triggering cell‐fate commitment in different PSCs 74. Nevertheless, despite these attempts to follow a single‐cell based differential network approach to study the differentiation process in heterogeneous PSCs populations, more sophisticated computational models are needed to address relevant questions in the field, such as the identification of signalling pathways and their downstream effects for the activation of master regulator transcription factors and epigenetics mechanisms determining cell fate decisions in different PSC sub‐populations; and which are the perturbations that prime cell subpopulations for differentiation into specific cell fates. These models should not only consider static differences in network topologies between different cell subpopulations, but also the dynamic nature of regulatory interactions at different regulatory layers during differentiation, defined by deterministic, as well as stochastic events. In addition, they should explore in more detail the bridge between network topology and dynamics to study response in different cell subpopulations to fluctuating environmental cues.

Specific gene regulatory network circuits regulate the balance between self‐renewal and differentiation

As discussed previously, the pluripotent state can be viewed as a functional cellular state formed by various combinations of interchangeable gene expression patterns of different cells within the population. Hence, combinations of interchangeable GRNs are responsible for stabilizing different states that characterize different cell subpopulations in PSCs. During stem/progenitor cell differentiation, the dynamic expression of key regulatory genes and their corresponding regulatory interactions in these subpopulations govern their lineage specification 77. Furthermore, similar to cellular reprogramming 78, during differentiation there seem to exist specific primed subpopulations of cells that are able to differentiate to specific cell types via pre‐existing pathways defined by their initial transcriptional and epigenetic states, as previously suggested 68. In addition, other cells outside of this primed subpopulation can also potentially enter these pathways as a result of stochastic events triggering transitions among cell subpopulations 5, 79, and thus, over time, reach a differentiated state. Consequently, a cell population shift can occur from the pluripotent state to the differentiated state.

Computational models that rely on single cell data have been implemented to study heterogeneity at the pluripotent state, and to explore how this heterogeneity is linked to cell fate commitment. For example, a data‐driven approach has been very useful for identifying switch‐like changes in expression of key regulatory factors, sequential waves of gene regulation, and expression of regulators triggering differentiation 80. A similar method has been developed for deriving single‐cell latent variable models (scLVM), which allows the identification of undetectable subpopulations of cells that correspond to different stages during the differentiation 81. Moreover, modeling gene expression changes in individual cells has proven to be effective for distinguishing cell subpopulations close to fate commitment, and for identifying putative regulators of commitment and probabilistic rules of transition between subpopulations 82. In general, these approaches rely on the identification of changes in the expression trends of individual genes in cell subpopulations profiled at the single cell level, thus they cannot account for the coordinated contribution of genes to trigger differentiation. In particular, these data‐driven approaches do not offer mechanistic insights into the interactions among lineage specifiers for triggering cell fate commitment, and cannot model the cross‐talk of different regulatory layers – i.e. transcriptional, epigenetics, signaling – during differentiation. Hence, integrative network modeling could constitute an interesting alternative approach for overcoming these limitations, and derive more realistic and comprehensive models of heterogeneity in the pluripotent state. Differentiation of primed subpopulations seems to be mediated by specific GRN motifs. For example, at the population level it has been shown that transcription factor cross‐repression plays a key role in the regulation of cellular differentiation. Indeed, a regulatory circuit composed of two TFs that inhibit each other and activate themselves constitutes a molecular mechanism – i.e. toggle switch – (Fig. 1, blue box) that has been shown to determine cellular commitment and provides stability to transcriptional programs mediating binary cell fate choices 83, 84, 85, 86, 87. This has been observed in other cellular systems, such as the common myeloid progenitor (mutual inhibition of Gata1 and PU.1) 88 and embryonic stem cells (mutual inhibition of Oct4 and Cdx2) 89. Therefore, toggle switches not only determine binary decisions in the cell fate tree (Fig. 1, blue box), but they seem to also play a key role in the stabilization of each possible cell fate and the maintenance of the progenitor cell state. More recently, a “seesaw model” has been proposed for explaining mesendodermal and ectodermal specification of ESCs. According to this model, the balanced expression of Oct4 and Sox2 mutually activate each other for maintaining the pluripotent state. Interestingly, pluripotency can be achieved by simultaneous upregulation of Oct4 and Sox2 or other lineage specifiers that are involved in mesendodermal and ectodermal specification 90, 91. On the other hand, deviation from this balanced equilibrium induces differentiation into mesodermal or ectodermal cell types. For example, increased expression of Gata3 has been shown to up‐regulate mesendodermal genes in mouse ESCs, thus replacing Oct4 as a lineage specifier 92. Hence, pluripotency seems to be a functional state maintained by a balance between differentiation forces exerted by antagonistic groups, or “teams,” of lineage specifier transcription factors, whereas, the over‐expression of members of each group directs cell fate commitment (Fig. 1, blue box). Transcription factors belonging to each group cooperate among themselves and compete with transcription factors of the other group (Fig. 1, blue box). These two types of interactions (cooperation and competition) are determined by the GRN topology. In fact, the previously described toggle switch circuit constitutes a clear example of this situation, where each interaction group is composed of one gene (Fig. 1). However, in general, more complex circuit architectures underlie this mechanism, especially when antagonistic groups contain more than one gene product. Positive feedback loops, known as positive circuits 93, generalize the role of toggle switches in reprogramming or differentiation, including several pluripotency regulators 4, 16, 22, which regulate the balance between self‐renewal and cell fate commitment. The presence of positive circuits in a gene regulatory network is a necessary condition for the existence of multiple attractors – i.e. multistability – in the Waddington's landscape 93, including those corresponding to the stem/progenitor and daughter cells. Hence, positive circuits guarantee the stability of attractors, and their perturbations can induce transitions among them.

Figure 1

Differential use of regulatory motifs in heterogenous PSCs subpopulations. Different subpopulations within the heterogeneous pluripotent state (red, orange, and blue), exhibit characteristic signaling, transcriptional and epigenetic patterns, caused by gene expression fluctuations. These different signaling, transcriptional, and epigenetics networks (green box), and the differential use of regulatory motifs (blue box), mediate different responses to niche's signals. Hence, due to different topologies of PGRNs, some specific subpopulations are insensitive to specific differentiation signals, due to a lack of expression of the specific receptors, or because the corresponding signaling pathways are not activated. As there is extensive cross‐talk between signaling, transcriptional, and epigenetics regulatory layers, the differences observed in different PSCs subpopulations determines dissimilar responses to the same environmental cue. This can be summarized by the use of different “Stability Determinants” with different topologies and involving different number of genes. (blue box). In these circuits, genes belonging to the same regulatory “team” (nodes surronded by circles of the same color) regulate the metastable equilibrium established in individual cells. Due to these different regulatory frameworks, environmental cues trigger different responses in each cell subpopulation, and the overexpression of regulators from the same team leads to the differentiation into one specific cell fate. However, cell fate commitment is not deterministic, and gene expression stochasticity could trigger differentiation to an alternative cell fate (gray box). In this regard, cells from the orange subpopulation have a high likelihood of remaining in a naïve pluripotent state, while the blue and red cell subpopulations are primed to differentiate to Astrocytes or Neurons, respectively. Due to the combination of deterministic (maintained by the regulatory framework) and stochastic (caused by gene expression heterogeneity and genes that can be expressed in bursts in some subpopulations) events, the propensity to commit to different cell fates can change depending on the environmental cues.

Modeling heterogeneity in the pluripotent state will be essential for devising novel strategies to improve the efficiency and fidelity of differentiation protocols

In heterogeneous stem/progenitor cell populations, distinct interconverting subpopulations are characterized by different gene regulatory networks, and therefore, positive circuits (Fig. 1) 12, 36, 94. We suggest that different subpopulations employ distinct combinations of positive circuits to stabilize the stem/progenitor state by maintaining the balance between antagonist groups of lineage‐specific regulators. Consequently, differential use of positive circuits in distinct subpopulations seems to determine their cell fate differentiation propensities and efficiencies. Moreover, loss of the expression balance between these antagonistic groups directs cell fate commitment. In this context, experimental results support the concept that up‐regulation of one of the lineage specifiers in PSCs subpopulations – i.e. which might result from stochastic fluctuations in the expression of some specific genes or in response to environmental cues 2, 3, 4, 5 – can disrupt this balance, and prime cells for differentiation into specific cell types. For example, up‐regulation of Gata2 in HSCs subpopulations 1 primes these cells for megakaryocytic and erythroid lineage commitment, and over‐expression of multiple genes involved in cardiac differentiation in mesenchymal progenitor subpopulations primes them for cardiomyocyte differentiation 95. It is worth noting that GRN dynamics, which are determined by stochastic and deterministic events acting on the initial epigenetic and transcriptional state, allows inter‐conversion between different subpopulations 5, 96. Additionally, different PSCs subpopulations could differentiate to a common cell fate following separate differentiation paths. Alternatively, cells in the same subpopulation, depending on the perturbation sensed from the microenvironment (Fig. 1), could differentiate into multiple different fates.

Intense research in this field has allowed the construction of PGRNs, and the identification of network motifs regulating the transitions among different cell subpopulations in the heterogeneous continuum expression landscape 4, 16, 22, 61, 62, 63. So far, the identification of simple network motifs – i.e. positive circuits – (Fig. 1) has been useful to explain the mechanisms controlling self‐renewal and differentiation in different PSCs models 4, 22, 61, 63. However, a systems level understanding of heterogeneity in the pluripotent state, that includes the involvement of more complete PGRNs and more complex regulatory motifs associated with different subpopulations, will be essential for identifying the checkpoints and peripheral regulators – i.e. those not located at the core of PGRNs motifs – fine‐tuning the transitions among them, modulating the response to differentiation signals for exiting the ground state to commit to a specific phenotype 4, 22. Moreover, in order to direct cell fate commitment, an appropriate strategy could rely on shifting the cell population distribution toward the subpopulation state primed for a specific cell fate. We propose that the generation of integrative GRNs models 72, 97, 98, gathering information at the epigenetics, transcriptional, and signaling levels, will be essential for devising novel strategies for increasing the efficiency and fidelity of differentiation. Computational network‐based approaches aiming at providing a mechanistic description of the regulation of heterogeneity in the pluripotent state, and cell fate commitment, have recently been developed 74, 99, 100, 101. However, these methods are not only scarce, but they significantly differ in the size of reconstructed networks, ranging from small circuits 100, 101 to genome‐wide networks 74, 99. Recently we have proposed a computational approach for the identification of network stability determinants in GRNs, mediating reprograming in different biological systems 73, 102, 103. This network‐based approach, which relies on population‐average gene expression data, allows the reconstruction of genome‐wide GRNs for different cellular phenotypes 103, 104. Differently to other network models, in which the direction – i.e. source and target genes – and effect – i.e. activation or inhibition – of the interactions are unknown or inferred computationally, our approach includes experimentally validated interactions obtained from Metacore™ database, which give us the opportunity to perform thorough computational modeling in the reconstructed networks. From this analysis, it is possible to identify the network stability determinants, which correspond to interconnected positive and negative circuits, stabilizing the phenotype‐specific gene expression pattern 73, 102, 103. This network‐based approach will significantly improve taking advantage of single cell expression data, allowing the analysis of homogeneous subpopulation‐specific gene expression data, for deriving the GRNs that determine fate commitment decisions in individual cells. Recently, we have extended this rationale to predict lineage specifiers performing a differential network analysis on single‐cell gene expression data in different biological systems 74. In contrast to previous computational methods for analyzing single cell data for identifying lineage specifiers 80, 81, 82, which rely on the identification of trends in the expression of individual genes in different subpopulations, in this study we have demonstrated that information of the interactions among genes in the subpopulation‐specific GRNs, is key for predicting the genes triggering differentiation reported experimentally 74. These network models 73, 74, 102, 103 relying solely on transcriptomics data could be significantly improved by overlaying epigenetics data 25, 32, 33, for contextualizing regulatory interactions at the transcriptional level depending on the chromatin state – e.g. open/close chromatin regions, and activating/inhibitory chromatin marks patterns 25, 26, 32, 33. Moreover, although it is not yet possible to perform single cell measurements of signaling pathways, the inclusion of phosphoproteomics data in these network models will include another layer of information for contextualizing network interactions in this layer, depending on the post‐translational modifications determining the activity of signalling proteins 105, 106.

More realistic integrative network models will allow the identification of the network stability determinants 4, 72, 75, 102, 103 sensing and transducing niche's signals, and how these signals trigger differentiation in specific PSCs subpopulations (Fig. 1). Thus, this computational approach will allow the identification of specific PSCs subpopulations, and those genes that could be perturbed to rationally shift the dynamic equilibrium in the pluripotent state, towards a primed PSC sub‐population for differentiation into a specific cell fate. Evidence exists for the role of oscillatory expression of multiple TFs (Ascl1, Olig2, and Hes1) for the regulation of self‐renewal in Neural Progenitor Cells 107, which is mediated by oscillations in the Notch signaling pathway 108. Moreover, it has been demonstrated that fluctuations of pluripotency and developmental regulators (Nanog, Gata4, Sox17, among others) during cell cycle are essential for regulating pluripotency and triggering cell‐fate commitment in human ESCs 109. However, the effect of spatiotemporal oscillations of signaling pathways, and how these fluctuations influence the epigenetics and transcriptional levels in heterogeneous PSCs populations has not been studied in detail. Interestingly, there are significant differences in the efficiency of protocols for deriving differentiated cell types. While dual SMAD inhibition is sufficient to induce rapid and complete neural conversion of human ESCs and iPSCs 110, there are limitations to efficiently differentiate PSCs to other cell types, such as cardiomyocytes 111, kidney cells 112, retinal pigment epithelium 113, or hepatocytes 114. In particular, the differentiated cells produced are largely immature, and resemble the fetal stages of development 111, 114, which hampers their ability to engraft and function in coordination with other cells of the tissue. In order to overcome these problems and provide more insights into cell‐fate commitment, computational modeling of single‐cell data will be essential for uncovering the mechanisms regulating self‐renewal, and triggering differentiation from a systems‐level perspective. However, the analysis of great quantities of data generated at the single cell level faces many challenges. Single‐cell measurements have a significant inherent noise 115, 116, as well as experimental drawbacks associated to batch effects and technical variability 117, 118, 119, and the computational methodologies for accounting for these statistical analysis limitations are still under development. Moreover, the reconstruction of GRNs from single‐cell data performed in several studies has yielded rather small networks, involving only a limited group of regulators 12, 99, 120. Despite the above‐mentioned limitations and the need for improvements in the statistical analysis of single‐cell data 118, 119, the generation of more comprehensive and accurate measurements of heterogeneity in the pluripotent state will allow the generation of integrative network approaches, such as the one proposed here. These approaches will be helpful for generating models for differentiation of PSCs into specific cell types, pinpointing the signaling pathways, transcriptional factors, epigenetic regulators, and the cross‐talk among them, that could be targeted for improving the efficiency of cell‐type specific differentiation protocols. In comparison to data‐driven computational models 80, 81, 82, integrative networks models will provide key mechanistic insights of the regulation of heterogeneity in the pluripotent state, and a systems‐level picture of the interactions among key regulators triggering differentiation to specific cell fates, which will be essential for devising more effective differentiation protocols.

Conclusions and outlook

In summary, heterogeneous PSCs populations comprise a continuum of distinct cellular subpopulations, characterized by different signaling, transcriptional, and epigenetics states, which underline their corresponding gene regulatory networks. These subpopulation‐specific networks display different topological features, and therefore include distinct combinations of network regulatory determinants involved in cell fate determination (Fig. 1). A differential use of these regulatory determinants can mediate differentiation to distinct cell fates. Moreover, experimental results indicate that PSCs may exit the pluripotent state via a continuum of intermediate states, which ultimately become primed for lineage specification 2, 79. In this regard, there is compelling evidence of the biological role of NSCs heterogeneity in vivo, which are activated in the brain in response to ischemia through interferon‐γ signalling activation of dormant NSCs subpopulations, to enter a primed state 56. Recent reports have shown the translational potential of PSCs for in vitro or ex vivo generation of differentiated cells that could be used for repairing damaged tissues 121. In this groundbreaking study, the authors proved that pluripotent cells cultured ex vivo can be successfully transplanted to form a corneal epithelium that recovers function in an experimentally induced animal model of corneal blindness. Thus, the study of the gene regulatory networks, signaling and epigenetic mechanisms controlling heterogeneity in the pluripotent state, and the generation of integrative models from single cell data will be key for getting a deeper understanding of self‐renewal and cell‐fate commitment. These integrative computational approaches will help devising novel strategies for increasing the efficiency and fidelity of differentiation, which will have a strong impact in regenerative medicine.

The authors have declared no conflicts of interest.