 Research article
 Open
 Published:
Classification of transient behaviours in a timedependent toggle switch model
BMC Systems Biologyvolume 8, Article number: 43 (2014)
Abstract
Background
Waddington’s epigenetic landscape is an intuitive metaphor for the developmental and evolutionary potential of biological regulatory processes. It emphasises timedependence and transient behaviour. Nowadays, we can derive this landscape by modelling a specific regulatory network as a dynamical system and calculating its socalled potential surface. In this sense, potential surfaces are the mathematical equivalent of the Waddingtonian landscape metaphor. In order to fully capture the timedependent (nonautonomous) transient behaviour of biological processes, we must be able to characterise potential landscapes and how they change over time. However, currently available mathematical tools focus on the asymptotic (steadystate) behaviour of autonomous dynamical systems, which restricts how biological systems are studied.
Results
We present a pragmatic first step towards a methodology for dealing with transient behaviours in nonautonomous systems. We propose a classification scheme for different kinds of such dynamics based on the simulation of a simple genetic toggleswitch model with timevariable parameters. For this lowdimensional system, we can calculate and explicitly visualise numerical approximations to the potential landscape. Focussing on transient dynamics in nonautonomous systems reveals a range of interesting and biologically relevant behaviours that would be missed in steadystate analyses of autonomous systems. Our simulationbased approach allows us to identify four qualitatively different kinds of dynamics: transitions, pursuits, and two kinds of captures. We describe these in detail, and illustrate the usefulness of our classification scheme by providing a number of examples that demonstrate how it can be employed to gain specific mechanistic insights into the dynamics of gene regulation.
Conclusions
The practical aim of our proposed classification scheme is to make the analysis of explicitly timedependent transient behaviour tractable, and to encourage the wider use of nonautonomous models in systems biology. Our method is applicable to a large class of biological processes.
Background
Development in wildtype organisms is robust to genetic and environmental variations. Conrad Hal Waddington introduced the notion of ‘canalisation’ to describe how developmental processes resist perturbations during embryogenesis [1–3]. The canalised nature of development explains, he argued, why most phenotypes are discrete and distinct. To illustrate these ideas, he developed the epigenetic landscape, one of his most wellknown concepts [3, 4].
In Waddington’s epigenetic landscape, the current state of a developing system is indicated by a ball on an undulated surface (Figure 1A, top panel) [3, 5, 6]. The topography of this landscape determines the developmental potential or repertoire of the system. The topmost edge of the surface shown in Figure 1A (top panel) represents the initial state of the system given by, for example, a particular set of initial protein concentrations in a cell. Valleys in the landscape symbolise the various differentiation pathways that are available. The landscape’s topography—together with the initial state—determine a developmental trajectory that follows a particular valley. The structure of the landscape is such that, if the system is slightly perturbed, the sloping valley walls will cause it to correct and readjust its trajectory. This behavour is called ‘homeorhesis’—the maintenance of a dynamic trajectory—in analogy to the more static concept of homeostasis—the maintenance of a (steady) state of the system [3]. The wider and deeper a valley is, the more canalised the developmental trajectory. Waddington named such canalised trajectories ‘chreodes’.
An additional feature of Waddington’s landscape is crucial in our context: it is the fact that its surface is not necessarily fixed, but can change over time due to the influence of genetic or environmental signals. Waddington represented this by pegs connected to the underside of the landscape by ropes (Figure 1A, lower panel) [3]. As the genetic or environmental context of the system changes, these ropes pull and stretch the surface, thereby changing its topography and hence its developmental potential.
Waddington’s epigenetic landscape was conceived as a conceptual tool to illustrate the nature of developmental robustness and its effect on evolutionary dynamics [1, 3–5, 7]. As such, it remained at a rather metaphorical level of explanation, since complex nonlinear biological processes are hard to formulate and analyse [5, 6, 8–12]. This is still the predominant way in which it is used in various reviews on contemporary stem cell research (see for example, [13–17]). However, in order to understand the precise nature of specific developmental trajectories or chreodes in real systems, we have to take Waddington’s landscape a step further: we have to calculate it based on experimental evidence, and use it to characterise the transient behaviours that govern the observed developmental dynamics [18–25].
The increasing availability of quantitative gene expression data renders this approach feasible. However, before we can successfully apply it, we also require new conceptual and mathematical tools to deal with the analysis of datadriven models that are formulated in terms of nonautonomous (i. e. explicitly timedependent) dynamical systems. Explicit timedependence is necessary to reproduce the changing topography of Waddington’s landscape. However, such systems are difficult to study in a rigorous mathematical manner, and few analysis tools exist at this point. In this work, we address this challenge by proposing a classification scheme for transient dynamic behaviours observed in a nonautonomous version of a simple gene network model. This scheme is meant to provide the foundation for the analysis of more complex timedependent models that reproduce the dynamics of specific, experimentally tractable, biological regulatory systems.
Dynamical systems
In this study, we focus on dynamical systems formulated in terms of ordinary differential equations, and illustrate how such models can help explain the function and potential of developmental gene regulatory networks in terms of their dynamical repertoires, that is, in terms of the set of behaviours that can be implemented by the system. A system’s behaviour is defined by its trajectories, which represent the change of the state of the system—e. g. consisting of a set of transcription factor concentrations—over time. The shape of these trajectories depends on the structure or organisation of the underlying regulatory network (see, for example, [3, 22–24, 26–31]). It is possible to gain a general qualitative understanding of what the system’s trajectories can and cannot achieve without the need to solve the dynamical system analytically [32, 33]. Most of this information comes from the geometrical analysis of phase space, i. e. an analysis of the number, nature and relative arrangement of the steady states of the system.
The phase space of a dynamical system is an abstract space, in which each dimension represents the value of a specific state variable. Here we use a well established doublerepressive feedback loop model, known as the toggle switch ([23], and references therein), to study transient dynamics in a twodimensional, timedependent gene regulatory network (see Figure 1B, panels 1 and 2, as well as ‘Model and methods’ below, for the full formulation). In this case, the state variables represent the concentrations of the transcription factors that constitute the network (denoted by X and Y; Figure 1B, panel 1). The graphical representation of phase space is called the phase portrait. It shows the rate of change of the system at any given state. This is known as the flow of the system. The flow of the toggle switch model is indicated by arrows of a given length and direction in the phase portrait shown in Figure 1B, panel 3. If we follow the flow from all possible initial states, we obtain the totality of possible dynamic trajectories.
It is evident from the inspection of the flow in Figure 1B (panel 3) that trajectories tend to converge to specific points in phase space: the steady states of the system. There are different kinds of steady states, those that are stable, and those that are unstable. The most simple stable steady state is an attractor point [32, 33]. Attractors, as their name implies, draw trajectories towards them. Furthermore, they have the special property that once a trajectory has reached an attractor, it will return to it if the system is slightly perturbed. An example of an unstable steady state is a saddle point (Figure 1B, panel 3). Saddle points attract trajectories from some directions, but repel them in others. Usually, the system will move away from a saddle upon perturbation, towards the nearest attractor. The repelling trajectory follows a structure in phase space called an unstable manifold. This manifold is defined by the the path that links a saddle with an attractor point. Note that unstable manifolds correspond to chreodes at the bottom of valleys in Waddington’s landscape.
In our phase portraits, we plot steady states as points coloured according to their stability: attractors in blue, and saddle points in red (Figure 1B, panel 3). The region of phase space around an attractor, from which all trajectories converge towards it, is called its basin of attraction. Curves known as separatrices set apart the different basins and their attractors (in the case of Figure 1B, panel 3, the separatrix is a straight line indicated in grey). Saddle points are always located on separatrices.
Attractors and saddles, with their associated basins and separatrices can be created or annihilated through the process of bifurcation [32, 34]. Bifurcations represent sudden qualitative changes in the structure of the phase portrait caused by small changes in the values of a given set of control parameters.
To date, the best example illustrating the importance and usefulness of geometric approaches for understanding the dynamics and function of specific, experimentally tractable, developmental regulatory systems comes from an analysis of the gap gene regulatory network involved in pattern formation during early embryogenesis of the vinegar fly Drosophila melanogaster. Manu and colleagues [35] identified features of phase space responsible for patterning and canalisation of spatiotemporal gene expression in the Drosophila blastoderm embryo. The analysis is based on lowdimensional projections of phase space to study the geometric arrangement of attractors in the fourdimensional system representing the change in protein concentration for four transcription factors encoded by the gap genes hunchback (hb), Krüppel (Kr), knirps (kni), and giant (gt). Gene expression domain boundaries that remain at a constant position over time could be attributed to movements of attractors in phase space, or nuclei switching between attractor basins depending on their position in the embryo. In contrast, more posterior expression boundaries, which keep shifting position over time, were associated with transient behaviours along a canalising unstable manifold, forming the equivalent of a valley in Waddington’s landscape [35].
Potential landscapes
Phase portraits of systems with two state variables can be visualised in terms of their associated potential landscape [32] (Figure 1B, panel 4). In this representation, the steepness of the potential landscape corresponds to the flow of the system (compare Figure 1B, panels 3 and 4). Trajectories travel downhill towards the attractor points. The topography of the potential landscape is therefore a direct result of the topology (number and nature of steady states) and the geometry (relative positions of the steady states and size of the flow) of the underlying phase portrait.
Potential surfaces can be thought of as mathematical representations of Waddington’s epigenetic landscape [18, 20, 24, 25, 36–40]. Attractors represent differentiated states (in agreement with earlier postulates, [1, 13, 14, 41–46]), separatrices form the ridges between the valleys, which are formed by unstable manifolds representing their associated chreodes (see Figure 1B, panels 3 and 4), and the flow is represented by the steepness of the slopes on the landscape. Canalisation is explained by the depth and width of each valley in the potential landscape. This provides a way of probing and understanding the features that confer robustness to the system (see, for example, [18, 19, 21, 24, 45, 47, 48]). Branching valleys can arise through particular types of bifurcation events. In particular, new attractor states can branch off from existing ones during development or stem cell differentiation through (supracritical) pitchfork bifurcations [19, 20, 25, 32, 47, 49–51].
Due to this analogy to Waddington’s epigenetic landscape, potential landscapes are becoming increasingly popular as explanatory tools in fields such as evodevo, developmental biology, and especially in stem cell research. In the case of stem cells, the positioning of the valleys in the landscape relative to each other explains which differentiation pathways can be reached by a given stem cell, and which cell fates can (or cannot) be transdifferentiated into each other [13–16, 36, 37, 40, 52–54]. This illustrates how specific biological meaning can be gained from studying the topography of potential landscapes, which ultimately, is nothing more than a visually accessible way of studying the underlying phase portrait.
Waddington meant his epigenetic landscape as a “diagrammatic representation” of development and warned explicitly against interpretations that were too rigorous or literal [2, 3]. Detailed topographical interpretations of Waddington’s landscape may be quite inaccurate and even misleading. Ferrell [25] points out that Waddington’s landscape, where valleys progressively split into an increasing number of branches, does not help explain many realistic cell differentiation processes. In particular, many inductive processes in development (e. g. vulval induction in the roundworm Caenorhabditis elegans, or mesoderm induction in vertebrates) involve saddlenode rather than pitchfork bifurcations, which correspond to the disappearance of valleys rather than to their creation by branching. The system shifts to a new attractor only once the old one has vanished. This is why it is important to move from metaphorical uses of Waddington’s epigenetic landscape to accurately calculated potential surfaces whenever possible, or more importantly, to the detailed analysis of the underlying phase portraits, which is where the dynamics of specific regulatory networks are determined.
Potential landscapes can only be calculated and visualised explicitly for a restricted range of dynamical systems, belonging to the class of gradient systems [32]. Note that the notion of a gradient system is defined mathematically by the absence of limit cycles or any other complex attractor structures in their phase portraits, and has nothing to do with biochemical or other biological gradients (see Model and methods for a detailed explanation). In cases where we do not know whether the system under study is a gradient system or not, we can approximate the actual potential using various numerical methods [19, 20, 38, 39, 48, 55]. The resulting approximations are called quasipotential landscapes.
However, even such quasipotential landscapes can only be visualised directly when the number of state variables of the system does not exceed two. This is not true for most biologically realistic systems. Nevertheless, (quasi)potential landscapes are still useful as conceptual tools for the analysis of higherdimensional regulatory networks. In some cases, it is possible to reduce a highdimensional system to a lowerdimensional one (see, for example, [56–58]). But even if this is not the case, the concept of the potential landscape provides two advantages. First, as discussed above, quasipotential surfaces link Waddington’s intuitively accessible concept of the epigenetic landscape to the biological interpretation of (highdimensional) phase space analysis. And second, potential surfaces are useful as visual guides to diagnose features of the underlying phase portrait of the system that are characteristic of specific dynamic behaviours of the regulatory network under study.
The importance of transient dynamics
Many models of biological systems are formulated with the assumption that the relevant dynamics occur near or at a steady state. For instance, Thom’s pioneering systematic and rigorous analysis of morphogenesis in terms of catastrophe theory explicitly and strongly relies on this assumption [8]. Similarly, work on robustness and evolvability, using ensembles of simulated networks, standardly assumes that the steady state pattern produced by a model can be taken as a satisfactory and realistic representation of the phenotype of the system (e. g. [59–62]). Furthermore, stem cell models that make use of potential landscapes are analysed with a strong focus on how their attractors govern system behaviour [19, 20, 47, 48, 56, 58].
In some cases, the steady state assumption makes obvious sense based on biological reasoning. One example is the analysis of the segment polarity network in Drosophila melanogaster, which amplifies and maintains a periodic input, and thus performs an intrinsically stabilising patterning function [63–65]. In most cases, however, biological pattern formation is highly dynamic and far from equilibrium, and the steady state assumption is justified based on methodological, rather than biological considerations. Focussing on asymptotic behaviour at or near steady state greatly simplifies the analysis of the system. First, it discretises and reduces phase (and hence phenotype) space into a small number of possible states—represented by the system’s attractors (see Figure 1C, panel 3). Second, it enables the powerful toolkit of linear stability analysis to be employed to examine the characteristic properties of system states [32, 33].
However, there are both theoretical and practical reasons indicating that steady state analysis misses essential and biologically relevant systems behaviours. One line of theoretical reasoning is provided by Waddington himself, who reminds us that “[i]n the study of development we are interested not only in the final state to which the system arrives, but also in the course by which it gets there”. [3]. His concept of homeorhesis and his representation of developmental trajectories or chreodes as descending valleys in the epigenetic landscape places the focus explicitly on the transient dynamics of cells on their way to their final, differentiated state. In the same spirit, other authors have suggested that phenotypes should be defined over developmental trajectories, rather than representing some sort of ‘final’ outcome or ‘end state’ of the system [22, 66]. To decide what a final pattern is, and when exactly it occurs, is always arbitrary to some degree, while transient features (such as intermediate stages and the timing of their transitions) are clearly important when considering the function and dynamics of developmental processes.
There are also practical reasons to consider transient dynamics explicitly. For developmental processes that consist of continuous transitions between patterns rather than the production of a final output, it is impossible to decide a priori whether the system is representing a nonautonomous succession of steady states (see below), or whether its behaviour is truly transient (i. e. far from steady state). In the case of gap domain shifts, we have evidence for the latter [35], although there is no reason to assume that the two situations need to be mutually exclusive. The gap gene model analysed by Manu and colleagues [35] exhibits boundary shifts that are caused by trajectories following a canalising unstable manifold. Assuming steady state dynamics would collapse the trajectories representing these shifts into a single attractor point at the final configuration of gene expression. Any such analysis would miss the relevant underlying features of phase space (the transient manifold), and therefore fail to provide a proper characterisation and explanation for the observed gene expression dynamics (the shift in domain position over time).
Other examples of developmental processes, where transient dynamics are clearly important are dorsoventral patterning of the vertebrate neural tube, which involves boundary shifts strikingly similar to those of the gap domains [67, 68], and vertebrate somitogenesis or shortgermband segmentation in arthropods where transient travelling waves of gene expression are an essential component of the underlying clockandwavefront patterning mechanism [69–75]. Incidentally, similar considerations can be made for models dealing with ecological networks, and several examples exist in the literature that consider transient dynamics explicitly (see models of coupled oscillating population dynamics between species [76, 77], or [78–80]).
Nonautonomy: explicit time dependence
Another important aspect of biological regulatory processes, which receives surprisingly little attention, is the explicit time dependence of these systems. As soon as we consider cellular dynamics, development, or evolution over a largeenough time span, the organisation of the underlying regulatory system starts to change. This affects the parameters—not just the state variables—of the system. Such explicitly timedependent dynamical systems are called nonautonomous [32, 81, 82]. Timedependent signalling cues and environmental conditions have long been known to shape many processes in the fields of evolutionary and developmental biology. Obvious examples of such phenomena are inductive processes or external (e. g. seasonal) cues that are essential to trigger many developmental pathways (as described in standard textbooks such as [83, 84]), or evolutionary dynamics driven by changing environmental conditions (examples, based on the simulation of gene regulatory network models, can be found in [85–88]).
Still, it is rare to find studies based on explicitly nonautonomous models in the literature, and most authors avoid the challenge of dealing with dynamical systems where the parameters representing external cues are timedependent. This is the case in the study of gap domain shifts by Manu and colleagues [35, 89], where maternal morphogen gradients providing regulatory input to the system were assumed to reach steady state before gap gene boundary positioning was analysed. Such simplifications can be risky, especially when describing biological phenomena where the time scales of change in parameters and state variables are of similar order. In such cases, time scales should not be separated, nor quasisteady states considered since it is easy for dynamical properties and behaviours of the system to be missed or misinterpreted under these conditions.
Recently there have been some attempts at including nonautonomy in biological models. Corson and Siggia [90], for example, offer an explanation for vulva development in C. elegans which explicitly considers temporal parameter changes due to inductive signalling cues in their model. Their model considers differentiation of cells into three differentiated states depending on inputs from two signalling pathways. Signalling inputs are encoded in the model by altering values of system parameters, which change and distort the geometry of the phase portrait by displacing separatrices from their original position (see Figures three–nine in [90]). When a cell receives a signal, its developmental trajectory comes to lie within another basin of attraction, inducing an alternative cell fate to that which would have been reached in the absence of the signal. This study illustrates the importance of nonautonomous dynamics in development. However, it remains somewhat limited in its implementation of explicit timedependence. Although parameter change is included in the model, signalling occurs before cells embark on their developmental paths, and trajectories develop purely autonomously thereafter. Therefore, this approach does not fully capture the transient dynamics of cell differentiation. In other words, although their model is able to offer an explanation for cell differentiation in vulval development, it cannot capture the full nonautonomous nature of the developmental process, since it does not reproduce developmental trajectories, chreodes, in an accurate and fully timedependent manner.
Similarly, most of the few other examples of nonautonomous models in the biological literature do not explicitly consider the effect of parameter changes on transient behaviour (e. g. [91–93]). This simplification may be justified in many cases and is necessary for any kind of rigorous analytical treatment of a model. In many situations, however, it fails to capture essential features of the system. For instance, a truly accurate analysis of gap gene regulatory dynamics would require the inclusion of both nonautonomy from the rapidly changing maternal morphogen gradients, and transient dynamics, which are known to underlie the temporal shifts in domain position. Before we can undertake such an analysis we must first build a conceptual toolkit for phase space analysis of transient, nonautonomous dynamics. Due to the limited amount of analysis possible in such systems, this toolkit will need to be developed in a pragmatic and empirical manner, using numerical simulation and exploration of a simple conceptual model as its basis.
In the following sections, we present such a simulationbased attempt at developing concepts to classify transient, nonautonomous behaviours. For this purpose, we use a simple twocomponent model of a genetic toggle switch, whose potential landscape can be explicitly visualised. We use time series of graphs and animations of systems dynamics on this potential to identify mechanisms leading to state transitions, and other forms of pattern formation. From this, we are able to identify four basic types of dynamical mechanism and behaviour—transitions, pursuits and two types of captures of trajectories—that can be used to classify and understand dynamical behaviour in more complicated and realistic models, such as a full nonautonomous version of the gap gene model. While this present paper provides the methodological foundation for such an analysis, detailed results for a realistic model of the gap genes will be presented elsewhere.
Model and methods
To develop our methodology for analysing transient behaviour in nonautonomous dynamical systems, we use a simple toggle switch model (see [22], and references therein) with timedependent parameters. We consider two interacting genes X and Y (Figure 1B, panel 1). Concentrations of the corresponding protein products are labelled x and y. X and Y mutually repress each other, are linearly activated by external signals and can autoactivate themselves (Figure 1B, panel 1). Protein products decay linearly dependent on their concentration. The mathematical formulation of our toggle switch model is thus given by
where parameters α _{ x } and α _{ y } represent the external activation on genes X and Y respectively. Sigmoid functions with Hill coefficients of 4 are used to represent autoactivation and mutual repression, where parameters a and c determine autoactivation thresholds, while b and d determine thresholds for mutual repression. Protein decay rates are represented by parameters λ _{ x } and λ _{ y }.
The toggle switch model (1) exhibits different dynamical regimes depending on the values of its parameters (Figure 2A–C). Its name derives from the fact that it can exhibit bistability over a wide range of parameters. When in this bistable region of parameter space, the underlying phase portrait has two attracting states and one saddle point (Figure 2B). All phase portraits associated with parameters in the bistable range are topologically equivalent to each other, meaning that they can be mapped onto each other by a continuous deformation of phase space called a homeomorphism [34].
The toggle switch model has two other dynamical regimes: monostable and tristable. Phase portraits associated with parameters in the monostable range have only one attractor point (Figure 2A), while those in the tristable range have three attractor states and two saddle points (Figure 2C). Again, phase portraits within each regime are topologically equivalent to each other. While phase space can be geometrically deformed within each regime (through movements of attractors or separatrices), its topology only changes when one regime transitions into another through different types of bifurcations [32, 34] (Figure 2). The transition from monostable to bistable is known to be governed by a supracritical pitchfork, the transition from bistable to tristable involves a subcritical pitchfork bifurcation, and the transition from tristable to monostable takes place through two simultaneous saddlenode bifurcations involving the two attractors labelled in darker blue in Figure 2C.
Definition of the potential landscape
Potential landscapes can only be calculated explicitly for the class of dynamical systems called gradient systems [32]. A twovariable gradient system is a dynamical system
which satisfies the following relationship between partial derivatives
For gradient systems, it is possible to calculate a closedform (explicit) potential function, V(x,y) such that
The local minima on the twodimensional potential surface given by V(x,y) correspond mathematically to the steady states of the system in (2) since, if (x ^{∗}, y ^{∗}) is such that
then
and, therefore, (x ^{∗}, y ^{∗}) is a steady state of (2).
Calculating quasipotential landscapes
Condition (3) will not always be met. In particular, dynamical systems representing gene interaction networks are not in general gradient systems, and therefore an associated potential function and landscape may not exist. In such cases, we can still take advantage of the visualisation power of potential landscapes by approximating the true potential using a numerical method. The numerical approximation method we adopt for our study was developed by Bhattacharya and colleagues [38] using a toggle switch model very similar to the one used here. This allows us to calculate a quasipotential landscape for any specific set of fixed parameter values.
The quasipotential, which we denote by V _{ q }, is defined to decrease along all trajectories of the dynamical system as they progress on the phase portrait over time
Δ x and Δ y are defined as smallenough increments along the trajectory such that d x/d t and d y/d t can be considered constant in the interval [(x,x+Δ x),(y,y+Δ y)]. In addition, $\mathrm{\Delta x}=\frac{\mathit{\text{dx}}}{\mathit{\text{dt}}}\mathrm{\Delta t}$ and $\mathrm{\Delta y}=\frac{\mathit{\text{dy}}}{\mathit{\text{dt}}}\mathrm{\Delta t}$, where Δ t is the time increment. Substituting into equation 7, we obtain
Δ V _{ q } has been formulated in such a way that, for positive time increments Δ t, Δ V _{ q } is always negative along the unfolding trajectory and is, in effect, a Lyapunov function of the twogene dynamical system [32]. This ensures that trajectories will always “roll” downhill on the quasipotential surface. Just as in the case of closedform potential, the steady states of the system (x ^{∗},y ^{∗}) correspond to the local minima on the quasipotential surface since Δ V _{ q }(x ^{∗},y ^{∗})=0.
We apply the numerical approximation method described above to trajectories with various initial points on the xy plane. This yields a sampled collection of trajectories with quasipotential values associated to every one of their points. Next, we apply the following two assumptions, in order to construct a continuous quasipotential surface from this sample of discrete trajectories [38]:

