Computational modelling of embryonic stem-cell fate control

Pluripotency is defined as the capacity of an individual cell to give rise to all lineages of the mature organism, including the germ line. This property is transient and restricted to a few epiblast cells in the early mammalian embryo, but it can be sustained in vitro through the establishment of embryonic stem cells (ESCs) (Evans and Kaufman, 1981; Martin, 1981). Mouse ESCs (mESCs) can be maintained in a naïve pluripotent state over many passages without losing their differentiation potential. Moreover, when injected into a host embryo, they return to normal development and contribute to chimera formation (Smith, 2001). Uncovering the molecular basis and understanding the design principles of pluripotency is, therefore, of utmost relevance for both stem cell research and developmental biology.

Like many other areas of modern life sciences, stem cell biology is increasingly influenced by novel technologies, such as molecular high-throughput techniques and (ultra-)high-resolution imaging. These new methods provide an unprecedented amount of information at levels previously inaccessible. Although we recognise the tremendous potential of this new level of quantification, it also generates new challenges, including new sources of error in experimental setup or analysis. Another challenge comes with the inherent complexity of biological systems: in the majority of cases, a functional and mechanistic understanding of the dynamic behaviour cannot be achieved by simply ‘adding up’ the individual effects observed for the different components of the system. As already Aristotle recognised, “the whole is more than the sum of its parts”, or, in more modern terms, the complex interplay of the individual components might lead to the emergence of unexpected behaviours.

It has become increasingly clear that a complementary application of methodologies from bioinformatics, statistics and dynamic modelling is required for an accurate description and a comprehension of complex interactions and emergent behaviours (Fig. 1). For example, computational approaches can define the structure and the dynamics of gene regulatory networks (GRNs, see Glossary, Box 1). However, care is required when interpreting such studies, as different structures potentially encode the same or similar biological functions. Moreover, the same network structure might generate different functional consequences. To distinguish between different potential explanations, the analysis of the network dynamics is as important as structure identification. Because such network dynamics might be non-intuitive, mathematical concepts and computational methods from dynamic systems theory are extremely helpful for their formal description and quantitative analysis. The fields of application of computational models are as versatile as the methodologies themselves (see Box 2).

In this Review, we will discuss strategies to build and understand computational models of GRNs in mESCs and in the blastocyst. First, we introduce the transcriptional core network of pluripotency. Then, we describe dynamic small-scale networks probing mechanisms that underlie transcription factor (TF) heterogeneity. We also summarise a selection of more complex pluripotency networks. Finally, we consider dynamic models of cellular reprogramming and lineage specification, with a brief outlook on the demand for multi-scale modelling. For a more detailed discussion on the biological aspects of mESC pluripotency and heterogeneity, we refer the reader to other, excellent reviews (Niwa, 2007; Chambers and Tomlinson, 2009; Martinez Arias and Brickman, 2011; Nichols and Smith, 2012; Torres-Padilla and Chambers, 2014).

Over the last decade, microarray and sequencing technologies have been used extensively to measure gene expression in mouse and human ESCs with the aim of identifying master regulators of pluripotency regulation and lineage commitment, and their interactions on a genome-wide scale (Ivanova et al., 2006; Loh et al., 2006; Matoba et al., 2006; Wang et al., 2006; Kim et al., 2008; Ding et al., 2009; Dowell et al., 2013). Although these large-scale studies uncovered an inherent complexity of pluripotency regulation, most attention has been paid to a core unit of three TFs, namely Oct4, Sox2 and Nanog. More than 15 years ago, Oct4 was shown to be essential for the establishment of pluripotency in the embryo and in cultured cells (Nichols et al., 1998). However, although a constant expression of Oct4 is required for self-renewal, the dosage also matters: both enhanced and reduced levels of Oct4 lead to lineage specification (Niwa et al., 2000). Oct4 cooperates with Sox2 to regulate the expression of multiple target genes (Chambers and Tomlinson, 2009), including fibroblast growth factor 4 (Fgf4) that initiates lineage specification (Kunath et al., 2007; Silva and Smith, 2008). The TF Nanog supports Oct4 and Sox2 in coordinating the activity of a number of target genes. Overexpression of Nanog alleviates the requirement for extrinsic culture medium supplements, such as leukaemia inhibitor factor (LIF) or the mitogen kinase inhibitor (MEKi), to prevent differentiation (Chambers et al., 2003; Mitsui et al., 2003; Ying et al., 2008). Although Nanog is essential to reach a pluripotent cell state, pluripotency can be maintained even in the absence of Nanog, albeit in an unstable fashion (Chambers et al., 2007; Silva et al., 2009).

Multiple experimental studies indicate that the core network constituted by Oct4, Sox2 and Nanog is characterised by a number of positive-feedback regulations (Chew et al., 2005; Rodda et al., 2005; Loh et al., 2006). The first ordinary differential equation (ODE) model (see Glossary, Box 1) of this core structure was published by Chickarmane et al. in 2006 (Fig. 2A) (Chickarmane et al., 2006). Here, the transcriptional regulation of the network components is mathematically described using Hill equations, which are widely used to model chemical reactions and approximate molecular interactions (Weiss, 1997). However, to simulate TF concentrations within this framework, parameters such as binding or transcription rates, as well as cooperativity parameters (that describe the stoichiometry of molecules) have to be defined. Unfortunately, a quantitative determination of biochemical parameters is challenging, if not impossible, and fitting them to experimental data often yields large uncertainties (Gutenkunst et al., 2007). Methodologies from dynamic systems theory, such as state space and bifurcation analysis (Box 3), are powerful tools to investigate the behaviour of a system qualitatively and quantitatively. Together with simulation studies and proper sensitivity analyses of model parameters, these approaches can identify whether the dynamics are consistent with the experimental observations and allow the effect of changes in network components or interaction rates to be predicted.

Chickarmane et al. investigated the dependency of stable Oct4 and Nanog expression states on external signals, demonstrating that the positive-feedback loops in the core network give rise to a system property termed bistability (Box 3) (Chickarmane et al., 2006). Bistability can be interpreted as the coexistence of two stable gene expression states associated with different developmental potentials of cells. The emergence and disappearance of stable expression levels as a function of a particular model parameter can be depicted by a bifurcation diagram, as shown in Fig. 2B. It reveals that the model illustrated in Fig. 2A can account for three distinct patterns of stable expression states, depending on the intensity of an input signal B, which negatively regulates Nanog transcription (Fig. 2B). The authors nominated bone morphogenetic protein 4 (BMP4) signalling and p53 as potential candidate factors for B. When signal B is weak, the pluripotency network is sustained and a stable expression state at high Nanog levels is established (termed ‘stem-cell state’, upper solid line in Fig. 2B). By contrast, if signal B is strong, Nanog autoregulation is suppressed and its transcription remains at a basal level. Thus, a second stable state at low Nanog levels is established (termed ‘differentiated state’, lower solid line in Fig. 2B). However, for an intermediate range of signal B, both states coexist, such that mESCs can attain either of the two expression states (region shaded in grey in Fig. 2B). Thus, a continuous increase of signal B intensity ‘moves’ cells from a monostable high into a monostable low state by transiently passing the bistable region. This process is reversible but the signal intensities at which cells change their state into one or the other direction (c₁ and c₂) are different. This property is called hysteresis and means that the behaviour of a system depends on its history. Due to the highly connected architecture of the underlying network, the model predicts the same qualitative behaviour for Oct4 and Sox2. Thus, in the model proposed by Chickarmane et al., all three TFs are either on or off, ensuring that the expression of common target stem cell and differentiation genes is mutually exclusive.

As shown above, the analysis of qualitative changes of gene expression patterns depending on model parameters provides insights into potential cellular states, their stability and robustness. However, to study dynamic properties of developing cell populations, the heterogeneity of gene expression and of cellular states has to be taken into account.

The maintenance of mESC pluripotency requires specific culture conditions. Standard mESC culture media contain serum factors and the cytokine LIF (Niwa et al., 1998), and promote a population-intrinsic heterogeneity reflected in the mosaic expression of pluripotency genes, including Nanog, Rex1 (also known as Zfp42 – Mouse Genome Informatics Database), Stella (also known as Dppa3 – Mouse Genome Informatics Database) and Esrrb (Chambers et al., 2007; Carter et al., 2008; Toyooka et al., 2008). Because of the key role of Nanog in the establishment of pluripotency, its heterogeneity has been studied extensively. Flow cytometry analyses of Nanog reporter cell lines show a stable and reproducible bimodal distribution of Nanog-high (NH) and Nanog-low (NL) cells (Chambers et al., 2007; Kalmar et al., 2009). Purified NH cells can give rise to NL cells and vice versa, revealing the dynamic capacity of mESCs to undergo state transitions (Box 3), and show a functional bias in their differentiation propensity (Singh et al., 2007; Marks et al., 2012). The biological relevance (Smith, 2013) and the molecular mechanisms underlying this intrinsic heterogeneity are still under debate (Miyanari and Torres-Padilla, 2012; Navarro et al., 2012; Karwacki-Neisius et al., 2013; Muñoz-Descalzo et al., 2013; Singer et al., 2014).

Complementary to experimental strategies, a few small-scale GRN models have been applied to improve the mechanistic understanding of dynamic heterogeneity in pluripotent mESCs (Kalmar et al., 2009; Glauche et al., 2010; Herberg et al., 2014). Three different mechanisms accounting for a bimodal Nanog distribution and reversible state transitions will be discussed here: first, an excitable model presuming monostable Nanog expression and extrinsic noise (Kalmar et al., 2009); second, an oscillation model based on a cyclic Nanog attractor (Glauche et al., 2010); and third, a fluctuation model assuming bistable Nanog expression and intrinsic noise (Glauche et al., 2010; Herberg et al., 2014). We review these models to evaluate their consistency with new experimental findings based on single-cell gene expression recently published by Singer et al. (2014).

