Pluripotency gene network dynamics: System views from parametric analysis

Multiple experimental data demonstrated that the core gene network orchestrating self-renewal and differentiation of mouse embryonic stem cells involves activity of Oct4, Sox2 and Nanog genes by means of a number of positive feedback loops among them. However, recent studies indicated that the architecture of the core gene network should also incorporate negative Nanog autoregulation and might not include positive feedbacks from Nanog to Oct4 and Sox2. Thorough parametric analysis of the mathematical model based on this revisited core regulatory circuit identified that there are substantial changes in model dynamics occurred depending on the strength of Oct4 and Sox2 activation and molecular complexity of Nanog autorepression. The analysis showed the existence of four dynamical domains with different numbers of stable and unstable steady states. We hypothesize that these domains can constitute the checkpoints in a developmental progression from naïve to primed pluripotency and vice versa. During this transition, parametric conditions exist, which generate an oscillatory behavior of the system explaining heterogeneity in expression of pluripotent and differentiation factors in serum ESC cultures. Eventually, simulations showed that addition of positive feedbacks from Nanog to Oct4 and Sox2 leads mainly to increase of the parametric space for the naïve ESC state, in which pluripotency factors are strongly expressed while differentiation ones are repressed.


Introduction
Pluripotency is a temporal state in embryogenesis artificially maintained in embryonic stem cells (ESCs) by specific components in the medium providing ESC self-renewal and inhibition of signaling pathways driving differentiation (reviewed in [1,2,3]). In mouse pluripotent cells there are three types of media used for this, 2i/LIF, serum/LIF and FGF2/Activin providing for naïve, formative and primed pluripotent states, respectively. These states differ in the expression of pluripotency and differentiation genes with naïve and primed states a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 biased correspondingly to pluripotency and differentiation. Somatic cells may be also driven to pluripotency on specific media, and obtained induced pluripotent stem cells (iPSCs) are very similar to ESCs in their molecular and functional traits [4], (reviewed in [5,6]).
Along with whole-genome technologies, high-resolution imaging and other experimental approaches, mathematical modeling becomes one of the tools for understanding of molecular mechanisms behind cell pluripotency and differentiation (reviewed in [7,8]). For simulations, a molecular genetic process represents as a regulatory gene network handling of input signals. Multiple experimental studies demonstrate that Oct4, Sox2 and Nanog (OSN) are the core factors in pluripotency gene network involved in induction, maintenance and loss of pluripotency (reviewed in [2,9]).
Chickarmane and coauthors developed the first model of this core circuit and follow-up series of extended models to investigate a landscape of possible embryonic stem cell (ESC) states [10,11,12]. Already analysis of the first model, where Oct4, Sox2 and Nanog (OSN) regulated by environmental signals activate each other, indicated that the positive-feedback loops in the OSN circuit give rise to bistability, which corresponds to existence of two stable cell states (self-renewal/pluripotency and differentiation) toggle-switched by external signals [10]. Mutual regulation within OSN core were withdrawn from whole genome experiments on identification of OSN targets in both human ESCs (hESCs) and mESCs [13,14]. Positive regulation of Nanog, Sox2 and Oct4 by Oct4/Sox2 heterodimer was also shown [15,16,17,18,19].
Meantime, regulatory links from Nanog to Oct4 and Sox2 appeared to be more complex. Nanog depletion resulted in Oct4 and Sox2 down regulation but Nanog overexpression increased Oct4 level to at least 150% whereas the Sox2 level remained unchanged [14]. Oct4 down regulation in response to Nanog knockdown was confirmed in [20], but Nanog overexpression in these experiments did not increase Oct4 concentration beyond the steady-state level. Furthermore, Nanog expression is low in cells with elevated Oct4 level and high in cells with low Oct4 [20,21]. In addition, it was reported that Nanog has minimal influence on Oct4 and Sox2 expression [22]. Besides OS activation by Nanog under the question, it was also shown that the main regulator of Nanog transcription is its autorepression. Thus, nowadays, the architecture of the OSN network should be refined to omit Nanog influence on OS expression and incorporate not positive but negative Nanog autoregulation. From this, the OSN gene network becomes a system where one of the genes is repressed by its own product. Understanding of such system behavior is not complete without consideration the transcriptional and translational delays in its functioning [23].
Nanog really plays a special role in OSN triad activity. Nanog overexpression in mESCs maintains pluripotency independently of LIF signaling [24,25], whereas Oct4 and Sox2 overexpression drives mESCs to mesendodermal and neural ectodermal differentiation, respectively [26,27,28]. However, Nanog was not included in the most efficient "cocktail" to induce pluripotency in mouse somatic cells [29] containing Oct4 and Sox2 along with Klf4 and cMyc. Nevertheless, Nanog overexpression accelerates reprogramming of somatic cells to a pluripotent state [30,31] and activation of the endogenous Nanog and Oct4 is one of the final events in this process [32]. Nanog is necessary for pre-iPSCs (dedifferentiated intermediates) to acquire ground state pluripotency [33]. Whereas Nanog autorepression is mediated by association of Nanog and Zfp281 proteins with NuRD (Nucleosome Remodeling and Deacetylase) complex and Zfp281 inhibits transition of pluripotent iPSCs to ground state by restricting activation of endogenous Nanog [34].
Transient up and downregulation of endogenous Nanog was recorded under serum/LIF conditions in individual mESCs predisposing them toward ground stage and differentiation, respectively [35,36,37,38]. Depending on culture conditions the percentage of Nanog-low cells may vary from 5% till 35% [36,39]. Mathematical modeling on Nanog heterogeneity demonstrated that it might arise from transcriptional noise [36,40]. Moreover, noise may be sufficient to trigger reactivation of the core pluripotency switch and reprogramming to a pluripotent state [41]. Nevertheless, Wu et al. investigated sensitivity and robustness of the Nanog gene network by means of a stochastic model and identified that the system dynamics is sensitive to positive regulation from the Oct4-Sox2 complex [42]. Before Nanog autorepression was shown experimentally by [22], in silico simulations demonstrated that Nanog autorepression even indirect provides for sustained oscillations in the system without noise [36,40]. Comparing sustained oscillations and stochastic fluctuations in an agent-based model suggested that the noise-driven model more consistently explain Nanog expression dynamics in mESC populations. Coupling in simulations the OSN network to extracellular signals provided evidence that stochastic autocrine feedback loops might also generate fluctuations in Nanog expression [43]. Meanwhile, a single-molecule fluorescent in situ hybridization used to detect Nanog mRNA in mESCs, cultured in serum/LIF or 2i/LIF media, offered evidence for the stochastic nature of Nanog expression in pluripotent mESCs, independently of the culture conditions implying that NANOG fluctuations are not dependent on autocrine ERK signaling mediated by FGF4 [44].
However, the mutual regulatory role of Nanog autorepression and the strength of the signals inducing Oct4 and Sox2 expression in self-renewal and differentiation of mESC as well as an independence/dependence of Oct4/Sox2 expression from Nanog activation have not been studied in silico yet. To sum it up so far, we propose the strength of Oct4 and Sox2 activation and complexity of the Nanog autorepression feedback as key forces governing pluripotency-differentiation switch (Fig 1B). The core transcriptional network of the factors orchestrating the pluripotency and differentiation genes (suggested by [10]). External A + and Bsignals activate and repress expression of Oct4, Sox2 and Nanog genes, correspondingly. Oct4 and Sox2 form a heterodimer, Oct4/Sox2, which positively regulates Oct4, Sox2 and Nanog expression. Nanog directly induces Oct4, Sox2 and its own expression. Oct4/Sox2 heterodimer and Nanog positively regulate pluripotency genes and repress differentiation genes. B: The revised core gene network suggested in this paper, in which transcription and translation processes were added; external signal B-is removed and positive signal A+ activates transcription of Oct4 Sox2 genes. Nanog represses its own transcription and does not influence on Oct4 and Sox2 expression.
Herein, we present results on modeling of the mouse ESC core regulatory circuit revisited according to the recent experimental data. The first model describing bistable behavior of the core circuit comprises very simple gene network-only core genes Oct4, Sox2 and Nanog enhancing activity of each other (Fig 1A). Here we revise this model [10] by adding new data on Nanog autoinhibition and taking off Nanog enhancement of Oct4 and Sox2 activity due to controversy records on this topic [22]. As shown in Fig 1B, unlike the initial model we also take into account compartmental localization of each gene transcription in the cell nucleus and maturation of functionally active proteins encoded by core pluripotency genes in the cytoplasm in order to developed model more precisely describes biomolecular events during differentiation and reprogramming processes.
Based on numerical simulations and the parameter continuation method [45] we identified patterns of the ESC states, which existence mainly depends on two following model parameters: the degree of Oct4 and Sox2 activation and molecular complexity of Nanog autorepression. We hypothesize that identified dynamical domains might correspond to naive and primed ES cells and intermediate states between them. To address heterogeneity in Nanog expression, we argue that Nanog autorepression and updated topology of the core network can be internal trigger for oscillations in Nanog expression levels. Eventually, we investigate extended model accounting for positive feedback loops, in which Nanog is an inducer of Oct4 and Sox2 transcription. Both model modifications have bistable, switch-like behavior in the same range of parameter values. Moreover, considering the additional positive feedback loops did not result in emergence of the new type of the model behavior but just enhanced the ESC naïve state traits.
It is worth to note that behavior of such complicated system as regulatory machinery of pluripotency\differentiation switch should not be constrained by consideration of only core gene network and extension of transcriptional networks incorporating epigenetic level of pluripotency regulation as well as regulatory circuits mediated via ncRNAs is needed [2,5,9]. However, obtained in silico results for the core regulatory circuit can not only provide the explanation and insight into dynamical behavior of the studied biological system, but serve as the basis or starting point to boost follow-up experimental-theoretical investigations.

