 Research Article
 Open Access
 Published:
Nanog induced intermediate state in regulating stem cell differentiation and reprogramming
BMC Systems Biology volume 12, Article number: 22 (2018)
Abstract
Background
Heterogeneous gene expressions of cells are widely observed in selfrenewing pluripotent stem cells, suggesting possible coexistence of multiple cellular states with distinct characteristics. Though the elements regulating cellular states have been identified, the underlying dynamic mechanisms and the significance of such cellular heterogeneity remain elusive.
Results
We present a gene regulatory network model to investigate the bimodal Nanog distribution in stem cells. Our model reveals a novel role of dynamic conversion between the cellular states of high and low Nanog levels. Model simulations demonstrate that the lowNanog state benefits cell differentiation through serving as an intermediate state to reduce the barrier of transition. Interestingly, the existence of lowNanog state dynamically slows down the reprogramming process, and additional Nanog activation is found to be essential to quickly attaining the fully reprogrammed cell state.
Conclusions
Nanog has been recognized as a critical pluripotency gene in stem cell regulation. Our modeling results quantitatively show a dual role of Nanog during stem cell differentiation and reprogramming, and the importance of the intermediate state during cell state transitions. Our approach offers a general method for analyzing key regulatory factors controlling cell differentiation and reprogramming.
Background
Phenotypic celltocell heterogeneity has been viewed as an important hallmark for both pluripotent embryonic stem (ES) cells [1,2,3,4] and multipotent adult stem cells [5,6,7]. Such nongenetic heterogeneity within stem cell populations suggests the coexistence of multiple cellular states manifested by transcriptomewide gene expression patterns [2, 5]. The dynamical interconversion among these gene expression states is driven by various types of fluctuations, including protein and RNA (mRNA, miRNA) expression noise in gene regulatory networks for cell fate determination [8,9,10], the noise of external input signals [11], the reestablishment of histone modification pattern on the chromatins of daughter cells after each cell division [12].
Nanog has been widely studied as one of the heterogeneously expressed genes in ES and induced pluripotent stem (iPS) cell populations [13, 14]. It plays a central role in maintaining the cellular selfrenewal and pluripotency [15,16,17]. In terms of the gene regulatory network defining cellular potency, Nanog has direct mutual interactions with two other core stem cell specific genes Oct4 and Sox2. They occupy the core position of the network, and link up most downstream genes controlling cellular stemness and cell fate specification [18,19,20,21].
Experimental data from mice ES cell lines indicate that the Nanog expression level in cell population shows a typical bimodal distribution: Approximately 80% cells show a relatively high level of Nanog, while the rest of 20% cells keep a low level of Nanog which is close to the Nanog expression level in differentiated cells [22,23,24]. This can be partly explained by cell differentiation due to some environmental factors in the experiments [22, 24]. However, after sorting and purifying multiple populations of ES cells with different Nanog expression levels, and the subsequent selfrenewal culture of these “subpopulations”, the same bimodal Nanog distribution can be eventually reestablished from all subpopulations with the scale of recovering time for about 10 days [22, 24]. These results suggest the coexistence of two cellular “substates” within the pluripotent stem cell state: highNanog substate and lowNanog substate, rather than one welldefined, homogenous pluripotent cellular state. Furthermore, the two substates spontaneously go through reversible interconversion toward a dynamical equilibrium.
Very few studies have addressed the functional roles of the bimodal heterogeneity of Nanog expression in the differentiation and specification process of stem cells. Previous experimental results have shown that the lowNanog stem cell subpopulations are more likely to differentiate than highNanog subpopulations, under the induction of differentiation signals [24,25,26,27]. One hypothesis based on these results is: During the stem cell differentiation process, the lowNanog state of stem cell functions as the “gatekeeper” state. The preparation of the “gatekeeper” state of stem cells and the presence of external signals are two necessary conditions for stem cell differentiation [10, 26, 27].
Some previous modeling works have attempted to study the Nanog bimodal distribution and the dynamical interconversion between high/lowNanog states based on network topology and gene regulatory interactions, including stochastic transition between the two states driven by expression noise [25] and gene expression oscillation [28, 29]. For instance, Glauche et al. modeled a threenode network including Oct4, Nanog and one unknown Nanog repressor based on bistable switch mechanism or oscillatory mechanism [25]. The stochastic bistable model often requires a high level of gene expression noise to fit the variation observed in the experiment, and lacks the ability to explain different variation levels for different genes. On the contrary, the excitable model can have relatively large excitable excursions in the phase space triggered by much smaller noise to avoid such disadvantages. Excitation mechanism has been utilized to explain bimodal distribution of gene expression and stochastic transitions among multiple metastable cell fates [30,31,32]. In particular, Kalmar et al. proposed a simplified twodimensional Oct4Nanog network based on excitation mechanism to model bimodal distribution of Nanog [24]. However, one limitation of the previous models for Nanog heterogeneity is that all lineage specification genes are neglected, thus the role of Nanog stochasticity in stem cell differentiation and iPS cell reprogramming cannot be quantitatively analyzed.
Here, we propose a mathematical model to investigate the contrary roles of Nanog in stem cell differentiation and iPS cell reprogramming. By introducing a fivenode stem cell gene regulatory network, we first verify the Nanog bimodal distribution and high/lowNanog state interconversion, and then investigate the functional roles of Nanog bimodality. For the role of Nanog in stem cell differentiation, we show that the lowNanog state serves as an intermediate state in the differentiation process, which is a stricter description than the “gatekeeper state”. The differentiation process through the lowNanog intermediate has a lower barrier on the potential landscape so that it is more prone to happen. While for the iPS cell reprogramming, we demonstrate the existence of lowNanog state dynamically leads to a negative effect, and additional Nanog activation is necessary for increasing reprogramming efficiency via preventing reprogrammed iPS cell from trapping in the lowNanog state by overactivated Oct4/Sox2.
Results
The fivenode stochastic gene regulatory network controlling cellular stemness
We first build a fivenode gene regulatory network model, including three pluripotent genes: Oct4, Sox2, Nanog, and two mutually antagonistic lineage specifiers: MEs (representing mesendoderm genes) and ECTs (representing ectoderm genes) (Fig. 1a). The modeling of interactions among Oct4, Sox2, MEs and ECTs is inherited from the previous fournode “Seesaw” model for the iPS cell reprogramming [33], including selfactivation of Oct4 and Sox2, mutual repression between MEs and ECTs, and multiple interactions among pluripotent genes and specifier genes. Both Oct4 & Sox2 and MEs & ECTs are assumed to be symmetric for simplicity. In the fivenode network, we introduce an additional stem cell gene: Nanog. The regulations among Oct4, Sox2 and Nanog consist of selfactivation of Nanog, activation of Nanog to Oct4 and Sox2 [18, 34], and one negative feedback of combined Oct4Sox2 to Nanog. For the negative feedback, low concentration of Oct4Sox2 weakly activates Nanog, while high concentration of Oct4Sox2 strongly represses Nanog, which is supported by some experimental evidence [18, 34]. Stochastic differential equations are used to formulate this fivenode network (see Methods and Materials for model details and Additional file 1: Table S1 for the choice of parameters).
Stochastic transition between highNanog state and lowNanog state
We performed simulations of the proposed model to examine the dynamics of Nanog gene expression. The computational results of trajectory density sampling in the phase space produced four regions with relatively high sampling probability, implying stable or metastable cellular states (Fig. 1b and c). Two of them represent the differentiated cellular states: ME state & ECT state, where one of the MEs or ECTs gene groups is highly expressed, and Oct4, Sox2, Nanog are low expressed for both two states (Additional file 2: Figure S1). The other two states correspond to the highNanog and lowNanog cell subpopulations within the heterogeneous pluripotent cellular state (both of them with high Oct4/Sox2 level) (Fig. 1d). The negative feedback loop between Nanog and Oct4/Sox2 plus the weak activation from Oct/Sox2 to Nanog constitute an excitable system, leading to a recurring stochastic transition between highNanog and lowNanog substates even with relatively low gene expression noise (Fig. 1d). In order to reveal the excitable mechanism in a more generalizable way, we further abstracted a minimal simplified twodimensional Oct4Nanog model and preformed the phase plane analysis as illustrated in Additional file 3: Figure S2A (see the formulation in Methods and Materials and the parameters in Additional file 4: Table S2).
The bimodal distribution of Nanog expression level within the pluripotent cell population under dynamical equilibrium can be naturally generated from our fivenode model (Fig. 2a and Additional file 3: Figure S2B), which is able to quantitatively fit the experimental data better than the original fitting by the Kalmar’s model [24]. Also, the simulation qualitatively matches the normalized results from two single cell RNAseq data sets (using mouse ES cells cultured in serum+LIF) [4, 35], which show a highlevel expression of Oct4 and a bimodal expression of Nanog in Fig. 2a, but a broader distribution of the Sox2 expression compared to the simulation. To verify that the distribution of high/lowNanog cell subpopulations are under interconversion toward dynamical equilibrium, we repeat the reestablishment process of Nanog expression bimodality from sorted stem cell subpopulations, initiating from all highNanog state cells, or all lowNanog state cells. Both two cell subpopulation groups could reestablish the Nanog bimodal distribution after reaching the dynamical equilibrium for ~ 10 days [18, 34], displayed in Fig. 2b. The distributions of the dynamical dwell time at the high/lowNanog states are shown in Fig. 2c. Note that they do not obey exponential distribution, which is one characteristic of stochastic transition between bistable switch systems, indicating an essential difference of excitable system and stochastic bistable switch.
Nanog as “catalyst”: Differentiation through the lowNanog state has a lower energy barrier
Now we consider possible functional roles of Nanog dynamics in the process of stem cell differentiation. In order to investigate the functional role of Nanog dynamics in the fivenode network model, it is necessary to introduce a control model that excludes the dynamics of Nanog. To build up the control model, the concentration value of Nanog in the fivenode model is fixed as a constant, and the steady state values of the other four genes remain unchanged (see details of the control model in Methods and Materials). Therefore, we are able to compare the parallel behaviors of the wildtype (WT) model (with Nanog dynamics) and the control model (without Nanog dynamics).
In order to quantitatively assess the tendency of cell state transition from the stem cell state (the highNanog substate in the WT model, or the homogeneous pluripotent state in the control model) toward one of the differentiated states driven by gene expression noise (rather than by parameter changing and bifurcation), we employ the minimum action path (MAP) to locate the transition path based on the WentzellFreidlin large deviation theory [36]. This approach allows us to calculate the most probable transition path between the initial and final states in a random dynamic system by minimizing the action functional, and has been successfully applied in many biological systems, such as gene switching in zebrafish hindbrain [37] and budding yeast cell cycle [38].
By calculating the transition paths (i.e., MAPs) of cellular differentiation for the WT model and the control model, the quantitative comparison of differentiation tendency is proceeded by two ways: (i) The value of the potential landscape along the transition path; (ii) The monotonous increasing action function along the transition path. For the first approach, we apply the landscape function for complex biological network dynamics, which often corresponds to the nonequilibrium systems, to serve as an analogy of the energy function in the equilibrium systems. There are three representative landscape theories in literatures: 1) the “potential landscape” from the steady state distribution of stochastic differential equations (SDEs) [39]; 2) the “quasipotential” from FreidlinWentzell’s large deviation theory [38]; 3) the potential function obtained through the SDE decomposition by the fluctuationdissipation theorem [40], which can be viewed as a mapping between a nonequilibrium dynamical system and a Hamiltonian system [41]. Here, we use the potential landscape defined by the negative logarithm of the steady state probability distribution P_{ ss } as the landscape function, i.e., − log(P_{ ss }), which is the same one used in [39].
The potential landscape along the transition path is similar to the free energy along the reaction coordinates, and it is very intuitive to illustrate the transition intermediates and the energy barriers. However, it should be noticed that the potential landscapes only provide rough sketches because the gene regulatory network systems are usually far from gradient systems. The transition path could be different from the reaction coordinates. Thus, we also apply the second approach to calculate the action function along the transition path for further validation. The action barrier is a much stricter physical quantity than the energy barrier to describe the transition tendency and estimate the transition rate [36].
We apply the minimum action method in [42] to compute the minimum action path of differentiation for the WT model (Fig. 3a) and the control model (Fig. 3b) (see numerical method in Methods and Materials). For the WT model, the most probable path (i.e., minimum action path) passes the lowNanog substate on the potential landscape, which is robust to the selection of initial path for the action minimization algorithm (Additional file 5: Figure S4). Also, it superimposes the first half of the excitable excursion of the typical ODE trajectory illustrated in green dotted line in Fig. 3a, showing the connection between the system excitability and the differentiation transition path. In the WT model, the energy along the transition path is a typical doublepeak plot (blue curve in Fig. 3c), and the action is a typical twostage stepwise function (blue curve in Fig. 3d), indicating the existence of one intermediate and two barriers along the transition path. The corresponding lowNanog substate on the landscape (Fig. 3a), the intermediate state of the potential landscape (blue curve in Fig. 3c), and the plateau part of the monotonous action (blue curve in Fig. 3d) all point to the same conclusion: the lowNanog state serves as the intermediate state of the differentiation process. For the transition path in the control model, the potential landscape (pink curve in Fig. 3c) has a single large peak, and the action is onestage stepwise (pink curve in Fig. 3d), indicating there is no intermediate state for the differentiation state transition process in the control model.
By comparing both barriers of potential landscape (Fig. 3c) and action values (Fig. 3d) of the entire transition path of the WT model and the control model, we arrive at an important conclusion that the existence of the lowNanog intermediate state during the differentiation process can significantly decrease the barrier of potential landscape and the action value along the transition path. This can be also understood qualitatively: low Nanog expression level leads to temporary low level of Oct4 and Sox2, which in turn release the inhibition of Nanog to the lineage specifier genes, making the high expression of one specifier gene and the differentiation process more prone to happen. Such a role of lowNanog intermediate state is intuitively similar to the enzyme in catalysis process. The enzymesubstrate complex is often an intermediate state between the substrate and the product that significantly decreases the energy barrier of the chemical reaction.
Intrinsic bimodal Nanog expression slows down iPS reprogramming, while exogenous induction of Nanog improves reprogramming efficiency
Next we investigate how bimodal heterogeneity of Nanog may affect induced cellular reprogramming process, i.e., the enforced transition process from somatic cell state (ME or ECT state) to the stem cell state. First, we notice that the minimum action path of reprogramming without additional exogenous Nanog induction does not pass through the Nanog intermediate state, very different from the path of differentiation (Additional file 6: Figure S5).
A common way of iPS cell reprogramming is to activate, or repress the expression of particular gene groups with various types of external inductions, such as virus transformation of vectors expressing target gene products, RNA interference, and small chemical molecule stimulation [13, 14]. Activating gene inductions can be modeled as constant terms indicating additional basal expression rates, noted as (C_{ O }, C_{ S }, C_{ N }, C_{ M }, C_{ E }) below, and repressive inductions can be represented by multipliers of regulated production rate, noted as (I_{ O }, I_{ S }, I_{ N }, I_{ M }, I_{ E }). Details of the mathematical formulations are described by Eqs. (3.1–3.5) in Methods and Materials.
One may adjust the induction parameters in different combinations to change the Waddington landscape of iPS induction to enable the induced somatic cells falling into the basin for stem the cell. At the end of the induction process, the induction parameters are reset to the original values, in order to model the retraction or degradation of the input vectors, RNA or other chemical molecules [43]. The original Waddington landscapes will be then restored and the successfully induced cells will be permanently switched to the stem cell state.
We notice that the pluripotent cellular state always splits up into two substates with high/low Nanog level due to spontaneous dynamical transition, i.e., all successfully reprogrammed cells have the heterogeneous stem cell state, with ~ 80% in high Nanog level and ~ 20% in low. However, the cells with high Oct4, Sox2 and low Nanog expression are usually regarded as the incomplete reprogrammed cells, as suggested by previous work that the high Nanog expression could activate the related downstream genes in the late stages of cell reprogramming process, to promote the incomplete reprogrammed cells towards the full pluripotent state [16, 17]. The time period of reprogramming usually lasts for tens of weeks [43]. If Nanog maintains a low expression level during the long induction process, the expression of other downstream genes of Nanog will remain repressed, leading to longer time to attain full pluripotency. Thus, the longtime existence of low Nanog level will increase the time of reprogramming [16, 17].
We next study the role of Nanog in the induction process by exploring four strategies previously studied in the “Seesaw” model of iPS reprogramming [33]. Specifically, the four strategies are: activating stemness genes Oct4 and Sox2; activating one stemness gene and another specifier gene group (Oct4 & ECT, or Sox2 & ME); activating one stemness gene while repressing another lineage specifier gene group (+Oct4 & ME, or +Sox2 & ECT); and equal activation of two lineage specifier genes ME and ECT. One important step among those strategies is the activation of Oct4 and Sox2, for which the feedbacks from Nanog to Oct4 and Sox2 provide additional regulations. When the activation strengths of Oct4 and Sox2 are mild (Fig. 4a & b), the highNanog cells dominate in the induction process. However, when the activation becomes stronger (Fig. 4c & d), the lowNanog intermediate state dominates. This can be intuitively understood as a result of the repression of Oct4/Sox2 on Nanog: when Oct4 and Sox2 are highly induced, Nanog will be repressed to a low level, leading to more incomplete reprogrammed cells and slower reprogramming process. We note that introducing Nanog activation in the induction process could rescue reprogramming by removing the intermediate Nanog state (Fig. 4e & f). This observation is consistent with the stable pluripotent state with highdoxycycline in D8H in [44], where Oct4 and Nanog have been rapidly upregulated to steady level, and remain after lowering the doxycycline concentration. The same mechanism can explain the two similar strategies of reprogramming involving the direct activation of Oct4 or Sox2: activating Oct4 and repressing MEs (Additional file 7: Figure S6AC), and activating both Sox2 and ECTs (Additional file 7: Figure S6DF), where additional Nanog activation could remove the lowNanog state. In the strategy of activating two lineage specifiers without Oct4 and Sox2, although the repression of Nanog due to Oct4/Sox2 overactivation may be less significant, additional Nanog activation may still benefit through removing the intermediate state (Additional file 7: Figure S6G & H).
Together, we suggest that additional Nanog activation can elevate Nanog level that is repressed by overactivated Oct4/Sox2 during induction, which can boost reprogramming efficiency by eliminating the Nanog intermediate state.
Discussion
Nanog has been recognized as a critical pluripotency gene in ES cell differentiation and iPS cell reprogramming [26, 27]. While the lowNanog state of stem cell has been speculated as a “gatekeeper” [25], the roles of the heterogeneous bimodal Nanog expression remain elusive. In this work, through analyzing a fivenode gene regulatory network, including both stemness genes Oct4, Sox2, Nanog and specifier genes MEs and ECTs, we investigated the role of Nanog in stem cell differentiation and iPS cell reprogramming. In differentiation, the lowNanog state is found to serve as an intermediate state in stem cell differentiation and to promote the differentiation process. On the other hand, the bimodal distribution of Nanog is found to have a negative effect on the iPS cell reprogramming. Additional Nanog activation is shown to prevent the reprogrammed iPS cell from trapping in the Nanog intermediate state due to overactivated Oct4/Sox2.
Through analyzing the Waddington landscape and calculating the minimum action path of differentiation, we suggest that the dynamic Nanog heterogeneity, induced by gene expression noise, can be utilized to maintain the stemness while, in the meantime, providing flexibility in cell fate decisions (Fig. 5). Differentiation drives cells in the highNanog state to overcome the barrier in the landscape. Once the cells arrive at the noiseinduced Nanog intermediate state, they become more prone to differentiate. The twostep process by going through the intermediate state allows “easier” differentiation than the onestep process bypassing the intermediate state. Our model shows that the negative feedback between Oct4Sox2 and Nanog is the key to introducing the intermediate state for cell differentiation. In principle, such intermediate state could be produced by other mechanisms. For instance, the histone modification in the chromatin region of Oct4 and Nanog is randomly reestablished during each cycle of cell division and chromosome duplication [12]. Consequently, stochastic inheritance of histone modifications may also control gene on/off switching, leading to the cellcell heterogeneity [45, 46] and potential intermediate states. Moreover, it will be interesting to develop a more comprehensive model to capture the fully dynamic processes of reprogramming by including the transcription factors (e.g. gene nodes of Klf4, Myc [47]) and specified ME/ECT lineage specifiers.
Recently, there is some controversial discussion that the genetic reporter for Nanog may cause a bifurcation in the underlying dynamics that leads to heterogeneous Nanog expression induced by the bistability [48]. Such problem of measurement in cells could induce the Nanog fluctuations. However, gene reporter strategy is generic, and may also cause the similar measurement problem of the other pluripotent genes such as Oct4 and Sox2 whose heterogeneity has not been observed in the literatures. Moreover, our present model shows the negative feedback between Nanog and Oct4Sox2 can constitute an excitable system, leading to a recurring stochastic transition between high and low Nanog states of pluripotent cells and the bimodal distribution of Nanog. Thus, the bistability required in [48] is not necessarily needed for the heterogeneous Nanog expression. Furthermore, perturbation of Nanog dynamics by the fluorenscent reporter system may not full account for the heterogeneity of Nanog expression under all culture conditions. The datasets using single cell RNASeq can support the bimodal Nanog expression for mouse embryonic stem cells cultured in serum and LIF [4, 35], where the genetic reporter is excluded.
As a strategy in regulating stem cell properties, multiple intermediate states have been discovered in epithelialmesenchymal transition (EMT) process to control cellular plasticity [49,50,51]. Similar to microRNA and Ovol2 in the EMT regulatory network, which is critical to inducing one or more intermediate (“hybrid”) states, Nanog is a key factor in the late stage of reprogramming process to promote the incomplete reprogrammed cells to attain the full pluripotent state. Identification of Nanogrelated external signals and genes will help to improve strategies for iPS cell reprogramming.
Conclusions
In summary, our modeling work has identified novel roles of Nanog in stem cell regulation. Such approach can include more biological details or other stem cell markers for more indepth exploration. Our work offers a general method for analyzing key regulatory factors controlling cell differentiation and reprogramming.
Methods
Stochastic differential equations and parameters of the fivenode model
The stochastic differential equations and parameters of the fivenode model are listed below, where [O], [S], [N], [M] and [E] represent the expression level of Oct4, Sox2, Nanog, MEs and ECTs respectively. The parameters of the model are listed in Additional file 1: Table S1.
The Oct4Sox2 protein complex will activate Nanog expression more sensitively, i.e. with a low critical value in the activation Hill function. While at the same time, Oct4Sox2 complex will repress Nanog with a high critical value in the repression Hill function. As the activation strength of Oct4Sox2 to Nanog (parameter β_{Na_OS} in the model) is much smaller than that of Nanog selfactivation (β_{ Na }), the Hillactivation term is simplified as a constant term (β_{Na_OS}) in the model. Finally, both MEs and ECTs repress the expression of Nanog. The downstream interactions of Nanog to lineage specifier genes are complicated, which are not considered in this model for simplicity. The noise in the system is added in the form of σ_{ X }[X] ∙ ξ(t), where σ_{ X } indicates the noise amplitude of each gene, and ξ(t) is the Gaussian white noise with zero mean and unit variance.
The setting of parameter values refers to the previous modeling works [2, 3], shown in Additional file 2: Figure S1A. Note that the copy numbers of Oct4 and Nanog in the model are adjusted to fit the results of mice E14IVc ES cell line cultured in serum and LIF [2], which can be adjusted and scaled for other stem cell lines. Sox2 is assumed to be symmetric to Oct4, so most parameters for Oct4 and Sox2 are set to be identical. Thus, the expression levels of both Oct4 and Sox2 are high in pluripotent cell states and low in differentiated cell states. If there is a strong initial imbalance of Oct4 / Sox2, such as high Oct4 and low Sox2, the system will lead to an imbalanced expression of lineage specifiers: high MEs and low ECTs. The high expression of MEs then represses Oct4, leading to differentiation toward the ME state, where both Oct4 and Sox2 are low in expression. On the other hand, if the initial imbalance of Oct4/Sox2 is lower than a threshold, which cannot trigger the differentiation process, the cell then returns to the pluripotent steady state, where both Oct4 and Sox2 are highly expressed. As a result, the imbalanced expression of Oct4 and Sox2 could only be a transient state, and such imbalance disappears after the cells either commit to differentiation or remain in the pluripotent state. The noise levels of each gene are adjusted to fit the previous experimental results of mice E14IVc ES cell line [2]. Note that the degradation rate 1 a.u.^{− 1} of the arbitrary timescale in the model is corresponding to 3 × 10^{−5} s^{− 1} under the real timescale. For the reestablishment time (Fig. 2b) and dwell time distribution (Fig. 2c), the timescale of the system is scaled from a.u. to day.
Simplified twodimensional Oct4Nanog model
The regulatory relationships between Oct4 and Nanog are extract out of the fivenode model (where the Sox2 is assumed totally symmetric to Oct4), and the concentration of MEs and ECTs are set up as the constant level under the pluripotent state. The stochastic differential equations are displayed below, where [O], [N] represent the level of Oct4 and Nanog respectively, and the parameters of the model are listed in Table S2. The behavior of excitable system can be illustrated on the twodimensional phase plane in terms of nullclines and vector fields (Additional file 3: Figure S2).
Parameter sensitivity analysis of the fivenode model
To test the robustness of the model, each parameter is individually increased/decreased by 20%, and the relative changes of the lowNanog distribution ratio (in the time course, blue bar), the average value of Oct4 expression level (green bar), and the average value of Nanog expression level within the highNanog state (red bar) is illustrated in Additional file 8: Figure S3. The parameters of the Oct4/Sox2 expression terms are more sensitive (the first panel from the left), and the parameters of the expression terms of Nanog, MEs and ECTs are more robust (the second, third panel from the left).
The control model excluding the dynamics of Nanog for MAM analysis
To verify the role of Nanog during the differentiation process, we set up a control model, where the dynamics and the regulatory relationships of Nanog are removed from the fivenode model. As the activation terms of Nanog to Oct4 and Sox2 (\( \frac{\left[N\right]}{\left[N\right]+{K}_{p\_ Na}} \) in the formulas 1.1 and 1.2) are the only input regulations from Nanog to the rest part of the network, the concentration value of Nanog in those two terms is set as the constant value of highly expressed steady state value of Nanog, so that the steady state values of the other four genes can remain unchanged at the same time.
The model with external induction input terms
In order to analyze the induced iPS reprogramming process, some constant input terms are added into the model. The input parameters for gene expression activation (C_{ O }, C_{ S }, C_{ N }, C_{ M }, C_{ E }) indicate the extra, constant protein production rates by exogeneous expression (by virus vector, for example) for each gene, which are set as 0 without any inductions. Note that the strengths of external induction terms are scaled (K_{ O } = K_{ S } = K_{ N } = K_{ M } = K_{ E } = 400). While the gene expression repressions due to specific miRNA molecules silencing the translation of certain mRNA molecules are modeled by the multiplier terms before the entire regulated expression rate terms for each gene (I_{ O }, I_{ S }, I_{ N }, I_{ M }, I_{ E }), which is set as 1 without any inductions. All other parameters remained unchanged.
The sampling method of the potential landscape
The landscape functions in the nonequilibrium systems have been widely used to model the gene switching and cell fate transition [52, 53], which serves as an analogy of the energy function in the equilibrium systems. However, several different landscape functions have been proposed by biophysicists and applied mathematicians [38,39,40]. The mathematical relations and differences between different landscapes are elucidated in [54]. Here we choose the simplest one, using the negative logarithm of the probability distribution in the phase space under the steady state,− log(P_{ ss }), as the definition of the potential landscape. Note that the gradient of − log(P_{ ss }) is no longer corresponding to the vector field direction of the dynamical system. And the landscape of− log(P_{ ss }) is actually, to some extent, limited as an intuitive sketch of cellular state distribution in the phase space. The probability distributions (P_{ ss }) in the phase space of the figures above are produced by simulating the system for more than 10^{7} time a.u., and the probability density of the trajectories are projected into the threedimensional space, with the x, y, z axis as \( \sqrt{\left[ Oct4\right]\bullet \left[ Sox2\right]} \), [Nanog] and [MEs] − [ECTs] (e.g. Fig. 1b); or twodimensional space in other figures, with the x, y axis as \( \sqrt{\left[ Oct4\right]\bullet \left[ Sox2\right]} \) and [Nanog] (e.g. Fig. 1c). The color scale of the potential landscape measures the energy value, indicating the probability density for the cell state to appear in that certain region.
The method of minimum action path
The WentzellFreidlin theory of large deviation gives an estimate of the probability of the paths in terms of an action functional. A key result of this theory is that the most probable path minimizes the action functional associated with the random dynamical system, i.e., the most probable path is the Minimum Action Path.
In order to find the MAP between two steady states, we follow the minimum action method in [42] to compute the numerical solutions with the time interval [0, 100]. We apply the BFGS algorithm for numerical optimization.
Abbreviations
 ECTs:

Ectoderm genes
 ES:

Embryonic stem
 iPS:

Induced pluripotent stem
 MAP:

Minimum Action Path
 MEs:

Mesendoderm genes
 WT:

Wildtype
References
 1.
Graf T, Stadtfeld M. Heterogeneity of embryonic and adult stem cells. Cell Stem Cell. 2008;3(5):480–3.
 2.
Hayashi K, Lopes SM, Tang F, Surani MA. Dynamic equilibrium and heterogeneity of mouse pluripotent stem cells with distinct functional and epigenetic states. Cell Stem Cell. 2008;3(4):391–401.
 3.
Martinez Arias A, Brickman JM. Gene expression heterogeneities in embryonic stem cell populations: origin and function. Curr Opin Cell Biol. 2011;23(6):650–6.
 4.
Kumar RM, Cahan P, Shalek AK, Satija R, DaleyKeyser AJ, Li H, Zhang J, Pardee K, Gennert D, Trombetta JJ, et al. Deconstructing transcriptional heterogeneity in pluripotent stem cells. Nature. 2014;516(7529):56–61.
 5.
Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S. Transcriptomewide noise controls lineage choice in mammalian progenitor cells. Nature. 2008;453(7194):544–7.
 6.
Pina C, Fugazza C, Tipping AJ, Brown J, Soneji S, Teles J, Peterson C, Enver T. Inferring rules of lineage commitment in haematopoiesis. Nat Cell Biol. 2012;14(3):287–94.
 7.
