Mapping, modeling, and reprogramming cell-fate decision making systems

Multistability is crucial for cells to adopt and maintain distinct fates. The compartmentalization of key signaling proteins within cellular structures, such as the cytoplasm and nucleus, creates spatial organization that supports this process. In eukaryotic cells, many signaling proteins actively shuttle between the cytoplasm, where they are activated through processes like phosphorylation, and the nucleus, where they can trigger transcriptional responses, before returning to the cytoplasm. This cyclical cytoplasm $\longrightarrow$ nucleus $\longrightarrow$ cytoplasm shuttling acts as a positive feedback loop that promotes multistability [40]. By compartmentalizing signaling pathways, cells can better process and retain information, enabling them to stabilize their state even amidst fluctuating signals. For example, MAPK proteins, which move between these compartments, have been shown to influence stem cell differentiation in vivo [41], highlighting that compartmentalized signaling is a key mechanism for establishing and maintaining cell-fates. This design principle was extracted by extrapolating the mathematical analysis of positive feedback systems to more general systems [42]

Robust perfect adaptation (RPA) is a phenomenon in which a stimulus triggers an initial response, often seen as an increase in the concentration of certain molecules, which then returns to baseline levels, even while the stimulus continues. This adaptive response allows cells to reset their signaling activity and avoid overstimulation, a crucial feature for maintaining stable cell states in the face of persistent signals. Aruajo and colleagues [43, 44] were able to develop precise algebraic criteria to define the network structures that can achieve RPA. There are only two core network configurations that are capable of producing this adaptation, meaning that all biological networks exhibiting RPA must adopt one of these designs. This insight into RPA’s design principles was achieved through abstraction, which allowed the identification of underlying structural requirements without exhaustive trial-and-error. In the context of cell-fate decision making, RPA provides cells with the ability to respond selectively to new signals while maintaining a stable fate, making it an essential signal-processing mechanism in both biological and technological systems.

Toggle switches allow cells to make binary decisions, often essential for specifying a cell-fate among competing options. Bistability, the ability to adopt one of two stable states, enables cells to lock into a specific fate when exposed to certain signals. Stochasticity is one of the hallmarks of many toggle switch architectures. Noise can help push the system between states, and many designs are only bistable in a stochastic regime. In development, toggle switches help cells decide between pathways, like becoming either a neuron or a glial cell. This principle is also seen in microbial population dynamics, where toggle switches provide population-level heterogeneity, supporting bet-hedging strategies that allow populations to survive in changing environments [45]. Efforts to identify the design principles of toggle switches have relied on a combination of modeling and statistical model selection. However, there remain many different designs that can exhibit toggle switch behavior [37]. This diversity suggests that toggle switches may represent a class of systems with so many possible implementations that it is difficult to define a single set of design principles [46, 37].

The final example we discuss here are Turing patterns, defined by Turing to be a purely chemical process by which stable spatial patterns can emerge. It is one of the fundamental patterning mechanisms in developmental biology, alongside positional information, rule-based mechanisms, and phase separation. Arguably positional information and Turing mechanisms have received the greatest attention [47, 48, 49, 50, 51], but were taken up primarily by the biological and mathematical communities, respectively. Over the past decade, however, there has been a convergence of these perspectives, with both mechanisms now widely recognized as working together to create the developmental diversity observed in nature [51, 31].

Turing’s and later work summarized the design principles of these patterning mechanisms as: “local activation and lateral inhibition” [52, 31]. In this process, a molecule $A$ in one cell increases through positive feedback, diffuses to neighboring cells, and activates $A$ in those cells as well. At the same time, it also activates a second molecule $B$, which diffuses and inhibits the production of $A$ in neighboring cells. This interplay of activation and inhibition enables stable spatial patterns to emerge across cell populations.

