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”, Fig1.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”, Fig1.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 (Fig1.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 in- vitro 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 for details).

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 possible configurations of binding of TFs and CREs. (Fig2.A). d₁ and d₂ are degradation rates of X and Y, respectively. Here, we considered a total of four CRE’s configurations as shown in Figure 2A (i.e., TFs bind to the corresponding CREs or not, 2²=4). Accordingly, depending on the transcription rates (i.e., r₀^x, r₁, r₂, r₃ in Eq1, similarly in Eq2) of each configuration, we can model the dynamics of TFs in the Shea-Ackers formalism [80, 81].

Figure 2.

Thus, the distinct logic operations (AND/OR) of two inputs (e.g., activation by X itself and inhibition by Y) can be further implemented by assigning corresponding profile of transcription rates in four configurations (Fig2.A). From the perspective of molecular biology, the regulatory logics embody the complicated nature of TF regulation that TFs function in a context-dependent manner. Considering the CIS network, when X and Y bind respective CREs concurrently, whether the expression of target gene is turned on or off depends on the different regulatory logics (specifically, off in the AND logic and on in the OR logic; Fig2.A). Notably, instead of exploring the different logics of one certain gene [44], we focus on different combinations of regulatory logics due to dynamics in cell fate decisions is generally orchestrated by GRN with multiple TFs.

Benchmarking the Boolean models with different logic motifs (Fig2.B; see Methods), we reproduced the geometry of the attractor basin in the continuous models resembling those represented by corresponding Boolean models (Fig2.C; see Methods). Under double AND and double OR motifs (termed as AND-AND motif and OR-OR motif, respectively), there are typically three stable steady states (SSSs) in the state spaces (Fig2.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).