In 2009, Kalmar et al. proposed a noise-driven, excitable model based on positive autoregulations and mutual activations of Oct4 and Nanog (Fig. 3A) (Kalmar et al., 2009). Additionally, the authors supposed that sufficiently high levels of Oct4 can repress Nanog expression, introducing a negative feedback. To formally describe these regulations, a stochastic differential equation (SDE) model (see Glossary, Box 1) has been used. Although the positive-feedback loops allow for bistability, parameters have been chosen such that only one stable NH state exists (compare with Fig. 2B: ‘Stem cell state’). Excursions from this state have been induced by a stochastic noise term (see ‘additive noise’ under the entry ‘Transcriptional noise’ in the Glossary, Box 1). Thus, mESCs stochastically escape from the stable expression state, which is characterised by high Oct4 and Nanog expression, towards a region of low Nanog and variable Oct4 expression. As no second (attracting) stable state exists, cells rapidly return to their origin and, thus, exhibit pulsing Nanog dynamics (Fig. 3A, bottom). This model predicted that excursions from the NH state are transient, providing a very short window of opportunity in which perturbations can become consolidated into a lineage commitment decision (Chambers et al., 2007).

Exploring the same topic, Glauche et al. suggested two alternative SDE models (Glauche et al., 2010). In the first model (model 1 in Fig. 3B) mESCs are bistable with respect to Nanog expression (compare with Fig. 2B, grey shaded region). State transitions are induced by a transcriptional noise (see ‘multiplicative noise’ under the entry ‘Transcriptional noise’ in the Glossary, Box 1) leading to stochastic fluctuations. In the second model (model 2 in Fig. 3B), a hypothetical factor X has been integrated such that Nanog is part of an activator-repressor system. The resulting negative feedback can generate a cyclic attractor and thus oscillations between high and low Nanog expression levels (Fig. 3B, bottom). Replacing the previously suggested direct activation of Oct4-Sox2 by Nanog with an indirect, double-negative feedback, Glauche et al. suggested a mechanistic explanation for the gatekeeper function of Nanog in the control of mESC differentiation (Silva and Smith, 2008). Indeed, assuming that Nanog prevents the effects of differentiation signals (represented by Y), only mESCs in the NL state can differentiate. Meanwhile, Muñoz-Descalzo et al. provided an alternative explanation, showing that Nanog buffers the differentiation-inducing activity of Oct4 through the formation of stable complexes (Muñoz-Descalzo et al., 2013) (see section below on lineage specification).

Recently, Herberg et al. put forward the fluctuation model to study whether different intensities of Fgf4/Erk signalling can account for culture-dependent differences in Nanog expression (Herberg et al., 2014). mESCs exposed to serum-free medium containing MEK/Erk and Gsk3β inhibitors, termed 2i, are captured in a pluripotent ground state characterised by a single, homogenous NH peak (Ying et al., 2008; Wray et al., 2010). Applying the model illustrated in Fig. 3C, the authors demonstrated that Fgf4/Erk signalling can regulate the existence of stable Nanog expression patterns by ‘shifting’ cells from a monostable high into a bistable regime as described above (compare with Fig. 2B, with signal B assumed to represent Fgf4/Erk). In particular, if Fgf4/Erk signalling is efficiently blocked (as in 2i), there is only one stable expression state at high Nanog levels. In this ground state, network-inherent stochasticity or (external) perturbations only induce micro-heterogeneity within an apparently uniform cell population. However, if Fgf4/Erk signalling is active, two stable expression states coexist (bistability), such that stochastic fluctuations generate functionally different subpopulations (macro-heterogeneity, Fig. 3C, bottom) (Huang, 2009; Herberg et al., 2014).

Although all three models consistently describe a stable, albeit dynamic, bimodal distribution of Nanog expression within a population of mESCs (compare with insets in Fig. 3, bottom row), they generate different hypotheses about the underlying mechanism and its dynamic outcome. These hypotheses guided experimental strategies to distinguish between the mechanisms (Box 2). For example, Glauche et al. demonstrated that stochastic fluctuations and deterministic oscillations are distinguishable based on the dynamics of the (re-)establishment of heterogeneity after cell sorting (Glauche et al., 2010). Furthermore, comparing the models revealed that measurements of the residence times of single cells in the NH and the NL state would allow the discrimination of mono- from bistable fluctuations. In a bistable system, residence times for both Nanog expression states would be prolonged, exceeding typical cell cycle times as predicted by Herberg et al. (2014), whereas a monostable system implies that unstable NL cells rapidly revert back into the stable NH state (Kalmar et al., 2009). Intriguingly, Singer et al. recently analysed single-cell expression dynamics, demonstrating that Nanog heterogeneity arises from a combination of stochastic transitions between coexisting gene expression states (i.e. macro-heterogeneity, bistability) and burst-like transcription within each state (i.e. micro-heterogeneity, noise) (Singer et al., 2014). They also showed that single mESCs remain in either one of the two Nanog states for multiple cell cycles. These findings clearly favour the fluctuation model suggested and described by Glauche et al. (2010) and Herberg et al. (2014).

