Dissecting Embryonic Stem Cell Self-Renewal and Differentiation Commitment from Quantitative Models

Quantitative model of ESC self-renewal and differentiation

The principles that guided us in constructing the model were based on an assumption about the establishment of mRNA abundance. That is, the expression profile or pattern of a gene is fully determined by transcription and degradation rates (Pérez-Ortín et al., 2007). Our model had four main components that determined the transcription and degradation rates: the first was the autoregulations and cross-regulations existing in ESC master regulators and lineage-specific pivotal genes, respectively, the second was the extrinsic cues, the third component was the switch effect of chromatin modifications, and the final component was the effect of antagonism between ESC master regulators and lineage-specific pivotal genes (Fig. 1). Thus, the key was to quantify the effects of these components on transcription and degradation rates and then integrate them into the differential equations of gene expression involved in cell fate choices.

OPEN IN VIEWER

FIG. 1.

Overview of quantitative model and approach. The quantitative model of ESC self-renewal and differentiation constructed by integrating the multifactors of the multilayered nature of regulatory mechanisms and -omics data sets. (a) Autoregulations and cross-regulations among ESC master regulators and extrinsic cues. (b) Chromatin switch effect on the gene expression. The epigenetic states of a gene are assumed to be activated, bivalent, and silenced. Moreover, the epigenetic state alteration or interconverting during cell-type conversions follows an ultrasensitivity and bistability law in response to chromatin regulators. (c) Antagonism between ESC master regulators and lineage genes. Antagonism is exerted by both mRNA degradation and epigenetic silence-regulating transcription abundance. In gray background, the dissociation constant of each gene is first computed from multifactors of the multilayered nature of regulatory mechanisms and multiple data sets. Note that the dissociation constants are time-varying in cell-type conversions to adapt to the new cell type and the variations follow ultrasensitivity and bistability law. By virtue of the dissociation constants, transcription rates are calculated from fractional occupancy of transcriptional factors binding target promoters. Degradation rates are obtained from the lineage antagonism and cross-regulations of transcriptional factors. Finally, the differential equations of gene expression involved in cell fate choices are established from the transcription and degradation rates according to kinetic strategies. For a detailed description, see main text and Supplementary Data (Supplementary Data are available online at www.liebertpub.com/dna). ESC, embryonic stem cell.

Autoregulations of ESC master regulators promote sustained ESC self-renewal (Niwa et al., 2000; Boyer et al., 2005; Chickarmane et al., 2006; Loh et al., 2006; Chickarmane and Peterson, 2008; Jaenisch and Young, 2008; MacArthur et al., 2009). The positive autoregulations contribute to transcriptional rates and are implemented by regulators binding to the promoters of their own genes and forming interconnected autoregulatory loops to regulate their own expression. In addition to autoregulations, the cross-regulations of ESC master regulators are also important for ESC self-renewal. For example, Oct4 and Sox2 form a heterodimer that positively regulates the expression of Pou5f1 (which encodes Oct4), Sox2, and Nanog (Fig. 1a). The experiment evidences suggest that the heterodimer of Oct4 and Sox2 regulates Nanog by epigenetic state alterations and by degradations as Oct4 overexpression, whereas the regulatory effect of Nanog on Oct4 is likely to be indirect (Chickarmane et al., 2006; Chickarmane and Peterson, 2008; Jaenisch and Young, 2008; MacArthur et al., 2009).

Since chromatin configuration determines the accessibility of promoters to transcription and the performance of the transcription machinery (Bibikova et al., 2008; Jaenisch and Young, 2008; Lu et al., 2009), we assumed that epigenetic regulations contributed to the dissociation constant alterations of transcriptional factors binding the promoters of target genes, which switch genes on or off. To model the epigenetic regulation effects, the epigenetic states of genes involved in cell fate choices were assumed to be activated, bivalent, and silenced. Moreover, the epigenetic state interconverting or dissociation constant variations during cell-type conversions follow an ultrasensitivity and bistability law in response to chromatin regulators (Dodd et al., 2007; Sneppen et al., 2008) (Fig. 1b).

Leukemia inhibitory factor (LIF) and bone morphogenetic protein (BMP) are the best-known factors to maintain the pluripotency of ESCs by upregulating Nanog and downregulating lineage-specific genes (Ying et al., 2003; Chen et al., 2008) (Fig. 1a). Thus, the extrinsic cues increase the transcriptional rate of Nanog, whereas they repress the lineage genes by increasing the degradation rates.

The experimental studies have indirectly revealed that upregulation of lineage-specific pivotal genes represses ESC master regulators (Chickarmane and Peterson, 2008; Mikkelsen et al., 2008; Enver et al., 2009). The negative effects of key regulators of a lineage on other lineages are referred to as antagonism (Fig. 1c). For instance, there is antagonism between Oct4 and Cdx2, and an excess of Cdx2 gives rise to the trophectoderm lineage, whereas an excess of Oct4 results in the stem cell state (Niwa et al., 2005). The endoderm lineage also arises from the balance between Nanog and Gata6 through mutual antagonism; an excess of Gata6 pushes the cell into the endoderm lineage, whereas an excess of Nanog results in the stem cell state. ShRNA loss-of-function techniques to downregulate both Nanog and Oct4 approaches suggest that most lineage-specific pivotal genes are silenced by ESC master regulators (Walker et al., 2007). There is evidence showing that lineage genes that are silenced or bivalent in ESCs, but are eventually expressed in differentiated cells, are targeted by various core PcG proteins and are marked by a modified chromatin structure in ESCs (Boyer et al., 2006; Jaenisch and Young, 2008). Combining the results, we suggested that the antagonism between lineage genes and ESC master regulators was exerted by both degradation and epigenetic silence-regulating antagonized genes.

Similarly, the transcription and degradation rates of lineage-specific pivotal genes were also established from the autoregulations, cross-regulations, epigenetic regulations, and especially the antagonism among lineages.

In addition to the known regulation relationship above, we also added tentative entries in the quantitative model to test and detect some indistinct regulatory relationship of two network components.

Integrating the contributions of multiple factors above into transcription and degradation rates, the differential equations of gene expression involved in cell fate choices were established according to the kinetic strategies (Pérez-Ortín et al., 2007). More details of the model are described in the Methods section.

Model validation

A critical test for our quantitative model was whether it could simulate cell-type maintenance and conversions and, moreover, model the expression patterns of Murine ESC fate choices. Since our model could describe ESCs, lineage-committed cells, and processes of ESCs differentiating into lineage-committed cells, we needed to assess them, respectively.