Evidently, the stem cell states exhibit different expression patterns between the two logic motifs. Stem cells in the AND-AND motif do not express X nor Y (Fig2.B top panel; express in low level in Fig2.C top panel), while in the OR-OR motif, stem cells express both lineage- specifying TFs (Fig2.B bottom panel; express in high level in Fig2.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 AND-AND motif ([0,0] state) needs to switch on lineage-specifying TFs to transit to downstream fates (Fig2.B top panel). Whereas in the OR-OR motif, fate transitions are subject to the switch-off of TF expression (Fig2.B bottom panel). Furthermore, we introduced the solution landscape method. Solution landscape is a pathway map consisting of all stationary points and their connections, which can describe different cell states and transfer paths of them [82–84]. From the perspective of the solution landscape, two logic motifs possess akin geometric topologies in their steady-state adjacencies (Fig2.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 (Fig2.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 (Fig2.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 AND-AND and OR-OR motifs under the noise- driven mode. Here, we assigned the stem cell state as the starting point in simulation. In biological systems, it is unlikely that the noise level of different genes is kept perfectly the same. Asymmetry of the noise levels was thus 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 AND-AND motif (Fig3.A-B), but toward LY in the OR-OR motif (Fig3.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 (Fig2.B-C). In the AND-AND 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 (Fig2.C top panel). Consequently, the fate decision of the stem cell population manifests a bias toward LX. Likewise, in the OR-OR motif, the stem cell population has a higher probability of entering LY basin following an increase in X’s noise (Fig2.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 (Fig3.E). Conversely, when two distinct logic motifs exhibiting the same fate-decision bias, cell populations need to employ opposite noise patterns (Fig3.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 AND- AND motif, transition of LY to S can be realized by increasing noise level of TF Y (FigS1.A). Meanwhile in the OR-OR 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 (Fig2.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 AND-AND motif, while in the OR-OR 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 [85], CHIR99021 in chemically induced reprogramming [86]. 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 AND-AND system decreases from three to two, where S attractor evaporates after a subcritical pitchfork bifurcation (Fig4.A, FigS2.A). Whereas in the OR-OR 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 (Fig4.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 AND-AND motif, the attractor basin of LX and LY started to adjoin and occupied the vanishing S attractor basin together (Fig4.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).

Figure 4.

Notably, in the AND-AND motif we observed a brief intermediated stage before S attractor disappears, where all three fates are directly interconnected (Fig4.C 2^nd panel and D 2^nd panel, Fig.4E). To manifest the generality, we globally screened 6,213 groups of parameter sets under the AND-AND motif, and this logic-dependent intermediated stage can be observed for 82.7% of them (see Methods; Table S1), indicating little dependence on particular parameter setting (1.8% in the OR-OR motif). Unlike the indirect attractor adjacency structure mediated by S attractor (Fig2.D), the solution landscape with fully-connected structure facilitates transitions between any two pairs of fates. Furthermore, this transitory fully-connected stage locates between the fate- undetermined stage (Fig4.C top panel) and fate-determined stage (Fig4.C 3^rd panel), comparable to the initiation (or activation) stage before the lineage commitment in experimental observations [87–89]. Therefore, we suspected that the robust fully-connected stage in the AND-AND motif may correspond to a specific period in cell fate decisions.

From the standpoint of reprogramming of differentiated cells back into progenitors, in the AND-AND motif, differentiated cells are more capable of maintaining their own fates during the symmetrical increase of the induction signals on both lineages (Fig4.A and C). Whereas in the OR-OR motif, the attractor basin of LX or LY is progressively occupied by the stem cell fate as u_xand u_y increase together (Fig4.F-G). In this scenario, the downstream fates are eventually reversed back to the undifferentiated state (Fig4.B and F). Namely, reprogramming engaged in the OR-OR 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 AND-AND motif, stem cells prefer to differentiate, while under the OR-OR 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, 91]. Take the lineage choices in hematopoiesis as an example, Some HSCs prefer myeloid over lymphoid [64, 92]. This fate- decision bias also further shifts along with aging and infection [88, 93]. 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 (Fig5.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 (Fig5.B-C, FigS3.A). However, the changes in the state space and the solution landscape follow different routes for two logic motifs. In the AND-AND motif, S attractor basin disappears at first, leaving a state space with two differentiated fates (Fig5.D-E). Then the basin of LY attractor shrinks and finally disappears (Fig5.D-E). Whereas in the OR-OR 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 (Fig5.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 AND-AND motif, the attractor basin of LX and LY adjoins (Fig5.D middle panel) when S attractor disappears due to the first saddle-node bifurcation (Fig5.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 AND-AND 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 OR-OR motif, the antagonistic fate, LY, disappears first (Fig5.C). The attractor basin of S and LX are adjacent in the state space (Fig5.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 AND-AND motifs and 6,634 groups of the OR- OR motifs; Table S1). We found that 96% AND-AND motifs and 70% OR-OR motifs 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 AND- AND motif, with the induction of X, LY directly transited into LX as the stem cell state disappears before LY (Fig5.D, FigS4.A). Intriguingly, for the OR-OR motif under the same induction, LY population first returned to the S state and then flows into LX (Fig5.F, FigS4.B). Namely, different logic motifs conduct distinct trajectories in response to identical induction in reprogramming. The AND-AND motif renders a one-step transition between downstream fates (FigS4.C). While in the OR-OR 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 [86] is more realizable in the OR-OR motif. Integrated with the foregoing symmetrical induction, we recapitulated that in the OR-OR 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 (Fig4.F, Fig5.F). Whereas in the AND-AND motif, it is substantially more difficult to achieve de- differentiation. This observation may explain why some cell types are not feasible to reprogram [94].

Section 5: The CIS network performs differently during hematopoiesis and embryogenesis

In prior sections, we systematically investigated 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 observations in computation can be mapped into real biological systems. And how to discern different logic motifs and driving modes is a prerequisite for answering this question.

To end this, we first evaluated the performance of different models, specifically in simulating the process of stem cells differentiating towards LX (Fig6.A). Under four models with different combinations of driving modes and logic motifs (FigS5.A-B), we assessed the expression level and expression variance (defined as the coefficient of variation) of TFs X and Y among the cell population over time in stochastic simulation. We observed that, under the same logic motifs, different driving modes change in the patterns of expression variance rather than expression levels (Fig6.B-C, FigS5.C-D). Overall, under the noise-driven differentiation from S to LX, the variance of expression exhibits a continuous and monotonic trend (FigS5.D) for both logic motifs. For different logic motifs, in the AND-AND motif, the expression variance of X (highly expressed in LX) declines (FigS5.D top panel). Whereas in the OR-OR motif, it is the expression variance of Y (low expressed in LX) displays a rising trend (FigS5.D bottom panel). Nevertheless, under the signal-driven mode, the expression variance increases and then decreases, exhibiting a non- monotonic transition due to signal-induced bifurcation. During S to LX differentiation, comparable to the noise-driven mode, it is the expression variance of TF X in the AND-AND motif and TF Y in the OR-OR motif display a nonmonotonic pattern.

Figure 6.

Such disparities between logic motifs originate from the location of S attractor (Fig6.D, Fig2.C). Although the target cell types are the same (LX), the AND-AND motif requests the expression of the TF X to be turned on, while the OR-OR motif requests the TF Y to be turned off. These key fate-transition genes, namely TF X in the AND-AND motif and TF Y in the OR-OR motif, both exhibit a sharp increase of variation in response to saddle-node bifurcation driven by u_x induction (1^st saddle node in Fig5.B; 2^nd saddle node in Fig5.C). Overall, these computational results suggest that we may be able to distinguish the two driving modes according to the expression variance 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 expression variance of the X gene exhibits a nonmonotonic 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 AND-AND-like motif.

To support our findings with real-world correspondence, we first focused on the differentiation of CMPs in hematopoiesis (Fig6.E). It is acknowledged that the transcriptional regulation of Gata1-PU.1 circuit dominates this cell fate decisions, which conforms to the CIS topology (Fig6.E). Mojtahedi et al. [85] 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 expression variance exhibits a nonmonotonic trend (Fig6.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 expression variance also presents a nonmonotonic pattern (Fig6.F bottom panel). The trends shown in the dataset resembles the signal-driven mode with the AND-AND 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 AND-AND motifs (Fig2.E, Fig6.D). Together, we suggested that the Gata1-PU.1 circuit performs in an AND-AND-like manner, and this differentiation system [85] is under the signal-driven mode.

Another paradigmatic model of fate decision is the differentiation of embryonic stem cells (ESC). Semrau et al. [95] found that under the retinoic acid (RA) exposure system in vitro, mouse embryonic stem cells (mESCs) 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 specification (Fig6.G). We observed that the expression variance in most of these fate-decision TFs (16/22 73%) are gradually increasing during time, and 14% (3/22) of them exhibit nonmonotonic behavior (FigS5.D), suggesting the process is more likely driven by noise (FigS6.E). Furthermore, we focused on potential key regulators: Gbx2 and Tbx3, the two likely targets of RA that are crucial for this fate decision [95]. The expression variances 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 OR-OR motif (Fig6.D and H). In short, we proposed that the mESCs differentiation system under RA exposure performs in an OR-OR- like manner, and its differentiation is under the noise-driven 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 OR-OR-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 [94, 96, 97]. 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. [98] recently achieved the direct chemical reprogramming of EBs to iMKs using a four-small-molecule cocktail (Fig7.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 [98]. We can observe the fate transition from the EB population (FLI1^low, KLF1^high) to the iMK population (FLI1^high, KLF1^low) (Fig7.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 [98]. Namely, the pattern of expression level is concordant with the OR-OR motif in our framework (Fig6.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 (Fig7.C). The bifurcation diagrams indicate that the signal-driven fate transition is mediated by the iPEM state (Fig7.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 OR-OR 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 expression variances during this fate transition by the model. Under the noise-driven mode, expression variances of FLI1 and KLF1 would gradually decrease and increase, respectively, until stabilizing (FigS7.B). Under the signal-driven mode, the expression variance of FLI1 would first decline and then remain nearly constant, while KLF1 would exhibit a nonmonotonic pattern (Fig7.E left panel). From the view of modeling, the nonmonotonic pattern presented by KLF1 originates from the rapid shut-off during the transition from iPEM state to iMK state (2^nd saddle node in Fig7.C).

Accordingly, we next quantified the expression variances in the real dataset over time. Impressively, the pattern emerging from the data accommodates the hypothesis of the signal- driven mode (Fig7.E). Altogether, we proposed that this reprogramming system [98] is the signal- driven process underlying the OR-OR-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 OR-OR 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. [98].

We then searched for genes with similar patterns of expression variance to those of KLF1 and FLI1, with a hypothesis that genes possessing comparable expression variance patterns, especially TFs, may synergistically perform fate-decision related functions. Thus, we applied the fuzzy c- means algorithm [99] to cluster genes based on their expression variances, 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 (Fig7.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 (Fig7.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 (Fig7.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 [98]). Of note, most of TFs filtered by expression variance patterns appear in the GRN constructed in the original article (8/10, 80%; Figure 6.G from [98]). 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 OR-OR-like motif, by comparing the expression and expression variance 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 expression variance, 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 [100, 101]. One of the most representative work is that Huang et al. [59] modeled the bifurcation in hematopoiesis to reveal the lineage commitment quantitatively. Compared to simply modularizing activation or inhibition effect by employing Hill function in previous work, our models reconsidered the multiple regulations from the level of TF-CRE binding.

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 variance 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 AND-AND and OR-OR motifs. Conversely, on the demand of the same bias, the progenitors with different logic motifs tend to employ a different profile of noise level (Fig3.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 [102–104]. It has been reported that aging HSCs show different level of DNA methylation and epigenetic histone modifications from young HSCs [92, 105, 106]. 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 variance has been utilized to locate the “critical transitions” in complex networks [107, 108], e.g., identify the critical transitions in diseases like lymphoma [109]. Instead of differentially expressed genes, Rosales-Alvarez et al. [110] harnessed differentially noisy genes to characterize the functional heterogeneity in HSC aging. Our work presents that the patterns of expression variance 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 expression variance patterns is generally not intuitively accessible [111]. 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. Notably, if the genes that constituting the CIS network are not specified, we can conversely leverage the patterns of temporal expression variance to nominate key regulators in a model- guided manner. Collectively, our framework provides a mechanistic explanation for expression variance patterns and qualitatively characterize key expression variance 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 AND-AND motif, it promotes differentiation, whereas the OR-OR 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 [112]. Our model suggests that restoring stemness by inducing two antagonistic lineage-specifying genes is more likely under the OR-OR-like motif (FigS2.C). This is in concert with our analysis that mESC differentiation system performs in an OR-OR-like manner. In addition, Mojtahedi et al. [85] 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 AND-AND motif, and we suggested that the core regulatory circuits in hematopoiesis performs in an AND-AND-like manner. More experimental validations would be needed to validate this hypothesis that “the seesaw model prefers the OR-OR 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 AND-AND-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 [82, 113, 114]. 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 AND-AND motif. We envisioned that this period corresponds to the priming stage of differentiation. Notably, this fully-connected stage was not found in the OR-OR 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 [115]. This is concordant with our result that Gata1-PU.1 circuit in hematopoiesis performs in an AND-AND-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 expression and expression variance. In hematopoiesis, we took fate choice between erythroid and myeloid as a paradigm, and assigned it an AND-AND-like motif under the signal-driven mode. In embryogenesis, we suggested the fate decision in RA exposure system is an OR-OR-like motif under the noise-driven mode. In reprogramming, the chemical-induced trans-differentiation is the signal-driven fate decisions incorporated an OR-OR-like motif. For simplicity and intuitiveness, we devised our model with two symmetrical combinations of regulatory logic (AND-AND/OR- OR). 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 [86] 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 AND-AND motif is more likely to enable direct conversion; model 4 mentioned in [86]). Recently, synthetic biology has realized the insertion of the CIS network in mammalian cells [22]. One of the prerequisites for recapitulating the complex dynamics of fate transitions in synthetic biology is systematical understanding of the role of GRNs and driving forces in differentiation. And the logic motifs are the essential and indispensable elements in GRNs. 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.

limitation of this study

Although our framework enables the investigation of more logic motifs, we chose two classical and symmetrical logic combinations for our analysis. Future work should involve more logic gates like XOR and explore asymmetrical logic motifs like AND-OR. The gene expression datasets analyzed here are only available for a limited number of time points. Though they meet the need for discerning trends, it is evident that the application to the datasets with more time points will yield clearer and less ambiguous changing trends to support the conclusions of this paper more generally. Notwithstanding the fact that the CIS network is prevalent in fate-decision programs, there are other topologies of networks that serve important roles in the cell-state transitions, like feed-forward loop, etc. The framework should further incorporate diverse network motifs in the future. In addition, for simplicity and intuition, we here considered signals as uncoupled and additive effects in ODE models, due to feasible mapping in real biological systems, such as ectopic overexpression.

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.