Taken together, the complementary use of computational and experimental strategies clearly indicates that Nanog heterogeneity is based on two coexisting expression states and stochastic fluctuations at the transcriptional level. Specifically, the presented mathematical models brought up hypotheses on potential regulatory mechanisms, evaluated their consistency with available experimental data and revealed experimental strategies required to distinguish between them (Box 2).

Although the GRN of mESC pluripotency is centred on Oct4, Sox2 and Nanog, numerous studies point to the involvement of additional transcriptional (co-) factors, such as Klf4, Esrrb, Dapp1, Zfp281, Tcl1 and Tbx3 (Ivanova et al., 2006; Matoba et al., 2006; Wang et al., 2006; Kim et al., 2008), to mention only a few. In this section, we introduce three studies in which statistical and bioinformatics methods have been applied to analyse more complex TF interactions under different conditions.

MacArthur et al. studied the role of Nanog heterogeneity in early fate decisions by mimicking its temporal fluctuations (MacArthur et al., 2012). In particular, they analysed global gene expression changes during downregulation of Nanog, with a focus on elements of a previously published GRN (Kim et al., 2008). All factors, including Oct4, Sox2 and Klf4, were downregulated considerably later than Nanog, providing a time window in which cells can be rescued from differentiation by reintroducing Nanog expression. Using clustering and support vector machine approaches (see Glossary, Box 1) for the analyses of single-cell expression profiles, the authors strengthened the view that mESCs exit pluripotency gradually and stochastically. They also established two classifiers that distinguish pluripotent and non-pluripotent cells by learning the expression profiles of cell samples taken from the Gene Expression Omnibus (GEO) database (Edgar et al., 2002). The application of these classifiers showed an intermediate, primed cell state in which differentiation and pluripotency genes are co-expressed. Finally, MacArthur et al. analysed the feedback structure of the GRN and found that fluctuations in Nanog expression transiently activate different sub-networks driving transitions between a self-sustaining state (feedback-rich) and a differentiation-sensitive state (feedback-sparse). Similar results have been reported by Trott et al., who analysed single-cell gene expression levels with the objective to decipher the activity of (sub-)networks in heterogeneous mESCs (Trott et al., 2012).

Recently, Xu et al. constructed a GRN consisting of 15 pluripotency regulators and 15 lineage markers from the ESCAPE database (Xu et al., 2013,, 2014). Based on single-cell gene expression data of mESCs cultured in LIF/serum and 2i/LIF conditions, the authors attempted to reveal differences in the network structure between the two conditions using a Boolean network model (see Glossary, Box 1) to learn the underlying network logic. In that model, continuous gene expression values had to be converted to binary values (0: factor is off, 1: factor is on), which was carried out using a k-means clustering step (see Glossary, Box 1). To derive the logical function of each network node, an exhaustive search with limitations on the operators to constrain the search space (e.g. nesting, XOR gates and autoregulatory interactions were neglected) was performed, yielding a group of network models consistent with the experimental data. Xu et al. used this set of models to predict several characteristics, including the most common feedback loops, the differences between network interactions in the two culture conditions and previously uncharacterised interactions. For example, the models indicate that the positive feedback between Oct4 and Sox2 is absent in 2i/LIF conditions, but present in LIF/serum. However, this and other predictions have yet to be validated experimentally.

In a methodologically similar study, Dunn et al. applied a data-constrained Boolean model to define a minimal GRN that allowed them to reproduce expression profiles of mESC populations in different culture conditions (Dunn et al., 2014). The authors started their analysis with a model composed of 17 pre-selected factors implicated in mESC maintenance. To reveal possible interactions, expression correlations between all pairs of TFs were quantified using Pearson's correlation coefficient (see Glossary, Box 1). Subsequently, a software tool was used to constrain the set of possible models to those able to reproduce the experimental data. Dunn et al. only kept the components and interactions that are highly correlated and present in all network structures. This selection resulted in 12 components and 11 interactions in addition to pre-specified target interactions. The authors applied the candidate models to investigate the consequences of different TF knockouts and tested their predictions experimentally. Because none of the models could explain the outcome of the Tbx3/Klf2 double-knockdown, the observed behaviour after knockdown was added as a new data constraint. Similar to the results of Xu et al. (2014), the positive link between Oct4 and Sox2 did not show up during the selection process but was integrated manually to fulfil all predefined constraints and predictions in 2i/LIF conditions. The authors claimed that the constructed GRN model constituted ‘the most parsimonious network’ that maintains naïve pluripotency. Although combining a Boolean network model with prior information on regulatory interactions is generally feasible, in our opinion this statement should be complemented with the limiting clause ‘within the given set of assumptions’.