1.
Two trajectories with different initial conditions that converge to the same steady state must also converge to the same final quasipotential level (normalisation within basins of attraction).

2.
Two adjacent trajectories that converge to different steady states will be taken to start from the same initial quasipotential level (normalisation between basins of attraction).
Finally, interpolation of all the normalised trajectories results in a continuous quasipotential landscape. Bhattacharya et al.[38] validated this approach by demonstrating that the quasipotential values of the steady states were inversely correlated with their probability of occurrence using a stochastic version of the toggle switch dynamical system.
Approximating nonautonomous trajectories
As we have argued in the Background Section, we cannot generally assume that parameter values remain constant over time when modelling biological processes. We take a stepwise approximation approach to the change in parameter values to address this problem (Figure 3). We chose a time increment (step size) as small as possible. Parameter values are kept constant for the duration of each time step. As a consequence, the associated phase portrait will also remain constant during this time interval, and is visualised for each step by calculating a quasipotential landscape as described in the previous section (Figure 3C, top row).
The smaller the time increments considered, the better we are able to approximate continuous changes in parameter values, and the consequent changes to the associated phase portrait and quasipotential landscape. Such accurate approximation allows us to faithfully reproduce nonautonomous trajectories produced by models with continuously timevariable parameters. This is done by integrating trajectories using constant parameters during each time step, and then using the resulting end position in phase space as the initial condition for the next time step. The resulting integrated trajectories can then be visualised by mapping them from the underlying phase plane onto the associated quasipotential landscape as described above. This allows us to track and analyse in detail how changes in the phase portrait and quasipotential landscape shape the trajectories as the values of the parameters are changing.
Results and discussion
Using the toggle switch defined in (1) as a conceptual model, we have set out to identify and catalogue mechanisms that affect the shape of trajectories in nonautonomous systems. In contrast to previous studies (see, for example, [19, 20]), we do not assume asymptotic behaviour of the system, but focus on transient dynamics—that is, the entire trajectory from initial to steady state.
In this context, we use the term ‘mechanism’ in an unusual, but precisely defined way. We do not mean it to imply any specific biochemical interactions—this is how most experimental biologists would use the term. Nor does our notion of ‘mechanism’ depend on any specific regulatory structure, such as network motifs (reviewed in [94]), or the different stripeforming mechanisms studied in [95]. Instead, we use the term ‘mechanism’ in a broader, more abstract sense: a mechanism is a causal explanation of dynamics in terms of the dependence of the flow on parameter changes, and how these affect the trajectories of the system. We will illustrate this very abstract definition using a number of concrete examples below.
Changing parameter values in nonautonomous systems can change the underlying phase portrait in two main ways: it can alter the geometry or the topology of phase space (or both) [32, 34]. On one hand, a change in geometry means that the phase portrait is still composed of the same basic elements—the same number, kind and relative placement of attractors and their associated basins, saddle points, and separatrices—but that the exact position, shape, and/or size of these elements has been altered. On the other hand, a topological alteration of phase space occurs when a change in parameter values cause the number or stability of steady states to change. For example, small changes in a parameter value can cause a pair of attractor and saddle points to appear or disappear through a saddlenode or fold bifurcation, or an attractor can be turned into a saddle point through a transcritical bifurcation. Topological changes to the phase portrait are therefore associated with bifurcations as studied in the context of autonomous dynamical systems [32–34].
In what follows, we examine how both geometrical and topological changes in phase space affect transient trajectories in nonautonomous systems. For this purpose, we use a numerical simulation approach with threedimensional visualisation of the quasipotential landscape associated with the phase portrait of our toggle switch model (see the Model and methods Section, for details). We identify four different basic types of nonautonomous mechanisms affecting dynamical behaviour: transitions, pursuits and two kinds of captures.
Transitions
Out of the four trajectoryshaping mechanisms that we have identified in this study, transitions are the most familiar. They have been found and described in previous approaches to the study of nonautonomous biological systems [78–80, 91, 92]. Transitions in our context are equivalent to Thom’s ‘catastrophes’ [8], and to what Scheffer has termed ‘critical transitions’, defined as a shift of the system from one attractor to another when it passes a given critical point in parameter space [92]. While Thom’s and Scheffer’s analyses focus on the resulting steady state behaviour of the system, we also consider the transient dynamics during a transition.
An example of a transition in the toggle switch model—occurring from the bistable via the tristable to the monostable regime—is shown in Figure 4 (see also Additional file 1, Supporting Movie S1). This transition is driven by an increase in the value of the autoactivation thresholds (a and c in equation 1). Changes in autoactivation have been previously used to explain the dynamics of stem cell differentiation [20, 21, 47, 48], and patterning by lateral inhibition [25].
In this case, we consider that the system is already at steady state at the outset, in its bistable regime with high x and low y (Figure 4A). In other words, we assume that the system has had time to converge, or that the initial condition coincides with the neighbourhood of this attractor. During the transition, two bifurcation events take place—a subcritical pitchfork bifurcation (Figure 4A to B), and two simultaneous saddlenode bifurcations (Figure 4D to E)—which make the initial steady state and its associated basin of attraction disappear while creating a new attractor at which both factors X and Y are present at low concentrations (Figure 4B to E).
Examining transient dynamics during the transition reveals important details that shape the trajectory and hence the behaviour of the system. First, we note that the trajectory initially at steady state does not necessarily remain anchored to its attractor (Figure 4C and D). As the change in parameter values causes the attractor state to move, the trajectory’s current state falls behind and reacts by travelling towards the moving attractor (see Pursuit below). Since the flow rate along the trajectory is smaller than the velocity of attractor movement, the system is not able to catch up with the moving steady state. Hence, it temporally reverts from asymptotic to transient behaviour. Obviously, such a reversal can never be observed if we only focus on steady state behaviour.
Second, we observe a delay between the subcritical pitchfork bifurcation creating the new attractor state (shown in light blue in Figure 4B) and the system switching into that basin of attraction (Figure 4E). This delay effect was already noted and described by Thom in his analysis [8]. The change from one basin to another (indicated by a change from dark to light blue on the trajectory shown in Figure 4E) only occurs once the tristable system undergoes two simultaneous saddlenode bifurcations, where the two saddles on the separatrices collide against the two outer attractors and annihilate each other (darker blue attractors in Figure 4D and E). A monostable system results from these bifurcations. The trajectory suddenly finds itself in a different basin, and eventually it will converge to this new attractor state (Figure 4E).
In this particular example, the shape of the observed transient trajectory does not differ much from its equivalent in an autonomous monostable system. This is because for most of the time the trajectory is in the basin of attraction of the attractor at low X and Y concentrations (light blue part of the trajectory in Figure 4E). Nevertheless, most developmental processes require tightly controlled timing and therefore the observed temporary reversal to transient behaviour and consequent delay in attractor switch can be significant. These aspects are therefore not negligible if we want to achieve a full understanding of the dynamics of the system. Steady state analyses, where more or less instant convergence to the new attractor state is assumed, will be limited in this regard since they are not able to address the question of developmental timing.
In addition to its effect on timing, the delay in switching basins of attraction also provides an explanation for the irreversibility of developmental pathways. This phenomenon has been studied in detail in the context of stem cell differentiation. Wang and colleagues [19, 20] use a nonautonomous toggle switch model similar to ours to study the transition from monostable (stem cell) to bistable behaviour (differentiated cells). The authors simulate this transition in both directions (from stem cell to differentiated state and back), and show that two different delay effects occur depending on the direction of the process. In other words, forward and backward pathways are very different, thus explaining why the system is unable to retrace its original differentiation pathway when the process is reversed.
Pursuit
Transient phenomena need to be taken into account to be able to fully understand the dynamical repertoire of a system. This becomes especially relevant when considering what we call pursuit mechanisms. We now consider what happens if the external activation on genes X and Y in our toggle switch model (α _{ x } and α _{ y } in 1) are altered over time. In biological terms, this could represent the system’s response to a changing external signal. This particular parameter change affects the geometry but not the topology of the phase portrait. The system remains in the bistable dynamical regime throughout the whole parameter range that we explored.
Let us first consider an increase in activation strength. As external production rates increase over time, the attractors move away from the origin along the direction of the x and y axes respectively (Figure 5, see also Additional file 2, Supporting Movie S2). This outward movement of the attractors leaves the position of the separatrix unchanged throughout the simulation, keeping the location and area of the basins of attraction constant over time (Figure 5). Since the separatrix does not move, trajectories are not able to change basins. Therefore, the initial conditions determine which basin of attraction the system will remain in, as they do in autonomous systems. However, the approach to the attractor is very different in this nonautonomous case, since trajectories are drawn towards a moving target. It is this continued pursuit that plays the main role in shaping them (Figure 5). Even though its associated attractor keeps moving over time, this sort of pursuit shows how the general direction of a trajectory can be stably maintained in a nonautonomous system.
Two specific pursuit scenarios can be distinguished. In the first, the trajectory eventually reaches the moving attractor. This occurs if the change in parameter values increases the flux around the attractor faster than it moves its position, or if further change in parameter values no longer alters the position of an attractor—in other words, the position of an attractor itself converges to a given location over time. Under these conditions, the trajectories of a nonautonomous system can show asymptotic behaviour and come to rest at steady state. However, it is not guaranteed that this should take place. There are many imaginable scenarios in which attractors will keep moving as parameters change and trajectories governed by pursuit mechanisms may never come to rest. It is probable that many developmental processes—at least to a certain degree—work in this regime, as external conditions and the cellular environment (signalling inputs or tissue context) constantly keep changing over time.
Just as in the case of transitions, it is obvious that pursuit mechanisms have a great impact on the timing of developmental dynamics. They can delay or even prevent the system from reaching steady state. In addition, attractor movements can also drastically influence and alter the shape of a trajectory. This becomes evident if we consider an alternative scenario, in which activating inputs decrease over time (Figure 6, see also Additional file 3, Supporting Movie S3). In this case, the attractors move towards the origin over time. This causes the attractor to ‘overtake’ the trajectory at a given point in time (Figure 6B), which induces a rather drastic change in the trajectory’s direction: initially it is travelling toward high x (Figure 6A), but later comes to approach the origin (low x) while pursuing the moving attractor (Figure 6B–D).
In summary, we have illustrated how pursuit mechanisms in nonautonomous systems can either lead to the stable maintenance of a transient trajectory (as in Figure 5), or to a drastic change in its direction (as in Figure 6) depending on the geometric arrangement of the attractor with regard to the converging trajectory. We have argued that pursuit mechanisms are likely to be almost ubiquitous in systems that are exposed to a changing environment or variable tissue context (e. g. signalling or external input by morphogen gradients). One specific example for pattern formation by pursuit are the transient gap domain shifts examined by Manu and colleagues [35, 89]. In this particular case, spatial shifts in domain position are caused by the timed up and later downregulation of a specific sequence of gap genes in nuclei located within the posterior region of the Drosophila blastoderm embryo (see [35], for details). Basically, each nucleus goes through a stereotypical series of expression pulses of different gap genes. Such pulses can easily be explained by the sort of trajectory observed in Figure 6, where an initial expression trend (upregulation of X) is later smoothly reversed (towards downregulation of the same gene).
Capture
In autonomous dynamical systems, trajectories never cross a separatrix [32]. Depending on the position of its initial conditions, a trajectory will find itself in a particular basin of attraction, that is, in the dynamical regime of a particular attractor. The trajectory will always remain in this basin and, given enough time, will eventually converge to the attractor.
The same is not true for nonautonomous systems. We know that the phase portrait is changing as the parameters change. This means that the position and number of steady states need not remain constant over time. We have seen how saddles and attractors can be created and annihilated by bifurcations, and how they can move through phase space in such systems. In particular, when saddle points change their position, they cause their associated separatrices to shift with them. Separatrices mark the boundaries between basins of attraction, and their movement will cause the size, shape and/or relative placement of these basins with regard to each other to change.
As trajectories progress on the changing phase portrait, it is possible for a moving separatrix to overtake them. We will call this event a ‘capture’. When a capture occurs, the affected trajectory is deviated, and often exhibits a sudden change in direction, since it is now attracted towards a different attractor at a different position of the phase plane than it was before. Captures will happen when the flow along the trajectory is smaller than the rate of change in separatrix position. The separatrix is recruiting points from one basin of attraction into another at a faster rate than the trajectory is travelling away from it. Therefore, the velocities of the trajectory and the separatrix relative to each other are determining whether a capture event is going to take place or not.
We identified two main ways by which capture events are brought about in the nonautonomous toggle switch model. In both cases, the mechanism shaping the trajectory is essentially the same: the trajectory suddenly finds itself converging towards a different attractor as a result of an abrupt change to the flow at its current position. The main difference between the two capture mechanisms resides in how the movement of the separatrix is caused.
In the first situation, the trajectory gets captured after a bifurcation event has lead to the creation of a new attractor state. This increases the number of attractor basins and, in this way, introduces new separatrices into the phase portrait. We simulate this situation using the bistabletotristable transition caused by an increase in the autoactivation threshold as described above (Figure 7, see also Additional file 4, Supporting Movie S4). In this example, a subcritical pitchfork bifurcation creates a new attractor and two associated saddles from a preexisting saddle point (Figure 7A,B). This results in a change in phase space topology. What used to be a single separatrix now ‘opens up’, giving rise to two different forked separatrices (Figure 7B,C). Further parameter changes then cause the new separatrices to move outward through phase space, catching up with, and overtaking, trajectories as they recruit points into the newly created and expanding basin of attraction (Figure 7B–D).
Additional file 4: Supporting Movie 4–Capture due to a change in the topology of the phase portrait. A capture results from a trajectory being recruited into a new basin of attraction due to the movement of a separatrix. In this example, the relevant separatrix is created and caused to move by a preceding bifurcation event, which leads to the appearance of a new attractor state, resulting in a change of phase space topology. This movie shows the simulation that Figure 7 is based on. It displays the changing potential landscape as the system transitions from the bistable to the tristable and then the monostable regime (see Figure 2) as autoactivation thresholds are increased. The corresponding phase portrait is shown in the inset on the upper right. Parameter values are displayed in the panel on the lower right. The current state of the system is shown as a black ball. Attractors are marked by blue dots, saddle nodes by red dots with associated separatrices shown as black lines. The history of the trajectory is shown on the phase portrait, colour coded according to which attractor the trajectory is converging towards at every individual step of the simulation. See main text and Figure 7 for details. (MP4 11 MB)
As with the transition mechanism, we notice a sometimes considerable time delay between the bifurcation and the capture event. The extent of this delay depends on how close a trajectory is to the new attractor at the time of bifurcation, and on the rate of change in separatrix position due to further parameter change. In our example, the delay effect can be clearly seen, as the trajectory first veers to the left, towards its original attractor, before being steered back towards the center of the phase portrait, where the new attractor lies. This introduces a marked bend into the trajectory (Figure 7D), which affects both time to convergence and the identity of transient states.
In the second scenario, we consider a situation where the capture event is not preceded by a bifurcation. In this case, a preexisting separatrix is moving through phase space due to parameter changes. This movement of the separatrix reconfigures the geometry of the basins of attraction without changing the topology of the phase portrait. We simulate this event in our toggle switch model by introducing asymmetric changes in the values of the thresholds determining the two mutually repressing regulatory interactions (i.e. b≠d in 1) (Figure 8, see also Additional file 5, Supporting Movie S5). This shifts the position of the separatrix in the bistable regime towards one of the two attractors without creating or annihilating any steady states.
Additional file 5: Supporting Movie 5—Capture due to a change in the geometry of the phase portrait. A capture results from a trajectory being recruited into a new basin of attraction due to the movement of a separatrix. In this example, the relevant separatrix is caused to move by a break in the symmetry of the repressive interactions between X and Y. No new attractor states or separatrices are created. Therefore, the change in landscape topography is purely geometrical and does not affect the topology of phase space. This movie shows the simulation that Figure 8 is based on. It displays the changing potential landscape as the threshold for repression of X by Y is raised, causing the system to be increasingly asymmetrical. The corresponding phase portrait is shown in the inset on the upper right. Parameter values are displayed in the panel on the lower right. The current state of the system is shown as a black ball. Attractors are marked by blue dots, the saddle node by a red dot with its associated separatrix shown as a black line. The history of the trajectory is shown on the phase portrait, colour coded according to which attractor the trajectory is converging towards at every individual step of the simulation. See main text and Figure 8 for details. (MP4 11 MB)
As in the case of capture with preceding bifurcation, the separatrix moves and recruits points from one basin into the other, thereby shrinking one basin of attraction while enlarging the other (Figure 8A–D). If the moving separatrix encounters a trajectory at one of the points that are being recruited into a different basin of attraction, a capture event will take place (Figure 8B). The affected trajectory will change direction as it now converges towards a different attractor than before the capture (Figure 8C,D).
Capture mechanisms illustrate how dramatic changes in direction and velocity of trajectories do not always require changes to the topology of the underlying phase portraits. In other words, they do not necessarily have to be associated with bifurcations, and even if they are, the observed sudden change in direction does not usually coincide with the bifurcation. Instead, it can be significantly delayed. This kind of phenomenon cannot be observed if only asymptotic behaviours are considered. In that case, only transitionlike events can occur since captures and pursuits crucially depend on transient dynamics.
There are several examples which illustrate how this might be important in practice. Corson and Siggia [90] considered an explicitly nonautonomous mechanism for determining vulval cell fates in their model. In their study, shifts in separatrices and captures happen instantaneously at the outset of the simulation. In other words, the initial condition of cells that are captured completely determines which basin of attraction they will find themselves in for the rest of the simulation. While this does not affect the analysis of the resulting differentiated (steady) states—as evidenced by the tight match between the model and the measured proportions of differentiated cells [90]—it could limit the model’s ability to reproduce the transient dynamics of cellular differentiation.
This becomes crucial when analysing more dynamic examples of gene expression, such as pattern formation driven by the gap gene system. In this example, it will be essential to study the interplay between changing maternal gradients and moving gap target domains. This was not really possible in the study of Manu and colleagues [35] since their model assumed the maternal gradients to remain constant. Our classification scheme of transient behaviour in nonautonomous systems should be useful to analyse such questions in more detail. For instance, capture events occurring at different times in different nuclei could explain how stable expression boundaries can be maintained in the presence of changing regulatory inputs. In contrast, moving domain boundaries in the posterior of the embryo are more likely to involve some type of pursuit mechanism (see above). Without doubt, the concepts developed in this paper will be useful to classify and characterise the boundaryforming and shifting mechanisms in a fully nonautonomous model of this system.
Conclusions
In this paper, we have argued that considering both transient dynamics and explicit timedependence (nonautonomy) of dynamical systems is essential for an accurate and complete quantitative understanding of many biological processes. In particular, we have argued that these two features are important if we are to fully capture the richness and explanatory power of Waddington’s metaphor of the epigenetic landscape [3]. Waddington’s concepts of developmental chreodes, and their canalised nature due to the regulatory process of homeorhesis, apply to transient or nonautonomous behaviour, not stable steady states of a system. Moreover, Waddington’s landscape changes its shape due to genetic or environmental signals as represented by the pegs and ropes in Figure 1A.
However, it remains difficult to analyse transient dynamics in nonautonomous systems. Where rigorous analysis is possible, it is mostly limited to special cases (e.g., systems with regular external forcing, [82]). No general theory or toolkit exists to deal with more complex nonautonomous systems, or the transient phenomena they exhibit.
In this paper, we have taken a first step towards overcoming this limitation. We have used a pragmatic simulationbased approach to classify different mechanisms that affect transient trajectories in nonautonomous systems. Our study is based on a simple conceptual model—a toggle switch with timedependent parameter values—to visualise the changing potential landscape which determines the dynamical repertoire of a regulatory circuit. We found four distinct mechanisms that affect trajectories on timevariable potential surfaces: transitions, pursuits and two kinds of captures. At this point it is impossible to say whether this list is complete. Other types of behaviour may be possible in higherdimensional phase spaces.
Despite its potential incompleteness, we believe that our preliminary classification scheme is useful, since it makes transient dynamics in nonautonomous systems more tractable in practice. The advantages of our framework are twofold: first, it allows us to describe and categorise patterning mechanisms. Each one of the mechanisms described above has its own characteristics, both in terms of their effect on the shortterm dynamics, as well as the longerterm evolutionary potential of the system. We have described specific examples of such characteristics in the Results and discussion Section. Second, it can be used as a diagnostic framework for the analysis of phase space, even in the case of higherdimensional systems where the potential surface can no longer be visualised explicitly. It shows us, for example, that we need not necessarily find a bifurcation associated with every observed critical transition of a system. Instead, drastic changes in systems behaviour can also be caused by pursuit or geometrical capture events as described in the Results and discussion Section.
Our framework is widely applicable to the analysis of realworld regulatory networks. In fact, it has been developed in the context of our analysis of a fully nonautonomous version of the gap gene model by Manu et al.[35]. Our classification scheme has allowed us to categorise different mechanisms for spatial pattern formation in this system, and to identify those features of phase space responsible for the spatiotemporal dynamics of gene expression. A full analysis of this system will be published elsewhere. We hope that our work will inspire further such analyses in systems where our framework could be useful, and where the required detailed and datadriven models are currently available. Examples such as vulval determination in C. elegans or neural patterning in vertebrates come to mind. We have no doubt that—considering the current rate of progress in systems biology—this list will be extended considerably in the very near future, leading to a more general investigation of transient dynamics in nonautonomous biological systems.
References
 1.
