Binary tree-like cell fate decisions are prevalent in biological systems (Zhou and Huang, 2011; Domcke and Shendure, 2023; Stadler et al., 2021; Macnair et al., 2019), orchestrated by a series of the CIS networks. Accordingly, we developed our ordinary differential equations (ODE) model based on this paradigmatic and representative topology (Equations 1 and 2; see ‘Materials and methods’ for details).

X and Y are TFs in the CIS network. n₁ and n₂ are the coefficients of molecular cooperation. k₁-k₃ in Equation 1 and k₄-k₆ in Equation 2 represent the relative probabilities for possible configurations of binding of TFs and CREs. (Figure 2A). 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 Equation 1, similarly in Equation 2) of each configuration, we can model the dynamics of TFs in the Shea-Ackers formalism (Shea and Ackers, 1985; Olariu and Peterson, 2019).

Figure 2

Download asset Open asset

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 the corresponding profile of transcription rates in four configurations (Figure 2A). 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; Figure 2A). Notably, instead of exploring the different logics of one certain gene (Kittisopikul and Süel, 2010), we focus on different combinations of regulatory logics due to dynamics in cell fate decisions is generally orchestrated by GRNs with multiple TFs.

Benchmarking the Boolean models with different logic motifs (Figure 2B; see ‘Materials and methods’), we reproduced the geometry of the attractor basin in the continuous models resembling those represented by corresponding Boolean models (Figure 2C; see ‘Materials and 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 (SSS) in the state spaces (Figure 2C): 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 (Figure 2B top panel; expressed in low level in Figure 2C top panel), while in the OR-OR motif, stem cells express both lineage-specifying TFs (Figure 2B bottom panel; expressed in high level in Figure 2C 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 (Loeffler and Schroeder, 2019; Palii et al., 2019). Intuitively, from the view of the 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 (Figure 2B top panel). Whereas in the OR-OR motif, fate transitions are subject to the switch-off of TF expression (Figure 2B 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 (Yin et al., 2020; Yin et al., 2021). From the perspective of the solution landscape, two logic motifs possess akin geometric topologies in their steady-state adjacencies (Figure 2D): 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 (Figure 2D). To investigate noise, we developed models with stochastic forms (see ‘Materials and methods’). Simulations display the primary distribution of cell populations, corresponding to SSSs in deterministic models (Figure 2E and F).

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 (Figure 3A and B), but toward LY in the OR-OR motif (Figure 3C and D). From the perspective of the state space, such properties intuitively originate from the distinctive status of stem cell attractors in two logic motifs (Figure 2B and 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 (Figure 2C 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 (Figure 2C bottom panel). Next, we simulated multiple sets of noise levels for X and Y. We quantified the 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 (Figure 3E). Conversely, when two distinct logic motifs exhibiting the same fate-decision bias, cell populations need to employ opposite noise patterns (Figure 3E and F). Collectively, if two of the three (noise profile, logic motif, fate-decision bias) are accessible, the last is inferential.

Figure 3 with 1 supplement see all

Download asset Open asset

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 the noise level of TF Y (Figure 3—figure supplement 1A). Meanwhile, in the OR-OR motif, it is the increased noise level of X that can drive the transition from LY to S (Figure 3—figure supplement 1B), which is also intuitive by viewing the basin geometry of the state space (Figure 2C). 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).

In addition to noise, cell fate decisions can also be driven by signals, e.g., GM-CSF in hematopoiesis (Mojtahedi et al., 2016), CHIR99021 in chemically induced reprogramming (Zhao, 2019). 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_x in Equation 3, u_y in Equation 4) to Equations 1; 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 (Figure 4A and Figure 4—figure supplement 1A). 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 (Figure 4B). 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 (Figure 4C and 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 (Figure 4—figure supplement 1B).

Figure 4 with 1 supplement see all

Download asset Open asset

Notably, in the AND-AND motif, we observed a brief intermediated stage before S attractor disappears, where all three fates are directly interconnected (Figure 4C second panel and Figure 4D second panel and Figure 4E). To manifest the generality, we globally screened 6231 groups of parameter sets under the AND-AND motif, and this logic-dependent intermediated stage can be observed for 82.7% of them (see ‘Materials and methods’; Supplementary file 1), indicating little dependence on particular parameter setting (1.8% in the OR-OR motif). Unlike the indirect attractor adjacency structure mediated by S attractor (Figure 2D), the solution landscape with fully-connected structure facilitates transitions between any two pairs of fates. Furthermore, this transitory fully-connected stage is located between the fate-undetermined stage (Figure 4C top panel) and fate-determined stage (Figure 4C third panel), comparable to the initiation (or activation) stage before the lineage commitment in experimental observations (Brand and Morrissey, 2020; Zhang et al., 2018; Arinobu et al., 2007). 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 (Figure 4A 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_x and u_y increase together (Figure 4F and G). In this scenario, the downstream fates are eventually reversed back to the undifferentiated state (Figure 4B 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.

According to experimental observations, the majority of fate decisions exhibit lineage preference, also known as ‘symmetry breaking of fate decisions’ (Stanoev and Koseska, 2022; Pei et al., 2020; Psaila et al., 2020; Chen et al., 2018). Take the lineage choices in hematopoiesis as an example, Some HSCs prefer myeloid over lymphoid (Pei et al., 2020; de Haan and Lazare, 2018). This fate-decision bias also further shifts along with aging and infection (Zhang et al., 2018; Reya et al., 2001). 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 (Figure 5A). 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 (Figure 5B and C and Figure 5—figure supplement 1A). 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 (Figure 5D and E). Then the basin of LY attractor shrinks and finally disappears (Figure 5D and 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 (Figure 5F and G).

Figure 5 with 2 supplements see all

Download asset Open asset

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 (Figure 5D middle panel) when S attractor disappears due to the first saddle-node bifurcation (Figure 5B). Notwithstanding the bias of differentiation toward LX, the initial population still possesses the possibility of transiting into LY (Figure 5—figure supplement 1B). 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 (Figure 5C). The attractor basin of S and LX are adjacent in the state space (Figure 5F). 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 (6207 groups of the AND-AND motifs and 6634 groups of the OR-OR motifs; Supplementary file 1). We found that 96% AND-AND motifs and 70% OR-OR motifs exhibit the same sequence of attractor vanishment mentioned above (Figure 5—figure supplement 1C and D and D; see ‘Materials and 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 (Figure 5—figure supplement 1E).

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 (Figure 5D and Figure 5—figure supplement 2A). Intriguingly, for the OR-OR motif under the same induction, LY population first returned to the S state and then flows into LX (Figure 5F and Figure 5—figure supplement 2B). 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 (Figure 5—figure supplement 2C). While in the OR-OR motif, it is a two-step transition mediated by the stem cell state (Figure 5—figure supplement 2D). This phenomenon suggests the observation that cells may be reprogrammed to distinct cell types depending on the induction dose (Zhao, 2019) 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 (Figure 4F and Figure 5F). 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 (Shi et al., 2017).

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 (Figure 6A). Under four models with different combinations of driving modes and logic motifs (Figure 6—figure supplement 1A and 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 (Figure 6B and C and Figure 6—figure supplement 1C and D). Overall, under the noise-driven differentiation from S to LX, the variance of expression exhibits a continuous and monotonic trend (Figure 6—figure supplement 1D) for both logic motifs. For different logic motifs, in the AND-AND motif, the expression variance of X (highly expressed in LX) declines (Figure 6—figure supplement 1D top panel). Whereas in the OR-OR motif, it is the expression variance of Y (low expressed in LX) displays a rising trend (Figure 6—figure supplement 1D 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 with 2 supplements see all

Download asset Open asset

Such disparities between logic motifs originate from the location of S attractor (Figure 6D and Figure 2C). 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 (first saddle node in Figure 5B; second saddle node in Figure 5C). 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 (Figure 6E). It is acknowledged that the transcriptional regulation of Gata1-PU.1 circuit dominates this cell fate decisions, which conforms to the CIS topology (Figure 6E). Mojtahedi et al., 2016 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 (Figure 6—figure supplement 2A), while the expression variance exhibits a nonmonotonic trend (Figure 6F top panel). Symmetrically, during GM-CSF/IL-3 induced differentiation toward GMPs, the expression level of PU.1 (highly expressed in GMPs) gradually increased (Figure 6—figure supplement 2B), while the expression variance also presents a nonmonotonic pattern (Figure 6F bottom panel). The trends shown in the dataset resemble 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 (Figure 6—figure supplement 2C; Wheat et al., 2020). They are at low levels in CMPs, corresponding to the expression patterns in the AND-AND motifs (Figure 2E and Figure 6D). Together, we suggested that the Gata1-PU.1 circuit performs in an AND-AND-like manner, and this differentiation system (Mojtahedi et al., 2016) is under the signal-driven mode.

Another paradigmatic model of fate decision is the differentiation of embryonic stem cells (ESC). Semrau et al., 2017 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 specification (Figure 6G). 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 (Figure 6—figure supplement 2D), suggesting the process is more likely driven by noise (Figure 6—figure supplement 2E). Furthermore, we focused on potential key regulators: Gbx2 and Tbx3, the two likely targets of RA that are crucial for this fate decision (Semrau et al., 2017). The expression variances over time of these two TFs are consistently increasing (Figure 6—figure supplement 2D). In addition, their initial expressions are at a high level, in agreement with that of the OR-OR motif (Figure 6D 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.

The foregoing cell fate decisions initiate from pluripotent cells in mice (mESC, CMP) 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 (Shi et al., 2017; Hamazaki et al., 2017; Abdallah and Del Vecchio, 2019). Likewise, trans-differentiations from one lineage to another have been realized by overexpression or chemical inductions (Xu et al., 2015), 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., 2022 recently achieved the direct chemical reprogramming of EBs to iMKs using a four-small-molecule cocktail (Figure 7A). Investigators presented that EBs underwent an induced bipotent precursor for erythrocytes and MKs (iPEM) to finally desired iMKs. It is acknowledged that the FLI1-KLF1 circuit with the CIS topology dominates this fate-decision process (Orkin and Zon, 2008; Palii et al., 2019). 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 (Qin et al., 2022). We can observe the fate transition from the EB population (FLI1^low, KLF1^high) to the iMK population (FLI1^high, KLF1^low) (Figure 7B). 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 (Qin et al., 2022). Namely, the pattern of expression level is concordant with the OR-OR motif in our framework (Figure 6D).

Figure 7 with 1 supplement see all

Download asset Open asset

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_y in model (Figure 7C). The bifurcation diagrams indicate that the signal-driven fate transition is mediated by the iPEM state (Figure 7C). In particular, overexpression of FLI1 renders sequential saddle-node bifurcations. Thus, EBs are converted to iPEMs before steering toward the terminal iMK state (Figure 7C 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 (Figure 7—figure supplement 1A). 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 (Figure 7—figure supplement 1B). 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 (Figure 7E 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 (second saddle node in Figure 7C).

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 (Figure 7E). Altogether, we proposed that this reprogramming system (Qin et al., 2022) 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., 2022.

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 (Kumar and Futschik, 2007) to cluster genes based on their expression variances, rather than expression levels. In total, we observed 12 distinct clusters of temporal patterns (Figure 7—figure supplement 1C; Supplementary file 2). We focused primarily on clusters 5 and 10, where KLF1 and FLI1 are found respectively (Figure 7F). To testify our hypothesis, we conducted an enrichment analysis using the gene set in clusters 5 and 10. Functions associated with specific cell types are significantly enriched (Figure 7G). 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 clusters 5 and 10 (11 in cluster 5; 4 in cluster 10). Next, by harnessing the TF interaction database (see ‘Materials and methods’; Supplementary file 3), we collected 21 regulons associated with 10 TFs to construct the TF regulatory Network (Figure 7H). 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 Supermodules 2 and 5, and their expression levels raise from low to high (Figure 6D from Qin et al., 2022). Of note, most of the TFs filtered by expression variance patterns appear in the GRN constructed in the original article (8/10, 80%; Figure 6G from Qin et al., 2022). 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 the 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.

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.

In GRN, each TFs is represented by a node and the edge between nodes represents the regulatory relationship. In our work, we considered a GRN comprised of 2 TFs (i.e., X, Y), formulated with 2 ODEs describing the change of each TF. Define ([X], [Y]) to be the concentration of TF X and Y, the ODE model is described as follows:

where H_X([X], [Y]) is the production rate of TF X that combines the effects from both activators and suppressors of the X and d_X is the decay rate from the X and Y, which integrate degradation and dilution.

The exact form of H_X([X], [Y]), H_Y([X], [Y]) are derived based on the simple GRN with nodes X and Y via a set of molecular interactions between these TFs themselves, genes that encode for them, and the mRNAs (Abdallah and Del Vecchio, 2019).

TFs in a GRN act as multimers to implement regulatory interaction. In our model, we treated TFs X and Y as acting in their homo-multimer forms, X_n1 and Y_n2, respectively (i.e., n₁ TF X monomers reversibly form an activated homo-multimerized form X_n1 and n₂ monomers of TF Y to reversibly form activated homo-multimer Y_n2). Of note, n₁ and n₂ here are able to be further generalized beyond the number of binding elements (Nam et al., 2022; Santillán, 2008). The multimerization biochemical reaction of TFs X and Y can be represented as follows,

There are many diverse mechanisms of regulation (Lambert et al., 2018; Balsalobre and Drouin, 2022). For simplicity, our model considered transcriptional regulation as a major mechanism since it is feasibly characterized in experiments and can well represent the interaction between genes in Cross-Inhibition with Self-activation (CIS) network. To activate downstream transcription, TFs bind with their Cis-Regulatory Elements (CREs). We denoted D_X for the no-bound CREs of TF X and D_Y for the TF Y in an independent manner. Here, we posited different TFs bind to exclusive, non-overlapping CREs to regulate target genes. Hence, there are totally eight binding patterns described as follows,

Next, we modeled the biochemical reactions of transcription and translation. Transcription is an elaborate process involving many steps, including initiation, elongation, and termination of an mRNA transcript. Likewise, translation includes peptide formation and elongation, followed by protein folding before function. However, to evaluate the behavior of GRN concisely, we treated transcription and translation as single-step reactions and then use lumped rates to encompass the time it takes for all steps in the elaborate machinery to complete (Abdallah and Del Vecchio, 2019). Hence, the transcription process can be represented as follows,

Here, we distinguished the basal rate of transcription of a gene without activation or repression (termed as constitutive transcription). The basal transcription rate is shown in the left panel in Equation 8.

Though we considered eight configurations, these ultimately led to two mRNA transcripts, m_X and m_Y. Translation can also be represented by a one-step process from these transcripts to their protein products. The one-step translation biochemical reactions are described as follows,

Next, we considered the decay of these proteins as well as their respective mRNA species. In general, these proteins and mRNA species will undergo decay to some extent which is a combination of both degradation and dilution.

here, the left and middle panels represent the degradation of TFs X and Y as well as mRNA transcripts m_X and m_Y. The right panel represents the dilution of TFs X and Y due to cell division. In general, the degradation rate γ can be determined using the following relations,

where t_1/2 represents the half-lives of each protein or mRNA transcript, t_doubling represents the cell division rate. For protein products X and Y, we used the notation d to illustrate the total rate of decay rather than two separate parameters in our model, i.e., d₁ = δ_x + β, d₂ = δ_y + β, respectively. In general, t_1/2 and t_doubling are robust and stationary, thus we treated these paramaters as constant in our model, i.e., d_i = Const, for i = 1, 2.

We have described the reactions that comprise the endogenous components of our GRN. In our work, fate decisions are classified into two modes, i.e., driven by the noise of intracellular gene expression or driven by extracellular signals. Under the noise-driven mode, a Gaussian white noise is added to the concentration of TFs to illustrate the random intracellular noise. Whereas if driven by signals, the expression of a gene will be affected by extracellular molecular cocktails, physical stimulation, etc. Therefore, we added the production rate of the TF’s mRNA from the ectopic DNA to generalize our model.

here u_x, u_y represent the additional mRNA species m_X and m_Y via ectopic overexpression at rates u_x and u_y, respectively. The biochemical reactions take place on a faster timescale. This enables us to derive the change in concentration of species from a biochemical reaction based on the law of mass action. For example, for multimerization reaction of X in Equation 6, once equilibrium, we can get

Also, for TF X bind with its CRE D_X, we can get

here, we define $K_1\triangleq \frac{p_1}{q_1}$ as a constant, i.e., chemical equilibrium constant. Furthermore, we define $K_1\triangleq \frac{p_1}{q_1},i=x,y,\mathrm{1,2},\cdots ,8$ for biochemical reactions in Equation 6 and Equation 7, respectively.

Then, we constructed an ODE model for the GRN. For a given species S, the principle to get a single equation is as follows:

To explain how this principle assists us to build the ODE model, we took species m_X as a running example. Summarizing all of these biochemical reaction rates involving m_X we can get the change in concentration [m_X],

Doing so for our biochemical reaction network yielded the equations of our ODE model,

The 14-dimension ODE model of the GRN contains overwhelming parameters that make it impractical to analyze for cell fate. To reduce the dimensionality of the model, we assumed the multimerization, DNA binding/unbinding and mRNA dynamic occurs sufficiently faster than protein production and decay, the temporal derivatives of the respective species can be set to 0, indicating that the species concentration reaches its quasi-steady state (Abdallah and Del Vecchio, 2019). Thus, we can get

To better demonstrate the deviation of the model, we take TF X as a running example. Since GRN is symmetric, we can deviate TF Y’s equation in the same way. Once the biochemical reaction reaches its equilibrium, using the law of mass action we can get

Noticing that $\left[\dot{m_X}\right]$ = 0 in Equation 18, the species [m_X] can be represented by CRE terms,

In our well-stirred system, for TF X’s CREs, the conservation law holds, i.e., the total CRE concentration of X, as a constant, equals the sum of CRE concentration that is bound by X, by Y, and by X and Y. By employing so and combining with Equation 19 we can get a representation of [D_X],

Substitute [D_X] in Equation 20 with Equation 21, we can get

Then finally we reach to the point that the deviation of [X]. In Equation 18 we substituted [m_X] with Equation 22, and got

here, we define ${\displaystyle {\gamma }_0^x={\alpha }_x\frac{D_{TX}{\kappa }_x}{{\eta }_x},{\gamma }_1={\alpha }_x^1\frac{D_{TX}{\kappa }_x}{{\eta }_x},{\gamma }_2={\alpha }_x^2\frac{D_{TX}{\kappa }_x}{{\eta }_x},{\gamma }_3={\alpha }_x^3\frac{D_{TX}{\kappa }_x}{{\eta }_x}}$ as relative protein generation rates, $k_1=K_1{K}_x,k_2=K_2{K}_y,k_3=K_x{K}_y{K}_1{K}_5$ as the weight of each TF binding patterns, as ectopic DNA (to simplify the notation, we may drop the tilde of ${\tilde{u}}_i,i=x,y$), then we can get

Similarly, for TF Y, we can get

Ultimately, by combining Equation 24 and Equation 25, we reduced our 14-dimension model to a 2-dimension dynamic system,

In Equation 22, parameters of the AND-AND motif (Figure 2C top panel) are: ($r_0^x$ , $r_0^y$ , r₁, r₂, r₃, r₄, r₅, r₆) = (0, 0, 0.95, 0, 0, 0.95, 0, 0), (k₁, k₂, k₃, k₄, k₅, k₆) = (1.5, 0.8, 1, 1.5, 0.8, 1) and (d₁, d₂) = (0.55, 0.55). Parameters of the OR-OR motif (Figure 2C bottom panel) are: ($r_0^x$ , $r_0^y$ , r₁, r₂, r₃, r₄, r₅, r₆) = (0.15, 0.15, 0.95, 0, 0.95, 0.95, 0, 0.95), (k₁, k₂, k₃, k₄, k₅, k₆) = (2.1, 1.7, 1.8, 2.1, 1.7, 1.8) and (d₁, d₂) = (0.55, 0.55). In both the AND-AND and OR-OR motifs, (n₁, n₂) = (2, 2). We simulated noise in the model using a Langevin equation (Foster et al., 2009):

Where $F_i\left(x\right)$ is the deterministic function of Equation 22. $W_i$ is a Wiener process which introduces additive noise with no dependence on the state x. We integrated the Wiener process over the interval which gives us $\sqrt{∆t}{\xi }_i$ and ${\xi }_i$ is a Gaussian white noise with zero mean and given variance (σ²). Thus we can compute the new state of the dynamical systems using the Runge–Kutta method.

We conducted global parameters screening of the fully connected stage and Progression-Accuracy trade-off in our models. First, to collect parameter sets with 3 SSSs, we used Latin hypercube sampling (LHS) to screen k-series parameters symmetrically (i.e., k₁ = k₄, k₂ = k₅, k₃ = k₆) ranging from 0.001 to 5 both in the AND-AND and OR-OR motifs. We ultimately collected 6,231 sets for the AND-AND motif and 6682 sets for the OR-OR motifs (Supplementary file 1).

To analyze the sequence of vanishing SSSs. We further filtered parameter sets with 2 SSSs remained as increasing ux (6207 sets for the AND-AND motif; 6634 sets for the OR-OR motif). For each SSS, we quantified the Minus of [X] and [Y] (Figure 5—figure supplement 1C, D).

Set an autonomous dynamical system (Yin et al., 2021)

where $F:R^n\to R^n\text{is a}{C}^r\left(r⩾2\right)$ function, and a point ${\displaystyle x^{\ast }\in R^n}$ is called a stationary point (or equilibrium solution) of Equation 28 if ${\displaystyle F\left(x^{\ast }\right)=0}$. Let $J\left(x\right)=\nabla F\left(x\right)$ denote the Jacobian of $F\left(x\right)$ . For a stationary point ${\displaystyle x^{\ast }}$ , taking ${\displaystyle x=x^{\ast }+y}$ in Equation 28, we have

where $\parallel *\parallel$ denotes the norm induced by the inner product. The associated linear system

is used to determine the stability of ${\displaystyle x^{\ast }}$ . Depending on the eigenvalues of $J\left(x\right)$ with positive, negative, and zero real parts, we can define unstable, stable, and center manifolds of the Jacobian $J\left(x\right)$ spanned by the corresponding eigenvectors, as ${\displaystyle W^u\left(x^{\ast }\right),W^s\left(x^{\ast }\right)\text{and}{W}^c\left(x^{\ast }\right)}$ .

From the primary decomposition theorem, ${\displaystyle {\mathbb{R}}^n}$ can be decomposed as a direct sum:

A hyperbolic stationary point is called a saddle if ${\displaystyle W^u\left(x^{\ast }\right)}$ and ${\displaystyle W^s\left(x^{\ast }\right)}$ are nontrivial. The hyperbolic stationary point ${\displaystyle x^{\ast }}$ is called a sink (source) if all the eigenvalues of ${\displaystyle J\left(x^{\ast }\right)}$ have negative (positive) real parts. The index of a stationary point ${\displaystyle x^{\ast }}$ is defined as the dimension of the unstable subspace ${\displaystyle W^u\left(x^{\ast }\right)}$ .

The HiOSD method (Yin et al., 2019) is designed for finding index-$k$ saddles of an energy function $E\left(x\right)$ . For gradient systems, $F\left(x\right)=-\nabla E\left(x\right)$ , and the Jacobian $J\left(x\right)=-{\nabla }^{2}E\left(x\right)=-G\left(x\right)$ , where $G\left(x\right)$ denotes the Hessian of $E\left(x\right)$ . The $k$-saddle ${\displaystyle x^{\ast }}$ is a local maximum on the linear manifold ${\displaystyle x^{\ast }+W^u\left(x^{\ast }\right)}$ and a local minimum on ${\displaystyle x^{\ast }+W^s\left(x^{\ast }\right)}$ . So, the high-index saddle dynamics (HiSD) for a index k saddle (k-saddle) is a transformed gradient flow

where $P_V$ denotes the orthogonal projection operator on a finite-dimensional subspace $V$. Here, $-P_{W^u\left(x\right)}F\left(x\right)$ is taken as an ascent direction on the subspace $W^u\left(x\right)$ and $F\left(x\right)-P_{W^u\left(x\right)}F\left(x\right)$ is a descent direction on the subspace $W^s\left(x\right)$ . The subspace $W^u\left(x\right)=span{\{v}_1,\dots ,v_k$ } where $v_i$ is the unit eigenvector corresponding to the smallest $i$-th eigenvalues, which can be obtained by many methods such as minimize the Rayleigh quotients. According to the above the HiSD for a k-saddle (k-HiSD) is:

which coupled with the initial condition

Where $I$ is the identity operator and β, γ>0 are relaxation parameters.

Similarly, the GHiSD (Yin et al., 2021) for a k-saddle (k-GHiSD) of the dynamical system (28) has the following form: