Novel Generic Models for Differentiating Stem Cells Reveal Oscillatory Mechanisms

Understanding cell fate selection remains a central challenge in developmental biology. We present a class of simple yet biologically-motivated mathematical models for cell differentiation that generically generate oscillations and hence suggest alternatives to the standard framework based on Waddington's epigenetic landscape. The models allow us to suggest two generic dynamical scenarios that describe the differentiation process. In the first scenario gradual variation of a single control parameter is responsible for both entering and exiting the oscillatory regime. In the second scenario two control parameters vary: one responsible for entering, and the other for exiting the oscillatory regime. We analyse the standard repressilator and four variants of it and show the dynamical behaviours associated with each scenario. We present a thorough analysis of the associated bifurcations and argue that gene regulatory networks with these repressilator-like characteristics are promising candidates to describe cell fate selection through an oscillatory process.


Introduction
In 1940, Waddington proposed representing the complex regulatory dynamics driving the process of cellular differentiation as an 'epigenetic landscape' [1] through which a single cell can be thought of as travelling. The differentiating cell encounters successive 'decision points' in the landscape morphology that correspond to differentiation events. These decision points emerge from dynamical changes in the underlying cellular gene regulatory network (GRN), and are mathematically described by variations in parameters that influence the GRN dynamics. For example in zebrafish it is known that Wnt signalling plays a fundamental role in specification and commitment of melanocytes [2][3][4], and the level of Wnt expression may therefore be considered as a parameter, driving the cellular GRN through one or more of the decision points hypothesised by Waddington. Despite the philosophical attractiveness of the epigenetic landscape metaphor, the details have remained unclear and no completely self-consistent mathematical description has emerged. For example there are debates regarding the types of bifurcations characterising differentiation: saddle-node bifurcations [5,6] have been proposed as being more consistent with the expected topologies of the core GRN responsible for differentiation than pitchfork bifurcations. Experimental support for these mathematical bifurcation phenomena, e.g. in embryonic stem cells (ESCs), is also lacking. For example, when progesterone is washed away from mature Xenopus oocytes, they do not dedifferentiate but remain mature [22], a finding that is more consistent with the saddle-node bifurcation than with Waddington's landscape.
Of course, the detailed biological understanding of cell differentiation dynamics has grown significantly recently and it is appreciated that qualitatively different mechanisms may be at work in different contexts. For instance, detailed study of the differentiation of ESCs into two cell types that form either embryonic or extraembryonic tissues indicates that this does not occur in a simple deterministic manner [7]. Instead, cells appear to pass transiently through intermediate states, each seemingly primed to differentiate into one of the possible cell fates [8,9]. Hence ESCs might maintain a collection of meta-stable states, and a better framework to understand and interpret experimental data is required [9].
Most fundamentally, the GRN topologies generating such dynamic stem cell states are poorly understood. In the present work we focus on the question of how broader multipotency (i.e. more than two fates) might generically be generated, in a way that allows differentiation to be biased towards locally-favourable outcomes. We note this is a topic of much current interest in systems biology, particularly noting the exploration of connections between GRN topology and oscillatory dynamics [12], generic links with differentiation dynamics [13][14][15], and e.g. reshaping the epigenetic landscape by rewiring the GRN using cross-repression among transcription factors (TFs) with self-regulation [14].
In this paper, we analyse and compare five minimally-constructed GRNs that admit cyclical dynamics and allow evolution towards fate commitment controlled by an external stimulus. Our GRNs are variants of the standard repressilator, a simple but fundamental circuit made up of three genes repressing each other. Theoretical modelling of this circuit and its subsequent engineering in E. coli have shown its capability to reproduce oscillatory behaviours under broad parameter ranges [17].
We are here interested in adopting this model and the variants proposed to mimick entry into the dynamical oscillatory state characterizing multipotency, and the exit from it, towards a fully differentiated cell state.
Therefore, besides exhibiting oscillations, the system must be equipped with a mechanism to exit the oscillatory regime and proceed towards differentiation. Hence the system must exhibit a sequence of at least two bifurcations, accounting for entering and exiting the oscillatory phase.
We identify two generic scenarios for entering and exiting the oscillatory regime, and we test both in our proposed GRN models. In the first scenario (S1), we hypothesise that a single parameter (i.e. one intervening signalling pathway) drives the system both into and out of the oscillatory regime. An example is Wnt signalling, required for the induction of neural crest, but also for specification of both melanocyte and sensory neuron fates through activation of key TFs (Mitf and Neurogenin, respectively) [3,4,[18][19][20][21][22][23], which we envisage as having this two-fold role. First, increasing Wnt signalling promotes neural crest cells to enter an oscillatory multipotent phase in which fate specification TFs are cyclically expressed; secondly, at higher signalling levels, it promotes the cell's exit from the oscillatory phase allowing specification of single cell fates through increased Mitf or Neurogenin expression which drive melanocyte or sensory neuron differentiation respectively.
In the second scenario (S2) we hypothesise again that a first bifurcation to oscillatory behaviour is driven by an external signalling pathway, but its influence is then blocked above a critical value. Exit from the oscillatory regime results from other signalling pathways, associated with other parameters. This would arise for a multipotent progenitor that in response to a specific signal starts oscillating, as in S1, driving expression of multiple fate specification TFs, each characteristic of a different 'sub-state'. Close to each transcriptional sub-state, the stem cell would express distinct cell signalling receptors; then sufficient activation of these receptors by environmental ligands could provide a mechanism that would move the cell out of the oscillatory phase, thus driving differentiation.
The paper is structured as follows. In Section 2, we consider the simplest GRN exhibiting oscillatory behaviour, the so-called repressilator [17]. We find that the system exhibits a parameter regime where trajectories tend towards an attracting limit cycle featuring a slowdown when passing near three sub-states; this appears as a rather weak effect but is numerically detectable. However, exit from the oscillatory regime leads to a single equilibrium, making this GRN unable to select between alternative differentiated states.
This result motivates Sections 3 and 4 where we extend the standard repressilator to include a second repressive circuit opposing the first; we therefore term this GRN the 'cross-repressilator'. Within this new GRN we consider whether the two TFs act on each gene as either an 'OR gate' or an 'AND gate', leading to two variants of this circuit. We find that for the 'OR gate', the circuit first transitions to oscillatory dynamics, as the standard repressilator, but then followed by a further transition to a stably differentiated cell state, in which only one of the genes is stably expressed. The 'AND gate' is even more interesting, since it allows simultaneous expression of two out of the three genes. In both cases our results appear to hold for a large class of modelling choices for these gene interactions.
In Section 5 we extend the 'AND gate' cross-repressilator circuit to four and five genes. We show that coexpression of up to three genes is possible depending on the topology of the GRN. These dynamics are compatible with both S1 and S2, as they provide an exit mechanism dependent on a single or multiple bifurcation parameters into a number of alternative equilibria. Furthermore, we note that in all the models, when a cell is in the oscillatory phase, scoring for gene expression of fate-specific TFs would result (in a snapshot view) in the cell being considered 'fate specified', despite the fact that it actually retains full multipotency.
A discussion of the biological implications of our models is presented in Section 6.

The standard repressilator
We first consider the repressilator circuit [17], whose GRN is shown in Fig. 1. Here each gene represents a master regulator transcription factor (TFs) of the differentiation process for a specific fate. Assuming fast (un)binding dynamics of the TFs to DNA, and fast mRNA dynamics, we can reduce the system by adiabatic elimination (see Supplementary Materials) [24], and describe its dynamics with ordinary differential equations (ODEs): Here, x i (i = 1, 2, 3) describe the concentrations of the proteins X i , b i are background expression rates, g i are maximal gene expression rates (which combine transcription, translation and mRNA degradation rates), α i are ratios of binding to unbinding rates (association constants) of TFs to DNA, h i are Hill coefficients describing possible cooperative effects (like multimerisation of the TFs [24]), and d i are TF degradation rates. Analysis (see e.g. [17]) shows that the repressilator supports oscillatory behaviour over wide regions of parameter space. Fig. 2 shows (from left) numerical solutions exhibiting the emergence of a family of stable periodic orbits through a Hopf bifurcation as the parameter g 1 increases, time series of sustained oscillations, and the limit cycle in phase space. Oscillatory behaviour exists when h > h * ≈ 2.20498, for the parameter values used in Fig. 2. Figure 2: a) Bifurcation diagram for the repressilator as g 1 varies. A family of stable periodic orbits (olive), shown with maxima and minima of x 1 , emanates from a supercritical Hopf bifurcation at g a and terminates at another supercritical Hopf bifurcation at g b . The stable (solid) branch of equilibria (black) becomes unstable (dashed) between these two bifurcations. The dotted vertical indicates the g 1 -value used for panels (b) and (c). b) Time courses of oscillatory concentrations x 1 , x 2 , x 3 when g 1 = 0.2, g 2 = g 3 = 0.05. c) The solution shown in panel (b) in the phase space (grey) converging to the periodic orbit (olive). The cross (×) is the unique saddle equilibrium, at x eq = (0.58469, 0.87298, 0.45578). Remaining parameter values are: b i = 10 −5 , α i = 50, d i = 0.01, h i = 3, for i = 1, 2, 3.
For h = 3, the oscillations exist between the two Hopf bifurcations at g a ≈ 0.00764 and g b ≈ 15.7919 suggesting that the system is an example of S1 -variations in a single parameter are sufficient to enter and also to exit the oscillatory regime. Our hypothesis for S1 is that external signalling might increase g 1 = g 1 (t) over time, letting the system sequentially explore non-oscillatory, oscillatory, and a new non-oscillatory regime. However, exit from the oscillatory regime at g 1 (t) = g b implies convergence onto the single available equilibrium, making the GRN unable to capture the idea of multiple equilibria necessary to describe stem cell differentiation.
Further, the period of the limit cycles remains bounded as g 1 varies. The shape of the limit cycle ( Fig. 2(c)) shows that the amplitudes of the oscillations for different protein concentrations may differ significantly, and generally appear to spike as the orbit approaches the coordinate axes; two proteins are expressed at relatively low levels while the third is higher. Our simulations show that even when the spiking of gene X 1 is more pronounced, with a huge increase in oscillation amplitude, only a logarithimic or algebraic dependency of the oscillation period appears (see Supplementary Materials), demonstrating the absence of any substantial slow-down around the genes maximal expression substates, and failing thereby to effectively describe any of the scenarios S1 or S2. Furthermore, no mechanism is provided in this model for the stem cell to differentiate into more than one cell type, as no alternative equilibria exist. Thus the standard repressilator does not capture the qualitative changes in behaviour associated with differentiation.