Waddington C: An introduction to modern genetics. 1939, London: Allen & Unwin
 2.
Waddington C: Canalisation of development and the inheritance of acquired characters. Nature. 1942, 150: 563564. 10.1038/150563a0.
 3.
Waddington C: The Strategy of the Genes. 1957, London: Allen & Unwin
 4.
Waddington C: Genes and Organisers. 1940, Cambridge: Cambridge University Press
 5.
Gilbert S: Epigenetic landscaping: Waddington’s use of cell fate bifurcation diagrams. Biol Philos. 1991, 6: 135154. 10.1007/BF02426835.
 6.
Fagan M: Waddington redux: models and explanation in stem cell and systems biology. Biol Philos. 2012, 27: 179213. 10.1007/s105390119294y.
 7.
Gilbert S: Diachronic biology meets evodevo. Am Zool. 2000, 40: 729737. 10.1668/00031569(2000)040[0729:DBMEDC]2.0.CO;2.
 8.
Thom R: Structural Stability and Morphogenesis. 1975, Boulder: CO: Westview Press
 9.
Thom R: An inventory of Wadington’s concepts. Theoretical Biology. Edited by: Goodwin B, Saunders P. 1989, Edinburgh: Edinburgh University Press,
 10.
Saunders P, Kubal C: Bifurcations and the epigenetic landscape. Theoretical Biology. Edited by: Goodwin B, Saunders P. 1989, Edinburgh: Edinburgh University Press,
 11.