We first tested the ability of ESC self-renewal and lineage-committed cell maintenance for various perturbations. In test model, the regulators, Nanog and Oct4, were assumed to be normally expressed, downregulated, and overexpressed. From model outputs, we observed that ESC master regulators maintained their expression patterns in culture, whereas LIF withdrawal induced strong downregulation of ESC master regulators in normal Nanog expression cells (Ying et al., 2003; Chen et al., 2008). In contrast, in cells of Nanog overexpression, ESC master regulators kept strong upregulation and lineage genes were effectively repressed (Fig. 2a). This is consistent with the observation that Nanog overexpression can bypass LIF (Chambers et al., 2007; Kalmar et al., 2009; Silva et al., 2009). In addition, model outputs showed that both downregulation and overexpression of Oct4 induced ESC differentiation (Niwa et al., 2000; Ivanova et al., 2006; Loh et al., 2006; Matoba et al., 2006; Lu et al., 2009) (Fig. 2b). The lineage-committed cells held a stable expression pattern under various perturbations. These results confirmed that our model could effectively simulate ESC self-renewal and lineage-committed cell maintenance.

OPEN IN VIEWER

FIG. 2.

Model validation. Model validation tests of whether it can simulate cell-type maintenance and conversions and, moreover, model the expression patterns of Murine ESC fate choices that were not used as input when fitting the parameters. (a) Differential equation outcomes of Oct4 expression levels at given Nanog expression and extrinsic cues. LIF, indicated by color, from gray (Oct4 low expression) to black (Oct4 high expression). In abscissa, 0 indicated normal Nanog expression and, Nanog expression levels are scaled as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \log _2} \left( {{{{m_N}} \mathord{ \left/ { \vphantom {{{m_N}} {{m_{N , 0}}}}} \right. \kern\nulldelimiterspace} {{m_{N , 0}}}}} \right)$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_{N , 0}}$$ \end{document} is normal expression abundance. (b) Stationary probability distribution of cells cumulated at endoderm lineage-specific gene Gata6 and trophectoderm lineage-specific gene Cdx2 from 10,000 cells (independent experiments). The initial conditions of cells are selected as ESCs in which Oct4 levels are randomly selected in equal probability from downregulation, normal expression, and overexpression, while other ESC master regulators are normally expressed and lineage genes are depressed. Stationary probability distribution of cells is calculated by averaging the differential equation outcomes at steady state (almost all genes expression levels approach to constants). Black hot spots indicate preferred genetic configurations at which cells will accumulate, and gray indicates low probability configurations; (c) and (d) are, respectively, the differential equation outcomes and measurements of ESC master regulators and lineage-specific genes for ESCs differentiating into embryoid bodies. Dashed lines are the differential equation outcomes. After 14D, the measurements are the extended data (see Supplementary Data). LIF, leukemia inhibitory factor.

As a further test of ESC differentiation, we compared the expression patterns computed from our model with the experiment measurements. We used two such data sets of ESCs differentiating into embryoid bodies (Boyer et al., 2006) and G9a^−/− and Dnmt3a/3b^−/− ESCs induced by retinoic acid (RA) (Feldman et al., 2006; Epsztejn-Litman et al., 2008) to test ESC differentiation. Our model, using the parameters trained on the original model, predicted expression of ESCs differentiating into embryoid bodies with mostly good accuracy (Fig. 2c, d and Supplementary Figs. S1–S3; Supplementary Data are available online at www.liebertpub.com/dna). In addition, we adjusted the dissociation constants in our model to simulate the epigenetic state alterations in experiment studies of ESCs, G9a^−/− and Dnmt3a/3b^−/− , treated by RA (Feldman et al., 2006; Epsztejn-Litman et al., 2008). The model output showed a consistent result with the experiment studies (Feldman et al., 2006; Epsztejn-Litman et al., 2008) and that a subpopulation of the differentiated cells restored into ESC expression profile after RA withdrawal.

Taken together, our validation tests provided the evidence that our model could capture some core principles governing pattern formation in the ESC fate choice network.

Attractor and basin analysis of quantitative regulation model

Having established the quantitative model, we sought novel insights into the cell fate choices from our model, which might be hard to find by experiment studies. The notion that cell types might correspond more generally to attractors of high-dimensional regulatory networks or differential equations was first proposed by Kauffman (1969) and has been examined extensively in the theoretical published literature since the late 1960s. Our model had two kinds of attractors: ESC self-renewal attractors and lineage-committed cell attractors.

Earlier works have shown that ESC self-renewal may be attributed to multiple factors, including extrinsic cues, ESC master regulator expression, and inhibition of lineage genes (Feldman et al., 2006; Bibikova et al., 2008; Chen et al., 2008; Karantzali et al., 2008). The results of our modeling supported these general notions and, moreover, provided novel insights into ESC self-renewal.

To investigate multiple factor effects on ESC self-renewal attractors, we designed combined perturbations to the model. First, from model outputs (Supplementary Table S1), we observed that the regulation effects of regulators, Oct4 and Nanog, on ESC self-renewal were consistent with the experiment studies (Niwa et al., 2000; Loh et al., 2006; Matoba et al., 2006; Chambers et al., 2007; Kalmar et al., 2009; Silva et al., 2009). In addition, we observed that in the absence of LIF, the model had attractors of differentiated cells, whereas in the presence of LIF, our model had two attractors, that is, a prominent peak in the region of high Nanog (HN) expression together with downregulation of lineage genes and a smaller peak of low Nanog (LN) expression accompanied with promiscuous expression of multiple lineage genes, for example, Gata4 and Gata6 (Fig. 3a, c and Supplementary Figs. S4–S7). Our results reinforced the previous studies showing that LIF sustains ESC self-renewal (Ying et al., 2003; Chen et al., 2008), and both ESC master regulators and lineage genes are expressed in a heterogeneous manner in cultures, for example, the two peaks of Nanog expression (Hayashi et al., 2008; Kalmar et al., 2009). Moreover, LN ESCs confirmed that Nanog is dispensable for ESC self-renewal (Chambers et al., 2007; Kalmar et al., 2009; Silva et al., 2009), and Oct4 being stably upregulated in both LN and Nanog-null ESCs also confirmed that Nanog has a small or indirect effect on Oct4 (dashed lines in Fig. 3a, b).

OPEN IN VIEWER

FIG. 3.