The cross-repressilator with an 'OR gate'
We now introduce a variant of the standard repressilator ('cross-repressilator'), where the inhibitory cycle is accompanied by a second inhibitory cycle arranged in the opposite direction, as depicted in Fig. 3. We assume that regulation of each gene follows an 'OR gate', i.e. each gene is active when both inhibitors are absent (see Supplementary Materials). X 1 X 3 X 2 Figure 3: Illustration of cross-repression among three genes X 1 , X 2 and X 3 formed with two copies of the circuit shown in Fig. 1, oriented in opposing directions.
Making the same assumptions as the standard repressilator, the dynamics of this GRN can be written as where the repression of each gene is described by the multiplication of two decreasing Hill functions to capture the assumed 'OR gate', (see Supplementary Materials). Like system (1), the last term in each equation represents protein degradation. For simplicity, and since it does not cause qualitative changes in the system dynamics, we remove basal expression here and from now throughout this paper. Also, in the numerical simulations we set parameter values equal for each gene, i.e. g i = g, α i = α, β i = β, d i = d, and h ij = h for i, j = 1, 2, 3; we return to this in Section 3.2. Numerical simulations of ODEs (2) are shown in Fig. 4. Figure 4(a) shows the bifurcation diagram for system (2) with respect to the common parameter g, all other parameters fixed. The (dashed) solid curves indicate (un)stable equilibria (black) or (maxima and minima of) periodic orbits (olive green). For g 0.4, system (2) has a single stable equilibrium which is (due to the equalities in parameter values) symmetric: x 1 = x 2 = x 3 . The cyclic symmetry apparent in (2) implies that the Jacobian matrix evaluated at this symmetric equilibrium has a pair of eigenvalues with equal real parts (see Supplementary Materials). When this complex conjugate pair crosses the imaginary axis, a supercritical Hopf bifurcation occurs generating a family of stable periodic orbits. As g increases further, the stable periodic orbits disappear at a global bifurcation involving three new equilibria created at saddle-node bifurcations. Figure 4: a) Bifurcation diagram with respect to g in log-scale for α = 9, β = 0.1, h = 3, and d = 0.2. b) The period of the periodic orbits shown in panel (a) increases to infinity while g tends to the SNIC bifurcation. c) Time courses of the expression level of genes X 1 (red), X 2 (green) and X 3 (blue) showing the sustained oscillations, for g = 1.5. d) For g = 1.6, the system stops oscillating and exits to a differentiated cell.
Moreover, there are two distinct global bifurcation mechanisms that remove the periodic orbits. For α = 9 (shown in Fig. 4), we find that the periodic orbits disappear at a SNIC (Saddle-Node on an Invariant Cycle) bifurcation where three saddle-node bifurcations, related by the cyclic symmetry, occur on the periodic orbit simultaneously. The SNIC bifurcation is well-known in simple dynamical systems, such as the nonlinear pendulum with forcing and damping [25] and in oscillators under external periodic perturbations [26,27]. For g > g SNIC ≈ 1.51449, system (2) has no periodic orbit but three equilibria related by the cyclic symmetry. For each of the equilibria, one of the x i is significantly larger than the other two. Initial conditions (ICs) determine which of these equilibria attract the trajectory. Figure 4(b) shows that, in contrast with the standard repressilator, the period of the periodic orbits diverges as g approaches g SNIC . Typical trajectories spend longer and longer in the vicinity of the three 'slow regions' of phase space where the saddle-node bifurcations are about to occur, and exhibit rapid transitions between these 'slow regions'. This phenomenon has been referred to as the 'ghost of SNIC' [30], and indicates the possibility of S2. We identify these slow regions with the sub-states discussed earlier.
In Fig. 4 (right), we show time courses of system (2) before (c) and after (d) the SNIC bifurcation with same ICs. Figure 4(c) shows the time evolution of the three genes for gg = 1.5 < g SNIC . Whenever one of the genes is highly expressed, the expression levels of the other two is low. In Fig. 4(d), we observe the exit from the oscillatory behaviour associated with a further increase of parameter g (g = 1.6). Here, gene X 1 is stably expressed at a high level, while genes X 2 and X 3 are stably expressed at very low levels.
In summary, for g increasing, the system shows three different behaviours: (i) convergence to a (unique) stable equilibrium; (ii) oscillations among three sub-states; and (iii) attraction to one of the stable equilibria. The natural interpretation, compatible with S1, is to see the first two cases as reflecting multipotency within the stem cell, while the third case represents attainment of one of multiple differentiated cell types.
The phenotype attained in the differentiated state depends on the IC used for the simulation and its intrinsic properties, here illustrated through different g-values, as shown in Fig. 5. Figure 5(a) colour-codes the (x 1 , x 2 )plane according to the final equilibrium for three different ICs for x 3 (t), and fixed g. Figure 5(b) shows the effect of varying parameter g on the same colourmap when x 3 (0) is fixed. The three coloured regions in the (x 1 , x 2 )-plane twists around a central point. Increasing x 3 (0) tightens the twists of these spiralling regions while the opposite occurs when g increases. Moreover, variations in g keep the central twisting point fixed at (x 1 (0), x 2 (0)) = (1, 1) whereas increasing x 3 (0) makes this point drift away. To obtain insight into the feasibility of S2 for this system, we explore the bifurcation behaviour in the (g, α)plane. Fig. 6(a) shows the continuation of the bifurcations identified for α = 9 in the (g, α)-plane. When α decreases, the location of the Hopf bifurcation curve (HB) moves to larger g whereas the curve of saddle-node bifurcations (SN) on which the SNIC bifurcation occurs moves to lower values of g. At α = α * ≈ 3, the SNIC regime terminates at a codimension-two bifurcation, indicated by a grey dot. Two new bifurcation curves emerge from this point, indicated by the labels HC (homoclinic) and SN (saddle-node). For α < α * the three saddle-node bifurcations (due to symmetry) occur away from the periodic orbit in phase space, resulting in a region of the (g, α)-plane in which these new equilibria (three stable and three saddle equilibria) coexist with the stable periodic orbits. The periodic orbits then collide with the three saddle equilibria at a global bifurcation and disappear. As in the SNIC bifurcation, the period of the periodic orbits goes to infinity as the orbits approach the global bifurcation, but the mathematical details of the behaviour differ from the SNIC case, as shown for example by the scaling of the period with distance to the bifurcation point being different (see Supplementary Materials).
We observe also that as α decreases further the heteroclinic bifurcation (HC) and SN curves move further apart from each other in the (g, α)-plane, with maximum separation at α ≈ 0.1.
For α < α * the HB (green) and HC curves (red) determine the left and right boundaries of the oscillatory region, respectively. Also note that the HB and HC curves meet and are tangent to each other at α = 0.1 which is not generic but arises here due to the additional symmetry existing when α = 0.1 = β. When α = β, ODEs (2) are symmetric under two independent symmetry operations: the cyclic permutation symmetry (x 1 , x 2 , x 3 ) → (x 2 , x 3 , x 1 ) and the interchange of any two of the three variables, e.g. (x 1 , x 2 , x 3 ) → (x 2 , x 1 , x 3 ). In group theory notation, the symmetry group is now the symmetries of an equilateral triangle D 3 rather than just the cyclic group of order three, Z 3 . This additional symmetry implies the tangency between the two curves. The analysis around this point will be discussed in detail in a separate paper [32] and has been observed in other contexts [28,29]. Figure 6(b) shows the bifurcation diagram, varying g, for α = 1 corresponding to the lower horizontal dashed line indicated in Fig. 6(a). Here we observe an interval of g-values in which the stable periodic orbit coexists with the three stable equilibria. Consequently we would expect that this transition to a stable equilibrium occurs rather SN HB HC Figure 6: a) Two-parameter continuation of saddle-node (black), Hopf (dark green) and heteroclinic (red) bifurcations, shown in Fig. 4(a), in the (g, α)-plane; vertical axis (α) is shown in log-scale. The dotted and dashed lines indicate the corresponding α values used in Fig. 4(a) (α = 9) and in panel (b) (α = 1), respectively. b) Bifurcation diagram of system (2) with respect to g in log-scale when α = 1. The diagram structure is mainly similar to the one in Fig. 4(a). However, here, a small region of multi-stability is created between the oscillatory regime and the stable equilibria such that the family of periodic orbits terminates at a heteroclinic bifurcation (instead of a SNIC).
smoothly since the sub-states already exist prior to the termination of the oscillations.
These observations support S2. Once the system has entered the oscillatory regime by tuning of g through a first signalling pathway, a second signalling pathway might intervene to lower α, stop the oscillations, and drive the system to differentiation. This scenario requires that the change of α happens faster than the permanence of the system close to a selected sub-state. This slowdown of the cycling dynamics is therefore essential for this mechanism.
Thus, after establishment of the stem cell, increasing the first signal initiates the cyclical expression of fatespecific TFs corresponding to the different differentiated cell-types. Subsequently, in each primed sub-state, as a result of TF specific activation of sensitivity to fate specification signals, the cell is sensitised to local signals. When a cell receives such a signal for sufficient time, this shifts α and removes the oscillations, thereby initiating differentiation.

Illustrating the two scenarios with time-dependent parameters
So far we have assumed that all the parameters of system (2) are fixed and do not vary in time. In this section, we explore more directly the two scenarios S1 and S2 by assuming that relevant parameters are regulated in time by an external signalling pathway. We will consider the case when either one or both of g and α are time dependent. In S1, parameter g only is responsible for shifting the cell from the multipotent regime into the oscillatory regime and then into the differentiated cell. Accordingly, we assume that g(t) follows the form of an increasing Hill function, from 0 and saturating atĝ > g SNIC as where we setĝ = 1.7, p = 1.25 and τ = 1000. The bottom panel of Fig. 7(a) illustrates the time-variation of g(t). The top panel in Fig. 7(a) shows how the system responds to g(t), moving through the three different regimes as g(t) increases and tends to the valueĝ (compare with Fig. 4(a)). In detail, the system first settles quickly from the IC to the single stable equilibrium and follows it as g(t) slowly approaches the Hopf bifurcation. After Hopf we observe the onset of rapid but finite-frequency oscillations. As g −→ g SNIC , the period of oscillations increases significantly, giving rise to dynamical phases close to each of the three sub-states. Finally, when g(t) crosses g SNIC , the oscillations disappear and the system settles at one of the three equilibria, each corresponding to a different expressed gene. For example, in Fig. 7, gene X 3 (blue) is expressed and the other two genes remain at a low level.
As discussed for S1, the precise state of the system when g goes through g SNIC , as well as the value of the saturation levelĝ play crucial roles in selecting cell fate. The shape of the profile of g(t) is also relevant for this selection.  Fig. 4(c). b) Development from a stable oscillatory regime into a differentiated cell according to S2, when the secondary parameter α starts decreasing att = 0 (shown in the lower panel). Here g = 1.5 is held constant, and all other parameters take the same values as in Fig. 4(c).
For S2, we hypothesise that two parameters, g and α, are involved in cell differentiation. As shown in Fig. 7(b), we first assume that parameter g shifts the system into the oscillatory regime while α remains fixed. At a certain time, we suppose a second pathway causes a decrease in α, modelled as where H is a Heaviside step function (i.e. H(s) = 0 when s < 0 and H(s) = 1 when s > 0),t is the time when the time dependency of α(t) starts, and the constants τ = 1000 and p = 1 control the profile shape.
S2 clearly allows more control over the final fate of the cell through the increased number of coefficients involved in specifying the time-dependency of the second bifurcation parameter. To understand these further dependencies, we analyse the effect of variations in the switch-on timet and the time scale parameter τ for the saturation of α(t). Shown in Fig. 8, we colour-code the parameter plane (τ,t) based on the expressed gene for g = 1.5 (a) and g = 1.51 (b). The (τ,t)-plane is divided into stripes which periodically alternate between red, green and blue. The lower value of g in Fig. 8(a) corresponds to narrower stripes as a result of the smaller oscillation period, with the system being more sensitive to changes in switch-on time and saturation time scale. Further, more horizontal stripes in Fig. 8(b) indicate that cell fates become less sensitive to τ for this larger value of g, at least in the regime close to the global bifurcation.

Breaking the symmetry
Until now, we assumed that system (2) is symmetric and all the association rates α i and β j , i, j = 1, 2, 3 for the TFs X 1 , X 2 and X 3 are equal. Here, we break the symmetry in the reciprocal repressing loops in Fig. 3 and consider the case that α i and β j take different values, as specified in Fig. 9. Note that we keep the remaining parameters g and d, unchanged. Figure 9(a) shows the bifurcation diagram with respect to g in log-scale. The structure remains qualitatively unchanged compared to the symmetric case. For g small, a unique stable equilibrium exists, which undergoes a Hopf bifurcation and becomes unstable. For g large, there are six branches of equilibria; the stable branches merge with the saddle-type ones and disappear as g decreases. The family of stable periodic orbits emerging from the Hopf bifurcation terminates at a SNIC bifurcation. However, here, only one equilibrium lies on the invariant cycle contrary to the symmetric case, with three equilibria on the invariant cycle. This happens because the saddlenode bifurcations, at which the six branches of equilibria appear, occur at different g-values (they are no longer symmetry-related), indicated by the vertical dotted lines.   Fig. 9(a). As before, each TF is expressed cyclically when the other two are at low levels. However, gene X 2 remains expressed for a significantly longer time than the other two, particularly X 3 . Due to the occurrence of the saddle-node bifurcation with high level of X 2 at smaller values of g, the slow region created around the bifurcation has a stronger effect which leads to a longer 'dwelling time' at expressed X 2 . In fact, the dwelling time for each expressed gene depends on the distance of g-value at which the corresponding saddle-node bifurcation (indicated with a dotted vertical line of the same colour) occurs.
Analysis of Fig. 9 determines that with the current parameters, the conditions favour gene X 2 ; therefore, when g crosses to the right of the SNIC bifurcation, the system settles at a stable equilibrium where X 2 is expressed. Changing this condition results in a shift in the location of the saddle-node bifurcations along the horizontal axis and might alter the order of saddle-node bifurcations; consequently, lengthening the dwelling time for another gene.

The cross-repressilator with an 'AND gate'
The third model we consider here is similar to system (2) with the difference that each gene is repressed only if all other genes are highly expressed (a so-called 'AND gate').

This model is described by
where again the decreasing Hill functions describe the effect of repressing genes (see Supplementary Materials). We also set g ij = g j for i, j = 1, . . . , 3 and d 1 = d 2 = d 3 = d. We also fix α i = β i = 1 and the Hill function exponents h ij = 3 for i, j = 1, 2, 3.
As previously, we explore system (3) numerically and compute the bifurcation diagram, varying g 1 while keeping g 2 and g 3 fixed (and unequal). Figure 10 (a) shows a bifurcation diagram for system (3) with respect to g 1 . For small values of g 1 , the system has a single stable equilibrium (black solid line) which undergoes a supercritical Hopf bifurcation and becomes unstable (dashed line). The stable periodic orbits (olive) are indicated by the maximum and minimum values of x 1 . The periodic orbit, similar to the 'OR gate' case, terminates at a SNIC bifurcation where it collides with a saddle-node bifurcation. Note that these branches are part of a single branch of equilibria terminating at another oscillatory regime at very large parameter values g 1 (not shown). Panels (b) and (c) in Fig. 10 show time courses of the expression level of each gene in system (3) for values of g 1 mentioned in the caption. Figure 10: : a) Bifurcation diagram of system (4) with respect to g 1 for g 2 = 0.76, g 3 = 1.9 and d = 0.3. b) Time courses of the expression level of genes X 1 (red), X 2 (green) and X 3 (blue) showing the emergence of an intermediate expression level in the oscillatory regime of system (3) when g 1 = 2.8. c) Gene expression levels converging to equilibrium in the tristable regime for g 1 = 3.
Clearly, the behaviour of the 'AND gate' model follows closely the analysis presented for system (2), but with one important and biologically relevant difference. Compared to system (2), in which only one gene could be highly expressed and the other two genes were necessarily at very low levels, in system (3), two genes can be expressed at high or moderate levels simultaneously. Of course, the precise levels depend on parameter values, but the 'AND gate' structure is crucial in admitting states with more than one high expression level. The form of the periodic orbit also changes and shows that it lingers in sub-states where two gene expression levels are nonzero. This notable feature is potentially useful in matching models with quantitative experimental data where it allows one to select a model that allows or prevents the existence of phenotypes dependent upon the expression of a single or a combination of master regulator TFs.

Extension to four and five genes
So far, we considered GRNs with only three genes. We now investigate the dynamics with more than three genes. For the same cross-repressilator in 'AND gate' configuration discussed in the previous section, an odd number of genes is necessary to create a negative feedback loop that supports sustained oscillations [31]. Hence a trivial extension of the cross-repressilator to four genes, preserving the cyclical topology, fails to reproduce oscillatory dynamics, and therefore we do not consider it here. Even numbers of genes in the cycle result in states where the gene expression levels alternate around the loop. However, adding cross connections to this four-gene network so as to build three-gene negative loops, embedded within the four-gene GRN as shown in Fig. 11(a), recovers the oscillatory dynamics. For five genes, the simple cyclical topology of the cross-repressilator, with only 'nearest neighbour' interactions and without additional cross-regulation, shown in Fig. 11(b), maintains oscillatory behaviour. Figure 11: a) Extension of the cross-repressilator GRN to four genes with cross connections that support oscillations. b) A cyclical arrangement of five genes, sufficient to drive oscillations without additional cross connections.