We use this to refer to the molecular census and structural arrangement of all the molecules comprising a cell. In practical terms, not everything will be experimentally accessible. In current applications in mammalian systems, including stem cell systems, we have access to mRNA measurements at the level of individual cells. We currently do not have the ability to sample protein abundances, epigenetic states, metabolome concentrations, and more, at the sufficient scale and quality. As a result, current Cell Atlas projects provide only a partial view of cell states.

Imaging data offers a complementary perspective by revealing cell morphology and the spatial organisation of structures like organelles and molecular complexes. Combined with molecular census data, this allows us to build a more integrated view of cell states.

Many definitions of cell states overlook or minimize the influence of environmental feedback, which crucially determines the cellular physiology. To address this limitation, we might consider defining a “renormalized” cell state, encompassing both the molecular inventory and structural arrangement within a cell, along with interactions with its environment.

A cell type can be usefully defined as an ensemble of cells that, from a physiological or biological perspective, can be considered as equivalent. In mathematical terms this could, for example, correspond to region in cellular phase space that, under suitable thermodynamic conditions, is enriched for experimentally observable cells.

Cell types are traditionally defined using morphology or cellular phenotypes. The extent to which these classical definitions of cell types align with molecular signatures remains an open question, one that we will revisit throughout this review. This question is also closely linked to understanding the sources and propagation of molecular noise within cells

We can usefully differentiate between fully differentiated cell types and cell types that have some differentiation and proliferation potential. The latter include all stem and progenitor cell types. The former include all cells which no longer undergo a cell cycle, and where the only possible fate is cell death.

This term refers to the cell type that can be reached from an initial cell type, with an inherent temporal aspect, representing a path (or set of paths) through cellular phase space between cell types. Along this path, individual cells can be observed moving in one or both directions.

Important questions include when is a cell’s fate determined at the molecular level versus the phenotypic or morphological level. Similarly important questions include how many potential fates or paths are open to a given cell type; which cell types can be reached from which others; whether cell-fates are reversible, and if they are, is the same path through cellular phase space traveled in both directions.

We use the term “cell-fate decision” as a shorthand for a complex set of processes, implying that cells possess some level of information-processing and computational ability, as well as agency: they respond and react to internal and external signals and cues and change their behavior accordingly. There are philosophical discussions to be had on this, we are sure, but we eschew these here.

To clarify the concepts of cell state, cell type and cell-fate developed above, we can adopt a mathematical perspective [53, 54]. Let $X$ denote the cell state, encompassing the abundances, and arrangements of all molecular components within the cell. We can define a suitable phase space, $\Omega_{X}$, in which $X$ resides,

$$X\in\Omega_{X}\subseteq\mathbb{R}_{0}^{+n}.$$

The temporal evolution of $X$ is described in terms of a stochastic differential equation (SDE),

$$dX=\overset{\text{\footnotesize\color[rgb]{0,0,1}Deterministic Dynamics}}{\overbrace{\color[rgb]{0,0,0}f(X;\theta,t)dt}}+\overset{\text{\footnotesize\color[rgb]{0,0,1}Stochastic Dynamics}}{\overbrace{\color[rgb]{0,0,0}g(X;\theta,t)dW_{t}}}+\overset{\text{\footnotesize\color[rgb]{0,0,1}External Influences}}{\overbrace{\color[rgb]{0,0,0}h(X;\Psi,t)dt}}.$$

(1)

Here $f(X;\theta,t)$ represents the deterministic dynamics and $g(X;\theta,t)$ the stochastic dynamics, where $W_{t}$ is a Wiener process capturing the inherent noise [55]. The term $h(X;\Psi,t)$ describes the (assumed to be deterministic – this can be relaxed) influences from the environment (e.g. a change of growth-medium). Parameters $\theta$ define the model, while $\Psi$ governs the environment. A cautionary note for readers with experience in dealing with SDEs: in chemical or molecular processes the functional forms of $f(X;\theta,t)$ and $g(X;\theta,t)$ are tightly linked; hence the noise depends explicitly on the system’s state, $X$. We therefore often refer to both terms jointly as the model $m$, which summarizes our understanding and hypotheses about the underlying processes in a biological system, here focusing on cell-fate decision making.