ESC self-renewal and differentiation. (a) Culture LIF inhibits lineage genes and counteracts the antagonism (dark gray lines) and thus sustains ESC self-renewal and (b) ESC differentiation in the absence of LIF. Mutual antagonism between ESC master regulators and lineage genes formed a positive feedback (dark gray lines) that enforced ESCs differentiating into lineage-committed cells. The dashed dark gray lines in (a and b) represent the regulation effect being likely to be indirect or indistinct. The temporal expression profiles of ESC master regulators and lineage-specific genes were scaled to one standard deviation and indicated by color, from gray (low expression) to dark gray (high expression). DC indicated the scaled dissociation constants, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \log _2} \left( {{ \rm{DC / D}}{{ \rm{C}}_{ \rm{0}}}} \right)$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \rm{D}}{{ \rm{C}}_{ \rm{0}}}$$ \end{document} is the dissociation constant corresponding to bivalent. (c) Stationary distribution of cells cumulated from 10,000 cells (10,000 independent experiments with initial conditions selected from uniformly random variables). Two small peaks are, respectively, LN and HN ESCs in culture and the prominent peak is the cell distribution after significantly reducing antagonism between ESC master regulators and lineage genes. (d) The temporal trajectory of differential equation outcomes evolves over time from 100 independent experiments (cells) under the condition of all lineage genes being silenced. The symbol “*” indicated the initial points (conditions) or perturbation of differential equations selected from the uniformly random variables and “O” the converged points at steady state (unless specified otherwise, these symbols elsewhere are the same as above). The expression levels of these regulators were scaled as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \log _2} \left( {{{ \rm{m}}_{ \rm{i}}}{ \rm{ / }}{{ \rm{m}}_{{ \rm{i , 0}}}}} \right)$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${{ \rm{m}}_{{ \rm{i , 0}}}}$$ \end{document} are, respectively, expression levels of Nanog and Oct4 in LN ESC. DC, dissociation constant; HN, high Nanog; LN, low Nanog.

To further investigate the mechanisms of ESC differentiation and heterogeneous expression in ESCs, we significantly reduced antagonism between ESC master regulators and lineage genes in our model. The model outputs showed that ESC master regulators were stably upregulated with HN expression and LN expression disappeared, meanwhile lineage genes were upregulated in both cultures, LIF and LIF withdrawal (Fig. 3c). This result suggested that the antagonism between ESC master regulators and lineage genes had a close connection with ESC differentiation and heterogeneous expression in ESCs.

On close examination of the quantitative model, we found that in the absence of LIF, the downregulation of ESC master regulators, especially Nanog, reduced the antagonizing effects on lineage genes and thus upregulated these genes that depressed further by ESC master regulators. The antagonism between ESC master regulators and lineage genes formed a positive feedback (the output was routed back to the noninverting [+] input) that enforced ESC differentiation into lineage-committed cells (Fig. 3b). In the presence of LIF, however, the antagonism effects were counteracted by LIF-inhibiting lineage genes (Fig. 3a). Culture LIF blockaded the further differentiation route of positive feedback and LN ESCs became trapped in an intermediate state between ESCs and differentiated cells, for example, epiblast cells, which corresponded to the LN attractor. On the other hand, upregulating Nanog or HN expression significantly inhibited lineage genes and constructed the HN attractor. Thus, the regulatory network had two attractors, that is, LN and HN ESCs. As perturbations to the network, antagonism between ESC master regulators and lineage genes functioned as a so-called bistable switch, orchestrating mutually exclusive outcomes, thus most cells cumulated at attractors of LN and HN ESCs and formed two peaks of Nanog expression (Fig. 3c and Supplementary Figs. S8 and S9).

Taken together, our findings suggested that ESC differentiation and heterogeneous expression in ESCs were mainly attributed to the antagonism between ESC master regulators and lineage-specific genes, for example, Gata6, Gata4, and Cdx2. Our model provided a new explanation why ESC master regulators and lineage genes were expressed in a heterogeneous manner in cultures.

When we withdrew LIF and inhibited Nanog, while significantly reducing the antagonizing effects in our model, the model output exhibited the ESC expression profile accompanied with promiscuous expression of multiple lineage genes, which was very similar to that of LIF sustaining ESC self-renewal with LN expression (correlation >0.85). This further reinforces aforementioned results showing that ESC differentiation was predominantly caused by antagonism between ESC master regulators and lineage-specific genes. We suggested that antagonism between ESC master regulators and lineage genes was a predominant reason for ESC differentiation, and the culture LIF and Nanog sustaining ESC self-renewal was mainly attributed to counteracting the antagonism. Moreover, the culture LIF and Nanog could replace each other in inhibiting lineage genes and reducing the antagonism for ESC self-renewal.

Our model outputs showed that the lineage genes with the bivalent chromatin state were not strictly silenced in ESCs and were easily triggered to be strongly upregulated by some perturbations (Supplementary Figs. S10 and S11). Upregulation of the lineage genes antagonized ESC master regulators and enforced ESC differentiation into lineage-committed cells. We speculated that ESC pluripotency might mainly be attributed to bivalent states of lineage development genes. Along with the differentiation, the ESC master regulators and other unrelated lineages were strictly silenced by epigenetic regulations. Our model analysis showed that the strict silence of unrelated lineages constructed a stable attractor with strong ability to suppress various perturbations (Supplementary Figs. S12 and S13). Therefore, after differentiation, it was not easy to trigger ESC master regulator upregulation again and the differentiated cells restored into ESCs with difficulty. Similarly, assuming that the lineage-specific genes were strictly silenced in ESCs, our model outputs also showed an ESC self-renewal attractor with a strong ability to suppress various perturbations, and the ESCs could extensively self-renew in the absence of LIF, Nanog-null, and Oct4 overexpression. That is, the basin of the attractor might have a larger size than LN and HN ESCs (Fig. 3d). These results implied that epigenetic states of unrelated lineages were important for cell-type maintenance.

To further test whether the strict silence of unrelated lineages was necessary for ESC self-renewal, we assumed that lineage-specific genes still kept bivalent states, whereas the expression of these genes was significantly depressed by increasing degradation rates. The model outputs indicated that ESC self-renewal could also be achieved in the absence of LIF, LN expression, Nanog-null, and Oct4 overexpression. Our findings suggested that the epigenetic regulation effects of lineage genes on ESC self-renewal were likely to be indirect by easily triggering upregulation of lineage genes. Upregulating lineage genes and antagonizing ESC master regulators enforced ESC differentiation.

