A stochastic and dynamical view of pluripotency in mouse embryonic stem cells

Pluripotent embryonic stem cells are of paramount importance for biomedical sciences because of their innate ability for self-renewal and differentiation into all major cell lines. The fateful decision to exit or remain in the pluripotent state is regulated by complex genetic regulatory networks. The rapid growth of single-cell sequencing data has greatly stimulated applications of statistical and machine learning methods for inferring topologies of pluripotency regulating genetic networks. The inferred network topologies, however, often only encode Boolean information while remaining silent about the roles of dynamics and molecular stochasticity inherent in gene expression. Herein we develop a framework for systematically extending Boolean-level network topologies into higher resolution models of networks which explicitly account for the promoter architectures and gene state switching dynamics. We show the framework to be useful for disentangling the various contributions that gene switching, external signaling, and network topology make to the global heterogeneity and dynamics of transcription factor populations. We find the pluripotent state of the network to be a steady state which is robust to global variations of gene switching rates which we argue are a good proxy for epigenetic states of individual promoters. The temporal dynamics of exiting the pluripotent state, on the other hand, is significantly influenced by the rates of genetic switching which makes cells more responsive to changes in extracellular signals.


I. INTRODUCTION
Embryonic stem cells derived from mammalian blastocyst are pluripotent: they show an indefinite capacity for self-renewal and the ability to differentiate into every cell type constituting an adult organism [1][2][3]. The development of healthy tissues hinges on the ability of these pluripotent stem cells to make critical decisions determining when and into which kind of cells to differentiate in response to their environment. Fates of embryonic cells are therefore decided through sophisticated biological computations orchestrated by a vast regulatory network consisting of genetic, epigenetic and signaling layers [3].
The dynamics of genetic expressions take place in molecular environments and are subject to intrinsic noise due to the discrete nature of the molecular copy numbers and to extrinsic noise from extracellular environment [4]. Thus, while at the level of the organism development is often predictable with a well-defined order of events, at the level of single cells fate determination is fundamentally stochastic [5,6]. An outstanding question involves understanding to what extent the inherent molecular stochasticity in embryonic cells is suppressed or exploited for functional purposes. Studies probing transcription in single cells have uncovered high variabilities of gene expression in embryonic stem cell populations [7][8][9]. Several hypotheses about the functional roles of stochasticity in embryonic stem cells have been put forward. Observations of dynamic and heterogeneous expression patterns of core transcription factors such as Nanog, Oct4 and Sox2 have promoted the view that stochastic excursions of these factors are governing the stability of pluripotent states [10,11]. A different hypothesis claims transcriptional noise to be advantageous by facilitating the exploration of the space of a gene network such that, at any instant, a subpopulation of cells is optimally primed to be responsive to differentiation signals [12]. The heterogeneity of populations of pluripotent cells has also raised some concerns that pluripotency is ill-defined on a single-cell level [12] and instead should be viewed as a statistical property of the macroscopic state emerging at the level of an ensemble of cells. A comprehensive physical picture of pluripotency at the single-cell level therefore still remains unclear.
In this regard, the roles of modeling and computational approaches are seen as especially important for bridging the gap between our understanding of molecular dynamics of regulatory networks and phenotypic outcomes. The rapid accumulation of data of gene expression through single-cell experiments has stimulated the techniques of statistical inference and computational systems biology for mapping embryonic stem cell networks [14][15][16]. In vitro studies of mouse embryonic stem cells (mESC) in different culture conditions, in particular, have become an ideal model system for exploring mechanistic issues surrounding pluripotency and lineage commitment in mammals [16]. In a tour de force study of mESC by Dunn et al. [13], regulatory relationships between transcription factors were uncovered through analysis of pairwise correlations in gene expression. Using known experimental data as constraints, a network topology with a minimal set of genes was derived with a high degree of predictive power with respect to perturbations of the network.  [13]. Each node corresponds to a given gene and their placement from left to right is chosen to indicate a general trend of downstreamness from the inputs. In our mechanistic model, each gene produces a unique transcription factor at a given rate, depending on the binding state of its promoter sites. These transcription factors go on to bind and activate (black arrow) or repress (red bar) other genes. The three nodes on left correspond to extra-cellular signals, which are either absent or present. The right panel shows our assumed molecular logic of transcriptional regulation. Shown are the outcomes of combinatorial binding of activators and repressors to promoter sites. Depending on the configuration of the promoter site, transcription factors are produced with rates 0,αm, andαmax, modeling the effects of cooperative binding.
Network topologies, however, remain silent about the roles of molecular noise and dynamics in stem cell differentiation governed by stochastic biochemical reactions. Furthermore, in order to validate that the inferred network reflects true microscopic reality of cell and is not a result of overfitting, one has to ultimately test the results using mass-action-based kinetics which integrate relevant molecular factors. The key challenge lies in finding the adequate resolution for the network which is able to be predictive and does not pose insurmountable computational burden.
In the present work, we outline a framework for extending boolean resolution networks into stochastic and dynamical models-with microscopic resolution of promoter architecture and individual gene switching events. The computational scheme utilizes static boolean information about the network topology and uses novel analytical model reduction to increase the computational efficiency, allowing for extensive searches in the space of microscopic reaction rates. This framework is successfully applied to the network topology inferred by Dunn et al. [13] in order to build a mass action based stochastic dynamic model which is capable of describing both the discrete states of all the genes and the populations of transcription factors. Starting from minimal assumptions about the rates of various reactions, we find a few distinct gene swithcing regimes where a remarkable agreement with the experimental gene expression profiles is achieved for all combinations of experimentally controlled external signals. We show that gene expression levels in complex regulatory networks are not a unique function of gene switching rates which cautions against over-interpreting boolean level networks and suggests strategies of inference which utilize higher moments in distribution of transcription factors. Using single cell experiments which have probed expression of pluripotency factors [7,8] we are able to argue that gene switching in pluripotent states happens primarily on the intermediate scale relative to the reactions of creation and degradation (dilution). This regime better agrees with the diverse set of experiments available and provides explanation for the multimodality in distributions of transcription factors and burst like expression dynamics for some genes.
In the second half of the paper, armed with a predictive and physically motivated model of pluripotency network, we explore the dynamics of lineage commitment driven by withdrawal of various well documented signals (LIF, 2i) for maintaining the naïve state of pluripotency. We find a number of non-trivial consequences of molecular noise and gene switching dynamics. In particular we show that intermediate gene switching regime generates higher sensitivity for the network when responding to external signals.

A. Molecular logic of transcriptional regulation in ESC
Here we outline the basic steps of constructing the microscopic genetic regulatory networks from the Boolean topology of genes inferred from experiments [13]. The network topology (Fig. 1) contains static information about the types of interactions between pairs of genes. The interactions are classified as being either repressing or activating.
To study the dynamics of complex genetic networks, one has to extend the Boolean level description to account for the molecular logic of gene regulation. This molecular logic specifies the precise relation between the binding of transcription factors to genomic sites and the regulatory outcome in terms of gene activation, repression, etc. In the case where the same sites can be bound to different transcription factors, the combinatorial nature of regulation can give rise to ambiguity in molecular logic. Such ambiguities can be resolved either by directly inferring regulatory logic or by simulating different combinatorial possibilities until sufficient agreement with experiments is reached. Once molecular logic is defined the topologic level description can be extended into a set of biochemical reactions.
The reactions in the present model include (i) the binding/unbinding of transcription factors to promoter sites of genes and (ii) the creation/degradation of populations of transcription factors. We adopt the network of mESCs inferred by Dunn et al. [13]. There are N = 12 genes in the network and we assign each gene G i to a discrete variable which describes the configuration of its promoter sites. The values of G i determine if the transcriptional activity is in the activated, repressed or neutral state. The n promoter sites of a particular gene allow binding of at most n transcription factors. Since a large majority of eukaryotic transcription factors in this network bind as dimers [17] we set n = 2 for all genes [18]. The binding reactions are First TF binding: Notice that we have adopted the uncooperative binding mechanism which assumes the transcription factors bind to any of the unbound promoter sites independently with a constant ratek on . Other cooperative binding mechanisms can be easily built into this framework. Each gene codes for its own transcription factor P i where indices i = 1, ..., N label both genes and the corresponding transcription factors. These promoters can act as repressors or activators of other genes: this information is contained in the experimentally-inferred topology. The molecular logic describing how the combinatorial binding of activators and repressors regulate each gene, however, is an assumption in the model.
When the rate of transcription is non-zero, genes produce copies of their respective transcription factors. We use a one step model of transcription [19] whereby transcription factors are produced via a unimolecular step. This approximation coarse grains several sequential intermediate reactions, such as mRNA production, into one effective step [20,21]. This approximation has become popular in models of protein production [22][23][24]. All of the transcription factors are assumed to have finite lifetimes set by the rate of degradation ensuring the existence of stable steady states with finite number of molecules. The reactions are TF production: We note that one can adopt a view with higher resolution of the network depending on the available experimental data and the type of question posed. For instance one may consider explicitly including the steps of cell cycle regulation, different epigenetic states, binding of non-coding RNAs etc. For simplicity and illustrative purpose, we consider the most simplified dynamical model which describes only the promoter configurations and the population dynamics of the transcription factors-this model aims to use the optimal resolution for capturing trends in gene expression while remaining feasible for efficient stochastic simulations capturing the variability in an ensemble of cells. After converting the network topology into a higher resolution network of biochemical reactions, our next goal is to exhaustively sample a vast space of parameter space to identify the optimal parameter regime with which the model best reproduces all of the experimental constraints.

