Introduction

Waddington’s epigenetic landscape is a fundamental and profound conceptualization of cell fate decisions [1]. Over decades, this insightful metaphor has facilitated researchers to distill a myriad of models regarding cell fate decisions [2–9]. While introducing various quantitative models and dissecting diverse fate-decision processes, researchers have further elaborated the Waddington landscape [10–15]. An outstanding question is whether the landscape is static or not, i.e., whether cell fate decisions are driven by noise or signal [14, 16, 17]. On one hand, some perspectives hold that cells reside in a stationary landscape, where decisions are made by switching through discrete valleys [18, 19], as a result of gene expression noise [20, 21]; termed as “noise-driven”, Fig 1.A). Meanwhile, some researchers argued that the epigenetic landscape is dynamic during fate decisions. That is, the distortion of the landscape orchestrates fate transitions [7, 22, 23] and is driven by extrinsic signals (termed as “signal-driven”, Fig 1.B).

Figure 1.

Under the noise-driven mode, the bias of cell fate decisions largely depends on the spontaneous heterogeneity of gene expressions in the cell population [24, 25]. Consequently, the initial cellular state predominantly impacts the direction of the fate decision. Chang et al. [20] uncovered that hematopoietic stem cell (HSC) population possesses intrinsic and robust heterogeneity of Scal-1 expression. Notably, populations with discrete expression levels of Scal-1 confer different propensities for downstream lineage commitment. Considering the signal-driven mode, cell fates are tightly steered by extrinsic signals (e.g., cytokines, chemical molecules, mechanical strength and genetic operations) that reshape the landscape (Fig 1.B). In this circumstance, the impact of the initial state on fate decisions is relevantly inconsequential. Additionally, due to the accessibility of signal manipulation, the signal-driven mode has been widely utilized for cell fate engineering [26, 27], leading to in-vitro induction systems centered on induced pluripotent stem cells (iPSCs) for obtaining desired cell types [28]. Recently, researchers reported a “fate-decision abduction” of erythroid-to-myeloid trans-differentiation induced by various type of cancer, which facilitates tumor escape from the individual’s immune system [29]. Collectively, driving forces couple the foundational and crucial features of fate decisions, serving as an essential basis for further decoding fate decisions and interpreting the development of organisms [16]. By examining the two driving modes, we can gain a better understanding and characterization of cell fate decisions, including in-vivo cell differentiation, oncogenesis, and invitro reprogramming systems.

Nevertheless, the driving forces that underlie fate decisions remain largely elusive. The intricate nature of gene regulatory networks (GRNs) presents a challenge in deciphering driving modes. It has been generally acknowledged that corresponding core GRNs orchestrate cell fate decisions [30, 31], where the lineage-specifying transcription factors (TFs) interact to implement fate-decision procedures. Furthermore, researchers transferred specific TFs into donor cells to reconfigure the intracellular GRNs for acquiring cell types of interest [32, 33]. Although some studies suggested that perturbation of a single TF is sufficient to transform certain cell fates [28], large number of TFs are inevitably involved in most differentiation/reprogramming processes [34]. In particular to orchestrate decisions among multiple cell fates, it is necessary for TFs to regulate target genes cooperatively [35]. As crucial determinants of cell fates, TFs function via binding to cis-regulatory elements (CREs, e.g., promoter and enhancer). CREs of a single gene in metazoans can simultaneously accommodate numerous TFs [36, 37]. While experimental protocols have been developed to assess TF binding and one-to-one up-or down-regulatory relationship, it is more challenging to quantify these combinatorial regulations. For instance, given two factors activate and inhibit the same target gene, respectively, does the target gene turn on or off when both factors are present in its CREs?

Computational approaches in systems biology can be utilized to tackle complex networks [38–42]. A concise GRN model typically entails the following two elements. The element 1 is the topology. Much research efforts have been devoted to investigating network topologies on cellular behaviors, e.g., toggle switch [31, 43], and feed-forward loop [44–46]. In particular, the Cross-Inhibition with Self-activation (CIS) network is one of the most studied two-node GRNs in cell fate decisions [47, 48], with examples found in GATA1-PU.1 and FLI1-KLF1 in hematopoiesis [49, 50], NANOG-GATA6 and OCT4-CDX2 in gastrulation [16, 51], and Sir2 and HAP in yeast aging [23]. In this topology, two lineage-specifying factors inhibit each other while active themselves. For example, in the well-known GATA1-PU.1 circuit, GATA1 directs fate of megakaryocyte-erythroid progenitor (MEP), and PU.1 (also known as SPI1) specifies the fate of granulocyte-monocyte progenitor (GMP) [47]. Namely, the antagonism of two TFs implicates two cell fates in competition with each other.