The model structure
To study regulatory mechanisms of the mESC maintenance and differentiation, the core circuit structure has been revisited. The previously developed model for interaction of core factors in mESC [10] as well as additional experimental data on Nanog expression regulation and its interplay with Oct4 and Sox2 [22] served as the basis for the new dynamical model (Fig 1; Materials and Methods). We have also taken into account occurrence of OSN transcription and following translation processes in separate cell compartments (nucleus and cytoplasm). It led to an addition of intracellular mRNA and protein transport simultaneously with the emergence of such parameters as basal transcription rates of the OSN genes, activation or inhibition thresholds, efficiencies of transcription factor binding to a promoter site, as an example. Values of the parameters and corresponding references are represented in the Materials and Methods (Table 1). Degradation rates for OSN mRNAs and proteins were taken from published data [27,46,47].
With the updates to the core network above, we added a set of differential equations (Materials and Methods) using generalized Hill functions [48]. The generalized Hill functions allow to elaborate both mathematical model construction and quantitative description of biological

Threshold in OS activation defining differentiation-pluripotency transition
OSN targets are subdivided into pluripotency target genes (PTGs) and differentiation target genes (DTGs) encoding for pluripotency and differentiation factors, respectively. To define a cell state, we considered the ratio between values of w 1 and w 2 denoting concentrations of proteins, encoding by PTGs and DTGs, respectively, at a stable steady state. If w 1 is considerably more than w 2 , the OSN gene network is in pluripotency state and if the ratio is opposite, it is in differentiation state. It is notable that unstable solution indicated by asterisks in Fig 2 and further, cannot be used for cell state determination and is referred only as unstable. The induction of pluripotency that corresponds to the switch from w 1 < w 2 to w 1 > w 2 is caused by A signal activating Oct4 and Sox2. The latter simulates either OS ectopic expression via transgenic factors or their activation by any other factors inducing pluripotency. The value of the parameter A equals to zero in differentiated cells before induction of pluripotency and during this process we increase it till 0.4. To study outcomes of above-mentioned structural reorganization of the    Nanog negative regulation of Nanog expression and A , the parameter determining a value of the signal A activating Oct4 and Sox2 expression. To consider a high-level complexity of Nanog autorepression, we initially constructed curves of w 1 and w 2 dependence on A values varying from 0 till 0.4 and the parameter h equaled to 6. As it follows from Fig 2A (initial condition-differentiated cells) and Fig 2B (initial condition-pluripotent cells), there are two following types of the steady state's dependencies on the parameter A, (i) stable differentiation, the steady state in the range 0 A < AÃ = 0.277 with turning point upon A = AÃ (Fig 2A) and (ii) stable pluripotency, the steady state at all A ! 0, including the range 0 A < AÃ (Fig 2B). Thus, depending on initial steady state values and at h = 6 generally there are three states in the range 0 A < AÃ, two steady states (pluripotency and differentiation) and one unstable state (differentiation/pluripotency), two of them starting from differentiated (Fig 2A) and one from pluripotent cells (Fig 2B). Whereas only one state exists when A > AÃ and this is the steady state from pluripotent condition. The existence of the unstable states (Fig 2A) coincide with the existence of the early phase of stochastic gene expression during induction of pluripotency [51].
If the parameter A value is higher than the threshold value AÃ = 0.277, pluripotency state was observed independently on initial values of the model variables. This simulation confirms experimental observations that once ground pluripotency is established, reduced Oct4 level in cells growing on 2i/LIF medium did not lead to loss of the pluripotent state [21,52]. These Oct4-low cells have nondynamic unimodal Nanog level [21] simulated in our model by the fixed value of the Hill coefficient for Nanog autorepression (h = 6). To consider the system behavior in more wide range of the h parameter we investigated the stability of stationary solutions under different values of the Hill coefficient for Nanog autorepression.

