Mechanistic models of blood cell fate decisions in the era of single-cell data.

Billions of functionally distinct blood cells emerge from a pool of hematopoietic stem cells in our bodies every day. This progressive differentiation process is hierarchically structured and remarkably robust. We provide an introductory review to mathematical approaches addressing the functional aspects of how lineage choice is potentially implemented on a molecular level. Emerging from studies on the mutual repression of key transcription factors, we illustrate how those simple concepts have been challenged in recent years and subsequently extended. Especially, the analysis of omics data on the single-cell level with computational tools provides descriptive insights on a yet unknown level, while their embedding into a consistent mechanistic and mathematical framework is still incomplete.

The production of blood cells as a dynamic equilibrium

The hematopoietic system is a primary example to illustrate how a population of resident tissue stem cells continuously contributes to the replenishment of a spectrum of mature cell types. While the dynamics and the particular contribution to the different blood lineages are tightly regulated processes, the overall theme, namely the differentiation of a pool of ‘undecided’, multipotent stem cells toward committed cell types, is remarkably robust. The accessibility but also the clinical relevance in the context of stem cell transplantation makes hematopoiesis probably the best-studied example of cellular decision-making in mammals.

On a more abstract level, the differentiation of hematopoietic stem cells into mature blood cells is the realization of a general principle that also governs the development of multicellular organisms, consisting of cells with distinct phenotypes and functions. During development, these different cell types arise from a common, totipotent cell (usually a zygote) in a sequence of decision steps. It is well established that a wealth of spatial, physical, and temporal factors influence this process, while the repertoire of developmental options is encoded in the genomic sequence. However, it is still one of the fundamental questions in biology, how this genotypic basis translates into a robust phenotypic decision process with adequate abundance, timing, and spatial orientation.

In the 1950s, Waddington [1] suggested a visualization of this decision-making process by comparing it to a marble rolling down a slope with emerging valleys and ridges and thus projecting the myriad of molecular, morphological, and functional features onto a three dimensional landscape. With the rise of molecular biology, this picture provided an interpretation of the valleys as cell types with distinct biochemical configurations, making Waddington's landscape a landmark visualization for stem cell differentiation. Already in the 1960s and 1970s, Kauffman [2] developed the formal concept that cell lineages can be viewed as attractors in gene regulatory networks and that cell differentiation represents a transition between those dynamic equilibria. However, understanding the nature of the branching point that initially separates two or more distinct lineages is still an ongoing struggle.

Stem cell decision-making and the concept of bistability

Mathematical modeling has long accompanied research on molecular decision-making [3, 4, 5, 6]. Corresponding models contributed particularly to the understanding of switch-like decisions in bacteria, as they commonly occur in response to changing environmental conditions [7]. In seminal work, Gardner et al [8] took a reductionist approach to demonstrate how a synthetic genetic toggle switch built on the mutual interaction of two repressible operons establishes bistability and exemplifies the minimal requirements of a binary molecular decision-making process, in which the choice of one option intrinsically represses the other. This simple but elegant example combined the mathematical concept with its experimental validation in a bacterial cell line.

Around the same time, efforts intensified to experimentally investigate the role of transcription factors in initiating and maintaining decision processes during hematopoietic differentiation [9,10]. Several studies revealed an antagonism between the transcription factor PU.1, which acts as a key regulator in myeloid and lymphoid cells and is constantly upregulated in granulocytic–monocytic progenitors, and the transcription factor Gata-1, which is required for the differentiation and maturation of megakaryocytic-erythroid cells (reviewed in the studies by Orkin [10] and Roeder and Glauche [10, 11]). The further discovery of direct, repressive interactions between the two factors PU.1 and Gata-1, as well as the identification of an uncommitted, upstream state, in which lineage-specific transcription factors are lowly expressed, led to a series of conceptual works to explore whether bistable switches can also serve as blueprint for decision-making in hematopoiesis [12]. Formalizing this concept in terms of ordinary differential equations, several models [11,13,14] predicted that the hematopoietic differentiation landscape qualitatively changes from a coexpression state of PU.1 and Gata-1 (usually referred to as priming) toward a bistable region in which the system converges to either of the two dominating states (indicating commitment, compare Figure 1). Such qualitative changes in the available system states (referred to as bifurcations) appear as a plausible conceptualization to account for progressive differentiation processes in tissue formation and maintenance and can also be facilitated by other molecular network motifs such as feed forward loops or positive autoregulation. The concept of direct interactions between potentially counteracting transcription factors has also been applied to related phenomena in hematology [15,16] and to other tissues [17,18].

Figure 1

While the PU.1-Gata-1 toggle switch model (left) has been successfully used to describe the concept of binary fate choice from a conceptual perspective, it evidently does not capture the full complexity of the early myeloid lineage decision. Single-cell gene expression analysis provides a better understanding of differentiation trajectories (often visualized as low-dimensional embeddings), while an explicit link between causative, mechanistic mathematical models and the high dimensional data is still an ongoing challenge.

Extensions to the toggle switch model

The initial ordinary differential equation–based models centered around the emergence of steady states through bifurcations and symmetry breaking, while sacrificing many of the molecular details. For the example of the PU.1-Gata-1 toggle switch, increasingly precise measurements revealed additional properties that had not been captured by the simple, conceptual models and fostered more sophisticated approaches. As first reported in the study by Warren et al. [19], only tens of PU.1 mRNAs are present in progenitor blood cells, prompting the development of stochastic toggle switch models [20], some of them also explicitly considering transcription and translation as two processes with distinct time scales [21,22]. Bayesian networks that included epigenetic and gene expression data successfully accounted for epigenetic changes during lineage specification [23]. Beyond the regulatory motif of only two transcription factors, Boolean networks can be used to represent major interaction axes and describe on/off gene expression states [24, 25, 26∗∗], including further upstream and downstream regulators. Modeling challenges remain, in particular the inherent nonstationarity of the system due to cell cycle that induces a continuous growth of cell volume and the number of cellular constituents, but also abrupt changes, for example, in the number of gene copies during DNA replication or during cell division.