The another element is the logic for regulatory functions [52, 53]. Exemplified by the CIS network, each node (e.g., X and Y) receives the activation by itself and inhibition by the counterpart. Hence there is naturally the logic function between these two inputs. Given the logic function is AND, in the context of biological mechanism of regulation by TFs, X gene expresses only when X itself is present in X’s CREs but Y is not. Researchers observed in the E. coli lac operon system that changes in one single base can shift the regulatory logic significantly, suggesting that logic functions of GRNs can be adapted on the demand of specific function in organisms [40, 54, 55]. Additionally, considering that the combination and cooperativity of TFs are of great significance in development [56–58], theoretical investigation of the logic underlying GRNs should be concerned in cell fate decisions. However, despite the existence of large number of mathematical models on fate decisions, the role played by the regulatory logic in cell fate decisions is still obscure. Some theoretical studies put emphasis on specific biological instances, adopting logic functions that best fit the observations derived from experiments [7, 39, 59]. As a result, the models incorporated different regulatory logic received limited attention. Other research assigned logic to large-scale multi-node GRNs, confining the interpretation of the role of logic [41, 60]. Collectively, the bridge between the logic of nodes in GRNs and cell fate decisions has not yet been elucidated systematically and adequately. Current research already encompassed a wealth of cell fate decisions: embryogenesis [61–63], lineage commitment [50, 64, 65], oncogenesis [66–68], in-vitro reprogramming [69–71], and large-scale perturbations [28, 72–74]. Analogous to the effort on taxonomy of cell types and tumors [75], how cell fate decisions can be classified and distilled to the common properties is a challenge for further exploring systematically and application on fate engineering [76, 77].

In this work, we integrated the fate-decision modes (noise-driven/signal-driven) and the classical logic operations (AND/OR) underlying GRNs in a continuous model. Based on our model, we investigated the impact of distinct logic operations on the nature of fate decisions with driving modes in consideration. Additionally, we extracted the difference in properties between the two driving modes. We unearthed that in the noise-driven stem cells, regulatory logic results in the opposite bias of fate decisions. We further distilled the relationship among noise profiles, logic motifs, and fate-decision bias, showing that knowledge of two of these allow inference of the third heuristically. Under the signal-driven mode, we identified two basic patterns of cell fate decisions: progression and accuracy. Moreover, based on our findings in silico, we characterized cell fate decisions in hematopoiesis and embryogenesis and unveiled their decision modes and logic motifs underlying GRNs. Ultimately, we applied our framework to a reprogramming system. We deciphered the driving force of this trans-differentiation, and utilized noise patterns for nominating key regulators. We underscored that clustering of gene noise patterns is an informative approach to investigate high-throughput dataset. Together, we underlined regulatory logic is of the significance in cell fate decisions. Our work presents a generalizable framework for classifying cell fate decisions and blueprint for circuit design in synthetic biology.

Result

Section 1 Mathematical model of the CIS network with logic motifs

Binary tree-like cell fate decisions are prevalent in biological systems [51, 75, 78, 79], orchestrated by a series of the CIS networks. Accordingly, we developed our ordinary differential equations (ODEs) model based on this paradigmatic and representative topology (Eq1, Eq2; see Methods).

X and Y are TFs in the CIS network. n₁ and n₂ are the coefficients of molecular cooperation. k₁-k₃in Eq1 and k₄-k₆ in Ep2 represent the relative probabilities for four configurations of binding of TFs and CREs. (Fig 2.A). d₁ and d₂ are degradation rates of X and Y, respectively. Of note, in Eq1 correspond to the transcription rate of gene X in four configurations (similarly in Eq2, Fig 2.A). Thus, the AND/OR logic configuration in GRN can be implemented by tuning the corresponding parameters above. With appropriate parameters, we are able to reproduce the Boolean-like attractor basin in the continuous models (Fig 2.B–C; see Methods). Under double AND and double OR motifs (termed as AA motif and OO motif, respectively), there are typically three stable steady states (SSSs) in the state spaces (Fig 2.C): two attractors near the axes representing the fate of lineage X (denoted as LX, X^highY^low) and the fate of lineage Y (denoted as LY, X^lowY^high), and the attractor in the center of the state space representing stem cell fate (denoted as S).

Figure 2.

Evidently, the stem cell states exhibit different expression patterns between the two logic motifs. Stem cells in the AA motif do not express X nor Y (Fig 2.B top panel; express in low level in Fig 2.C top panel), while in the OO motif, stem cells express both lineage-specifying TFs (Fig 2.B bottom panel; express in high level in Fig 2.C bottom panel). The difference in the status of S attractors relates to the co-expression level of lineage-specifying TFs in stem cells in real biological systems [11, 50]. Intuitively, from the view of Boolean model, stem cell state in the AA motif ([0,0] state) needs to switch on lineage-specifying TFs to transit to downstream fates (Fig 2.B top panel). Whereas in the OO motif, fate transitions are subject to the switch-off of TF expression (Fig 2.B bottom panel). Furthermore, we introduced the solution landscape as a pathway map consisting of all stationary points and their connections [80–82]. From the perspective of the solution landscape, two logic motifs possess akin geometric topologies in their steady-state adjacencies (Fig 2.D): when there are three fates coexisting in the state space, S attractor resides in the middle of LX and LY as the possible pivot for fate transitions (Fig 2.D). To investigate noise, we developed models with stochastic forms (see Methods). Simulations display the primary distribution of cell populations, corresponding to SSSs in deterministic models (Fig 2.E–F).