The presented studies are examples of computational methods suitable to evaluate gene expression data and to infer (TF- or culture-dependent) interactions between larger sets of genes. However, they also make clear that modelling the dynamics of more complex GRNs requires an even more simplified description of the system. Boolean model frameworks are commonly used to describe GRNs, but they present limitations (e.g. discretised expression levels) that have to be considered when interpreting the results (see Discussion). Applied in a proper way, these models are extremely useful to predict the consequences of perturbations (Box 2) and to identify candidate genes that play decisive roles in the loss or regain of pluripotency, as discussed in the next section.

In 2006, Takahashi and Yamanaka evaluated a number of candidate genes with respect to their potential to induce pluripotency in somatic cells (Takahashi and Yamanaka, 2006). The discovery that the overexpression of four factors (Oct4, Sox2, Klf4 and c-Myc) can direct the reprogramming of somatic cells into induced pluripotent stem cells (iPSCs) constituted a paradigm shift in developmental biology. However, details and mechanisms of the reprogramming process during iPSC generation still need to be fully elucidated (Takahashi and Yamanaka, 2013).

A current limitation of the cellular reprogramming process is its low efficiency, for which there are two conceptual explanatory models: the deterministic model, in which only some cells have the potential to generate iPSCs within a fixed and uniform time period (latency), and the stochastic model, in which most or even all cells are competent for reprogramming, but the latency differs (Hanna et al., 2009; Yamanaka, 2009). Hanna et al. combined time-series experiments with computational modelling to gain insights into the nature of reprogramming (Hanna et al., 2009). In their experiments, they monitored the proportions of iPSCs generated and found that reprogrammed cells appear in most clonal populations, if cultured for a sufficiently long time. However, the time to conversion apparently differs between clones, supporting the view that reprogramming is a continuous, stochastic process. Accordingly, the authors developed a mathematical model that considers reprogramming as a one-step stochastic process with a constant cell-intrinsic rate. Fitting the model to their experimental data, this intrinsic rate was estimated for different experimental settings. The authors argued that Nanog overexpression accelerates the reprogramming kinetics through a mechanism most likely independent of the cell proliferation rate. However, in that study, the molecular changes during reprogramming were analysed in heterogeneous cell populations, making sequential events occurring in single cells inaccessible. Combining single-cell expression analysis with a Bayesian network model (see Glossary, Box 1), Buganim et al. demonstrated that the reprogramming process can best be described by two phases: an early stochastic phase with high variation in gene expression and a subsequent, more hierarchical phase of gene activation (Buganim et al., 2012). The network model predicted that the activation of Sox2 initiates a number of consecutive steps that finally lead to fully reprogrammed iPSCs. The authors applied the hierarchical model to predict sets of TFs capable to induce pluripotency. Subsequently, these sets were tested experimentally, demonstrating that all of them facilitate iPSC generation but with different efficiencies, ranging from 0.2% (a combination without Oct4, Klf4, c-Myc) to 22.2% (using Oct4, Nanog, Esrrb, Klf4 and c-Myc). Interestingly, there is limited correlation between those genes that can facilitate efficient reprogramming and those of which the expression is predictive of future iPSC generation. For example, the endogenous expression of Oct4 does not necessarily predict iPSC generation, whereas the expression of Utf1 and Esrrb does.

Given the importance of both the genetic and epigenetic control of pluripotency, Artyomov et al. developed a computational model of pluripotency induction that couples both layers of regulation (Artyomov et al., 2010). In this study, all genes responsible for a particular cellular identity (e.g. Oct4, Sox2 and Nanog for pluripotency) were described as a single ensemble module. Moreover, cellular states, arranged in a hierarchical tree-like structure, were described not only by the expression levels of master genes, but also by the state of their epigenome. Artyomov et al. defined a set of rules on how protein expression can modify the epigenetic state and vice versa during cell cycle progression. Simulating thousands of independent reprogramming experiments, the authors found rare pathways leading to successful reprogramming.

As briefly illustrated by the studies discussed above, a combination of dynamic single-cell expression data and computational models provides insights into the different phases of the reprogramming process. Although an optimisation of the process seems hardly possible within the stochastic phase, targeted activation of pathways or genes in the hierarchical phase can enhance the generation of fully reprogrammed iPSCs. Selecting cells or cell colonies that express predictive markers might help to further increase the proportion of iPSCs, especially because environmental cues and cell-cell communication, which can be manipulated experimentally, play pivotal roles (compare with section on spatio-temporal dynamics).