Saunders P: The organism as a dynamical system. Thinking about Biology. Edited by: Varela F, Stein W. 1993, Reading, MA: Addison Wesley,
 12.
Slack J: Conrad Hal Waddington: The last renaissance biologist?. Nat Rev Genet. 2002, 3: 889895.
 13.
Enver T, Pera M, Peterson C, Andrews P: Stem cell states, fates, and the rules of attraction. Cell Stem Cell. 2009, 4: 387397. 10.1016/j.stem.2009.04.011.
 14.
Graf T, Enver T: Forcing cells to change lineages. Nature. 2009, 462: 587594. 10.1038/nature08533.
 15.
Hochedlinger K, Plath K: Epigenetic reprogramming and induced pluripotency. Development. 2009, 136: 509523. 10.1242/dev.020867.
 16.
Yamanaka S: Elite and stochastic models for induced pluripotent stem cell generation. Nature. 2009, 460: 4952. 10.1038/nature08180.
 17.
Ladewig J, Koch P, Brüstle O: Leveling Waddington: the emergence of direct programming and the loss of cell fate hierarchies. Nat Rev Mol Cell Biol. 2013, 14: 225236. 10.1038/nrm3543.
 18.
Huang S: Reprogramming cell fates: reconciling rarity with robustness. Bioessays. 2009, 31: 546560. 10.1002/bies.200800189.
 19.