Section 2 Two logic motifs exhibit opposite bias of fate decisions under the noise-driven mode

We first investigated the difference between the AA and OO motifs under the noise-driven mode. To mimic differentiation, we assigned the stem cell state as the starting point in simulation. Asymmetry of the noise levels was introduced. First, we set the noise level of TF X higher than that of Y (σ_x=0.18, σ_y= 0.12). Under this asymmetric noise, we observed that stem cells shifted toward LX in the AA motif (Fig 3.A–B), but toward LY in the OO motif (Fig 3.C–D). From the perspective of the state space, such properties intuitively originate from the distinctive status of stem cell attractors in two logic motifs (Fig 2.B–C). In the AA motif, the stem cell state resides at the origin of coordinates. Thus, with increasing X’s noise level, the stem cell population crossed the boundary between S and LX basins with a rising probability (Fig 2.C top panel). Consequently, the fate decision of the stem cell population manifests a bias toward LX. Likewise, in the OO motif, the stem cell population has a higher probability of entering LY basin following an increase in X’s noise (Fig 2.C bottom panel). Next, we simulated multiple sets of noise levels for X and Y. We quantified distribution of cell types, which was determined by the basin in which the final state of each round of stochastic simulation ended up. We observed that stem cell population displays almost opposite differentiation preference under identical noise levels but distinct logic motifs (Fig 3.E). Conversely, when two distinct logic motifs exhibiting the same fate-decision bias, cell populations need to employ opposite noise patterns (Fig 3.E–F). Collectively, if two of the three (noise profiles, logic motif, fate-decision bias) are accessible, the last is inferential.

Figure 3.

Next, we wondered whether noise could act as a driving force for reprogramming (e.g., from LY to S). We assigned LY state as the starting cell type in simulation. Apparently, in the AA motif, transition of LY to S can be realized by increasing noise level of TF Y (FigS1.A). Meanwhile in the OO motif, it is the increased noise level of X that can drive the transition from LY to S (FigS1.B), which is also intuitive by viewing the basin geometry of the state space (Fig 2.C). These observations suggested that under the noise-driven mode, experimental reprogramming strategy need to take consideration of the regulatory logic (e.g., in reprogramming of LY to S, perturb the high expression TF of LY in the AA motif, while in the OO motif, perturb the low expression TF).

Section 3 Two logic motifs decide oppositely between differentiation and maintenance under the signal-driven mode

In addition to noise, cell fate decisions can also be driven by signals, e.g., GM-CSF in hematopoiesis [83], CHIR99021 in chemically induced reprogramming [84]. The change conducted by signals corresponds to the distortions of the cell fate landscape. To simulate the signal-driven mode, we focused on the effect of parameters in the mathematical models on the system’s dynamical properties. To simulate models feasibly and orthogonally, we added parameters u (u_xin Eq3, u_y in Eq4) to Eq1-2:

The increase of u represents an elevation in the basal expression level of lineage-specifying TFs, reflecting an induction signal from the extracellular environment. From an experimental standpoint, this signal can be the induction of small molecules or overexpression by gene manipulations, such as the transfection of cells with expression vectors containing specific genes.

We first explored the impact on the system when the two induction parameters are changed symmetrically (u=u_x=u_y). As the increase of u, the number of SSS in the AA system decreases from three to two, where S attractor evaporates after a subcritical pitchfork bifurcation (Fig 4.A, FigS2.A). Whereas in the OO motif, after the increase of u, LX and LY attractors disappear with saddle-node bifurcations respectively. Only the SSS representing stem cell fate is retained in the state space (Fig 4.B). We then portrayed all the topology of the steady-state adjacency that accompanied the increase of u, from the perspective of the solution landscape. In the AA motif, the attractor basin of LX and LY started to adjoin and occupied the vanishing S attractor basin together (Fig 4.C–D). Accordingly, the stem cells cannot maintain themselves and decided to differentiate into either one of the lineages. Moreover, if the cell population possesses the same noise levels in both X and Y, then the fate decisions are unbiased (FigS2.B). Notably, only in the AA motif we observed an intermediated stage during the increase of u, where all three fates are directly interconnected (Fig 4.C 2^nd panel and D 2^nd panel, Fig.4E). We globally screened 6,213 groups of parameter sets under the AA motif, and this temporal fully-connected stage can be observed for 82.7% of them (see Methods; Table S1), indicating little parameter sensitivity. Unlike the indirect attractor adjacency structure mediated by S attractor (Fig 2.D), the fully-connected solution landscape facilitates transitions between any two pairs of fates. Furthermore, this fully-connected stage locates between the fate-undetermined stage (Fig 4.C top panel) and fate-determined stage (Fig 4.C 3^rd panel), comparable to the initiation (or activation) stage before the lineage commitment in experimental observations [85–87]. Therefore, we suspected that the robust fully-connected stage in the AA motif may correspond to a specific priming period in cell fate decisions.

Figure 4.