We assumed that in lineage-committed or differentiated cells, ESC master regulators and other lineage genes were strictly silenced. Our model outputs showed that lineage-committed cells could maintain extensively their expression patterns under various perturbations and the basin of attractor had a large size (Supplementary Fig. S13). To further investigate the effect of antagonism between ESC master regulators and lineage genes on lineage cell attractors, the ESC master regulators were reset as bivalent or activated state in differentiated cells. From the model outputs, we found that the ESC master regulators were easily triggered to be upregulated, which antagonized the lineage genes and enforced the differentiated cells to restore into ESCs in the presence of LIF.

Thus, the regulatory model had two attractors, that is, the differentiated cells and ESCs. ESC master regulator upregulation reduced the basin size of lineage-committed cell attractor (Supplementary Figs. S14 and S15). This was very similar to upregulation effects of lineage genes on ESC self-renewal. The results were in agreement with the epigenetic state alterations and differentiated cells restoring into ESCs in the experiment studies of ESCs treated by histone deacetylase inhibitor, Trichostatin A (Karantzali et al., 2008). Our findings suggested that ESC master regulators antagonized lineage-specific genes and upregulation reprogrammed differentiated cells into ESCs. The maintenance of lineage-committed cell types required other lineage genes, including ESC master regulators, being strictly silenced. This further reinforced aforementioned results showing that antagonism among multiple lineages largely influenced the cell-type maintenance.

Our model analysis showed that metazoans, such as Murine, had two types of attractors that characterize different cell types. One was the attractors constructed by strictly silencing all unrelated lineages, for example, lineage-committed or differentiated cells, and the other was the attractor of two or multiple lineages being bivalent or activated simultaneously, for example, LN and HN ESCs. To distinguish these attractors, we referred to the former as the type I attractor and the latter as the type II attractor. Comparatively, the type I attractor had a larger-sized basin and was more robust against perturbations since it was more difficult to overcome epigenetic barriers and convert into other attractors or lineages. In contrast, the type II attractor had a small basin and was metastable and sensitive to perturbations since LIF withdrawal, Oct4 overexpression, and downregulation upset the balance of stable state and the type II attractor disappeared. The attractor basin size largely depended on whether the antagonized lineages were upregulated. It had a large size if antagonized lineages were downregulated, otherwise a reduced size (Supplementary Figs. S10, S11, and S16).

In summary, antagonism between ESC master regulators and lineage-specific genes was a predominant reason for ESCs being prone to differentiation. Counteracting the antagonism by either epigenetic silence or downregulating lineage-specific genes, ESCs could self-renew extensively in the absence of LIF, LN expression, or Nanog-null and Oct4 overexpression. On the other hand, if the ESC master regulators were activated or bivalent, they largely influenced the maintenance of differentiated cells by antagonism as well. Antagonism loss or counteraction by extrinsic cues, for example, culture LIF, allowed multiple lineage genes, including ESC master regulators, being simultaneously activated, thus ESCs showed promiscuous and heterogeneous expression of multiple lineage genes.

A critical role of lineage antagonism in cell plasticity

Contrary to cell-type maintenance, a cell can convert into another cell type either spontaneously by extrinsic cues or by gene perturbation experiments. The ability of a cell to convert into another cell type is referred to as cell plasticity. Earlier work on cell plasticity has shown that cell plasticity largely depends on the epigenetic distance between the starting and target cells (Eminli et al., 2009; Hochedlinger and Jaenisch, 2009; Selvaraj et al., 2010). Our model supported these general notions and found further new factors that affected cell plasticity, in particular antagonism between target and starting lineages, and hence ultimately mapping to cell fate decisions.

To investigate cell plasticity by quantitative model, we first tested reprogramming lineage-committed cells into ESC-like cells. We set the starting cells as lineage-committed cells in which ESC master regulators were strictly silenced. By adding extrinsic cues, for example, viral integration of four genes Oct4, Sox2, cMyc, and Klf4, our model outputs indicated that it was difficult for most differentiated cells to overcome epigenetic barriers and convert into ESCs. This was consistent with low efficiency of induced pluripotent stem (iPS) cells (Mikkelsen et al., 2008; Eminli et al., 2009; Hochedlinger and Jaenisch, 2009; Selvaraj et al., 2010). We then assumed that in starting cells, the target lineage genes were activated or bivalent. It is a frequent event in lineage reprogramming since extrinsic cues or perturbations may trigger two or multiple lineages to be activated simultaneously (Mikkelsen et al., 2008). Our model outputs exhibited two attractors (mutually exclusive outcomes) corresponding to both starting and target cells and a large proportion (40–50%) of starting cells were converted into target cells by upregulating the target lineage genes (gene expression belongs to the attractor basin of target lineage). Although mutually exclusive outcomes caused by antagonism were similar to LN and HN ESCs above, the epigenetic states were rather different from them. Along with the reprogramming, the starting lineage genes were gradually silenced and the reprogrammed cell became a type I attractor (Fig. 4b, c and Supplementary Fig. 17). Taken together, our results reinforced previous studies showing that cell plasticity largely depends on the epigenetic distance between the starting and target cells (Eminli et al., 2009; Hochedlinger and Jaenisch, 2009; Selvaraj et al., 2010).

OPEN IN VIEWER

FIG. 4.

Cell plasticity and antagonism. (a) Top panel: closed-loop network of mutually exclusive lineages (dark gray lines), and bottom panel: closed-loop network of viral integration, the starting and the unrelated lineage genes in induced pluripotent stem cells (dark gray lines and gray lines). (b) The temporal outcome trajectory of the closed-loop network of top panel figure in (a). The attractor basin of target lineage gene is the upper triangle area. If the gene expression belongs to this area, starting cells can be converted into target cells. (c) Stationary distribution of cells cumulated from 10,000 cells (independent experiments). Two small peaks are mutually exclusive outcomes caused by lineage antagonism of top panel figure in (a). The prominent peak is the cell distribution after counteracting antagonism between starting and target lineages. (d) The differential equation outcome trajectory of the closed-loop network in bottom panel figure in (a). Evidently, except for the mutually exclusive outcomes, the new converged set or attractor (5.2, 2.8) shows that starting and the unrelated lineage genes are simultaneously upregulated, which represents the partially reprogrammed or intermediate state of starting and unrelated lineages. In (b) and (d), 100 initial points (conditions) or perturbations are selected from uniformly random variables.