Wang J, Xu L, Wang E, Huang S: The potential landscape of genetic circuits imposes the arrow of time in stem cell differentiation. Biophys J. 2010, 99: 2939. 10.1016/j.bpj.2010.03.058.
 20.
Wang J, Zhang K, Xu L, Wang E: Quantifying the Waddington landscape and biological paths for development and differentiation. Proc Nat Acad Sci. 2011, 108: 82578262. 10.1073/pnas.1017017108.
 21.
Zhou JX, Huang S: Understanding gene circuits at cellfate branch points for rational cell reprogramming. Trends Genet. 2011, 27: 5562.
 22.
Jaeger J, Irons D, Monk N: The inheritance of process: a dynamical systems approach. J Exp Zool (Mol Dev Evol). 2012, 318B: 591612.
 23.
Jaeger J, Crombach A: Life’s Attractors. Evolutionary Systems Biology. Edited by: Soyer O. 2012, New York: Springer,
 24.
Huang S: The molecular and mathematical basis of Waddington’s epigenetic landscape: a framework for postDarwinian biology?. BioEssays. 2012, 34: 149157. 10.1002/bies.201100031.
 25.
Ferrell J: Bistability, bifurcations, and Waddington’s epigenetic Landscape. Curr Biol. 2012, 22: 458466. 10.1016/j.cub.2012.03.045.
 26.