Whereas stimulating particular components of the pluripotency GRN is crucial for the re-acquisition of a pluripotent cell state, limiting its activity is essential for its maintenance, as shown by studies on mESC differentiation, which we will cover in the following section.

Mimicking the first lineage specification observed in the embryo, mESCs can differentiate into epiblast cells (Epi; giving rise to embryo proper), primitive endoderm (PrE; generating extra-embryonic membranes) and trophectoderm (TE; giving rise to extra-embryonic tissues), depending on the culture environment.

The differentiation into PrE and TE can be forced by ectopic expression of the lineage determinants Gata6 for PrE (Fujikura et al., 2002) and Cdx2 for TE (Beck et al., 1995; Strumpf et al., 2005), but also through changes in the expression of Oct4. Whereas an acute repression of Oct4 leads to differentiation into TE, an increase causes differentiation into PrE (Niwa et al., 2000). Therefore, Oct4 exerts a dose-dependent dual function as both a key factor for pluripotency and a signal driving cells toward lineage specification. However, how these dual functions of Oct4 are regulated remains to be elucidated.

To describe the dose-dependent function of Oct4, Chickarmane and Peterson extended their bistable core network (ODE model; see above) by integrating two mutual antagonisms: one between Oct4 and Cdx2, and one between Nanog and Gata6 (Chickarmane and Peterson, 2008). The simplified GRN, outlined in Fig. 4A, is based on a set of experimental findings reviewed by Niwa (2007). Additionally, the authors suggested a novel motif to account for the non-linear (‘bell-shaped’) relation between stable Nanog states and expression of Oct4, as illustrated in Fig. 4B and described by Matoba et al. (2006). In this motif, Oct4 activates Gata6, and both factors cooperatively suppress Nanog transcription (dashed lines in Fig. 4A). With this assumption, the network model can give rise to three different dynamics, depending on Oct4 expression, which is assumed to be controlled by an external signal A (e.g. Bmp or Wnt/β-catenin or signalling). If Oct4 expression is low, the stem cell switch including Nanog is off and Cdx2 is constantly expressed at high levels (TE-like state; Fig. 4B). Enhanced Oct4 expression initially leads to an increase of Nanog and, thus, to a stable stem cell state. However, at a certain threshold, Oct4 induces a low-level expression of Gata6 and both factors cooperatively repress Nanog transcription, leading to bistability (Fig. 4B, region shaded in grey). In this range, the fate of a cell is determined by the initial concentrations of the pluripotency factors. If Oct4 is further increased or overexpressed, the autoregulatory capacity of Gata6 is switched on and Nanog transcription is efficiently suppressed (PrE-like state; Fig. 4B,C). According to the model prediction that the stem cell state is most ‘stable’ for an intermediate range of Oct4, cellular reprogramming has also been suggested to be most efficient when Oct4 is ‘overexpressed’ within a specific range (Chickarmane et al., 2012). Otherwise, Oct4 expression is either too low to reactivate the pluripotency network or sufficiently high to activate Gata6, which, in both cases, leads to a stable differentiation state. Consequently, the question arises as to how the concentration range of Oct4 is or can be narrowed to efficiently promote pluripotency.

Combining quantitative single-cell analysis and mathematical modelling, Muñoz-Descalzo et al. pointed out that post-translational interactions allow the buffering of the differentiation-inducing activity of Oct4 (Muñoz-Descalzo et al., 2013). They demonstrated that, under self-renewing conditions, Nanog and β-catenin proteins form stable complexes with Oct4, and, thus, restrict its levels to a range that facilitates pluripotency. However, under LIF/serum conditions, Oct4 levels are high enough to cause Nanog bistability, potentially due to the activation of Fgf4/Erk or Gata6. According to the authors, Nanog fluctuations impact the amount of unbound (‘active’) Oct4 proteins that trigger differentiation. In line with this, lowering Oct4 levels in pluripotent mESCs has been shown to reduce Nanog heterogeneity and to establish a more robust pluripotency state (Karwacki-Neisius et al., 2013).

Modelling the mutual transcriptional and post-transcriptional regulation of pluripotency and lineage factors provided possible mechanistic explanations for the dose-dependent activity of Oct4 and its association with the instability of the pluripotent cell state when Nanog expression is low.

So far, we have focused on the analysis of temporal dynamics of intracellular GRNs. However, mESCs do not exist in isolation. Similar to their in vivo counterparts, they self-organise as spheroid aggregates, termed embryoid bodies (EBs), or arrange in colony structures. Embedded in these structures, cells send and receive signals which might affect gene expression and, thus, establish spatial patterns or cell type arrangements.