B. Multi-scale simulation of complex genetic networks
Biochemical reactions in gene networks are of fundamentally probabilistic nature. The most rigorous way to simulate a network of reactions is by numerically solving the master equation which accounts for all possible states of the network down to the level of single molecules [4,19]. However in high dimensions, i.e., when the number of species is large, this approach is not computationally efficient. Instead, kinetic Monte Carlo algorithms are the most straightforward way to generate sample paths of the stochastic process. Fully individual-based models, however, still suffer from a steep scaling of computational time with the number of components. This fact renders them inefficient for simulating large gene networks, especially when it comes to scanning or exploring the parameter space.
A wide range of approximate schemes have been employed to simulate large scale gene regulatory networks. Most conventional approximations so far have been the size-expansion methods which are known for being problematic when the molecular noise induced by genetic switching becomes non-negligible. On the other hand, for embryonic stem cell networks it is essential to account for the stochastic nature of genetic switches which give rise to multiple attractor states corresponding to various phenotypes. A network with N genes can theoretically have ∼ 2 N attractors. Even if populations of all the other species are present in large quantities, the stochastic fluctuations caused by the genetic switches (due to stochastic binding-unbinding events of the transcriptional factors) between ON and OFF states cannot be ignored, unless the switching is operating in the extremely fast limit compared to any other reactions, known as the "adiabatic regime" [25][26][27]. In the other cases-the non-adiabatic regime-gene switching can completely dominate the dynamics in the network [25,28]. Eukaryotic gene regulatory networks are often found in the intermediate regime where gene switching events are dynamically interwoven with the rest of the reactions in the network and cannot be ignored [5,29]. Single-cell studies of mESC in particular have shown bursty behavior in gene expression with sudden jumps in levels of proteins resulting in multi-modal distributions of core transcription factors [7].
Many of the early computational models of embryonic stem cells have focused on small fragments of pluripotency network, typically involving bistable switches [10,30]. These early pioneering studies have yielded many insights on stochastic decision making in regards to fate determination and self-renewal. A few studied have looked at larger portions of the regulatory network [31][32][33][34] but at the cost of making near-adiabatic switching approximations thereby ignoring explicit stochastic treatment of genetic switching dynamics. Lastly, with very few exceptions, the kinetic parameters in those models were not throughly explored or data-informed and had to be picked relying on physical intuition alone. All of these limitations have been completely overcome in the present work. A series of recent studies [22-25, 35, 36] on gene expression dynamics have developed a novel computational framework utilizing piecewise-deterministic Markov process (PDMP) [37] as approximations to the fully individual-based model. This approach treats genetic switching events exactly, while assuming noise due to the finite nature of populations of transcription factors to be small in comparison. As will soon be seen, the assumptions underlying the PDMP approach turn out to be sound as we go on to obtain a nearly perfect quantitative agreement with full blown kinetic Monte Carlo schemes even for the case of the complex networks of ESC operating in the intermediate gene switching regime.
The PDMP simulations carried out in the present study show nearly O(10 3 ) fold faster generating stochastic trajectories compared to conventional individual-based kinetic Monte Carlo techniques (the algorithm for generating PDMP sample paths is provided in the Supplementary Information). This rigorous and rapid sampling of gene switching events has not only allowed us to investigate the stochastic dynamics of the regulatory network at a longer timescale compared to conventional kinetic Monte Carlo methods, but also enabled us to explore a vast parameter space efficiently. We have used the obtained information to derive microscopic resolution models of ESC. The mathematical details of the model described in greater detail can be found in the SI.