The four gene network
Here, we describe the dynamics of a GRN with four genes, as illustrated in Fig. 11(a). The network is fully connected such that each gene is repressed by the other three genes, operating in 'AND gate' configuration. The results in (see Supplementary Materials for details): For simplicity, we again set g ij = g j and d i = d for i, j = 1, . . . , 4. We also keep fixed the values g j = 0.5 for j = 5, . . . , 7, α i = β i = γ i = 1 and h ij = 3 for i, j = 1, . . . , 4 (i = j). Figure 12(a) shows the bifurcation diagram of system (4) with respect to g 1 , together with time courses of the expression levels of all genes in the system for two specific values of parameter g 1 in the multipotent (b) and differentiated (c) regimes. In Fig. 12(a), a single branch of stable equilibria (solid; black) undergoes a Hopf bifurcation and becomes unstable/saddle (dashed). The envelope of stable periodic orbits (olive) emanating from the Hopf bifurcation terminates at a SNIC bifurcation where new branches of equilibria appear. Four of these new branches are stable, and collide with the other four unstable branches, disappearing at saddle-node bifurcations.
Panels (b) and (c) in Fig. 12 show time series of the expression levels of the genes in the oscillatory and differentiated regimes for g 1 = 3 and g 1 = 3.6, respectively. One notable, and biologically relevant, feature of this model is that all four genes can be expressed, with one gene highly expressed, two genes expressed at intermediate levels, and one gene only slightly expressed. This scenario is consistent with differentiation reflecting a combinatorial role for multiple TFs. The relative concentrations of the expressed genes would open up possibilities for differentiation to multiple cell types. Figure 12: a) Bifurcation diagram of system (4) with respect to g 1 for g 2 = 1.6, g 3 = 0.5, g 4 = 1.5 and d = 0.4. (b) and (c) Time courses of the expression level of genes X 1 (red), X 2 (green), X 3 (blue) and X 4 (magenta) with two different behaviours: in the multipotent oscillatory regime of system (4) when g 1 = 3 (b) and in the multistable regime for g 1 = 3.6 (c).

The five gene network
We now consider the five-gene network depicted in Fig. 11(b). The network comprises two inhibitory loops arranged in opposing directions so that gene X i is repressed by genes X i−1 and X i+1 , again in 'AND gate' configuration. Therefore, similar to the three-gene case, the dynamical equations are where the maximal expression rates imposed by the anticlockwise and the clockwise loops respectively equal to g 1 , g 2 and d i to d. We also fix the values α i = β i = γ i = 1 and h ij = 3 for i, j = 1, . . . , 5 (i = j). Figure 13(a) shows the bifurcation diagram remains unchanged with respect to g 1 except for the number of stable (solid; black) and saddle (dashed; black) branches in the multistable regime. Figure 13(b) exhibits time courses of gene expression levels when g 1 = 2. Figure 13(c) shows gene concentrations for g 1 = 2.5 when the cell differentiates. In the differentiated regime, the number of expressed genes increases to three compared with the three-gene case. We conjecture that more generally for a cyclic network of n genes, n+1 2 genes are expressed at high or moderate (but unequal) levels and n−1 2 are either not expressed at all or are only weakly expressed. Biologically, the stable generation of cells expressing mixtures of key fate specification TFs shown by these fourand five-gene models matches closely the concept of combinatorial TF expression defining similar cell-types. This is most keenly demonstrated in the peripheral and central nervous systems where distinct neuron types express mixtures of TFs that together specify their individual fate [33]. It may be that the fully differentiated phenotype represents the impact of a mix of TFs expressed at different levels, even where we tend to think of a single dominant master regulator, such as Mitf in melanocyte development.