This model $m$ captures cell-fate dynamics and is often associated with Waddington’s epigenetic landscape [56, 57, 55, 58]. These models range from simple, qualitative “back-of-the-envelope” types to more comprehensive frameworks. They may be designed to align with experimental data [59, 58, 60] or to explore broad phenomenological principles [61, 62, 57]. At the most detailed level, a whole-cell model, termed a ”CellMap” by Sydney Brenner [63], provides mechanistic, quantitative, and testable insights into cellular processes.

While we focus here on SDE models of the form in Eqn. (1), other mechanistic modeling frameworks can provide valuable perspectives. These include Boolean and generalized logical networks, stochastic Petri nets, agent-based models, and interacting spin-systems. Much of what we discuss here can be applied, with minor adjustments, to these other modeling frameworks. What they share is a foundational mechanistic structure that enables us to test and refine our understanding of the underlying mechanisms driving cellular behavior.

Mathematical models enable us to explore fundamental questions in cell-fate decision making systems such as:

1. 1.

   How many cell types exist?

2. 2.

   How heterogeneous can cells of the same type be? Or in other words, how many states and how varied can these states be within a single cell type?

3. 3.

   How stable are cell types?

4. 4.

   What causes a cell to change its type and make a cell-fate decision?

We cannot answer these questions without a tight integration of experiment, observation, and theoretical analysis.

We now have the DNA sequences of many species, and for several, we also have high-resolution data on genomic diversity, revealing genetic variability within populations. Evolutionary and comparative genomics, supported by bioinformatics tools, enable us to annotate and interpret these genomes, especially for newly sequenced organisms.

In the context of cell-fate decisions, the heritable genome plays a lesser role, except when comparing different species, such as embryonic development across mammals or varying phenotypes in yeast species, where even closely related species can exhibit very different phenotypes. Heritable genomic elements—such as gene sequences, regulatory regions, and arrangements—shape many molecular and phenotypic differences, influencing factors like gene activation timing, mRNA processing, and protein isoform diversity. However, not all sequence-level differences are shaped by selection pressures [68], and some aspects, like the full diversity of protein isoforms, remain largely speculative [69].

Since the first draft of the human genome was published, our understanding of genome content and the diverse functions of non-coding DNA has transformed. We now recognise, though still incompletely, that the genome is regulated by a complex epigenetic machinery that actively opens and closes DNA sequences, modifying access to coding regions [70].

These epigenetic mechanisms, while reversible, are sensitive to environmental factors and influence which genes are accessible for transcription, thereby affecting cell fate.These include DNA methylation and demethylation, histone modifications, and chromatin remodeling. Methylation is the addition of a methyl group (CH₃) to DNA (or proteins), often suppressing gene expression by limiting transcriptional access, playing a crucial role in maintaining or altering cell states. Histones are the proteins around which DNA is wrapped and the tightness of this wrapping also leads to decreased accessibility of DNA to the transcriptional machinery. Chromatin remodellining adjusts the structure of chromatin to reveal or hide specific DNA sequences, directly impacting which genes are active in determining cell identity.

A raft of experimental techniques are available to study these epigenetic and dynamic aspects of the genome and its in vivo geometry. ATAC-seq is used to understand patterns of DNA methylation and to analyze chromatin remodeling. Histone modifications are assessed through CHiP-seq (Chromatin immunoprecipitation) and similar assays. Hi-C (chromsome conformation capture) is widely used to assess the 3D structure of the genome, but other methods, such as genome architecture mapping and related approaches give a more direct representation of 3D genome structure, including contacts between chromosomes. Real-time methods like live-cell imaging and single-molecule tracking operate at a lower throughput compared to other techniques. Unlike most methods, they allow us to monitor the same cells over time, though they are typically limited to tracking a small number of molecules. Collectively, these tools help us connect the physical and dynamic organisation of the genome to its functional role in cell fate decisions, enhancing our ability to map and interpret complex epigenetic landscapes.