From the standpoint of reprogramming of differentiated cells back into progenitors, in the AA motif, differentiated cells are more capable of maintaining their own fates during the symmetrical increase of the induction signals on both lineages (Fig 4.A and C). Whereas in the OO motif, the attractor basin of LX or LY is progressively occupied by the stem cell fate as u_xand u_y increase together (Fig 4.F–G). In this scenario, the downstream fates are eventually reversed back to the undifferentiated state (Fig 4.B and F). Namely, reprogramming engaged in the OO motif can be accomplished by bi-directional induction of downstream antagonistic fates. In sum, we found that under symmetrical signal induction, the behavior of stem cells is subject to core GRN’s logic motifs. In the AA motif, stem cells prefer to differentiate, while under the OO motif, the stem cell population inclines to maintain its undifferentiated state.

Section 4 The trade-off between progression and accuracy of cell fate decisions under the signal-driven mode

According to experimental observations, the majority of fate decisions exhibit lineage preference, also known as “symmetry breaking of fate decisions” [14, 64, 65, 89]. Take the lineage choices in hematopoiesis as an example, Some HSCs prefer myeloid over lymphoid [64, 90]. This fate-decision bias also further shifts along with aging and infection [86, 91]. In studying this preference in fate decisions, we broke the symmetry in the signal-driven models, by solely increasing u_x while keeping u_y =0 (Fig 5.A). First, it is apparent that the fate decision will significantly steer toward LX along with the increase of u_x, regardless of the logic motifs. Ultimately the state spaces contain only LX attractor when u_x is sufficiently high (Fig 5.B–C, FigS3.A). However, the changes in the state space and the solution landscape follow different routes for two logic motifs. In the AA motif, S attractor basin disappears at first, leaving a state space with two differentiated fates (Fig 5.D–E). Then the basin of LY attractor shrinks and finally disappears (Fig 5.D–E). Whereas in the OO motif, LY attractor disappears first. Then S attractor, with an enlarged basin, shares the state space with LX attractor. Finally, S attractor basin abruptly disappears by a saddle-node bifurcation (Fig 5.F–G).

Figure 5.

The distinct sequences of attractor basin disappearance as u_x increasing can be viewed as a trade-off between progression and accuracy. In the AA motif, the attractor basin of LX and LY adjoins (Fig 5.D middle panel) when S attractor disappears due to the first saddle-node bifurcation (Fig 5.B). Notwithstanding the bias of differentiation toward LX, the initial population still possesses the possibility of transiting into LY (FigS3.B). That is, in the AA motif, as the increase of induction signal u_x, the “gate” for the stem cell renewal is closed first. Stem cells are immediately compelled to make fate decisions toward either LX or LY, with a bias toward LX but a nonignorable probability of entering LY. Albeit the accuracy of differentiation is therefore compromised, the overall progression of differentiation is ensured (i.e., all stem cells have to make the fate decisions downward. This causes the pool of stem cells to be exhausted rapidly). Whereas in the OO motif, the antagonistic fate, LY, disappears first (Fig 5.C). The attractor basin of S and LX are adjacent in the state space (Fig 5.F). In this case, the orientation of the fate decisions is generally unambiguous since the stem cell population can only shift to LX, ensuring the accuracy of differentiation. Next, to check if the observed sequences of basin disappearance are artifacts of specific parameter choice, we randomly sampled parameter sets to check the sequence of attractor changes in their state spaces (6,207 groups of the AA motifs and 6,634 groups of the OO motifs; Table S1). We found that 96% AAs and 70% OOs exhibit the same sequence of attractor vanishment mentioned above (FigS3.C-D; see Methods). These results of the global screen demonstrated that the sequence of attractor vanishment is robust to parameter settings. In sum, we proposed that logic motifs couple the trade-off between progression and accuracy as a general phenomenon in the signal-driven asymmetrical fate decisions (FigS3.E).

Next, we examined the trans-differentiation from LY into LX by increasing u_x. In the AA motif, with the induction of X, LY directly transited into LX as the stem cell state disappears before LY (Fig 5.D, FigS4.A). Intriguingly, for the OO motif under the same induction, LY population first returned to the S state and then flows into LX (Fig 5.F, FigS4.B). Namely, different logic motifs conduct distinct trajectories in response to identical induction in reprogramming. The AA motif renders a one-step transition between downstream fates (FigS4.C). While in the OO motif, it is a two-step transition mediated by the stem cell state (FigS4.D). This phenomenon suggests the observation that cells may be reprogrammed to distinct cell types depending on the induction dose [84] is more realizable in the OO motif. Integrated with the foregoing symmetrical induction, we recapitulated that in the OO motif, the bi-directional induction or a unilateral induction from a counterpart (e.g., solely induced Y to realize reprogramming of LX to S) confer downstream cell fates to return to the undifferentiated state (Fig 4.F, Fig 5.F). Whereas in the AA motif, it is substantially more difficult to achieve de-differentiation. This observation may explain why some cell types are not feasible to reprogram [92].

In prior sections, we evaluated two logic motifs under the noise-and signal-driven modes in silico. With various combinations of logic motifs and driving forces, features about fate-decision behaviors were characterized by computational models. Next, we questioned whether and how fate-decision behaviors in real biological systems correspond to computational characterizations.