Bonzanni N, Garg A, Feenstra KA, Schutte J, Kinston S, MirandaSaavedra D, Heringa J, Xenarios I, Gottgens B. Hardwired heterogeneity in blood stem cells revealed using a dynamic regulatory network model. Bioinformatics. 2013;29(13):i80–8.
 8.
Graf T, Enver T. Forcing cells to change lineages. Nature. 2009;462(7273):587–94.
 9.
Eldar A, Elowitz MB. Functional roles for noise in genetic circuits. Nature. 2010;467(7312):167–73.
 10.
Balazsi G, van Oudenaarden A, Collins JJ. Cellular decision making and biological noise: from microbes to mammals. Cell. 2011;144(6):910–25.
 11.
Sokolik C, Liu Y, Bauer D, McPherson J, Broeker M, Heimberg G, Qi LS, Sivak DA, Thomson M. Transcription factor competition allows embryonic stem cells to distinguish authentic signals from noise. Cell Syst. 2015;1(2):117–29.
 12.
Sasai M, Kawabata Y, Makishi K, Itoh K, Terada TP. Time scales in epigenetic dynamics and phenotypic heterogeneity of embryonic stem cells. PLoS Comput Biol. 2013;9(12):e1003380.
 13.
Mitsui K, Tokuzawa Y, Itoh H, Segawa K, Murakami M, Takahashi K, Maruyama M, Maeda M, Yamanaka S. The homeoprotein Nanog is required for maintenance of pluripotency in mouse epiblast and ES cells. Cell. 2003;113(5):631–42.
 14.