One example of spatio-temporal patterning can be found in the early phases of mouse development when the inner cell mass (ICM) segregates into Nanog⁺ epiblast cells (Epi) and Gata6⁺ PrE cells [reviewed by Artus and Hadjantonakis (2012)]. After an initial phase of co-expression, both TFs display a mutually exclusive salt-and-pepper pattern (Chazaud et al., 2006) and finally arrange into two layers, with PrE cells occupying the outer layer that faces the blastocoel (Rossant and Tam, 2009). This patterning depends on Fgf4 signalling, which is induced by cells that express Nanog (Yamanaka et al., 2010; Frankenberg et al., 2011).

To reveal a potential mechanism for the establishment of the salt-and-pepper pattern, Bessonnard et al. developed an ODE model describing the interactions between Nanog, Gata6 and the Fgf4/Erk signalling (Fig. 5A) (Bessonnard et al., 2014). First, they demonstrated that the model can account for the coexistence of three stable expression states within a small range of Fgf4 activity (Fig. 5B). These stable states correspond to an ICM-like state (co-expression of Nanog and Gata6), an Epi-like state (Nanog expression) and a PrE-like cell state (Gata6 expression, Fig. 5C). Subsequently, 25 ‘cells’ (i.e. GRNs) were arranged in silico on a squared grid to specify neighbour relationships between individual cells. Cell-cell communication was modelled by averaging the concentrations of Fgf4 produced by a cell itself and by its closest neighbours. This average concentration determined the intracellular activity of Erk signalling. The model assumes an initial heterogeneous distribution of Fgf4, and this leads to a random Nanog versus Gata6 expression pattern, consistent with experimental findings (Chazaud et al., 2006; Plusa et al., 2008). Thus, the model provides one potential mechanistic explanation for the salt-and-pepper cell pattern observed in the ICM. However, it involves a rather large number of arbitrary parameters that had to be confined to a specific range in order to capture the presented dynamics (i.e. tristability) (the so-called ‘over-fitting problem’, discussed further below). Another issue comes with the static nature of the spatial dimension. The ICM is subject to dynamic changes and re-arrangements (e.g. due to cell proliferation or mechanical forces), which impact the emergent pattern. However, these changes have been neglected in this study. Krupinski et al. considered these dynamic properties and combined a three-dimensional (3D) cell-based model with GRNs to simulate both the spatial structure and gene expression levels of a growing blastocyst (Krupinski et al., 2011). Herein, each cell is represented by an incompressible ellipsoid that interacts with its local environment through elastic adhesion forces. Furthermore, each cell is characterised by TF concentrations [determined according to the GRN developed by Chickarmane and Peterson (2008), described above], cell cycle and polarity. By coupling the intrinsic properties of ICM cells with the mechanical forces present in their environment, this model framework allows us to address general, but very important, questions, such as: Does cell position determine gene expression? Do geometric constraints affect spatial expression patterns? Are differential adhesion strengths sufficient to facilitate cell sorting? The ability to assess these questions in silico is particularly significant, as this type of question is still challenging to address experimentally although technologies for in vivo imaging and quantification are constantly improving. This is partly due to the inaccessibility of emerging embryonic structures but also due to the lack of quantitative measures of single-cell properties and spatial patterns.

Recent findings on the impact of cell density and clustering on the self-renewal and differentiation of mESCs (Peerani et al., 2009; van den Brink et al., 2014; Warmflash et al., 2014) emphasise once more that gene expression levels, cellular properties and spatial effects have to be considered in order to achieve a true systemic understanding of stem cell fate decisions. To link these different regulatory layers and to comprehensively analyse the emergent dynamic behaviour, multi-scale mathematical modelling is needed.

At present, live-cell imaging of fluorescently labelled cells is the most suitable method to acquire structural and transcriptional information of cells within an ‘undisturbed’ environment. However, not all cell properties influencing inter- and intracellular regulatory processes, e.g. cell adhesions or mechanics, can be measured or derived directly from imaging data. In this case, multi-scale models can be used to explore systematically the consequences that functional changes at the cellular level might have at the population level. For example, the multi-scale blastocyst-model of Krupinski et al. demonstrates how gene expression, cell proliferation and mechanical properties jointly structure an early embryo. In particular, it revealed that the formation of TE is most robust when assuming a position-based regulation of Cdx2 expression, and that differential adhesion strengths are crucial for the sorting of Epi and PrE cells (Krupinski et al., 2011). Similar relationships between cell position and TF expression have been reported by Herberg et al. (2015). Comparing spatial gene expression patterns of in vitro and in silico mESC colonies on the basis of quantitative measures revealed that cells with a high self-renewing capacity are located in the interior of a colony structure due to TF-related differences in proliferation and adhesion (Herberg et al., 2015). By simulating spatio-temporal patterns of Oct4 expression in EBs, White et al. showed that competing influences between Oct4⁺ and Oct4^– neighbours can account for experimentally observed clusters during pluripotency loss, independently of EB structure, size or cell division (White et al., 2013).