Section 5 The CIS network performs differently during hematopoiesis and embryogenesis

To answer that question, we first computationally evaluated the performance of different models, specifically in simulating the process of stem cells differentiating towards LX (Fig 6.A). Under four models with different combinations of driving modes and logic motifs (FigS5.A-B), we assessed the expressions and noise levels of TFs X and Y among the cell population over tim e in stochastic simulation. We observed that, under the same logic motifs, different driving modes change in the patterns of noise rather than expression levels (Fig 6.B–C, FigS5.C-D). Overall, under the noise-driven differentiation from S to LX, the level of noise exhibits a continuous and monotonic trend (FigS5.D) for both logic motifs. For different logic motifs, in the AA motif, the noise level of X (highly expressed in LX) declines (FigS5.D top panel). Whereas in the OO motif, it is the noise level of Y (low expressed in LX) displays a rising trend (FigS5.D bottom panel). Nevertheless, under the signal-driven mode, the noise level increases and then decreases, exhibiting a non-monotonic, “Pulse-like” behavior due to signal-induced bifurcation. During S to LX differentiation, comparable to the noise-driven mode, it is the noise level of TF X in the AA motif and TF Y in the OO motif display a “Pulse-like” pattern.

Figure 6.

Such disparities between logic motifs originate from the location of S attractor (Fig 6.D, Fig 2.C). Although the target cell types are the same (LX), the AA motif requests the expression of the TF X to be turned on, while the OO motif requests the TF Y to be turned off. These key fate-transition genes, namely TF X in the AA motif and TF Y in the OO motif, both exhibit a sharp increase of variation in response to saddle-node bifurcation driven by u_x induction (1^st saddle node in Fig 5.B; 2^nd saddle node in Fig 5.C). Overall, these computational results suggest that we may be able to distinguish the two driving modes according to the noise level over time series, then logic motifs can be correspondingly assigned by the expression level of the genes in the target cell types. For instance, if the noise level of the X gene exhibits a “Pulse-like” pattern and X is highly expressed in target cell types, then this cell fate decision can be assigned as the signal-driven fate decision in an AA-like motif.

To support our findings with real-world correspondence, we first focused on the differentiation of CMPs in hematopoiesis (Fig 6.E). It is acknowledged that the transcriptional regulation of Gata1-PU.1 circuit dominates this cell fate decisions, which conforms to the CIS topology (Fig 6.E). Mojtahedi et al. [83] stimulated murine multipotent hematopoietic precursor cell line EML with erythropoietin (EPO) or granulocyte macrophage colony-stimulating factor/interleukin 3 (GM-CSF/IL-3) to examine the commitment into an erythroid or a myeloid fate, respectively. Based their curated dataset, we found that the expression level of Gata1 (highly expressed in MEPs) gradually increased during EPO induction (FigS6.A), while the noise level exhibits a “Pulse-like” trend (Fig 6.F top panel). Symmetrically, during GM-CSF/IL-3 induced differentiation toward GMPs, the expression level of PU.1 (highly expressed in GMPs) gradually increased (FigS6.B), while the noise level also presents a “Pulse-like” pattern (Fig 6.F bottom panel). The trends shown in the dataset resembles the signal-driven mode with the AA motif. In addition, we quantified the expression of Gata1 and PU.1 via the single molecule FISH dataset (FigS6.C) [24]. They are at low levels in CMPs, corresponding to the expression patterns in the AA motifs (Fig 2.E, Fig 6.D). Together, we suggested that the Gata1-PU.1 circuit performs in an AA-like manner, and this differentiation system [83] is under the signal-driven mode.

Another paradigmatic model of fate decision is the differentiation of embryonic stem cells (ESC). [93] found that under the retinoic acid (RA) exposure system in vitro, mouse embryonic stem cells (mESC) differentiated into two lineages: extraembryonic endoderm (XEN)-like and ectoderm-like. The investigators recapitulated that two clusters of TFs with the CIS topology determined this lineage decision-making (Fig 6.G). We observed that the noise levels in most of these fate-decision TFs (16/22 73%) are gradually increasing during time, and 14% (3/22) of them exhibit “Pulse-like” behavior (FigS5.D), suggesting the process is more in the noise-driven mode (FigS6.E). Furthermore, we focused on Gbx2 and Tbx3, the two likely targets of RA that are crucial for this fate decision [93]. The noise levels over time of these two TFs are consistently increasing (FigS6.D). In addition, their initial expressions are at high level, in agreement with that of the OO motif (Fig 6.D and H). In short, we proposed that the mESCs differentiation system under RA exposure performs in an OO-like manner, and its differentiation is under the noisedriven mode in this experimental setting.

Section 6 The chemical-induced reprogramming of human erythroblasts (EBs) to induced megakaryocytes (iMKs) is the signal-driven fate decisions with an OO-like motif