Chambers I, Colby D, Robertson M, Nichols J, Lee S, Tweedie S, Smith A. Functional expression cloning of Nanog, a pluripotency sustaining factor in embryonic stem cells. Cell. 2003;113(5):643–55.
 15.
Ivanova N, Dobrin R, Lu R, Kotenko I, Levorse J, DeCoste C, Schafer X, Lun Y, Lemischka IR. Dissecting selfrenewal in stem cells with RNA interference. Nature. 2006;442(7102):533–8.
 16.
Silva J, Nichols J, Theunissen TW, Guo G, van Oosten AL, Barrandon O, Wray J, Yamanaka S, Chambers I, Smith A. Nanog is the gateway to the pluripotent ground state. Cell. 2009;138(4):722–37.
 17.
Schwarz BA, BarNur O, Silva JC, Hochedlinger K. Nanog is dispensable for the generation of induced pluripotent stem cells. Curr Biol. 2014;24(3):347–50.
 18.
Boyer LA, Lee TI, Cole MF, Johnstone SE, Levine SS, Zucker JP, Guenther MG, Kumar RM, Murray HL, Jenner RG, et al. Core transcriptional regulatory circuitry in human embryonic stem cells. Cell. 2005;122(6):947–56.
 19.
Chickarmane V, Troein C, Nuber UA, Sauro HM, Peterson C. Transcriptional dynamics of the embryonic stem cell switch. PLoS Comput Biol. 2006;2(9):1080–92.
 20.