When the antagonism between starting and target lineages was significantly reduced in our model, the model outputs showed rather different outcomes from the mutually exclusive outcomes above, and the genes from both the target and starting lineages were upregulated simultaneously (prominent peak in Fig. 4c). This meant that there was no epigenetic state difference or gap between starting and target lineages, yet the starting cells could not be converted into target cells if there was no antagonism of starting and target lineages. In effect, when two or multiple lineages are simultaneously activated, antagonism among lineages functions as a so-called bistable switch, orchestrating mutually exclusive outcomes. Thus, each cell's transcriptional profile converged to one of starting and target lineages (two small peaks in Fig. 4c). Reducing antagonism, however, blocked further cell-type conversions and cells became trapped in intermediate states of starting and target lineages being activated simultaneously.

During lineage reprogramming, gene perturbations, for example, viral integration of four genes in iPS, may trigger strong upregulation of genes from unrelated lineages, for example, primer endoderm marker Gata6 (Mikkelsen et al., 2008). Our model showed that upregulating unrelated lineage genes might extensively affect iPS. On close examination of two regulatory networks (Fig. 4a), we found that the viral integration was to counteract the antagonism effects and establish a balance between the starting and unrelated lineages by upregulating the unrelated lineage genes. Counteracting the antagonism by viral integration blocked further reprogramming and thus cells became trapped in the partially reprogrammed state that was similar to MCV6 cells (Mikkelsen et al., 2008) (Fig. 4a, d). This further reinforced aforementioned results showing that cell-type conversions largely depended on lineage antagonism. Our model provided an explanation of why the partially reprogrammed cells became trapped in intermediate states. Since multiple lineages were activated, the stable state was evidently a type II attractor. By increasing the degradation rates of the unrelated lineage genes and reducing the dissociation constants of ESC master regulators, our model output gradually converged to the ESC expression profile. This meant that the balance between the starting and unrelated lineages was upset and the attractor disappears by repressing the unrelated lineage genes and thus the intermediate state cells could be further reprogrammed into ESC-like cells. This is supported by siRNA against unrelated lineage genes, and then ESC-like cells appear at a significant frequency in all examined populations of MCV6 cells (Mikkelsen et al., 2008).

Taken together, we concluded that in addition to the epigenetic distance, antagonism between starting and target lineages had a critical role in cell-type conversions. Counteracted or reduced antagonism, a cell could not be converted into other cell types and become trapped in intermediate states of two or multiple lineages being activated simultaneously.

Ordered expression alterations and cell fates

Since multiple lineage development genes are bivalent, ESCs thus have multiple differentiation directions into lineage-committed cells. ESCs differentiating into a desired lineage or homogeneous differentiations are very important for metazoan development. Our model provided more knowledge of how gene expression changes over time during ESC differentiation, hence ultimately mapping to cell fates. In this study, we provided a novel insight into dynamic regulation effects of ESC master regulators on differentiation commitment from the quantitative model.

To investigate the relationship between differentiated lineages and ESC master regulator expression alterations over time, we designed perturbations to the regulators, Oct4 and Nanog. When Oct4 was overexpressed, the model outputs showed that Nanog was immediately repressed. Oct4 overexpression prevented differentiation into trophectodermal lineage (Ivanova et al., 2006; Loh et al., 2006; Lu et al., 2009), while significant repression of Nanog released the barrier of ESCs differentiating into ensoderm and primer endoderm lineages, the model outputs thus exhibited the expression patterns of primer endoderm and ensoderm lineages. On the other hand, as Oct4 was depressed, Nanog was postponed repressed. Repressing Oct4 canceled the barrier of ESCs differentiating into trophectodermal lineage, while the persisting expression of Nanog kept the obstacle of ESCs differentiating into endoderm and primer endoderm cells, thus the trophectodermal lineage genes exhibited a strong upregulation in our model. When Nanog was downregulated, we observed that Oct4 was postponed repressed in model output. The maintenance of Oct4 upregulation blocked ESCs differentiating into trophectodermal lineage, while downregulating Nanog released the barrier of ESCs differentiating into endoderm cells. Finally, ESCs differentiated into endoderm lineage (Fig. 5a and Supplementary Figs. S18–S20). The model outputs and lineage commitments were consistent with the experiment studies (Ivanova et al., 2006; Loh et al., 2006; Chambers et al., 2007; Lu et al., 2009).

OPEN IN VIEWER

FIG. 5.

Cell fate decisions. (a) Oct4 and Nanog are ordered and repressed; the quantitative model has only a single attractor. (b) When Oct4 and Nanog are downregulated (disordered) simultaneously, the model exhibited at least three attractors. In (a) and (b), 100 initial points (conditions) or perturbations are selected from uniformly random variables. (c) Ordered expression alterations of ESC master regulators ensure that the populations of ESCs differentiate into homogeneous lineage, otherwise promiscuous cell types of multiple lineages. (d) Antagonism among lineages functions as a so-called bistable switch, orchestrating mutually exclusive outcomes. Antagonism loss or counteracting antagonism leads to cells becoming trapped in intermediate states of multiple lineages being upregulated simultaneously. Bicolor cells indicated intermediate state cells of two lineages.

Our findings suggested that the ordered repression of Nanog and Oct4 had a close connection with ESC differentiation fates. Similar results were also found in other lineage commitment approaches, for example, in hematopoietic lineages, the order of expression of key transcription factors is critical for their interplay to selectively drive lineage specification programs, by which stem cells can generate multiple lineage-committed cells in a hierarchical manner (Iwasaki et al., 2006).

When Oct4 and Nanog were depressed simultaneously, we observed that there was competition or antagonism among multiple lineage genes, and the regulatory model had multiple attractors or stable states. Evidently, ESCs could not differentiate into a single (desired) lineage, for instance, the trophectodermal lineage by Oct4 repression above. The differentiated cells exhibited promiscuous cell types of multiple lineages (Fig. 5b, c). This was consistent with the experiment studies (Walker et al., 2007), in which LIF withdrawal caused simultaneous downregulation of ESC master regulators and triggered upregulation of genes from multiple lineages. The counterexample further confirmed that homogeneous differentiation requires an ordered expression alteration of ESC master regulators. We also found that both LN and Nanog-null ESCs were prone to differentiate into epiblast as Oct4 downregulation. This was not consistent with the differentiation commitment of HN ESCs and further confirms the effect of Nanog-postponed repression on differentiation commitments.