Conclusions
We have presented a detailed bifurcation analysis of four variants of the standard repressilator model, and discussed their suitability as novel generic models for the process of cellular differentiation. While the original repressilator fails to provide convincing dynamical behaviour, we find that the variants proposed here exhibit biologically relevant features. Our modelling moves beyond Waddington's vision of the epigenetic landscape since such a picture cannot capture highly dynamical internal states of cells which our model places centre stage.
Bifurcation analysis of the systems using [16] reveals that the oscillatory behaviour emerges from a Hopf bifurcation. At the late stages of multipotency, the oscillating cell visits a set of sub-states where the relevant fate-specific transcription factor(s) are highly expressed relative to the other gene(s). The oscillations terminate at a global bifurcation and the cell settles at a state resembling that of a committed cell-type. From a biological perspective, we can think of the individual fate-specific transcription factors (TFs) as defining individual cell-fates. Where they oscillate, the stem cell retains multipotency (it sequentially transits nearby all states), but at the same time it becomes periodically biased to adopt individual fates due, we propose, to changing exposure to fate-specification signals. When such a signal is received at a sufficiently high level, it may force oscillations to cease and a single TF becomes stably expressed at the expense of the others --cell differentiation has begun. We have focused on two possible scenarios for driving the oscillatory behaviour. In the first scenario, the same pathway is responsible for initiating and terminating the oscillations, by the increase of the chosen parameter. In the second scenario, a first pathway brings the system in the oscillatory region, and a second one intervenes to drive the exit from it. Our results show that the second scenario, however, gives more control over cell fate. Moreover, the timing of controlling signals is essential in our framework, since this combines with ICs to determine the specific TF that remains stably expressed, and thus the differentiated cell-type.
Depending on the type of gates assumed in the network, one or two genes (or more in networks of higher dimensionality) can be co-expressed during both the oscillatory phase and in the selected differentiated cell type. During the oscillations, combinations of genes may be expressed transiently --these might correspond to partiallyrestricted cell intermediates (expressing a subset of fate-specific TFs); however, we emphasise that in our scenarios these retain oscillations and hence must be thought of as fully multipotent, despite their 'snapshot appearance'. In the selected differentiated cell-type, the combination and the levels of expressed TFs vary with ICs and timing of the fate-specification factor; thus different cell-types, with specific quantitatively distinguishable combinations of TFs, are generated. Such a model now provides a view of how multipotency might be exhibited in stem cells, in contrast to the series of bipotential intermediates that underpin standard thinking in development. Note that the external signals drive a fully multipotent cell to a specific committed state directly, without further intermediates with restricted potency. Finally, the case with which stable (committed) states consisting of mixed expression of two or more TFs can be generated may well relate to the mechanism producing related cell-types, e.g. neuron types in both the central and peripheral nervous system.
There is a wider context for the results presented here; the identification of dynamical scenarios typical for repressilator-type GRNs, particularly in view of their proposed role in clock circuits underlying global patterning processes in development [34]. Although the developmental context considered there is different, the formulation and behaviour of the mathematical models has some points of similarity with our analysis, such as the occurrence of Hopf and SNIC bifurcations, which are generic bifurcations involving time-periodic oscillations.
A key difference with our work is ordering of the different regimes. Model 1 of [34] shows a progression from oscillations to a single stable equilibrium, and then a saddle-node bifurcation that creates multiple equilibria. Model 2 of [34] has no Hopf bifurcation but instead a SNIC separating oscillatory from multistable regime; there is no single equilibrium. In contrast, our model robustly indicates a progression from a single equilibrium to oscillations via a Hopf bifurcation, with the oscillations then terminating at a global bifurcation where multiple equilibria emerge. We note that the distinction between the SNIC mechanism and the combination of the global and saddle-node bifurcations is likely to be extremely difficult to detect experimentally due to the high level of resolution required in time, and the effects of noise. It is however mathematically important in terms of comparisons between models.
In terms of model construction a key point of difference is that both Model 1 and Model 2 in [34] are constructed by interpolating, rather artificially, between two mechanisms potentially driven by distinct and unrelated morphogens. Our dynamical model is structurally more robust and parsimonious as our work has shown that the transition from equilibrium to oscillations whose frequency drops to zero, to 'differentiated' equilibria, is a natural result of a single GRN mechanism. Naturally, spatial patterning is not part of our model since we focus on the dynamics within a single cell, but establishing connections between temporal dynamics and spatial patterning is an intriguing research direction that we will explore in a future paper.

Novel Generic Oscillatory Mechanisms in Models for
Differentiating Stem Cells

S1 Model Derivation
Following [1], we here provide a derivation of the Hill functions dynamics utilised in the construction of the models in this paper.

S1.1 Regulation by repressive transcription factors
Many processes are involved in the regulation of a protein. Here, we consider only few of them and assume that the other processes such as transcription and mRNA translation are implicit. The model includes (un)binding of the transcription factor (TF) R, degradation of the protein P, and synthesis of P when gene G is in the active state, denoted G * , as follows.
where e R and r R represent the 'empty' and 'regulated' regulatory elements of G, respectively. The parameters over/under the arrows indicate the rate of the processes. If [·] denote the number of a chemical per cell, the processes in (S1) are described with Binding and unbinding of R to the regulatory element is significantly faster than the protein regulation and degradation (α 0 , α 1 ≫ g, d P ). By quasi-steady-state approximation d[r R ]/dt = 0, we have where α = α o /α 1 .
Here, the TF R regulates the protein P as a monomer giving rise to a Hill coefficient one. Hill functions with higher coefficients can be derived when regulation proceeds through binding/unbinding of complexes (dimers, trimers, etc).

S1.2 Regulation by two repressors arranged in an 'OR gate'
Now, we consider regulation of gene G by two independent and repressive TFs, A and B. In an OR gate configuration, the presence of one TF is sufficient for the repression of G. Similar to the previous section, we have with steady states and . (S5) where α = α 0 /α 1 and β = β 0 /β 1 . Because A and B are repressing and act independently, gene G remains active when both binding sites are empty. Accordingly, Therefore, the dynamics of the protein P is described by

S1.3 Regulation by two repressors arranged in an 'AND gate'
In an AND gate configuration, the presence of both TFs A and B is necessary for the repression of G. Accordingly, gene G is active when at least one of the binding sites is empty. In this case, the average number of activated genes is where [r A ] and [r B ] are given by equation (S5). For generality, we here assume that the production rate of protein P depends on the configuration of bound TF [2]. Therefore, by assuming production rates g 1 , g 2 and g 3 for [G * 1 ], [G * 2 ] and [G * 3 ], respectively, the dynamics of protein P is described by This can be rewritten as which is the form used in (3) and (5).

S1.4 Regulation by three repressors arranged in an 'AND gate'
Similarly to the AND gate with two repressors, derived above, the average number of activated genes for the AND gate with three repressors is which is the form used in (4).

S2 Analysis of the oscillation period in the standard repressilator
Here we show a numerical exploration of the dependency of the period of oscillations of the standard repressilator on the two significant parameters in the model. We restrict attention to the fully-symmetric version for simplicity: g 1 = g 2 = g 3 , d 1 = d 2 = d 3 , α 1 = α 2 = α 3 and b 1 = b 2 = b 3 in all cases. We consider separately the variation in the period T of the oscillation with changes in g 1 and in d 1 ; these parameters are as defined in (1). Our simulations are presented in Fig. S1.  Figure S1: Fit of the oscillations period dependency on g 1 and d 1 parameters in the repressilator model. In both panels the points (blue squares) represent the result of the simulations for d 1 = d 2 = d 3 = 0.01 fixed (left) and g 1 = g 2 = g 3 = 0.05 fixed (right), while the solid lines correspond to best fitting curves as detailed in the text. The resulting fitting parameters are a 0 = 1300, a 1 = 202, a 2 = 1.2, and β = 4/3. In all cases b 1 = b 2 = b 3 = 10 −5 and The fitting functions assumed in these cases are as follows: T(g 1 ) = a 0 + a 1 log(g 1 ) for d 1 fixed (S13) and for g 1 fixed (S14) The simulations show a very weak dependency of the period on the parameter g 1 -over seven orders of magnitude in g 1 the period increases by a factor of only around seven. As d 1 increases the period drops much more rapidly, in line with the power-law scaling T ∝ d −4/3