Chang R, Shoemaker R, Wang W. Systematic search for recipes to generate induced pluripotent stem cells. PLoS Comput Biol. 2011;7(12):e1002300.
 21.
Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS Comput Biol. 2013;9(8):e1003165.
 22.
Chambers I, Silva J, Colby D, Nichols J, Nijmeijer B, Robertson M, Vrana J, Jones K, Grotewold L, Smith A. Nanog safeguards pluripotency and mediates germline development. Nature. 2007;450(7173):1230–4.
 23.
Singh AM, Hamazaki T, Hankowski KE, Terada N. A heterogeneous expression pattern for Nanog in embryonic stem cells. Stem Cells. 2007;25(10):2534–42.
 24.
Kalmar T, Lim C, Hayward P, MunozDescalzo S, Nichols J, GarciaOjalvo J, Martinez Arias A. Regulated fluctuations in nanog expression mediate cell fate decisions in embryonic stem cells. PLoS Biol. 2009;7(7):e1000149.
 25.
Glauche I, Herberg M, Roeder I. Nanog variability and pluripotency regulation of embryonic stem cells  insights from a mathematical model analysis. PLoS One. 2010;5(6):e11238.
 26.
Niwa H. How is pluripotency determined and maintained? Development. 2007;134(4):635–46.
 27.