III. RESULTS
A. In the pluripotent state the mean levels of gene expression are a hardwired feature of network architecture whereas heterogeneity is shaped significantly by stochastic dynamics of genetic switches.
In this section we first identify certain parameters which provide results which are consistent experiment. Following this the model becomes predictive about dynamical features which are not contained in the topology.
Having reduced the individual-based model to a computationally tractable PDMP, we are at a point where we can tune the model parameters to match experimental data. The experimental data used for constraining the rate coefficients are the binarized gene expression levels of transcription factors under well-defined culture conditions consisting of different combinations of leukemia inhibitory factor (LIF), glycogen synthase kinase 3 (CH), and mitogen-activated protein kinase (PD). Identical culture conditions were used by Dunn et al. [13] to infer the original Boolean topology are LIF, 2i (same as PD+CH), LIF+2i, LIF+PD, LIF+CH. We employ the Hamming distance (the number of discrepancies between the simulated and experimental profiles) as a cost-function for optimization. The Hamming distance is minimized through multiple rounds of simulations where we vary the three free parameters: the rates of promoter binding k on and dissociation k off , and the basal transcription rate α m . Through this procedure we find that, in certain parameter regimes, our assumed molecular logic closely reproduces the gene expression profiles from experiments for different combinations of external signals. The dimeric binding-modeled in our framework as two binding sites (n = 2)-of TF to promoter sites in particular leads to a more sensitive switching and emerges as an important feature of the network which we show to be advantageous for reproducing experimental data [38]. We find localized parameter sets in three distinct regimes of gene switching: slow, intermediate and fast; all of which successfully reproduce the experimental gene expression patterns (Fig. 2) corresponding to pluripotent and lineage committed cells. The fact that the rate parameters of the network occupy finite regions and cover different regimes suggests that distinct gene expression profiles can tolerate fluctuations in reaction rates. Such rate fluctuations, reflecting the effect of extrinsic noise, are inevitable in dynamic cellular environments of embryonic cells which experience frequent epigenetic and extracellular perturbations [16,39,40]. In a way changes in gene switching rates can be seen as a proxy for how global epigenetic changes govern the rates of transcription factor binding to target genomic regions.
From the methodology point of view the absence of a unique regime of rates implies the following: inferring networks using only mean levels of gene expression (as is done for Boolean networks) may lead to the loss of valuable information contained in higher moments of distribution. Thus new approaches of inference need to be developed in order to account for broad distributions of transcription factors. For this reason, we look beyond comparisons of mean expression levels and turn to comparing stationary distributions of transcription factors observed in experiments with the computationally generated distributions in three chosen parameter regimes identified in Fig. 2. The pluripotent state of mESCs has been the subject of intense investigations by nucleic acid-based single-cell techniques such as RNAseq, sm-FISH, qPCR and there is now extensive data on the steady state distributions of RNA and transcription factors maintained under pluripotency-favoring culture conditions [7,8,41]. These experiments have revealed a heterogeneous nature of gene expression with many core pluripotency factors, such as Nanog, having long tailed or bimodal distributions. Additionally, the single-cell stochastic trajectories of TFs have shown sharp, bursty transitions Although the PDMP approach accurately captures the effects of genetic switching, it assumes demographic noise arising from finite populations to be negligible. To test the validity of this assumption and assess the contribution of different sources of noise in establishing the steady state distribution of the pluripotent state, we carry out individualbased simulations for the intermediate switching regime of the network, Fig 4. In the individual-based model all reactions are treated stochastically thereby accounting for all of the sources of noise in the system. The resulting gene expression profiles follow closely those obtained by PDMP simulations, showing that the noise arising from stochastic switching events of promoter configuration accounts for the significant part of overall variability in the network. Trajectories of individual transcription factors show that indeed most of the variance in the the molecular distributions are generated by genetic switching events which appear as abrupt stochastic jumps. Consistent with experiments, we also find that under LIF gene expression is more heterogeneous than under 2i, but also that the overall levels of expression are higher [7,8]. In all three regimes of genetic switching supporting pluripotency (LIF+2i,LIF+PD, LIF+CH, LIF and 2i), core transcriptions factors such as Nanog, Oct4, Sox2 are highly expressed. The same factors are also repressed in conditions favoring differentiation (CH, PD and none) irrespective of gene switching regime. This shows that the pluripotent and differentiated states, as determined by the pattern of gene expression, are hardwired in the architecture of the genetic network independent of genetic switching rates. Nevertheless, as we will show in the next section, the routes and dynamics of lineage commitment from pluripotent states strongly depend on the level of molecular noise generated in different gene switching regimes.  [13] Leukemia factor LIF based signaling have been shown to provide a stable environment for maintaining pluripotency of stem cells in vitro [3,13]. Conversely, withdrawal of either LIF or 2i leads to irreversible lineage commitment after a 24-hour period. Despite a similar ability to guard pluripotent cells against lineage commitment, these factors deploy different regulatory mechanisms reflected in distinct distributions of pluripotency factors. As a result, stem cell differentiation by withdrawal of different signals proceeds via different routes. To gain a mechanistic understanding of how the interplay of signaling, molecular noise and network architecture give rise to the steady state expression profiles, we study the dynamics of transitioning between the pluripotent and lineage committed states induced by rapid initiation and withdrawal of signaling conditions (LIF, CH, PD).
The temporal evolution of distributions of the TFs exiting (LIF/2i withdrawal) and entering (LIF/2i immersion) pluripotent states reveals rich dynamical signatures of these transitions (Fig 5). To reveal the role of stochasticity in these transitions, we compare the intermediate regime-which is dominated by genetic switching-to the fast regimein which transitions are largely governed by the network topology and the fluctuation of promoter configuration is almost-completely ignored. The irreversible nature of transitions manifests clearly in different routes exiting and entering the pluripotent state (transitions to and from the None state in Fig 5).
In the intermediate regime of genetic switching rates, expression noise greatly facilitates transitions out of a pluripotent state by making the network more responsive to changes in environmental signaling. In contrast, in the fast switching regime, upon withdrawal of pluripotency signals the downstream regulation happens on a much slower time-scale with some factors remaining virtually unresponsive to changes of signaling. This signaling enhancement in the intermediate regime reveals the importance of molecular noise in making pluripotent states more sensitive to environmental conditions. There are qualitatively different patterns of re-entrance into pluripotent states upon LIF vs 2i addition, with LIF being much more efficient at reversing pluripotency compared to the 2i. The different potential of signaling culture for pluripotency reversal has been established in experiments [3] which have shown that in the later stages of commitment only the LIF is able to reverse lineage-primed cells back to their naive pluripotent states. The exit and re-entrance from pluripotency upon withdrawal/addition of LIF shows complex signatures of hysteresis and bifurcations. This suggests that there can be multiple pathways of entering or exiting pluripotency. The transition times for all signaling-induced changes of the steady states of the network are visualized on a kinetic diagram (Fig 6). The kinetic diagram shows an underlying structure to these transitions where conversion among pluripotent states takes place with lower "activation barriers" compared to transitions accompanying loss of pluripotency. In the intermediate switching regime there is a clear time-scale separation between transitions which keep cells in pluripotent states and transitions out of pluripotent states. One may argue that such a time-scale separation between signals inducing differentiation and pluripotency allows embryonic cells to execute developmental decisions more faithfully. In the fast switching regime this clear time-scale separation is partially lost where only the loss of all three signals is separated from the rest of the transitions.
Detailed analysis of individual distributions and trajectories of transcription factors can be very informative due to their high information content. It is, however, not immediately clear how changes in the expression of individual genes contribute to global changes corresponding to different phenotypic transitions. To reveal such global changes, we project stochastic trajectories of all transcription factors onto the first two eigenvectors obtained by principal component analysis of the reference pluripotent steady state (LIF+2i). Most of the variance of transcription factors is well captured by the few principal components. The high-dimensional steady state of the cellular network can thus be conveniently projected onto a 2-dimensional subspace, allowing us to visualize the attractor states of the network as probability landscapes π(P C 1 , P C 2 ) which are often masked by heterogeneous distributions.
Comparing these probability landscapes with different gene switching regimes reveals the distinct roles played by gene switching-induced molecular noise and the deterministic network topology in guiding the transition out of the pluripotent states (Fig 7). The intermediate gene switching regime, once again, appears to be the more viable regime underlying pluripotent states since the probability landscape shows up as a broad attractor with interconnected states. Going towards the limit of slow switching results in the fragmentation of the landscape into states separated by high barriers. This gene switching-induced remodeling of attractors shows the potential for regulation via global epigenetic changes which are purported to act via silencing or activating entire sets of genes at once. Thus one may view gene switching rates as a proxy for genome wide acetylation/methylation patterns which can dramatically alter the access of transcription factors to key target genomic sites. The sequential removal of pluripotency inducing signals reveals  6. Calculations of transition times between the stationary distributions of different external conditions. A larger, darker arrow indicates that a given transition takes a longer time to converge to its stationary state. This time-scale is measured by simulating a large ensemble (10 5 ) of PDMP sample paths to provide a simulated probability density, and finding the Jensen-Shannon divergence [42,43] between the instantaneous distribution of each TF and its final stationary distribution. The time for the each divergence to fall below a threshold (:= 0.3) is recorded, and we choose the largest of these as a quantification of the transition time.
a consistent change in the size of the attractor towards occupying smaller regions on the landscape. This argues for the physical state of network corresponding to pluripotent states to be the one with maximal variance of regulatory transcription factors where lineage commitment is accompanied by their gradual constraining and repression. A similar idea which views pluripotency as a macrostate emerging from an ensemble of cells which try to maximizes the information entropy with respect to regulatory transcription factors has been postulated before [12,44]. The analysis of steady state stochastic dynamics of pluripotency network in this work appears in agreement with his view. Furthermore, we are able to reveal the microscopic origin of this entropic paradigm. By analyzing pairwise correlation among different transcription factors we find that signals like LIF/2i create greater independence between the expression of core transcription factors leading them to explore larger range of values. Hence, removal of these signals leads to more constrained and inter-dependent patterns of gene expression for the same transcription factors which greatly diminishes overall variance. (Fig S1).