The strength of Oct4 and Sox2 activation and non-linearity of Nanog autorepression determine the choice between pluripotency and differentiation states and their stability
The investigation of the system stability showed that there is a domain defined by h values above threshold, hÃ = 8.68, where system behavior is practically independent from A value A ! 0. If we consider the range A ! 0, h > hÃ all states of the system will be unstable and represented only by an oscillatory mode.
Describing the dependence of steady state stability on parameters A and h, we found four following domains on the plane (A, h) separated by the lines A = AÃ = 0.277 and h = hÃ = 8.68 (Fig 3): Note that the minimal value of the Hill coefficient, which equals 2, is explained by the fact that Nanog homodimer is the core protein complex in ESCs [53,54]. It was experimentally demonstrated, that Nanog protein dimerization is vital for stem cell self-renewal and pluripotency. Higher level of the parameter value can be interpreted by means of accounting for larger protein complexes with Nanog participation and/or extended gene regulatory circuits for negative feedbacks regulating Nanog expression [55,56,57]. As can be seen from Fig 2, D 4 domain contains three states. The stationary solution diagrams of these states reproduce analogous diagrams and stability properties represented in Fig 2. It means that two solutions in the domain are stable and correspond to differentiation and pluripotency, while the third solution is unstable and tend to the differentiation according to the model analysis (Fig 2A). There is the single and stable stationary solution in D 1 domain, which corresponds to the pluripotent cell and illustrated at h = 6, Fig 2B. Stationary solutions in D 2 and D 3 domains comply with the stationary solutions in D 1 and D 4 domains, correspondingly, but the qualitative difference is that all stationary solutions are unstable in D 2 and D 3 domains. Model solutions in D 2 and D 3 domains exhibit oscillations for any initial conditions. Altogether D 2 and D 3 domains may constitute the formative pluripotency suggested by Smith (2017) as intermediate phase between naïve and primed pluripotency [1]. Thus, the formative phase might be associated with increased complexity of Nanog repression simulated via h in our model. Coinciding with this Nanog is absent from immediate post-implantation epiblast [58,59] corresponding to the formative phase of ESCs [1], Nanog downregulation is necessary for initiation of lineage specification [27,35] and Nanog overexpressing ESCs resist differentiation to EpiSCs [60].
It is noteworthy that the number of domains and their properties match to observed ESC developmental progression from naive to primed pluripotency (reviewed in [61]). It was also demonstrated that both bFGF/Activin A and L-proline induce naïve to prime transition, but after L-proline treatment ESCs mainly reached a fully reversible early phase of this transition [62]. Naïve to prime transition is reversible until complete dissipation of ground state factors occurs [59,61]. This denotes reaching the transitional stage at which neither ground state factors (e.g. Nanog, Essrb and Tcp3l1) not lineage markers (e.g. Bry, Sox17 and Brn2) are expressed, whereas Oct4 and Sox2 are present. Cells at the transitional stage become competent for lineage specification. Nanog is high at the ESC ground state [21,59]. Nanog ; D 4 domain contains three states, from which two (pluripotency and differentiation) are stable and one (transition between these states is unstable (according to Fig 2A). Domains (a) predicted by the model and (b) their correspondence to developmental progression of ESCs from the naïve pluripotency (the ground state) to lineage commitment according to [61]. The initial phase of exit from the ground state is asynchronous in the cell population and reversible until the complete dissipation of naïve state factors (reviewed in [59,61,62]). Cells reaching a transitional point after 2i withdrawal are competent for lineage specification and characterized by the absence of both groups, naïve factors and lineage markers. The late phase of pluripotency (primed pluripotency) is characterized by nascent expression of lineage specification factors. The "clock model" was proposed as a route of consistent transitions with the dual mechanism of hour hand movement depending on the initial cell state: pluripotent (counter-clockwise movement of black solid arrows) or differentiated (clockwise movement of black solid arrows). Red arrows, in turn, reflects directions from naïve to reversetransition-primed stages (developmental progression during differentiation) or from primed to transition-reverse-naïve states (developmental progression during reprogramming into pluripotent state), while dotted black arrows were added to underline intermediate reverse and transitional states to which domains D 2 and D 3 correspond, respectively.
Revealed in silico boundaries of the parameter values that constrain functional ESC states and transitions between them, in turn, can be represented by "clock model" with the dual mechanism of hour hand movement depending on the initial cell state: pluripotent or differentiated and the road map depending on h value (Fig 3A and 3B). At h = 6 corresponding approximately intermediate Nanog level (Fig 2B) we fail to reach differentiation state by decreasing Oct4 and Sox2. This again emphasizes the key role of Nanog downregulation in differentiation suggesting its occurrence via increasing complexity (nonlinearity) of Nanog negative direct and indirect feedbacks. After increase in the Nanog level by decreasing its downregulation feedbacks, the resulting diagrams of stationary solutions at h = 2 (Fig 2C)  Theunissen and Silva (2012) suggested that the key role of high endogenous Nanog in reprogramming might closely resembles its role in establishment of the naive pluripotent epiblast in early mouse embryogenesis [63]. At least initially, a weak expression of Nanog is observed and its level is variable between blastomers [64], which correspond to D 2 and D 3 domains in our diagrams. This is followed by general upregulation of Nanog expression with further differentiation of primitive endoderm and preimplantation epiblast. The latter one corresponds to the ground ESC state (reviewed in [65]).
Moreover, the graph in Fig 4 indicates that intersection point A 0 = 0.173 at h = 2 raises within the bounds 0.16 A 0 < 0.277 simultaneously with the increase of the parameter h from 2 until 3.5. But starting from h = 3.5 and higher A 0 value keeps constant and equals to A 0 = AÃ = 0.277 making in this area the boundary condition between differentiation and pluripotent states not dependent on non-linearity (complexity) of the Nanog negative feedback.
Whereas in the region of weak activation of OS genes and low values of the Hill coefficient, which reflects low nonlinear effect of Nanog on its own expression, the dependency curve grows gradually. Thus, a weak activation of OS genes and low values of the Hill coefficient is sufficient to pre-iPSC/iPSC transition upon the pluripotency induction, while a high-order nonlinearity of Nanog repression (more complex, involving more number of regulatory interactions and factors) needs the higher value of OS expression. Increase in OS expression is limited and at the highest value is independent on values of Hill coefficient (Fig 4).
This simulation fits several experimental observations. Reprogramming occurs gradually and increase in endogenous expression of core genes also occurs gradually starting from very low doses [32,66]. The observation corresponds to the graph part, where value of the Hill coefficient increases simultaneously with A 0 value. Further increase of h parameter results in repression of Nanog level after a certain threshold, while an activation of Oct4 and Sox2 transcription reaches the own threshold. It is known, that ESC self-renewal requires Oct4 and Sox2 maintaining within narrow limits exceeding which lead to ESC differentiation [26,67]. Our model confirmed the existence of these thresholds and showed that during reprogramming increase in OS levels is accompanied by increase in complexity of Nanog autoregulation. Whereas on serum/LIF media Nanog-positive and Nanog-negative mESCs were recorded, under 2i/LIF conditions where main mESC signaling via MAPK/ERK1/2 and GSK3 is strongly inhibited, the mESCs uniformly express pluripotency markers including Nanog [68,69]. Nanog autoinhibition persists in 2i/LIF [22], but Nanog-negative cells were not recorded. It is obvious that signaling pathways in serum/LIF are more complex than in 2i/LIF. It allows to suggest that regulation of Nanog autoinhibition in serum/LIF is more complex and order of the nonlinearity is higher. Strikingly, the upper boundary of the parameter range (h = 6) quantitatively corresponds the value used in Miyamoto et al. (2015) study, in which they examined the gene expression dynamics model with epigenetic feedback regulation to show that differentiation with the loss of pluripotency progresses from the embryonic stem cell state with oscillatory expression of pluripotent genes [70]. Apparently, the high value of Hill parameter in both models reflects the complex structural organization of Nanog expression regulation.
Decreasing Oct4 to the level present in Oct4+/-mESC shifts Nanog to the higher level detected in wild type mESCs due to reducing the proportion of Nanog-low or negative cells [21]. The increase in Oct4 level results in establishing of Nanog heterogeneity, in other words increase in Oct4 activation is accompanied by increase in non-linearity of Nanog autoinhibition. The 2i/LIF conditions may correspond to pluripotent stem cells from the lower part of the curve in Fig 4, whereas serum/LIF corresponds to its upper part. In other words, in mouse ES cells under 2i/LIF conditions the values of both A 0 and h are lower than in serum/LIF medium. Pluripotency gene network dynamics: System views from parametric analysis It was also experimentally found that using 2i/LIF media at transition from pre-iPSCs to iPSCs is more preferable than serum/LIF due to increase in the efficiency of reprogramming [71].

Nanog expression heterogeneity in D 2 and D 3 domains
In the proposed model the parameter A characterizes OS activation by the signal A. OS expression dynamics linearly depends on the A value, while dynamics of Nanog expression is not so intuitive clear.
Expression of Nanog and some other genes including Nanog targets is heterogeneous on serum/LIF media [35, 36, 68, 69, 72, 73; 74, 75]. The heterogeneity is dynamically maintained, with individual cells exhibiting transient changes in expression levels. Nanog heterogeneity has been widely studied by mathematical modeling approach and it was shown that the phenomenon may be inferred from properties of OSN interaction circuit, activity of signal transduction pathways or transcriptional stochasticity induced transitions [12,36,40,43,76]. Our model also describes this behavior in domains D 2 and D 3 , in which all revealed states are unstable ( Fig 3A). In particularly, both Nanog mRNA concentration and protein levels in nuclear and cytoplasmic compartments have oscillation dynamics (Fig 5).
The emergence of these attractors occurs at h > 8. 68 and as shown in Fig 5A and 5B, both Nanog and pluripotent/differentiation factors (Fig 5C and 5D) have expression oscillations and it does not depend on the level of OS activation (Fig 5A and 5C: A < 0.277, Fig 5B and 5D: A > 0.277). In dynamics of pluripotency (w 1 ) and differentiation (w 2 ) factors A values influence the proneness of the whole system to differentiation or pluripotency: at A < 0.277, it biased to differentiation and at A > 0.277 it has pluripotent tendencies.
Negative feedback is a general requirement for oscillatory behavior [77]. For Nanog regulation, the shortest feedback is Nanog autoinhibition. This may have different number of intermediate steps due to involvement of histone acetylation [78], histone methylation [79], DNA methylation [80,81] or additional regulators [55]. The diversity of additional steps in Nanog autorepression can explain diverse unstable states represented by oscillations in the domain D 3 .
Therefore, we suppose that obtained unstable states may represent different routes to and types of pre-iPSC states or other incomplete reprogramming derivates. At least oscillations in expression were indicated by real-time imaging for several differentiation genes in partly differentiated cells (for example, for neural progenitors) where these oscillations serve for maintenance of multipotency before the sustained expression of these genes in their lineages after cell fate decisions [82,83]. Whereas sustained expression of the key differentiation gene in the lineage is the trait of the lineage differentiation steps, simultaneous oscillation expression of the key differentiation genes of several lineages marks the multipotent state of the common proliferating progenitor. Thus, during developmental transitions, the oscillatory expression of several fate determination factors were recorded at the multipotent state, whereas after choosing cell fate only one of them was sustainably expressed in the course of chosen differentiation. Fluctuations and oscillations in gene expression were suggested as a basic character of stemness, potentiality both to proliferate and differentiate [84]. While a loss of oscillatory dynamics leads to differentiation and the loss of stemness, reactivating of oscillatory expression of several key lineage specific genes was predicted to restore pluripotency. It was also shown that counteracting lineage specifiers synergistically induced pluripotency in the absence of both Oct4 and Sox2 [85,86]. The malignant transformation also looks like differentiation-stemness transition. We can only speculate about the role of increase in non-linearity of Nanog autorepression in this process. Nevertheless, Zfp281 mediates Nanog autorepression by means of the NuRD complex recruitment and inhibits somatic cell reprogramming by repressing Nanog activity [34]. The ZNF281 increased expression played an important role in tumor cell formation and this ZNF281 function was only recently discussed and reviewed [87]. Thus, we hypothesize that in our simulations the border between D 3 and D 4 domains may describe some experimentally observed  = 0.2 (a, c) and A = 0.3 (b, d). The other parameters were fixed. c: The pluripotent factors w 1 were suppressed and the differentiation factors w 2 were expressed. This state corresponds to differentiation. d: Pluripotent factors were highly expressed, and differentiation factors were suppressed. This state corresponds to pluripotency. transitions from one steady state to another (differentiation-pluripotency or pluripotency-differentiation). In the next section we shall consider in details conditions for this bistability.

Conditions for differentiation/pluripotency bistability in the D 4 domain
As we pointed out in the previous section, the oscillations of differentiation factors were recorded at the multipotency state before differentiation [34,83]. In mESC cultivated on 2i media where differentiation driving signals (MAPK/ERK1/2 and GSK3) are inhibited, there is uniform expression of pluripotency genes and low to absent expression of differentiation genes [68,69]. Whereas on LIF/serum media mESCs have fluctuations in expression of Nanog, Pecam1, Rex1(Zfp42), Dppa3 (Stella), Tbx3, Klf4, Esrrb Tcl1, Fgf5, Bry and Dnmt3b [35,36,68,69,72,73,74,88]. There are both pluripotency and lineage specific genes among them. Due to this, mESC on LIF/serum media are more prone to differentiation than those on 2i/LIF [68,69]. Between naïve pluripotency and differentiation states there is a transition stage ( Fig  3A; [61]) with fluctuations in expression of both pluripotency and differentiation genes.
To find out model stationary solutions for pluriponency/differentiation and differentiation/pluripotency transitions with pronounced transition zone we investigated how the model parameters determine the dynamics of the PTGs and DTGs. More precisely, we studied the cell state dependence on parameters, a 3 and a 7 , related to the free energies of A signaling molecule binding to the promotor regions of regulated genes in terms of the statistical mechanics approach [10] . The analysis demonstrated that there is a range of the parameters a 3 and a 7 , (Fig 6) for which the system has a bistable switch-like behavior (Fig 7).
As can be seen from Fig 6, one observes the parameters range marked by blue color in which the switch-like behavior is demonstrated. It was numerically found that a straight line a 3 = a 7 in the plane (a 3 , a 7 ) splits the wedge-shaped region into two areas: the first one is an approximate boundary between a 3 < a 7 , in which the differentiation state only exists at A ! 0, and the second one is an area a 3 > a 7 , in which both pluripotency and differentiation criteria can be performed upon a certain value of the external A signal.
In During ESC differentiation, the dynamics is more complicated: Oct4 level increased but Sox2 level was repressed at initiation of mesondoderm differentiation and in opposite low Oct4 and high Sox2 levels characterized the beginning of neuroectoderm differentiation [27]. Cell fate selection by decreasing of one of OS factors and increasing another one was also proved by the facts that Oct4 or Sox2 overexpression led to specific differentiation.
Oct4 overexpression resulted in differentiation to mesendoderm precursors in the presence of LIF and to ectoderm in the absence of LIF and BMPs [26,89]. Increase in Sox2 expression triggered mESCs to neuroectoderm and mesoderm [67].

Positive feedback loops from Nanog to Oct 4 and Sox 2 drive the system towards the ESC naïve stage
Lakatos with coauthors [43] considered in their simulations different variants of OSN transcriptional regulatory circuit, in which Nanog activates or not activates Oct4 expression; Oct4 represses Nanog and Nanog autorepression and combination of two latest regulations. The model in all these modifications demonstrated bistable, switch-like behavior. To verify the significance of the feedback loops from Nanog to Oct4 and Sox2 in the revised network we also considered the model extended by an addition of positive feedback loops that led to the update of the system of differential equations (Fig 8; Materials and Methods).
The system also exhibits bistability for the same range of parameter values as not extended variant. Moreover, accounting for the additional positive feedback loops did not result in emergence of the new type of steady states (subsection "Threshold in OS activation defining differentiation-pluripotency transition"). However, it led to significant qualitative change in  depending on (a 3 , a 7 ) parameters and at h = 6. Highlighted region is the range of parameter values, having which the system has switch-like behavior. Furthermore, the analysis indicated that a straight line a 3 = a 7 divides the plane (a 3 , a 7 ) it into two areas. When a 3 < a 7 , the cell has differentiated state at all values A ! 0. When a 3 > a 7 , there will be some A 0 , that upon A > A 0 the cell is pluripotent, while at A < A 0 the cell will differentiate.
https://doi.org/10.1371/journal.pone.0194464.g006 functional cell state on the plane (a 3 , a 7 ) (Fig 8). The numerical analysis showed that the straight line a 3 = a 7 that belong to the wedge-shaped region in the plane (a 3 , a 7 ) is an approximate boundary between an area a 3 < a 7 , in which both pluripotency and differentiation criteria can be performed upon a certain value of the external A signal and an area a 3 > a 7 , in which the pluripotent state only exists at A ! 0.
Therefore, additional activation of Oct4 and Sox2 expressions through Nanog positive feedback gives rise to much stronger enhance of the PTG expression simultaneously with increasing of repression effect on DTG expression. This corresponds to naïve pluripotency where differentiation genes are strongly repressed and expression of pluripotency genes is stable and homogeneous [68,69].

Concluding remarks
In this work, a new kinetic model of revisited minimal regulatory circuit for mouse pluripotent cell induction, self-renewal and differentiation was proposed. We have conducted a thorough parametric analysis of the developed model. Numerical simulations suggest that the system dynamics is mainly sensitive to variations of two model parameters: h , which is a Hill coefficient determining the nonlinear effect of Nanog autorepression and A parameter, which is a value of the external signal A for activation of OS expression. The model predicts for a mouse pluripotent cell the existence of four dynamical domains with different numbers of stable and unstable steady states, which, as we suppose, can present a developmental progression from ground state ESC to lineage commitment and vice versa. Taken together, the computational study indicates that molecular mechanisms of Nanog regulation and OS activation are the most essential for differentiation/pluripotency transition.

Dynamical model of the revised core pluripotency network
We considered an autonomous system of differential equations representing a dynamical model of the revised core pluripotency network (Fig 1B; S1 File). The model constitutes of three groups of equations with parameters a i , i = 1,2,. . .,29, A for the first group, parameters b i , i = 1,2,. . .,14, h for the second group and parameters c i , i = 1,2,. . .,10, for the third group. Initial values of the parameters and corresponding references are represented below (Table 1).  depending on (a 3 , a 7 ) parameters and at h = 6. Highlighted region is the parameters range, for which the switch-like behavior has existed. Furthermore, the analysis indicated that a straight line a 3 = a 7 in the plane (a 3 , a 7 ) divides it into two areas. When a 3 < a 7 , there will be some A 0 , that upon A > A 0 the cell is pluripotent, while at A < A 0 the cell will differentiate. When a 3 > a 7 , the cell has pluripotent state at all values A ! 0. https://doi.org/10.1371/journal.pone.0194464.g008 The first group of equations reflects concentration changes for both mRNAs and proteins of Oct4 and Sox2, considers changes in Oct4-Sox2 heterodimer level and includes 7 differential equations for u i , i = 1,2,. . .,7 This group is independent on other following two groups of equations, but defines dynamical behaviors of variables in these groups. The equation system for the first group is written as: where u 1 -Oct4 mRNA concentration, u 2 -Sox2 mRNA concentration, u 3 -Oct4-Sox2 heterodimer concentration, u 4, u 6 -Oct4 and Sox protein concentrations in the nucleus, correspondingly, u 5, u 7 -Oct4 and Sox protein concentrations in the cytoplasm, correspondingly. The second group representing concentration changes for both Nanog mRNA and protein (in the cytoplasm as well as in the nucleus) is described by the system of equations for v 1, v 2, v 3 variables: where v 1 -Nanog mRNA concentration, v 2 -Nanog protein concentration in the nucleus, v 3 -Nanog protein concentration in the cytoplasm.
The last group contains equations describing the dynamic of concentration changes for pluripotency and differentiation factors, w 1 and w 2 , correspondingly: The numerical analysis of the steady states based on a method of solution continuation with respect to a parameter allowed for the investigation of functional cell state's dependence on the external signal A+ concentration. To define cellular attractors we employed diagrams of stationary solutions that represent dependence of w 1 and w 2 concentrations on parameter A and by means of the next criterions: w 1 > w 2 -criterion for pluripotent state, w 1 < w 2 -criterion for differentiate state.

The model extension by accounting of several positive feedback loops in the revised core network
Chickarmane and coauthors [10] showed that positive feedback loops between Oct4-Sox2-Nanog factors in the core transcriptional network give rise to bistable switch-like behavior. To verify the role of the feedback loops in the revised network we considered a model extended by an addition of the positive feedback loops (S1 File). For this the system of differential Eq 1 was written as: where q 1 , q 2 , q 3 , q 4 , p 1 , p 2 , h 1 , h 2 are additional to a i , b j , c k , A, h parameters introduced for taking into account the positive feedback loops. Equations in the systems (2)-(3) remain unchanged. It should be noted that feedback loops impose interdependencies between the systems (4) and (2) for the extended model.