Goodwin B: Development and evolution. J Theor Biol. 1982, 97: 4355. 10.1016/00225193(82)902752.
 27.
Alberch P: Developmental Constraints in Evolutionary Processes. Evolution & Development. Edited by: Bonner JT. 1982, Berlin/Heidelberg: Springer,
 28.
Oster G, Alberch P: Evolution and bifurcation of developmental programs. Evolution. 1982, 36: 444459. 10.2307/2408093.
 29.
Alberch P: From genes to phenotype: dynamical systems and evolvability. Genetica. 1991, 84: 511. 10.1007/BF00123979.
 30.
Webster G: Goodwin B: Form and Transformation: Generative and Relational Principles in Biology. Cambridge. 1996, : UK: Cambridge University Press
 31.
François P, Siggia E: Phenotypic models of evolution and development: geometry as destiny. Curr Opin Genet Dev. 2012, 22: 627633. 10.1016/j.gde.2012.09.001.
 32.
Strogatz S: Nonlinear dynamics and chaos. Boulder. 1994, CO: Westview Press
 33.
Hirsch M, Smale S, Devaney R: Differential Equations, Dynamical Systems, and an Introduction to Chaos. 2013, Amsterdam: Elsevier
 34.
Kuznetsov Y: Elements of applied bifurcation theory. 2004, New York: Springer
 35.