We refer to [72] for a recent review of descriptive methods for single-cell analysis. But we like to issue a note of caution: the high dimensionality, (tens of) thousands of genes across (tens of) thousands of cells has given rise to a suite of dimensionality reduction techniques. Arguably their popularity is tied to our preference to visualize data and this is only possible in low dimension. Techniques such as t-SNE and uMap are producing deceptively pretty plots and it can be tempting to apply quantitative interpretations to such visualizations, but this should be avoided as the purposes of analysis and visualization of high-dimensional data diverge and can be incompatible.

Networks, as discussed above, offer us the opportunity to study interacting sets of genes/proteins and gain some systems-level insights. Single-cell mRNA has given rise to a new set of network inference methods. They take the matrix of gene $\times$ cell gene expression measurements and seek to identify sets or pairs of genes that are statistically informative about each other’s expression levels. Network inference has been a notoriously challenging problem in systems biology [87, 88, 89]. At the core sits the large-$p$-small-$N$ problem, i.e. the fact that we have more potential hypotheses of the presence of pairwise interactions, than measurements $N$. Single-cell data with its abundances has overcome much of this.

Two types of approaches have been successful: approaches that are robust to noise and which can capture non-linear dependencies [90]; and flexible statistical or machine learning approaches, especially if combined with additional data [91]. The most promising methods still have appreciable error rates but serve as interpretative tools or hypotheses-generation methods for further testing. Crucially these methods, especially for carefully designed experiments, allow us to infer networks corresponding to different cell types or developmental stages. Capturing the complex temporal dependencies among cells (using e.g. geometric perspectives) allows us even to generate cell specific networks [92]. This field is still in its infancy, but will likely result in more detailed testable hypotheses about the workings of biological processes.

There have been mathematical approaches to cell cycle, cell signaling, and gene expression among others. One of the recurring themes has been the origin and maintenance of cell-to-cell variability and its impact on cell-fate decision making. We distinguish between two types of noise: intrinsic noise is due to the random effects affecting molecular reactions and is captured by stochastic modeling approaches; the second term in Eqn. (1) aims to capture these effects. Extrinsic noise, the second component we consider refers to the variability that is due to factors that are not explicitly captured in our model. This can, for example, include the cellular environment but also cell cycle stage, differences in ribosome abundance, or mitochondrial activity between different cells [93]. Capturing both types of noise is important; dissecting sources of noise is, however, not always possible, though some experimental designs allow this in principle and in practice.

In addition to SDEs of the form  (1), we can also employ other modeling frameworks, both exact and approximate, to simulate and analyze molecular processes. For a limited number of cases exact solutions or statements can be derived, but simulation is currently the only route for many situations. We are getting better at describing the design principles for noise propagation and attenuation, but both further experimental and theoretical work is required. For integrating transcriptomic and protein levels stochastic modeling can already help a lot and we can explore different alterations and refinements to the model and its dynamics, for example, by including opening and closing of chromatin, degradation of mRNA and protein, or threshold effects [94].

To model developmental processes biologists and mathematicians often turn to Waddington’s epigenetic landscape [53, 3] – a metaphor in which cell differentiation is visualized as a journey through a landscape of valleys and hills. Cells start as less differentiated at the top of the landscape, move downhill, and encounter branching points along the way. These branches represent cell-fate choices, with each path leading to a distinct valley corresponding to a stable cell type. This metaphor, introduced nearly a century ago, has since inspired quantitative and computational frameworks to describe key qualitative dynamics of developmental systems, such as the stability of cell types, bifurcations leading to different cell-fates, and the role of intermediate cell states. Over time, approaches to modeling the developmental landscape have evolved into two broad categories: “gene-centric”[59, 58, 60] and “gene-free models”[61, 62, 57], each offering distinct perspectives, advantages and disadvantages.