Silva J, Smith A. Capturing pluripotency. Cell. 2008;132(4):532–6.
 28.
Suzuki N, Furusawa C, Kaneko K. Oscillatory protein expression dynamics endows stem cells with robust differentiation potential. PLoS One. 2011;6(11):e27232.
 29.
Rabajante JF, Babierra AL. Branching and oscillations in the epigenetic landscape of cellfate determination. Prog Biophys Mol Biol. 2015;117(2–3):240–9.
 30.
Lindner B. Effects of noise in excitable systems. Phys Rep. 2004;392(6):321–424.
 31.
Suel GM, GarciaOjalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006;440(7083):545–50.
 32.
Maamar H, Raj A, Dubnau D. Noise in gene expression determines cell fate in Bacillus subtilis. Science. 2007;317(5837):526–9.
 33.
Shu J, Wu C, Wu Y, Li Z, Shao S, Zhao W, Tang X, Yang H, Shen L, Zuo X, et al. Induction of pluripotency in mouse somatic cells with lineage specifiers. Cell. 2013;153(5):963–75.
 34.
Loh YH, Wu Q, Chew JL, Vega VB, Zhang W, Chen X, Bourque G, George J, Leong B, Liu J, et al. The Oct4 and Nanog transcription network regulates pluripotency in mouse embryonic stem cells. Nat Genet. 2006;38(4):431–40.
 35.