The foregoing cell fate decisions initiate from pluripotent cells in mice (mESCs, CMPs) corresponding to a typical “downhill” in Waddington’s metaphor. In 2006, Yamanaka et al. accomplished the reprogramming from mouse fibroblasts into iPSC state via the noted “OSKM” factors, representing “uphill” in Waddington’s metaphor [92, 94, 95]. Likewise, trans-differentiations from one lineage to another have been realized by overexpression or chemical inductions [26], whether they correspond to direct “trespassing” of the ridge or an “up-and-down” through the peak in Waddington landscape are still elusive. We then applied our models to reprogramming systems, with a primary focus on hematopoiesis. Qin et al. [96] recently achieved the direct chemical reprogramming of EBs to iMKs using a four-small-molecule cocktail (Fig 7.A). Investigators presented that EBs underwent an induced bipotent precursor for erythrocytes and MKs (iPEMs) to finally desired iMKs. It is acknowledged that the FLI1-KLF1 circuit with the CIS topology dominates this fate-decision process [47, 50]. To deduce the logic motif of the FLI1-KLF1 circuit, we quantified the expression patterns of FLI1 and KLF1 based on published single-cell RNA-seq data [96]. We can observe the fate transition from the EB population (FLI1^low, KLF1^high) to the iMK population (FLI1^high, KLF1^low) (Fig 7.B). According to their expression level, the cell populations can be primarily classified into three clusters. In addition, both FLI1 and KLF1 are highly expressed in the intermediate cell population suspected to be the progenitors of iMKs and EBs [96]. Namely, the pattern of expression level is concordant with the OO motif in our framework (Fig 6.D).

Figure 7.

Next, to investigate the driving force of this reprogramming system, we simulated the fate transition from EB (corresponding to LX, blue) to iMK (LY, purple) under both driving modes. Under the noise-driven mode, we assumed that the reprogramming system facilitated the noise levels of both TFs. For simplicity, the starting cell population (EBs) was assigned symmetrical high noise levels. While under the signal-driven mode, we assumed that the four-small-molecule cocktail upregulated the expression of FLI1 (highly expressed in iMKs). Then, we simulated the transition from EBs to iMKs by lifting the basal expression level of FLI1, corresponding to parameter u_yin model (Fig 7.C). The bifurcation diagrams indicate that the signal-driven fate transition is mediated by the iPEM state (Fig 7.C). In particular, overexpression of FLI1 renders sequential saddle-node bifurcations. Thus, EBs are converted to iPEMs before steering toward the terminal iMK state (Fig7.C and D), which is consistent with the experiment’s findings.

Furthermore, we next assessed the dynamic behaviors in models of this system under different driven modes with the OO logic. As discussed before, there is no discernible difference in the expression levels between the two driving modes. Overall, the common tendency is an up-regulation in KLF1 and a down-regulation in FLI1 (FigS7.A). We then quantified the noise levels during this fate transition by the model. Under the noise-driven mode, noise levels of FLI1 and KLF1 would gradually decrease and increase, respectively, until stabilizing (FigS7.B). Under the signal-driven mode, the noise level of FLI1 would first decline and then remain nearly constant, while KLF1 would exhibit a “Pulse-like” pattern (Fig 7.E left panel). From the view of modeling, the “Pulse-like” pattern presented by KLF1 originates from the rapid shut-off during the transition from iPEM state to iMK state (2^nd saddle node in Fig 7.C).

Accordingly, we next quantified the noise levels in the real dataset over time. Impressively, the pattern emerging from the data accommodates the hypothesis of the signal-driven mode (Fig 7.E). Altogether, we proposed that this reprogramming system [96] is the signal-driven process underlying the OO-like motif. Moreover, the high expression level of FLI1 induced by small molecules is the key driving force of the fate transition, as suggested by the properties of the OO motif. We underlined that it is a classical two-step fate-decision process mediated by the upstream progenitor state, which is in agreement with the phenomenon articulated by Qin et al. [96].

We then searched for genes with similar noise patterns to KLF1 and FLI1, with a hypothesis that genes possessing comparable noise patterns, especially TFs, may synergistically perform fate-decision related functions. Thus, we applied the fuzzy c-means algorithm [97] to cluster genes based on their noise levels, rather than expression levels. In total, we observed 12 distinct clusters of temporal patterns (FigS7.C; Table S2). We focused primarily on clusters 5 and 10, where KLF1 and FLI1 are found respectively (Fig 7.F). To testify our hypothesis, we conducted an enrichment analysis using the gene set in cluster 5 and 10. Functions associated with specific cell types are significantly enriched (Fig 7.G). In particular, cluster 5 with FLI1 is largely related to platelet-related functions, whereas cluster 10 with KLF1 is largely related to the energy metabolism of blood cells. Furthermore, we filtered out 15 TFs in cluster 5 and 10 (11 in cluster 5; 4 in cluster 10). Next, by harnessing the TF interaction database (see Methods; Table S3), we collected 21 regulons associated with 10 TFs to construct TF regulatory Network (Fig 7.H). Intriguingly, in the original article, genes were classified into six “Supermodules” according to the patterns of expression levels. Genes in cluster 5 are located in Supermodules 1 and 3, representing a decreasing tendency for expression level. Meanwhile, genes in cluster 10 are distributed in Supermodule 2 and 5, and their expression levels raise from low to high (Figure 6.D from [96]). Of note, most of TFs filtered by noise patterns appear in the GRN constructed in the original article (8/10, 80%; Figure 6.G from [96]). As a result, we underscored that EGR3 and ETS1, which did not appear in the GRN of the original article, have been suggested to play important roles in the trans-differentiation. In addition, we observed that MYC and FOS possess the largest connectivity in our 10-node GRN, suggesting that these two TFs as hubs are essential regulators of this reprogramming system.