Although there are currently only very few multi-scale models of mESC differentiation available, their insights already indicate the importance of considering not only transcriptional but also cellular interactions in the regulation of pluripotency.

Here, we have discussed a number of different studies, which all apply computational approaches to quantitatively analyse regulatory processes that underlie fate decisions in mESCs and in the blastocyst. Although all these in silico approaches share this aim, the specific questions they aim to answer and the methods applied to do so are quite different. From a general perspective, one can categorise these questions into those aiming to establish a quantitative description of a biological process (i.e. asking “What is going on?”) and those geared towards understanding the underlying mechanisms (i.e. asking “Why does it happen?”) (Mainzer, 2014). Answering the first type of questions calls for a detailed description of the behaviour of the system, without necessarily requiring a deep understanding of the underlying mechanisms. By contrast, addressing the mechanisms that determine the behaviour of the system requires the integration of different levels of information, potentially at the expense of a detailed description of the system.

Although there are methodological overlaps, these two perspectives are usually linked to specific computational approaches (Fig. 6). Typical tools to address ‘What?’-questions are bioinformatics methodologies. The availability of public databases such as GEO (Edgar et al., 2002), FunGenES (Schulz et al., 2009) or ESCAPE (Xu et al., 2013), together with the application of sufficient computing resources, allows us to very quickly answer scientific questions such as: Which genes are differentially expressed across tissue samples? What is the effect of activating or repressing a gene in a defined context? However, one problem of purely data-driven approaches is the sheer abundance of data. Even with efficient high-performance computers, it will (presumably for a long time) not be possible to process and analyse all possible configurations of a biological system. This means that we still need to guide these algorithms by formulating the right questions, implying that we still need to investigate the underlying design principles. To do so, one needs to study interaction rules and regulatory mechanisms that cause dynamic, emergent properties and functions. In terms of Aristotle's characterisation of complexity (“the whole is more than the sum of its parts”), we need to focus on the ‘more’. To this aim, the theory of dynamic systems, including its mathematical framework, is extremely helpful to conceptualise and to describe formally dynamic phenomena. As described above, these techniques require a simplified description of the systems. Whereas these simplifications disregard various details, they nevertheless facilitate the discovery of regulatory principles and mechanisms.

From this, it becomes clear that a systemic, mechanistic understanding of biological systems ultimately requires the combination of a detailed description of the components of the system and of the regulatory principles that rule system dynamics. Thus, data-driven bioinformatics approaches and hypothesis-driven mechanistic modelling approaches should not be seen as competing but as complementary computational disciplines.

However, any kind of mathematical model, regardless of its complexity or formalism, is based on a set of assumptions. Thus, conclusions drawn from a particular model are constrained by these assumptions, and any interpretations neglecting them can, therefore, be misleading. Importantly, this statement is not restricted to mathematical models, but holds true for experimental models, too. It is therefore important to consider the limitations and pitfalls (which are often inherent to the applied modelling strategy) to correctly interpret the results obtained from both experimental and modelling experiments.

Gene expression is often modelled as a binary (i.e. on/off) process, using Boolean network approaches. These models are extremely useful to capture the essence of the regulatory relationships, but they might also miss important details, such as the fact that relationships between TFs and target genes are dose-dependent rather than binary. Closely related to this issue, when interpreting the relationships between network components, one should also consider the properties of the particular measure used to describe the links. For example, the widely used Pearson's correlation coefficient is constructed to detect linear, but not (in general) non-linear relationships. Thus, GRNs with links that have been deduced using this metric might potentially disregard any kind of non-linear interactions.

Another problem, which might be related to the process of adapting a computational model to experimental data, is the so-called over-fitting. It occurs if a model (i.e. its structure and/or its parameters) is adapted to one particular dataset. Indeed, a ‘perfect’ model fit to a particular set of data can always be achieved by a sufficiently large number of model parameters. Thus, it must not be interpreted as the only explanatory model but as one consistent explanation. The reliability of a model increases with the number of independent datasets or experimental conditions that it consistently explains. As different models might potentially be consistent with a certain dataset, it is suggested, according to the parsimony principle [also denoted as ‘Occam's razor’, after William of Ockham (ca. 1287-1347): “pluralitas non est ponenda sine necessitate” (“pluralism must not be presumed except if it is necessary”)], to prefer the simplest one with the fewest assumptions. By contrast, identifying a condition that cannot be described by a given model will disprove the model in its current formulation. This highlights a very important requirement for any model: it has to provide experimentally testable predictions, because model validation strictly depends on the comparison with data.

These remarks do not imply that computational approaches are less reliable than experimental studies. On the contrary, we would like to emphasise explicitly their strengths and potentials. However, to take advantage of their application in an optimal way and to avoid potential misinterpretations, any modelling work should be accompanied by a detailed description and discussion of the assumptions that went into the analysis. We would also like to highlight the necessity of independent confirmations of model-derived hypotheses and predictions.