Kolodziejczyk AA, Kim JK, Tsang JCH, Ilicic T, Henriksson J, Natarajan KN, Tuck AC, Gao XF, Buhler M, Liu PT, et al. Single cell RNAsequencing of pluripotent states unlocks modular transcriptional variation. Cell Stem Cell. 2015;17(4):471–85.
 36.
Freidlin MI, Wentzell AD: Random perturbations of dynamical systems. 2nd edition. 1998.
 37.
Zhang L, Radtke K, Zheng L, Cai AQ, Schilling TF, Nie Q. Noise drives sharpening of gene expression boundaries in the zebrafish hindbrain. Mol Syst Biol. 2012;8:613.
 38.
Lv C, Li XG, Li FT, Li TJ. Energy landscape reveals that the budding yeast cell cycle is a robust and adaptive multistage process. PLoS Comput Biol. 2015;11(3):e1004156.
 39.
Wang J, Xu L, Wang EK. Potential landscape and flux framework of nonequilibrium networks: robustness, dissipation, and coherence of biochemical oscillations. P Natl Acad Sci USA. 2008;105(34):12271–6.
 40.
Ao P. Potential in stochastic differential equations: novel construction. J Phys aMath Gen. 2004;37(3):L25–30.
 41.
Xing JH. Mapping between dissipative and Hamiltonian systems. J Phys aMath Theor. 2010;43(37):375003.
 42.
Weinan E, Ren W, VandenEijnden E. Minimum action method for the study of rare events. Commun Pure Appl Math. 2004;57(5):637–56.
 43.
Hanna J, Saha K, Pando B, van Zon J, Lengner CJ, Creyghton MP, van Oudenaarden A, Jaenisch R. Direct cell reprogramming is a stochastic process amenable to acceleration. Nature. 2009;462(7273):595–601.
 44.
Hussein SMI, Puri MC, Tonge PD, Benevento M, Corso AJ, Clancy JL, Mosbergen R, Li M, Lee DS, Cloonan N, et al. Genomewide characterization of the routes to pluripotency. Nature. 2014;516(7530):198.
 45.
Rando OJ, Verstrepen KJ. Timescales of genetic and epigenetic inheritance. Cell. 2007;128(4):655–68.
 46.