Together, in mapping our framework to the real reprogramming system, we assigned the cell fate decisions of EBs to iMKs to the signal-driven mode incorporated the OO-like motif, by comparing the noise and expression patterns measured from real dataset with pseudo-data produced by our models in a “top-down” fashion. According to the model, the reprogramming is primarily driven by induced up-regulation of FLI1. Additionally, from the view of noise patterns, we recapitulated a concise 10-node GRN and identified some TF nodes not previously recognized as major regulators, like EGR3 and ETS1.

Discussion

Comprehending the driving forces behind cell fate decisions is crucial for both fundamental scientific research and biomedical engineering. Recent advances in data collection and statistical methods have greatly enhanced our understanding of the mechanisms that regulate cell fate. However, despite the widespread use of the Waddington landscape as a metaphor in experiments, there has been little examination into whether the fate decision observed in a particular experiment corresponds to a stochastic shift from one attractor to another on the landscape (i.e., noise-driven) or an overall distortion of the landscape (i.e., signal-driven). The application of appropriate computational models in systems biology can aid in uncovering the underlying mechanisms [98, 99]. Our computational investigations have emphasized the importance of the combinatorial logic in connecting gene expression patterns with the driving forces underlying cell fate decisions. Utilizing a representative network topology known as the CIS, our analysis demonstrated how both driving forces and regulatory logic jointly shape expression patterns during fate transitions. In turn, mean and noise in gene expression patterns can reveal logic motifs and driving forces. Our analytical framework promotes the interpretability of fate decisions and can be employed to speculate on the driving factors of fate decisions using a “top-down” approach, thereby providing a reference for investigating the causality of fate decisions and experimental validation.

The roles of noise as a possible driver of fate transitions are intriguing. By our models, the relationship among noise configuration, logic motifs, and fate-decision bias has been unveiled, and we noticed the opposite fate bias for the AA and OO motifs. Conversely, on the demand of the same bias, the progenitors with different logic motifs tend to employ a different profile of noise level (Fig 3.E–F). Therefore, we suggested that under the noise-driven mode, the logic motif works like a “broker” to shape the fate preferences. Based on the assumption that the preference of fate decisions is the result of evolution and adaptation, we posited that if an organism has a functional demand for a specific bias, the noise profile will be iterated via logic motifs. In turn, changes in noise levels will be mediated by logic motifs to shape the differentiation bias. One intuitive example is the significant shift in differentiation preferences of HSCs over aging [100–102]. It has been reported that aging HSCs show different level of DNA methylation and epigenetic histone modifications from young HSCs [90, 103, 104]. How changes in epigenetic level shape the noise profile of the cell population and further affect the shift of fate-decision preference is a fascinating question. We underscored that the dissection of logic motifs underlying associated GRNs is a prerequisite for answering that question.

With the ever-expanded of single-cell sequencing data, characterizing genes by their mean expressions does not make the most of high-throughput datasets. As an intrinsic characteristic in central dogma, expression noise has been utilized to locate the “critical transitions” in complex networks [105, 106], e.g., identify the critical transitions in diseases like lymphoma [107]. Instead of differentially expressed genes, Rosales-Alvarez et al. [108] harnessed differentially noisy genes to characterize the functional heterogeneity in HSC aging. Our work presents that the patterns of noise can also be used to indicate driving forces and key regulators during the fate-decision processes. Nevertheless, compared to traditional gene expression, the interpretation of noise patterns is generally not intuitively accessible [109]. Additionally, extra a priori knowledge is needed to filter out the cluster of interest. To this end, our framework enables researchers to locate functional clusters via mathematic model based on appropriate hypothesis. Collectively, our framework provides a mechanistic explanation for noise patterns and qualitatively characterize key noise patterns to locate core regulators of fate decisions without reliance on a priori knowledge.

Comparing to tuning noise, altering signals is a more accessible approach for experimentally manipulating cell fates. When the basal expressions of two lineage-specifying genes grow symmetrically, we have shown that opposite trends of fate transitions occur under two logic motifs: In the AA motif, it promotes differentiation, whereas the OO motif stabilizes stem cell fates. This is reminiscent of the “seesaw” model where maintenance of stemness can be achieved by overexpression of antagonistic lineage-specifying genes [110]. Our model suggests that restoring stemness by inducing two antagonistic lineage-specifying genes is more likely under the OO-like motif (FigS2.C). This is in concert with our analysis that mESC differentiation system performs in an OO-like manner. In addition, Mojtahedi et al. [83] found that, under simultaneous induction of two antagonistic fates (EPO and GM-CSF/IL-3), although differentiation was delayed, it eventually occurred. This observation is consistent with models in the AA motif, and we suggested that the core regulatory circuits in hematopoiesis performs in an AA-like manner. More experimental validations would be needed to validate this hypothesis that “the seesaw model prefers the OO motif”. Conversely, insight from the “seesaw” model also provides a candidate approach to further testify logic motifs underlying GRNs in experiments. For instance, stem cells with an AA-like circuit are expected to display differentiation rather not maintenance in response to bidirectional induction.