Examining the quantitative model, we found that in the beginning of ESC differentiation, the pivotal genes from multiple lineages were bivalent and subsequently were gradually silenced along with ESC differentiation and were strictly silenced until the desired lineage genes were stably upregulated (Supplementary Figs. S18–S20). The strict silence of unrelated lineage genes was implemented by upregulating the desired lineage genes. Thus, if there was no ordered repression of ESC master regulators, the unrelated lineage genes with bivalent state might possibly be strongly upregulated before the stable upregulation of desired lineage genes, which antagonized the genes from the desired lineage and resulted in ESCs differentiating into the undesired lineages.

From the differentiated cell fates in ESCs and hematopoietic stem cells (Iwasaki et al., 2006), we suggested that the ordered expression alterations of key or master regulators played a central role in controlling differentiated cell fates of ESCs and adult stem cells that had multiple lineages being activated or bivalent and thus multiple differentiation directions.

Model description

We developed a unified computational model for gene expression involved in cell fate choices. The genes involved in the quantitative model were ESC master regulators Oct4, Sox2, Nanog, Esrrb, Tbx3, and Tcl1, etc.; lineage-specified pivotal genes, including trophectoderm lineage: Hand1, Eomes, Cdx2, Elf5, and Tead4; ectoderm lineage: Lhx5, Otx1, and Hoxb1; endoderm lineage: Gata6, Gata4, Foxa2, Gcnf, and Hex1; mesoderm lineage: Myf5, T, and Gsc; and other lineages, such as neural crest lineages.

Transcription and degradation rates were computed as a function of intrinsic transcription factors, extrinsic cues, epigenetic modifications, and antagonism between ESC master regulators and lineage-specific pivotal genes. Positive autoregulations and cross-regulations of ESC master regulators and lineage-specific pivotal genes contributed to transcriptional rates. The extrinsic cue effects were to enhance the transcriptional rates of ESC master regulators and inhibit lineage genes by increasing degradation rates. Antagonism between ESC and lineage as well as different lineages was exerted by both degradation and epigenetic regulation alterations. Epigenetic state change of ESC master regulators and lineage genes during ESC differentiation and lineage reprogramming contributed to dissociation constant alterations of transcriptional factors binding the promoters of target genes.

Combining above description, the differential equations of gene expression were established according to kinetic strategies. We denoted the mRNA levels of ESC master regulators Oct4, Sox2, and Nanog by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_O}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_S}$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_N}$$ \end{document} and the j-th genes of the i-th lineage-specified pivotal genes by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_{i , j}}$$ \end{document} . Note that these variables were the functions of time variable; for simplifying, we omitted the time variable. According to previous studies, the levels of genes were determined by real-valued, nonlinear regulation functions that took the Michaelis–Menten and Hill form (Turcotte et al., 2008; Kalmar et al., 2009). Thus, the differential equations of ESC master regulators were \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { d { m_N } } { dt } } &= T { R_N } - { D_N } { m_N } \\&= { \alpha _N } + { \frac { { \beta _ { N , N } } { { \left( { { m_N } } \right) } ^n } + { \beta _ { A , N } } { A^n } + { \beta _{ O , N } } { { \left( { { m_O } } \right) } ^n } + { \beta _ { OS, N } } { { \left( { { m_ { OS } } } \right) } ^n } } { K_N^ { } \left( t \right) + { { \left( { { m_N } } \right) } ^n } + { { \left( A \right) } ^n } + {\rho} _ { O , N } { { \left( { { m_O }} \right) } ^n } + { \rho _ { OS , N } } { { \left( { { m_ { OS }} } \right) } ^n } } } \\ &\quad - { d_N } { m_N } - { F_1 } \left( { { m_O } } \right) { m_N } - \mathop \sum \limits_ { i , j} { { F_2 } \left( { { m_ { i , j } } } \right) { m_N }} \tag {1}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { d { m_O } } { dt } } &= T { R_O } - { D_O } { m_O } \\ &= { \alpha _O } + { \frac { { \beta _ { O , O } } { { \left( { { m_O } } \right) } ^n } + { \beta _ { N , O } } { { \left( { { m_N} } \right) } ^n } + { \beta _ { OS , O } } { { \left( { { m_ { OS} } } \right) } ^n } } { K_O^ { } \left( t \right) + \left( { { m_O } } \right) ^n } + { \rho _ { N , O } } { { \left( { { m_N } } \right) } ^n } + { \rho _ { OS , O } } { { \left( { { m_ { OS } }} \right) } ^n } } \\ & \quad- { d_O } { m_O } - \mathop \sum \limits_ { i , j } { { F_3 } \left( { { m_ { i , j } } } \right) { m_O }} \tag { 2 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _N} , { \alpha _O}$$ \end{document} were the basal transcriptions of Nanog and Oct4, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_{OS}}$$ \end{document} was mRNA abundance of the heterodimer of Oct4 and Sox2, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$K_N^{} ( t )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$K_O^{} ( t )$$ \end{document} were the dissociation constants, A was the extrinsic signal strength, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{N , N}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{O , O}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{O , N}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{O , O}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{A , N}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{OS , N}}$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{OS , O}}$$ \end{document} were the transcription strength of autoregulations and cross-regulations between Nanog and Oct4, as well as the heterodimer of Oct4 and Sox2. If Nanog did not positively regulate Oct4 or the regulation effects of Nanog on Oct4 was very weak, the transcriptional strength \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{N , O}}$$ \end{document} and the coefficient \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \rho _{N , O}}$$ \end{document} would disappear or set to zero by optimization trainings from experiment data. The degradation rates of Nanog and Oct4 should comprise entries \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_2} \left( {{m_{i , j}}} \right) {m_N}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_3} \left( {{m_{i , j}}} \right) {m_O}$$ \end{document} , which were the functions of lineage-specific gene expression levels. In addition, the dissociation constants of Nanog and Oct4, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$K_N^{} ( t )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$K_O^{} ( t )$$ \end{document} , were the functions of lineage-specific gene expression levels, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${d_N}{m_N}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${d_O}{m_O}$$ \end{document} were the basal degradations. The superscript, n, was the cooperativity exponent of Hill functions. In general, n = 4.

If the free parameters in above equations identified by the optimization training were significantly smaller and the contributions were also smaller compared with others, we neglected these components or regulation relationship. This meant that there were no direct regulation relationships between two network components. Thus, from our model, we could capture some unknown regulation relationships. Motivated by the idea, we added a tentative entry to test whether Nanog had a positive regulation effect on Oct4. In addition, we could also test other regulation relationships, for example, the regulation relationship between Nanog, Oct4, and lineage-specific genes, by this method.