Surkova A, Spirov A, Gursky V, Janssens H, Kim A, Radulascu O, VanarioAlonso C, Sharp D, Samsonova M, Reinitz J, Manu: Canalization of gene expression and domain shifts in the drosophila blastoderm by dynamical attractors. Plos Comput Biol. 2009, 5: e100030310.1371/journal.pcbi.1000303.
 36.
Macarthur B, Ma’ayan A, Lemischka I: Toward stem cell systems biology: from molecules to networks and landscapes. Cold Spring Harb Symp Quant Biol. 2008, 73: 21110.1101/sqb.2008.73.061.
 37.
Macarthur B, Ma’ayan A, Lemischka I: Systems biology of stem cell fate and cellular reprogramming. Nat Rev Mol Cell Biol. 2009, 10: 672681.
 38.
Bhattacharya S, Zhang Q, Andersen M: A deterministic map of Waddington’s epigenetic landscape for cell fate specification. BMC Syst Biol. 2011, 5: 8510.1186/17520509585.
 39.
Zhou J, Aliyu M, Aurell E, Huang S: Quasipotential landscape in complex multistable systems. J R Soc Interface. 2012, 9: 35393553. 10.1098/rsif.2012.0434.
 40.
Furusawa C, Kaneko K: A dynamicalsystems view of stem cell biology. Science. 2012, 338: 215217. 10.1126/science.1224311.
 41.
Delbrück M: Discussion. Unités biologiques douées de continuité génétique. 1949, Paris: Editions du Centre National de la Recherche Scientifique,
 42.
Jacob F, Monod J: Genetic regulatory mechanisms in the synthesis of proteins. J Mol Biol. 1961, 3: 318356. 10.1016/S00222836(61)800727.
 43.
Kauffmann S: Homeostasis and differentiation in random genetic control networks. Nature. 1969, 224: 177178. 10.1038/224177a0.
 44.
Kauffmann S: The Origins of Order: SelfOrganization and Selection in Evolution. 1993, Oxford: Oxford University Press
 45.
Huang S, Eichler G, BarYam Y, Ingber D: Cell fates as highdimensional attractor states of a complex gene regulatory network. Phys Rev Lett. 2005, 94: 128701128705.
 46.
Chang H, Hember M, Barahona M, Ingber D, Huang S: Transcriptomewide noise controls lineage choice in mammalian progenitor cells. Nature. 2008, 453: 544547. 10.1038/nature06965.
 47.
Huang S, Guo Y, May G, Enver T: Bifurcation dynamics in lineagecommitment in bipotent progenitor cells. Dev Biol. 2007, 305: 695713. 10.1016/j.ydbio.2007.02.036.
 48.
Kim K, Wang J: Potential energy landscape and robustness of a gene regulatory network: toggle switch. PLoS Comput Biol. 2007, 3: e6010.1371/journal.pcbi.0030060.
 49.
Roeder I, Horn M, Glauche I, Hochhaus A, Mueller M, Loeffler M: Dynamic modeling of imatinibtreated chronic myeloid leukemia: functional insights and clinical implications. Nat Med. 2006, 12: 11811184. 10.1038/nm1487.
 50.
Guantes R, Poyatos J: Multistable decision switches for flexible control of epigenetic differentiation. PLoS Comput Biol. 2008, 4: e100023510.1371/journal.pcbi.1000235.
 51.
Andrecut M, Halley J, Winkler D, Huang S: A general model for binary cell fate decision gene circuits with degeneracy: indeterminacy and switch behavior in the absence of cooperativity. PLoS One. 2011, 6: e1935810.1371/journal.pone.0019358.
 52.
Goldberg A, Allis C, Bernstein E: Epigenetics: a landscape takes shape. Cell. 2007, 128: 635638. 10.1016/j.cell.2007.02.006.
 53.
Hemberger M, Dean W, Reik W: Epigenetics dynamics of stem cells and cell lineage commitment: digging Waddington’s canal. Nat Rev Mol Cell Biol. 2009, 10: 526537. 10.1038/nrm2727.
 54.
Furusawa C, Kaneko K: Chaotic expression dynamics implies pluripotency: when theory and experiment meet. Biol Direct. 2009, 4: 1110.1186/17456150411.
 55.
Wang J, Xu L, Wang E: Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations. PNAS. 2008, 105: 1227112276. 10.1073/pnas.0800579105.
 56.
Choi M, Shi J, Jung SH, Chen X, Cho KH: Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage. Sci Signal. 2012, 5: ra83
 57.
Qiu X, Ding S, Shi T: From understanding the development landscape of the canonical fateswitch pair to constructing a dynamic landscape for twostep neural differentiation. PloS One. 2012, 7: e49271
 58.
Li C, Wang J: Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS Comput Biol. 2013, 9: e100316510.1371/journal.pcbi.1003165.
 59.
Wagner A: Does evolutionary plasticity evolve?. Evolution. 1996, 50: 10081023. 10.2307/2410642.
 60.
Siegal M, Bergman A: Waddington’s canalization revisited: developmental stability and evolution. Proc Nat Acad Sci. 2002, 99: 1052810532. 10.1073/pnas.102303999.
 61.
Wagner A: Robustness and Evolvability in Living Systems. Princeton. 2005, NJ: Princeton University Press
 62.
Wagner A: The Origins of Evolutionary Innovations: A Theory of Transformative Change in Living Systems. 2011, Oxford: Oxford University Press
 63.
von Dassow G, Meir E, Munro E, Odell G: The segment polarity network is a robust developmental module. Nature. 2000, 406: 188192. 10.1038/35018085.
 64.
Albert R, Othmer HG: The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003, 223: 118. 10.1016/S00225193(03)000353.
 65.
Ingolia N: Topology and robustness in the Drosophila segment polarity network. PLoS Biol. 2004, 2: e12310.1371/journal.pbio.0020123.
 66.
Fusco G: How many processes are responsible for phenotypic evolution?. Evol Dev. 2001, 3: 279286. 10.1046/j.1525142x.2001.003004279.x.
 67.
Dessaud E, Yang L, Hill K, Cox B, Ulloa F, Ribeiro A, Mynett A, NovitchG B, Briscoe J: Interpretation of the sonic hedgehog morphogen gradient by a temporal adaptation mechanism. Nature. 2007, 450: 717720. 10.1038/nature06347.
 68.
Balaskas N, Ribeiro A, Panovska J, Dessaud E, Sasai N, Page K, Briscoe J, Ribes V: Gene regulatory logic for reading the Sonic Hedgehog signaling gradient in the vertebrate neural tube. Cell. 2012, 148: 273284. 10.1016/j.cell.2011.10.047.
 69.