1
. This is to be expected since the degradation rate d 1 is a key timescale in the problem and the decay rate of the concentrations of TFs is essential to the nature of the 'repressilator' dynamics. In neither case does the qualitative nature of the periodic oscillation change over these ranges of parameters. This allows us to conclude that the parameter values used, for example in figure 2 in the main text, are entirely typical of the repressilator's dynamical behaviour. Finally we note that the dependence of the period and shape of the oscillation on the production rates b 1 , b 2 and b 3 is very weak in this parameter range where typical concentrations of TFs are maintained well above each of the b i by the nonlinear prodution terms.

S3 Analysis of mathematical concepts in the cross-repressilator
Three mathematical concepts that need to be explained and analysed more closely are symmetric structure of systems (2) and (3), the growth rate of period of limit cycles close to SNIC/heteroclinic bifurcation and determining factors in the order of expressed genes when the system oscillates.

S3.1 Analysis of the three-gene Models in the symmetric case
In the case that the models (2) and (3) are cyclically symmetric under the permutation (x 1 , x 2 , x 3 ) → (x 2 , x 3 , x 1 ), i.e. when the parameters g, α, β, k and h do not take different values in each equation, the symmetric structure of the problem enables some semi-analytic progress to be made since the Jacobian matrix at the symmetric equilibrium point (x, x, x) is constrained to always take the form Figure S2: Growth scale of periodic orbits. Period of the family of periodic orbits ending at a SNIC bifurcation for α = 9 (left) and heteroclinic bifurcation for α = 1 (right). The rest of parameters of the system are the same as Fig. 4(b).
the saddle point equilibrium that collides with the periodic orbit (in the case that the periodic orbit is stable, which it is here). The factor of 3 arises due to the Z 3 symmetry and the three distinct saddle points that all collide with the periodic orbit simultaneously. Numerically, we find that λ + ≈ 0.05922 (5 dp); therefore 1/λ + = 16.88619 (5 dp). Hence the theoretical estimate for the slope of the straight line is 3 × 16.88619 = 50.65857 (5 dp). From the right panel of Fig. S2 we estimate the slope of the curve to be 50.17112 based on the values of the period T at |g − g HC | = 3 × 10 −7 and |g − g HC | ≈ 3 × 10 −4 . This therefore provides quantitative confirmation of the bifurcation behaviour in the differential equations.

S3.3 Ordering of the gene expression in the oscillations
As discussed in the main text, the genes are expressed in a particular order in the oscillatory regimes. For instance, in the three-gene case with the OR gate configuration, the gene expression occurs cyclically clockwise. That is, specifically, first gene X 3 is expressed and then genes X 2 and finally X 1 . This order of gene expression is due to the values of the parameters of (2). In (2), we chose α = 9 > β = 0.1 which makes the inhibitions of the inner loop in Fig. 3 stronger; this results in the genes being expressed in clockwise order. However, if we set α < β, the order of gene expression in the oscillatory regime is reversed.
Mathematically the oscillatory direction of the family of periodic orbits emanating from a (supercritical) Hopf bifurcation is determined by the relative signs of the real and imaginary parts of the eigenvectors corresponding to the pair of complex conjugate eigenvalues (λ 2,3 = Re(λ 2,3 ) ± Im(λ 2,3 )). Since these eigenvectors remain non-zero and vary continuously along the Hopf bifurcation curve on which these eigenvalues have zero real part, we can conclude that the oscillation changes direction when the imaginary parts of the eigenvalues pass through zero. Figure S3 shows the imaginary parts of λ 2,3 and its conjugate −λ 2,3 along the curve of Hopf bifurcations (dark green) in Fig. 6 projected on g and α (in log scale) axes. For small values of g (large values of α), let assume that Im(λ 2,3 ) is positive (blue); therefore, its conjugate, −Im(λ 2,3 ), is negative (red). As g increases (α decreases), Im(λ 2,3 ) and its conjugate, −Im(λ 2,3 ), cross each other. Thus, the direction of oscillation of periodic orbits switches. The crossing point of Im(λ 2,3 ) and −Im(λ 2,3 ) at 0 corresponds to the point that the Hopf curve (HB) touches the homoclinic curve (HC, red). Precisely at the parameter value where the complex conjugate pair of eigenvalues are zero, the Hopf bifurcation theorem breaks down and further analysis is required; this discussion will be presented in future work [3].

S4 Numerical Methods
We have employed different tools and software packages for computing solutions, bifurcation diagrams and colourmaps. All the time series have been integrated in XPPAUT, freely available at http://www.math.pitt.edu/~bard/xpp/xp Figure S3: Imaginary parts of eigenvalues along the Hopf curve. The variation in the imaginary parts of the complex conjugate eigenvalues along the Hopf bifurcation curve (HB, dark green) in Fig. 6 shown versus g (left panel) and α in log scale (right panel). Other parameter values are as in Fig. 6.
The bifurcation diagrams are computed using pseudo-arclength continuation methods in software package AUTO [4,5].
Finally, the colourmaps are computed in MATLAB ® (Mathworks Inc.) starting from different initial conditions. In Fig. 5, the initial-condition planes are colour-coded depending on the steady state (one of the three) that the system converges to. Similarly, in Fig. 8, the final state of the system when g increases and saturates at a value larger than g SNIC determines the colour of the initial conditions.