Similar to Nanog and Oct4, we could also establish the differential equations of other ESC master regulators, such as the heterodimer of Oct4 and Sox2, Sox2, Esrrb, Tbx3, and Tcl1.

We could also obtain the differential equations of lineage-specific gene expression: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \frac { d { m_ { i , j } } } { dt } } &= T { R_ { i , j } } - { k_ { i , j } } { m_ { i , j } } \\ &= { \alpha _ { i, j } } + { \frac { { \beta _ { i , j } } { { \left( { { m_ { i , j } } } \right) } ^n } + { \beta _ { O , i , j } } { { \left( { { m_O } } \right) } ^n } + \mathop \sum \limits_ { q = 1 } ^Q { { \beta _ { i , q } } { { \left( { { m_ { i , q } } } \right) } ^n }} } { { K_ { i , j } } \left( t \right) + { { \left( { { m_ { i , j } } } \right) } ^n } + { \rho _ { i , j } } { { \left( { { m_O }} \right) } ^n } + \mathop \sum \limits_ { q = 1 } ^Q { { \beta _{ i , q } } { { \left( { { m_ { i , q } } } \right) } ^n } } } } \\ \quad &- { d_ { i , j } } { m_ { i , j } } - { F_7 } \left( A \right) { m_ { i , j } } - { F_5 } \left( { { m_N } , { m_O } } \right) { m_ { i , j } } \\\quad &- \mathop \sum \limits_ { l = 1, l \ne i } ^L { { F_6 } \left( { { m_l } } \right) {m_ {i , j}}}, \qquad {i = 1} , {\cdots} {j = 1} , {\cdots} \tag {3} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \beta _{i , j}}$$ \end{document} was the transcriptional strength of autoregulation and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${m_{i , q}} , {q = 1} , \cdots \ {,q} {\ne j}$$ \end{document} were the other genes specified by the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i$$ \end{document} th lineage, which were the cross-regulations of the genes specified by the same lineage.

Since we considered only the lineage-specific pivotal genes that were located almost at the top layer of the hierarchical regulatory network, there were relatively fewer cross-regulations. In addition, our investigation and other experiment studies (Niwa et al., 2005) also found that the effects of these cross-regulations on ESC self-renewal, differentiation, and lineage cell-type maintenance were smaller than the autoregulations. Thus, these cross-regulations might either be neglected or taken as tentative entry similar to Equations (1) and (2).

In cell-type maintenance, only the desired lineage genes were activated and other lineage-specific genes, including ESC master regulators, were silenced, which was unlike ESC self-renewal. Thus, the quantitative model would only have the desired lineage genes on and other lineage genes off. In lineage reprogramming, however, the extrinsic cues, for example, viral integration of four genes Oct4, Sox2, cMyc, and Klf4, might trigger upregulation of unrelated lineage genes. In ESC self-renewal, LN ESCs also exhibited multiple lineage-specific genes, which were simultaneously activated. Thus, in these cases, there were multiple activated lineages in the quantitative model.

Dissociation constants and degradation functions

In this study, we required further quantifying the relationships between the dissociation constants and the associated chromatin regulators and the time-varying law of epigenetic state alterations during lineage reprogramming and cell-type conversions for ESC master regulators and lineage-specific genes.

The activated level state of a gene, for example, H3K4me3, was denoted by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$a$$ \end{document} and the repressed level, for example, H3K27Me3, was denoted by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$r$$ \end{document} . By the chromatin regulatory factors, which recruited specific histone acetyltransferases (HATs), histone demethylases (HDMs), and histone deacetylases (HDACs) to the promoter of the gene, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$a$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$r$$ \end{document} could be described as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { da } { dt } } = T { R_a } - { k_a } a - k { \prime _a } r \tag { 4 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { dr } { dt } } = T { R_r } - { k_r } r - k { \prime _r } a \tag { 5 } \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$T{R_a} , {k_a} , T{R_r} , {k_r}$$ \end{document} were, respectively, the contributions of the regulatory factors, which recruited the specific HATs, HDMs, and HDACs to the promoter of the gene, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$k{ \prime _a} , k{ \prime _r}$$ \end{document} were the antagonism between the activated and silenced regulators.

Take Oct4 for example, as mentioned above, the ESC master regulators regulated the epigenetic states of their own genes. Thus, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$T{R_a}$$ \end{document} was a function of Oct4 and Nanog mRNA levels. In addition, the lineage-specific gene upregulation silenced the ESC master regulators and then contributed to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${k_a}$$ \end{document} . So, we had \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { d { a_O } } { dt } } = T { R_a } \left( { { m_O } , { m_N } } \right) - { k_a } \left( { { m_ { i , j } } } \right) { a_O } - k { \prime _a } { r_O } \tag { 6 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ \frac { d { r_O } } { dt } } = T { R_r } \left( { { m_ { i , j } } } \right) - { k_r } \left( { { m_O } , { m_N } } \right) r - k { \prime _r } { a_O } \tag { 7 } \end{align*} \end{document}

If \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${a_O} \to c > 0$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${r_O} \to 0$$ \end{document} , then Oct4 was activated. Otherwise, this gene was strictly silenced. If \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${a_O} \to 0 < {c_1} < c$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${r_O} \to {c_2} > 0$$ \end{document} , Oct4 might be bivalent. Since chromatin configuration determined the accessibility of promoters to transcription and the performance of the transcription machinery, the dissociation constants were the functions of the activated or silenced states of a gene, which might be inversely proportional to activated level \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$a$$ \end{document} and silenced level \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$r$$ \end{document} , that is, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${K_O} ( t ) = {k_O} \left( {{a_O} , {r_O}} \right)$$ \end{document} .

Certainly, the differential Equations (4)–(7) were much more complicated if we substituted these dissociation constants into the quantitative model. For simplifying, we alternately took Hill functions of multiple variables to approximate the dissociation constants, that is, the dissociation constants of Nanog as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ K_N } ( t ) &= { G_N } \left( { { m_ { OS } } , { m_N } , { m_ { i , j } } } \right) = \left( { { \alpha _ { N , N }} - { \frac { { \alpha _ { N , N } } { { \left( { { m_N } } \right) } ^n } } { { c_ { N , N } } + { { \left( { { m_N } } \right) } ^n } } } } \right) \\ \quad &\times \left( { { \alpha _{ N , O } } - { \frac { { \alpha _ { N , O } } { { \left( { { m_ { OS } } } \right) } ^n } } { { c_ { N , O } } + { { \left( { { m_{ OS } } } \right) } ^n } } } } \right) \times \mathop \sum \limits_ { i = 1 } ^I { { \frac { { \alpha _N } { { \left( { { m_{ i , j } } } \right) } ^n } } { { c_N } + { { \left( { { m_ { i, j } } } \right) } ^n } } } } \tag { 8 } \end{align*} \end{document}