Cooke J, Zeeman E: A clock and wavefront model for control of the number of repeated structures during animal morphogenesis. J Theor Biol. 1976, 58: 455476. 10.1016/S00225193(76)801312.
 70.
Cooke J: A gene that resuscitates a theory–somitogenesis and a molecular oscillator. Trends Genet. 1998, 14: 8588.
 71.
Pourquié O: The segmentation clock: converting embryonic time into spatial pattern. Science. 2003, 301: 328330. 10.1126/science.1085887.
 72.
Chipman A, Akam M: The segmentation cascade in the centipede Strigamia maritima: involvement of the Notch pathway and pairrule gene homologues. Dev Biol. 2008, 319: 160169. 10.1016/j.ydbio.2008.02.038.
 73.
Pueyo J, Lanfear R, Couso J: Ancestral Notchmediated segmentation revealed in the cockroach Periplaneta americana. Proc Nat Acad Sci. 2008, 105: 1661416619. 10.1073/pnas.0804093105.
 74.
Sarrazin A, Peel A, Averof M: A segmentation clock with twosegment periodicity in insects. Science. 2012, 336: 338341. 10.1126/science.1218256.
 75.
Oates AO, Morelli L, Ares S: Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock. Development. 2012, 139: 625639. 10.1242/dev.063735.
 76.
Hastings A: Transient dynamics and persistence of ecological systems. Ecol Lett. 2001, 4: 215220. 10.1046/j.14610248.2001.00220.x.
 77.
Hastings A, Higgins K: Persistence of transients in spatially structured ecological models. Science. 1994, 263: 11331136. 10.1126/science.263.5150.1133.
 78.
Veraart A, Fassen E, Dakos V, van Nes E, Lürling M, Scheffer M: Recovery rates reflect distance to a tipping point in a living system. Nature. 2012, 481: 357359.
 79.
Hirota M, Holmgren M, van Nes E, Scheffer M: Global resilience of tropical forest and savanna to critical transitions. Science. 2011, 334: 232235. 10.1126/science.1210657.
 80.
Barnosky A, Hadly E, Bascompte J, Berlow E, Brown J, Fortelius M, Getz W, Harte J, Marquet P, Martinez N, Mooers A, Roopnarine P, Vermeij G, Williams J, Gillespie R, abd C Marshall JK, Matzke N, Mindell D, Revilla E, Smith A, a Hastings: Approaching a state shift in the Earth’s biosphere. Nature. 2012, 486: 5258. 10.1038/nature11018.
 81.
Kloeden P: Nonautonomous Dynamical Systems. Providence. 2011, RI: American Mathematical Society
 82.
Rasmussen M: Attractivity and Bifurcation for Nonautonomouos Dynamical Systems. 2007, Heidelberg: Springer
 83.
Gilbert S: Epel D: Ecological Developmental Biology: Integrating Epigenetics, Medicine, and Evolution. Sunderland. 2009, MA: Sinauer Associates
 84.
Gilbert S: Developmental Biology. Sunderland. 2013, Sinauer Associates: MA
 85.
Kashtan N, Alon U: Spontaneous evolution of modularity and network motifs. Proc Nat Acad Sci. 2005, 102: 1377313778. 10.1073/pnas.0503610102.
 86.
Kashtan N, Noor E, Alon U: Varying environments can speed up evolution. Proc Nat Acad Sci. 2007, 104: 1371113716. 10.1073/pnas.0611630104.
 87.
Crombach A, Hogeweg P: Evolution of evolvability in gene regulatory networks. PLoS Comput Biol. 2008, 4: e100011210.1371/journal.pcbi.1000112.
 88.
Crombach A, Hogeweg P: Evolution of resource cycling in ecosystems and individuals. BMC Evol Biol. 2009, 9: 12210.1186/147121489122.
 89.
Surkova A, Spirov A, Gursky V, Janssens H, Kim A, Radulascu O, VanarioAlonso C, Sharp D, Samsonova M, ReinitzJ J, Manu: Canalization of gene expression in the drosophila blastoderm by gap gene cross regulation. Plos Biol. 2009, 7: e1000049
 90.
Corson F, Siggia E: Geometry, epistasis, and developmental patterning. PNAS. 2012, 109: 55685575. 10.1073/pnas.1201505109.
 91.
Scheffer M, Carpenter S, Lenton T, Bascompte J, Brock W, Dakos V, van de Koppel J, van de leemput I, Levin S, van Nes E, Pascual M, Vandermeer J: Anticipating critical transitions. Science. 2012, 338: 344348. 10.1126/science.1225244.
 92.
Scheffer M: Critical Transitions in Nature and Society. Princeton. 2009, NJ: Princeton University Press
 93.
Wang R, Dearing J, Langdon P, Zhang E, Yang X, Dakos V, Scheffer M: Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature. 2012, 492: 419422. 10.1038/nature11655.
 94.
Alon U: Network motifs: theory and experimental approaches. Nat Rev Genet. 2007, 8: 450461. 10.1038/nrg2102.
 95.
Cotterell J, Sharpe J: An atlas of gene regulatory networks reveals multiple threegene mechanisms for interpreting morphogen gradients. Mol Syst Biol. 2010, 6: 425
Acknowledgements
We thank Gregory Boyle for his help with rendering the supplementary movies, and Damjan CicinSain for programming support. We would also like to thank Nick Monk, and the members of the Jaeger Lab for critical discussions and reading of the manuscript. This work was supported by a ’la Caixa’ fellowship awarded to BV. AC was supported by the BioPreDyn consortium, funded by European Commission grant FP7KBBE20115/289434. The laboratory of JJ is funded by the MECEMBL agreement for the EMBL/CRG Research Unit in Systems Biology. Additional financial support was provided by SGR Grant 406 from the Catalan funding agency AGAUR, and by grants BFU2009 10184 and BFU2012337758 from the Spanish Ministerio de Economia y Competitividad (MINECO).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
BV, AC and JJ conceived and designed the study. BV implemented computational tools and performed computer simulations. BV, AC and JJ analysed results. BV and JJ wrote the paper. All authors read and approved the manuscript.
Authors’ original submitted files for images
Additional file 1: Supporting Movie S1–Transition. A transition indicates the switch of a nonautonomous system from one attractor state to another. This movie shows the simulation that Figure 4 is based on. It displays the changing potential landscape as the system transitions from the bistable to the tristable and then the monostable regime (see Figure 2) as autoactivation thresholds are increased. The corresponding phase portrait is shown in the inset on the upper right. Parameter values are displayed in the panel on the lower right. The current state of the system is shown as a black ball. Attractors are marked by blue dots, saddle nodes by red dots with associated separatrices shown as black lines. The history of the trajectory is shown on the phase portrait, colour coded according to which attractor the trajectory is converging towards at every individual step of the simulation. See main text and Figure 4 for details. (MP4 15 MB)
Additional file 2: Supporting Movie S2–Pursuit stabilising the direction of a trajectory. Pursuit behaviour results from the movement of an attractor. In this case, both attractor and trajectory progress in the same direction, which is therefore stabilised even before the system enters any asymptotic regime. This movie shows the simulation that Figure 5 is based on. It displays the changing potential landscape as activation strength is increased. The corresponding phase portrait is shown in the inset on the upper right. Parameter values are displayed in the panel on the lower right. The current state of the system is shown as a black ball. Attractors are marked by blue dots, the saddle node by a red dot with its associated separatrix shown as a black line. The history of the trajectory is shown on the phase portrait. See main text and Figure 5 for details. (MP4 11 MB)
Additional file 3: Supporting Movie 3–Pursuit altering the direction of a trajectory. Pursuit behaviour results from the movement of an attractor. In this case, the direction of the trajectory is altered since the attractor moves against the flow. This movie shows the simulation that Figure 6 is based on. It displays the changing potential landscape as activation strength is decreased. The corresponding phase portrait is shown in the inset on the upper right. Parameter values are displayed in the panel on the lower right. The current state of the system is shown as a black ball. Attractors are marked by blue dots, the saddle node by a red dot with its associated separatrix shown as a black line. The history of the trajectory is shown on the phase portrait. See main text and Figure 6 for details. (MP4 9 MB)
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Saddle Point
 Phase Portrait
 Unstable Manifold
 Transient Behaviour
 Capture Event