IV. DISCUSSION
Recent studies of ESC have increasingly emphasized systems level perspective on pluripotency and fate determination as outcomes of complex biological computations orchestrated by non-random networks of genes [45,46]. Rapid growth of gene expression data collected from controllable in vitro experiments on mESCs has made it possible to make reliable inferences of gene regulatory networks which govern the state of pluripotency [13]. These inferences show a highly interconnected nature of regulatory networks where signaling molecules, genomic binding, epigenetics and biochemical feedback act in a concerted manner in balancing pluripotency and decision making. Single-cell experiments have also revealed significant dynamic heterogeneity of gene expression in the population of cells [7,8], suggesting the important roles played by molecular noise and non-equilibrium processes. Schematic network topology models are no longer sufficient for rationalizing the level of detailed quantitative information obtained in single-cell studies. In the present work we have developed a multi-scale computational scheme for converting experimentally-inferred boolean topologies into quantitative and predictive models of networks with microscopic resolution of gene expression dynamics. The employed computational model is based on previously proposed hybrid stochastic approaches [22][23][24][25] in which the switching dynamics of individual genes are considered exactly while the rest of the biochemical reactions are approximated as deterministic. This hybrid stochastic approach is approximately a thousand fold faster than conventional kinetic Monte Carlo methods. This allows us to simulate large scale gene regulatory networks of ESC under different culture conditions and gene switching regimes. To infer the parameters in the model from the experimental data, we use hybrid simulations to exhaustively sample the space of rates and identify the parameter regime in which the model prediction best matches experimental data. The approximation using the hybrid scheme is validated by carrying out fully stochastic simulations of the network for the identified parameter sets. This agreement also shows that the switching events of genes-due to stochastic TFs binding to the promoter sites-to be a significant source of variance in the ESC networks. Furthermore we see the changes in gene switching rates as a proxy for global epigenetic modifications that can alter the rates of access of transcription factors to sites buried under chromatin structures. Thus the significant remodeling of steady state landscape that we see by varying the global gene switching rates gives us a glimpse of potentially powerful leverage that epigenetic changes of chromatin have over stability of pluripotent states. We find that the intermediate regime, in which gene switching rate is comparable to the other reaction rates in the network, is most consistent with single-cell measurements [7][8][9]. In this regime transcription factors show bursty dynamics which lead to heterogeneous distributions with some showing long tailed and bimodal features. We find signaling by LIF and 2i to be a major driving force maintaining the stability of pluripotent states. Withdrawal of either LIF/2i initiates lineage commitment via a robust pattern of reduced expression of Nanog/Oct4/Sox2 triad. To characterize the dynamics of lineage commitment, we compute transition times from pluripotent to differentiated steady states. We find that higher levels of molecular noise generated by slower gene switching make network more responsive to changes in signaling condition. Next, by carrying out principal component analysis on ensembles of gene expression profiles, we find a much simpler description of pluripotency and lineage commitment in terms of effective probability landscapes. As the signals safeguarding pluripotency are removed, these landscapes reveal a gradual narrowing of the steady state attractor explored by the network. Thus, we see a hierarchical organization of differentiation landscapes where pluripotent states pose the largest attractor which is maintained through the extracellular signals and the molecular noise of gene switching.
At last we believe the computational framework developed here should also be useful in organize and interpreting experimental data of other complex gene regulatory networks within a coherent and unified computational models thereby making it easier to conceive and vet physical hypothesis in a systematic way. Given the rapid rise of information from high throughput single-cell nucleic acid based techniques (RNA-seq, RNA-FISH, qPCR, etc), we expect such microscopic resolution models to play important roles for bridging the systems-level behavior of genetic networks with the underlying molecular-level processes of binding and regulation at the genomic sites.