and the dissociation constant of regulator Oct4 as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { K_O } ( t ) &= { G_O } \left( { { m_ { OS } } , { m_ { i , j } } } \right) \\ &= \left( { { \alpha _ { O , O } } - { \frac { { \alpha _ { O , O } } { { \left( { { m_O } } \right) } ^n} } { { c_ { O , O } } + { { \left( { { m_O } } \right) } ^n } }} } \right) \mathop \sum \limits_ { i = 1 } ^I { { \frac { { \alpha _O } { { \left( { { m_ { i , j } } } \right) } ^n } } { { c_O } + { { \left( { { m_ { i , j } } } \right) } ^n } } } } \tag{9} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _{N , O}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _{N , N}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _N}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${c_N} , {c_{N , N}} , {c_{N , O}}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _{O , N}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \alpha _O}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${c_O} , {c_{O , O}}$$ \end{document} were the free parameters, which were unknown and needed to be trained or estimated by the learning algorithm.

Similarly, dissociation constants of lineage-specific pivotal genes were also the functions of ESC master regulators and other lineage genes. Thus, we had \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { K_ { i , j } } \left( t \right) &= { G_ { i , j }} \left( { { m_ { OS } } , { m_N } , { m_ { l , q } } } \right) \\ &= { \frac { { \alpha _ { N , i , j } } { { \left( { { m_N } } \right) } ^n } } { { c_ { N , i , j } } + { { \left( { { m_N } } \right) } ^n } } } { \frac { { \alpha _ { O , i , j } } { { \left({ { m_ { OS } } } \right) } ^n } } { { c_ { O , i , j } } + { { \left( { { m_ { OS } } } \right) } ^n } } } \\ \quad &\times \mathop \sum \limits_ { i = 1 } ^I { \left( { { \alpha _ { i , j }} - { \frac { { \alpha _ { i , j } } { { \left( { { m_ { l , i } }} \right) } ^n } } { { c_ { i , j } } + { { \left( { { m_ { l , i} } } \right) } ^n } } } } \right) } \\ \tag { 10 } \end{align*} \end{document}

Evidently, the dissociation constants in Equations (4)–(10) remained constant in ESC self-renewal and cell-type maintenance and changed to adapt to the new cell types in ESC differentiation and lineage reprogramming.

We observed that as the ESC master regulators were upregulated and lineage genes were downregulated, the dissociation constants were very small, which represents the activated state of ESC master regulators. Otherwise, at the silenced state of ESC master regulators, the mRNA abundance of the ESC master regulators was very low and lineage-specific genes were upregulated, and the dissociation constants would become significantly great constants. Therefore, the increment or decrement of the dissociation constants adjusted the accessibility of promoters to transcription and the performance of the transcription machinery, which could effectively represent the activated, bivalent, and silenced states of genes.

In addition to determining the dissociation constants, we also needed to determine the functions in degradation rates, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_i} \left( x \right)$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$i = 1 , 2 , \ldots\,$$ \end{document} . In this study, we assumed that the function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_i} \left( x \right)$$ \end{document} could be approximated by Taylor series expansion as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{F_i} \left( x \right) = {f_{i , 0}} + {f_{i , 1}}x + {f_{i , 2}}{x^2} + \cdots \tag{11} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${f_{i , l}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$l = 0 , 1 , \ldots$$ \end{document} were the unknown constants needed to be estimated from optimization training and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${f_{i , 0}}$$ \end{document} could be combined into the basal degradations. In the approximation procedures, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_i} \left( x \right)$$ \end{document} could be approximated by the first-order Taylor series expansion, that is, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${F_i} \left( x \right) \approx {f_{i , 1}}x$$ \end{document} .

Parameter selection and adjustment

The model had three types of free parameters for each gene, representing values that were typically unknown: (1) the strengths of autoregulations, cross-regulations, and extrinsic cues; (2) the dissociation constants and their variations during cell-type conversions; and (3) the degradation parameters.

Unified parameters were fitted to minimize the error between the measured and model outputs in each cell state, including ESC self-renewal and lineage cell maintenance, ESC differentiation, and lineage reprogramming. The full mathematical details and fitting procedures of our model are described in the Supplementary Data.

We took Nanog, Oct4, Sox2, Cdx2, and Gata6 as examples to test the estimation of the free parameters. The results showed that it could model the expression patterns of each cell state with good accuracy. The parameter set for these pivotal genes is given in Supplementary Table S2.

Model analysis

To quantitatively investigate the mechanisms of cell-type maintenance attractors, cell plasticity, ESC differentiation, and lineage reprogramming, we designed perturbation vectors (Oct4-E, Nanog-E, LC-E, LIF, ESC-Ep, LC-Ep, EL-A, LL-A) to ESC master regulators and lineage genes, including downregulation and overexpression, withdrawal and addition of extrinsic cues, and alteration of epigenetic states. In the perturbation vectors, -E represented gene expression levels that were respectively downregulation, normal, and overexpression; LIF had two states of withdrawal and addition; -Ep represented the epigenetic states of activated, bivalent, and silenced; LC indicated the lineage-committed cell; EL-A and LL-A were the antagonism between ESC master regulators and lineage genes, and antagonism among lineages, respectively, which were chosen to be either zero or values from optimal training.

In addition to perturbation vectors above, we also assumed that the transcriptional rates, degradation rates, and dissociation constants followed a Gaussian distribution with the optimal training values as mean values to simulate the intrinsic and extrinsic noises of gene expression. Thus, the model outputs were scattered in a limited region, instead of converging or focalizing at a point.

As perturbations to our model, the mechanism of cell fate choices was found from the model output profiles of numerous independent experiments, for example, 10,000 independent experiments. Each independent experiment might be equivalent to a single cell, thus the mechanism or result of cell fate choices was an averaged result over numerous cells. The results of cell-type maintenance and conversion under various perturbations are shown in Supplementary Table S1.