Dodd IB, Micheelsen MA, Sneppen K, Thon G. Theoretical analysis of epigenetic cell memory by nucleosome modification. Cell. 2007;129(4):813–22.
 47.
Polo JM, Anderssen E, Walsh RM, Schwarz BA, Nefzger CM, Lim SM, Borkent M, Apostolou E, Alaei S, Cloutier J, et al. A molecular roadmap of reprogramming somatic cells into iPS cells. Cell. 2012;151(7):1617–32.
 48.
Smith RCG, Stumpf PS, Ridden SJ, Sim A, Filippi S, Harrington HA, MacArthur B. Nanog fluctuations in embryonic stem cells highlight the problem of measurement in cell Biology. Biophys J. 2017;112(12):2641–52.
 49.
Tian XJ, Zhang H, Xing J. Coupled reversible and irreversible bistable switches underlying TGFbetainduced epithelial to mesenchymal transition. Biophys J. 2013;105(4):1079–89.
 50.
Lu M, Jolly MK, Levine H, Onuchic JN, BenJacob E. MicroRNAbased regulation of epithelialhybridmesenchymal fate determination. Proc Natl Acad Sci U S A. 2013;110(45):18144–9.
 51.
Hong T, Watanabe K, Ta CH, VillarrealPonce A, Nie Q, Dai X. An Ovol2Zeb1 mutual inhibitory circuit governs bidirectional and multistep transition between epithelial and mesenchymal states. PLoS Comput Biol. 2015;11(11):e1004569.
 52.
Zhou JX, Aliyu MDS, Aurell E, Huang S. Quasipotential landscape in complex multistable systems. J R Soc Interface. 2012;9(77):3539–53.
 53.
Ge H, Qian H. Landscapes of nongradient dynamics without detailed balance: stable limit cycles and multiple attractors. Chaos. 2012;22(2):023140.
 54.
Zhou PJ, Li TJ. Construction of the landscape for multistable systems: potential landscape, quasipotential, Atype integral and beyond. J Chem Phys. 2016;144(9):094109.
Acknowledgements
The authors are grateful to Tiejun Li for helpful discussions.
Funding
LZ was partially supported by the National Natural Science Foundation of China 11622102, 91430217, and 11421110001. CT was partially supported by the National Natural Science Foundation of China 91430217 and the Chinese Ministry of Science and Technology 2015CB910300. QN was partially supported by NSF grant DMS1562176, and NIH grants P50GM076516, R01GM107264, and R01NS095355.
Availability of data and materials
The codes for generating the data are available from the corresponding author on reasonable request.
Author information
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1:
Table S1. Parameters used in Eq. (1) for the fivenode model. (DOCX 50 kb)
Additional file 2:
Figure S1. Typical temporal trajectories of stochastic gene expressions at the ME differentiated cell state. ME state is a stable state, and the noisedriven transition from differentiated states (low Oct4, Sox2 and Nanog) to pluripotent states (high Oct4 and Sox2, low MEs and ECTs) cannot occur spontaneously. (TIFF 1916 kb)
Additional file 3:
Figure S2. The simplified twodimensional Oct4Nanog model on the phase plate and the distribution of Oct4. (A)The nullclines and the vector field of the simplified twodimensional Oct4Nanog model on the phase plate. A typical trajectory is illustrated to indicate the excitable mechanism of the model. (d[Oct4]/dt = 0: Red line; d[Nanog]/dt = 0: Blue line.) (B) Distributions of Sox2 level within simulated cell population (N = 10,000). (PDF 102 kb)
Additional file 4:
Table S2. Parameters used in Eq. (2) for the simplified Oct4Nanog model. (DOCX 42 kb)
Additional file 5:
Figure S4. The MAPs of the differentiation process with two different initial paths in the WT model. The MAPs (white curves) starting from the pluripotent state (the green point) to the ME differentiated state (the blue point) are insensitive to different initial conditions (purple curves): (A) a smooth curve passing by the lowNanog state; (B) a smooth curve far from lowNanog state. (PDF 614 kb)
Additional file 6:
Figure S5. The MAP of the reprogramming process in the WT model. The MAP (white curve) starting from the ME differentiated state (the blue point) to the pluripotent state (the green point) is different from that of differentiation process (Fig. 3A). The green dotted line is the ODE trajectory to compare with the MAP. (PDF 3338 kb)
Additional file 7:
Figure S6. Three different strategies of reprogramming demonstrate additional Nanog activation is necessary to maintain the high Nanog level and promote the efficient cell reprogramming. (AC) Strategy by of activating Oct4 and repressing MEs. (A) C_{0} = I_{ m } = 0.3; (B) C_{0} = I_{ m } = 0.5; (C) C_{0} = I_{ m } = C_{ n } = 0.5; (DF) Strategy of activating Sox2 and ECTs. (D) C_{ m } = 0.3, C_{ s } = 0.06; (E) C_{ m } = 0.5, C_{ S } = 0.1; (F) C_{ m } = 0.5, C_{ S } = 0.1, C_{ n } = 0.5; (GH) Strategy of activating MEs and ECTs. (G) C_{ m } = C_{ e } = 0.3; (H) C_{ m } = C_{ e } = C_{ n } = 0.3. (PDF 2322 kb)
Additional file 8:
Figure S3. Parameter sensitivity analysis for the model. Illustration of the relative changes of the lowNanog distribution ratio (blue bar), the average Oct4 level (green bar), and the average Nanog level of highNanog population (red bar). (TIFF 699 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Gene network
 Stem cells
 Cell differentiation
 iPS cell reprogramming
 Intermediate cellular state
 Nanog