In quantifying the signal-driven landscape changes, the solution landscape enables intuitive interpretation even in high-dimension GRNs [80, 111, 112]. In this work, we used both the state space and the solution landscape, in order to relate them for further investigations involving more than two TFs. Interestingly, from the perspective of the solution landscape, we found a robust fully-connected stage in the AA motif. We envisioned that this period corresponds to the priming stage of differentiation. Notably, this fully-connected stage was not found in the OO motif, suggesting that the necessity for priming during differentiation may be subject to the logic motifs of core GRNs.

Actual cell fate decisions are seldom purely unbiased. Under the asymmetrical signal-driven mode, we summarized the progression-accuracy trade-off in cell fate decisions: If a large number of cells are ensured to differentiate, then concessions have to be made in the accuracy of differentiation, and vice versa (FigS3.E). An intuitive example is the large-scale apoptosis occurs daily in hematopoiesis. Hence, maintaining homeostasis in vivo inevitably requests cells to respond rapidly to differentiation. A recent study reported that in response to systemic inflammation by polymicrobial sepsis, pool of CMPs is rapidly depleted to accelerate the production of downstream cell fates [113]. This is concordant with our result that Gata1-PU.1 circuit in hematopoiesis performs in an AA-like manner. On the other hand, in embryonic development like C.elegans, the accuracy of cell fate decisions is considerably emphasized. How this nature of differentiation has been adapted in evolution is an interesting question. Our work highlights that this property is associated with the logic motifs of GRNs, suggesting the emphasis on progression or accuracy may be embedded in the logic motifs of core GRNs.

We classified three examples of cell fate decisions based on patterns of noise and expression. In hematopoiesis, we took fate choice between erythroid and myeloid as a paradigm, and assigned it an AA-like motif under the signal-driven mode. In embryogenesis, we suggested the fate decision in RA exposure system is an OO-like motif under the noise-driven mode. In reprogramming, the chemical-induced trans-differentiation is the signal-driven fate decisions incorporated an OO-like motif. For simplicity and intuitiveness, we devised our model with two symmetrical combinations of regulatory logic (AA/OO). Albeit there are merely four types of cell fate decisions in consideration, our framework enables to be generalized and expanded to accommodate multi-node GRNs and complex logic combinations. Plenty of studies zoomed in one particular fate-decision events. However, from the standpoint of systems biology, we underlined that classification of fate decisions is a vital step for further investigation, as is the case for the typing of cells and tumors. Theoretically, appropriate classification of fate-decision systems enables the enrichment of common properties. So, accumulated knowledge can be inherited to new fate-decision cases. Taking reprogramming as an example, Zhao [84] recapitulated five kinds of trajectories in chemical-induced reprogramming. We suggested that the reprogramming trajectory is coupled with the logic motifs. On one hand, it is possible to answer why a certain reprogramming system exhibit a particular trajectory. On the other hand, it is possible to postulate achievable reprogramming according to the logic motifs of core GRNs (e.g., the AA motif is more likely to enable direct conversion; model 4 mentioned in [84]). Recently, synthetic biology has realized the insertion of the CIS network in mammalian cells [22]. Our work also provides a blueprint for designing logic motifs with particular functions. We are also interested in validating the conclusions drawn from our models in a synthetic biology system.

Supplemental information

Supplemental information includes seven figures and three tables and can be found with this article online.

Acknowledgements

We thank Y. Zhao for advice of model application on real biological process; D. Grün for insightful and generous feedback on noise interpretation; J. Zhang, Y. Li, and X. Pei for providing original raw-data processing R scripts; Z. Zhou for illustrations and useful discussion; and the entire Zhiyuan laboratory for support and advice.

Funding

This work was supported by the National Key Research and Development Program of China (No. 2021YFF1200500, 2021YFA0910700), and National Natural Science Foundation of China (No.12225102, 12050002, and 12226316). It is also supported by grants from Peking-Tsinghua Center for Life Sciences.

Author contributions

Conceptualization, G.X. and Z.L.; Methodology, G.X., Xiaoyi Z. and Z.Z.; Software, G.X. and Xiaoyi Z.; Formal Analysis, G.X.; Investigation, G.X.; Resources, G.X. and D.Z.; Data Curation, G.X. and W.L.; Writing – Original Draft, G.X.; Writing – Review & Editing, G.X., Xiaoyi Z., W.L., Lu Z., Z.Z., Xiaolin Z., Lei Z and Z.L.; Visualization, G.X., Xiaoyi Z. and Z.Z.; Supervision, Lei Z. and Z.L.; Project Administration, Lei Z. and Z.L.; Funding Acquisition, Lei Z. and Z.L.

Declaration of interest

The authors declare no competing interests.