Single-cell omics profiling of lineage decisions in blood

The single-cell omics revolution that we currently witness provides an unprecedentedly wide and detailed view on cellular differentiation (reviewed in the study by Watcham et al. [27]). After pioneer RNA sequencing of thousands of genes in thousands of single blood cells in mouse embryos [25] and adult humans [28], it is now recognized that differentiation shows many features of a continuous process instead of a sequential transition through clearly separable states [29,30]. To arrive at this conclusion and harness the large, high-dimensional single-cell data, we rely on computational methods. A variety of tools has been established to allow trajectory inference and pseudotime estimation in high-dimensional state spaces (compared in the study by Saelens et al. [31]), recently also exploiting splicing gradients in single cells to predict future gene expression through RNA velocity [32,33]. Nowadays, multiple tools are even bundled in convenient packages [34,35]. Single-cell RNA sequencing has recently also been used to link molecular profiles to stem cell functionality in transplantation assays [36] and to identify hierarchical mixed-lineage states preceding differentiation [37]. The combination of single-cell transcriptomics and lineage tracing methods can further reveal to which extent differentiation preferences are clonally fixed or subject to external influences [38]. Beyond gene expression, chromatin accessibility via single-cell ATAC sequencing in defined hematopoietic populations revealed heterogeneities in human progenitor populations and allowed analyzing regulation and expression of transcription factors [39], while CITE-seq allows us to measure both protein abundance and gene expression [40]. Notably, recently launched commercial kits to simultaneously profile gene expression and chromatin state or surface markers in the same cell will allow for increasingly precise analyses of cell identity. The combination of these two data types has also been used for the assembly of regulatory networks from gene–gene correlations [41, 42, 43∗] and points toward a more comprehensive approach beyond mRNA abundance that also integrates information of the epigenetic level.

Tracking hematopoietic decision-making with real-time imaging

In contrast to the snapshot data obtained from the simultaneous analysis of thousands of cells, quantitative live-cell imaging techniques with fluorescent reporters allow us to continuously track the abundance of a few essential regulators over time [44]. While RNA sequencing is powerful to identify differential expression patterns along a pseudotime trajectory, image-based tracking can measure gene expression in real time. Here, however, the analysis is limited by the number of labeled factors and separable fluorescence channels. The resulting data catch the dynamics of potentially interacting proteins within the same cell and provide essential information about timing, coexpression, and inheritance of protein expression but also about their variability between cells. While fluorescence labeling of critical transcription factors during hematopoietic differentiation yields great potential for dynamical assessment of the interaction network [45], few works took live-cell imaging data to the point at which mechanistic mathematical models could be evaluated [46]. In a recent work, single-molecule RNA imaging of PU.1, Gata-1, and Gata-2 was combined with time-lapse microscopy [47]. Stochastic modeling suggests that differentiation is preceded by a reversible transition between different coexpression states that most likely results from the intrinsic stochasticity of gene expression.

Linking single-cell data with mechanistic models

The notion of a continuous differentiation process together with the appearance of intrinsic molecular heterogeneity makes the identification and understanding of molecular principles governing hematopoietic decisions a fascinating and challenging endeavor. The conceptual notion of differentiation as a transition between dynamical attractors in larger gene regulatory networks has been conceptually discussed since the late 1960s [2,48]. The single-cell omics revolution largely confirms this interpretation and describes differentiation as a continuous trajectory in a high-dimensional state space. While this topological embedding is increasingly accepted, we need to further identify the ‘driving forces’ that direct cells along those trajectories and execute decision processes at the branching points.

Several works aim to link the quasi-potential Waddington landscape with regulatory networks in more formal [49,50] or data-driven [51] manifestations, but up to now, mathematical models are rarely applied to the rapidly increasing volume of single-cell data to foster a mechanistic understanding of the decision processes (Figure 1). While population approaches can be fitted to the overall gene expression trajectory of a differentiation process [52], a stochastic transition model (such as in Wheat et al. [47]) captures the sequence of coexpression states preceding differentiation onset. In contrast to descriptive approaches, a general framework to evaluate mechanistic interaction models with single-cell data faces a number of challenges: How can we identify the crucial regulators initiating and directing differentiation? How to infer the true differentiation time from snapshot data? And which model classes can sufficiently represent the expression and temporal heterogeneity observed in single-cell data? Modern machine learning methods, such as informed neural networks [53], graph neural networks, or autoencoders [54] might help to reduce the number of genes and features to a minimum and merge different data modalities and prior knowledge.

More generally, we may need to answer the question whether the low-dimensional representation of a differentiation process as proposed by Waddington and propagated by modern embeddings truly emerges from binary decision processes that are captured with the mathematics of bifurcation theory. There is no inherent reason why gene regulatory networks should favor the cross regulation of only two factors, especially if molecular redundancies can stabilize a decision process and make it more robust to variation. The theory of gene regulatory networks is a powerful theoretical approach to explain cellular plasticity and decision-making, but it needs to be extended to include other data types and to deal with variability, cellular heterogeneity, and higher-dimensional regulations. However, it remains to be shown how the phenotypic but also molecular data obtained from advanced sequencing and imaging techniques can truly be embedded into such an extended mathematical framework describing mutual interactions of key transcription factors on the mechanistic and biochemical level and whether our basic understanding of molecular regulations and interactions provides the right framework to account for the apparent heterogeneity and robustness of lineage specification processes.

Conflict of interest statement

Nothing declared.