Gene-centric models are typically rooted in SDEs of the form (1) that describe the evolution of gene expression profiles over time. By presuming a network of gene interactions, these models seek to understand how individual genes and their interactions shape the developmental landscape. Specifically, the probability $P(X)$ of a cell’s gene expression vector $X$ is used to determine valleys (stable states) and paths (transitional states) within the so-called “quasi-potential” landscape:

$$U(X)\propto-\log(p(X)).$$

(3)

The simple picture conceived by Waddington and analyzed by many since then, has provided valuable insights into the molecular mechanisms underlying stability and transition points within the developmental landscape. These models are effective at capturing the influence of specific gene interactions on cell behavior, which has been important for understanding the roles of critical genes in driving cell-fate choices and maintaining stable cellular states [95, 96].

As we connect this metaphor to the underlying mathematics and biological data, however, certain limitations surface. For example, these models typically rely on the a gradient assumption, where the drift component of the underlying SDE is derived from the gradient of a potential function. The system’s dynamics are thus guided by an “energy-minimizing” behavior, with cells tending to move “downhill” in this landscape to reach stable states. This gradient-based perspective has limitations, especially in biological systems with complex, oscillatory, or non-conservative behaviors. Such processes do not strictly follow energy-minimizing paths, which means that they cannot be fully captured by a gradient potential alone. Moreover, the gradient assumption is typically applied within a quasi-steady-state context, where gene expression dynamics are presumed to stabilise over time. This assumption may not hold in developmental systems that undergo active transitions, bifurcations, or shifts between cell states [97].

Stochastic dynamics, intrinsic to cellular behavior, can reshape the landscape both quantitatively and qualitatively, introducing elements such as noise-induced transitions. These can drive transitions across the landscape in a way that diverges from simple “downhill” behavior. Differential equation models, however, typically use deterministic frameworks or assume additive noise. Gene-centric models also overlook cell-cell interactions and spatial or temporal factors, which are critical in many developmental processes but are challenging to integrate into a fixed gene network framework.

In parallel to connecting models of landscapes to gene regulatory systems, recent approaches explore the structure of these landscapes without relying on specific assumptions about gene networks. This shift, pioneered in part by the work of Briscoe, Siggia, and Rand, [61, 62, 57] leverages the work of renowned mathematician René Thom, who had a keen interest in applying catastrophe theory to biological systems. This “gene-free” perspective diverges from traditional gene-centric models, instead focusing on the global structural features that characterize developmental processes. These approaches seek to describe the behavior of the system as a whole, using concepts like stability points and saddle nodes to study the dynamics of cell state transitions. Rather than anchoring models in predefined gene networks, this framework captures the system’s global behavior, focusing on the topological and geometrical features that shape developmental landscapes.

Despite their innovative approach, these methods rely on the same gradient-based assumptions as gene-centric models and are thus subject to similar limitations when connecting these with data. However, they do offer several advantages. By not relying on specific gene interactions, gene-free approaches provide flexibility in capturing the complex and adaptive behaviors of developmental systems, accommodating phenomena that are difficult to model within a gene-centric framework. The structural features identified through catastrophe theory, such as bifurcations and attractors, offer robust insights into cell state stability and the conditions that lead to cell-fate transitions. Such methods lend themselves well to high-dimensional single-cell data, where the explicit modeling of all gene interactions is computationally prohibitive and often impractical due to incomplete information on gene regulatory networks. As the field progresses, integrating these approaches with empirical data may illuminate previously uncharacterized aspects of developmental biology, such as transient cell states and non-equilibrium dynamics, offering a pathway toward comprehensive, system-level models of development.

Ultimately, combining gene-centric and gene-free models may yield a more comprehensive and nuanced understanding of developmental processes, accommodating both gene-specific interactions and broader system-level behaviors. This combined framework has the potential to capture a richer, multidimensional picture of development, opening new pathways for exploring cell-fate determination and dynamic regulatory landscapes.