Parameter asymmetry and time-scale separation in core genetic commitment circuits

Abstract

Theory allows studying why Evolution might select core genetic commitment circuit topologies over alternatives. The nonlinear dynamics of the underlying gene regulation together with the unescapable subtle interplay of intrinsic biochemical noise impact the range of possible evolutionary choices. The question of why certain genetic regulation circuits might present robustness to phenotype-delivery breaking over others, is therefore of high interest. Here, the behavior of systematically more complex commitment circuits is studied, in the presence of intrinsic noise, with a focus on two aspects relevant to biology: parameter asymmetry and time-scale separation. We show that phenotype delivery is broken in simple two- and three-gene circuits. In the two-gene circuit, we show how stochastic potential wells of different depths break commitment. In the three-gene circuit, we show that the onset of oscillations breaks the commitment phenotype in a systematic way. Finally, we also show that higher dimensional circuits (four-gene and five-gene circuits) may be intrinsically more robust.

Graphical abstract

Keywords

systems biology / theoretical biology / gene regulation / nonlinear dynamics / stochasticity

INTRODUCTION

Mathematics sheds light on commitment and evolution

The realization that bifurcation theory underlies Conrad Waddington’s historical metaphor for cellular differentiation [ 1] adds mathematical substance to the original insight; see Ferrell in [ 2]. There, the author stresses the need to identify what type(s) of bifurcation(s) underlie commitment. We use the word “commitment” for the notion that a dynamical system permits, affects and maintains a chosen state from a set of alternates. At issue is whether alternate fates are still available to a committed cell beyond a bifurcation. In the case of a supercritical pitchfork bifurcation [ 3], alternative states persist past the decision point (Waddington’s original conjecture). However, in the case of a saddle/node bifurcation, the alternative state disappears. Ferrell argues for the latter [ 2]: Nature picks a framework.

The related notion that dynamical systems theory is also offering a plausible unifying framework for understanding the evolution of gene regulatory systems seems very deep. This aspect was recently reviewed by Jaeger and Monk [ 4]. Thus, mathematics underpins the process of not only understanding “how” living systems work, but much more interestingly, “why” they do so in particular ways, rather than others. Hence, insight is to be gained through the study of not only the genetic circuits implemented by Nature, but from others as well, so as to address the “why”.

Biochemical noise contributes in subtle ways

In recent years, we have learned that biochemical fluctuations (i.e., “noise”) are an integral part of how living systems not only mechanistically work, but also why Evolution picked certain regulation schemes over others [ 5– 12]. Therefore, it is not only continuous nonlinear dynamics that is relevant, but it is also stochasticity induced by the paucity of some key molecular regulator(s) in the system. Thus, the notion of stochastic potential is of paramount importance (see Wang et al. [ 13– 23]). The role of energy in the establishment and maintenance of living systems at steady-states far from equilibrium is related, but is not addressed here (see Qian et al. [ 24– 29]).

Robustness to parameter asymmetry and time-scale separation induced fragility

In this study, we investigate core commitment gene regulation circuit topologies and the relationship to the commitment phenotype. Recently, much understanding on the relationship of circuit topology and parameter space to delivered phenotype, was acquired via massive computerized search and analysis that, in essence, can be argued to recapitulate Evolution in-silico [ 30, 31]. Herein however, the focus is on a few key circuit topologies rather than massive exploration, relying heavily on bifurcation analysis for effective parameter space exploration. Moreover, the central driving idea of the present work is the realization that Nature does not appear to fine tune parameters of gene regulation circuits to achieve phenotype. Therefore, some degree of robustness can be expected to exist in order to effectively fend off phenotype-spoiling effects induced by unavoidable parameter asymmetry. Parameter asymmetry is understood to be the level of inequality of strengths in otherwise corresponding but directionally opposite regulation branches (edges linking nodes) of a commitment circuit. The current work thus attempts to parcel out how circuit topology may assist in providing robustness against asymmetry and what the role of noise might be, in the light of time-scale separation between translation and transcription, the two basic underlying manifolds of the non-linear dynamical systems involved.

RESULTS

The core genetic commitment circuits focused on are shown on Figure 1. The circuit shown on Figure 1A, which will be referred to as DS1, is a canonical mutually repressive positively auto-regulating commitment switch. The circuit shown on Figure 1B, DS2, is the generalization of DS1 to three genes. The circuits on Figure 1C, DS3, and Figure 1D, DS4, are progressive generalizations of the same. For all circuits, mRNA transcription regulation is assumed to be additive. This means the effect of two different transcription factors is modeled by the addition of their respective terms. Why focus on these circuits? DS1 is a simple two mutually exclusive gene circuit with positive parity that yields bi-stability [ 32]. The concept of bi-stability underlies β-galactosidase induction in bacteria, lysis-lysogeny decision in bacteriophage λ and maturation of Xenopus oocytes, and cell differentiation [ 2] and references within. DS2, DS3 and DS4 are the simplest tractable systematic generalizations of DS1 that possess the same requirements for multi-stability. Since the purpose of the study is theoretical, no requirement is placed at the outset that these circuits must be found implemented exactly as such in Nature or not, to justify investigating what features in their dynamics could make them suitable or unsuitable. Discovering reasons supporting either outcome is interesting.

The regulation of DS1

Parameter-asymmetry breaks the commitment phenotype

Figure 2A and 2B shows a bifurcation diagram of “2D” DS1 under two basic modes of regulation. Here “2D” refers to the infinite time-scale separation version of the circuit. In the 2D circuit, there are only two state variables, the protein levels X₁ and X₂. The first mode of regulation is one in which there is no parameter asymmetry in the circuit (bmod=1, continuous lines). The second mode of regulation is one where there is a modest 10% parameter asymmetry in the circuit (bmod=0.9, dash-dot lines). Here, “bmod” is a helper multiplicative parameter that acts in only one branch of the circuit, and not the other. Thus, it acts asymmetrically by multiplying the mutual repression parameter “b” of that branch only. This kind of parameter is a convenient device to make the analysis more straightforward. Herein it is referred to as “asymmetry provider” parameter.

It is evident from Figure 2 (A and B) that just a modest amount of parameter asymmetry (here, 10%) completely spoils the delivery of the phenotype. This core regulatory circuit is well-known [ 2]. In the parameter-symmetric version (bmod=1) of DS1, the circuit delivers bi-modal commitment beyond a bifurcation, here located at b~0.3. Beyond the bifurcation, stability of the steady state switches from stable to unstable, and the dynamics opens up to two alternative stable states: one at high protein expression and the other one at low protein expression. In terms of the underlying stochastic potential, out of one valley are born two valleys separated by an intervening ridge appearing at the bifurcation.

However, in the parameter-asymmetric version of DS1 (bmod=0.9), the delivery of the phenotype is clearly broken. As can be seen on the diagram, over a wide range of the repression parameter “b” up to about b~.38, only one high X₁ (or one low X₂) steady state is available in the dynamics. But further above b~.38, beyond a saddle/node bifurcation topologically disconnected from the main branch, the dynamics reverts to one similar to the parameter symmetric case: high and low steady states separated by an unstable one. On Figure 2 (C and D), the equivalent “4D” DS1 model bifurcation behavior is shown. In the 4D version of DS1, there are 4 state variables: the protein levels X₁, and X₂, and the messenger RNA levels mRNAX₁ and mRNAX₂. Apart from the occurrence of unphysical states (paired negative-concentration unstable steady states) that can be ignored since they have no biological meaning, the dynamical behavior is the same as the 2D version.

Time-scale separation breaks the commitment phenotype

Figure 3 shows the phase diagram of DS1. Tracks are numerically integrated starting from different initial conditions S1 and S2 on the phase plane. From location S2, both the 2D and 4D tracks head for stable fixed point #3. From S1 however, the 2D track reaches fixed point #1, but the corresponding 4D track reaches fixed point #3. The simulation was performed for a parameter asymmetry of 0.9. Compared to the 4D version of DS1, the 2D version has infinite intrinsic time-scale separation between the mRNA and protein manifolds; the 2D version is obtained by setting the mRNA dynamics at rest in the 4D version. These results reveal a key biologically significant point: intrinsic time-scale separation between translation and transcription is a handle on the commitment phenotype.

Noise and time-scale separation break the commitment phenotype

Figure 4 shows how intrinsic noise impacts commitment in DS1. The noise is sourced only in biochemical fluctuations induced by the small size of molecular populations in the system. Figure 4 (A and B), each show the average residency metric from one hundred independent stochastic simulations. Residency measures the probability of the stochastic system to be at a certain value of X₁ and X₂ (see Methods). The higher the residency, the higher the probability and consequently, the lower the stochastic potential well depth. Plots are base 10 logarithmic in residency. DS1 is simulated in low time-scale separation (TSS) (configuration #5) on Figure 4A, and in high time-scale separation on Figure 4B (configuration #5.2). Simulations are all started near the upper fixed point. The study reveals that, in low time-scale separation, the depths of the stochastic potential wells underlying the two fixed points differ drastically. As shown on Figure 4A, the upper well is much deeper than the lower well. Thus, in low time-scale separation, the upper fixed point dominates the commitment phenotype. In contrast, in high time-scale separation (Figure 4B), the alternative fate is almost as likely. Although the simulations are still all started near the upper fixed point, in high time-scale separation, the DS1 dynamics populates each well more equally. Hence, in high time-scale separation, both fixed points have similar well depths and similar residency. In Figure 4A, the location of the 2D and 4D deterministic separatrix (see Methods) differs markedly but on panel B they are seen to overlap. The 2D system has, by construction, infinite time-scale separation compared to the 4D system (see Methods).

Thus this study reveals an unintuitive aspect of stochastic dynamics absent from deterministic dynamics: the underlying stochastic potential well depth associated with fixed points of the dynamics depends strongly on the time-scale separation present in the system. Time-scale separation is a handle on the commitment through the depths of the stochastic potential wells. The depth of a stochastic potential well greatly impacts the trapping ability of a realistic biological system in any dynamical state (here static but elsewhere, dynamic). This aspect of stochastic dynamics was encountered previously, albeit in a completely different context, as described in previous work [ 11, 12].

A moving separatrix breaks the commitment phenotype

Figure 5 summarizes the effect of changing the mutual repression parameter “b” and “bmod” in the system. On the four panels, the results of a numerical search algorithm (see Methods) used to locate the position of the deterministic separatrix (2D and 4D) as function of “b” and asymmetry parameter “bmod” are reported. In the case of the parameter-symmetric circuit (bmod=1, panel B and D), no matter the value of “b”, the 4D and 2D separatrix always overlap. Furthermore, the separatrix remains immobile, always bisecting the phase plane at π/4.

However, in the case of the parameter-asymmetric circuit (bmod=0.9, panel A and C), particularly at lower values of “b”, the 4D and 2D separatrix differ significantly. Specifically, we see that the separatrix is no-longer immobile with “b” changing; it no-longer simply bisects the phase plane at π/4 for all “b”. Instead, Figure 5 shows that as “b” increases, the separatrix location moves significantly over the phase plane. Why is this significant? The reason is because the separatrix is the boundary between the two basins of attraction in the system: a moving separatrix will result in a switchover between alternate fates, as the controlling parameter(s) (“b”, “bmod” and time-scale separation) is (are) changed, with all other conditions remaining unchanged. This was verified to be the case (data not shown for brevity). Note that the stochastic system is intrinsically 4D, and thus the 4D separatrix compares to the stochastic system, and does so exactly in the limit of infinite number of molecules. These studies demonstrate that both parameter asymmetry and time-scale separation together impact DS1 commitment in un-intuitive ways.

The regulation of DS2

The parameter-symmetric DS2

New simulation and analysis software necessary to conduct studies of circuit DS2 (and others) was built (see Methods). Figure 6 shows the X₁, X₂, X₃ phase volume of DS2. The surfaces shown are three null-surfaces of the X₁, X₂ and X₃ manifolds, respectively. A stochastic simulation in the volume is also shown. The infinite time-scale separation version of DS2 has three state variables: X₁, X₂ and X₃, thus it is three dimensional (3D). In contrast, the 6D version has both transcription and translation manifolds and six state variables (X₁, X₂, X₃, mRNAX₁, mRNAX₂ and mRNAX₃) thus it is referred to as the 6D model. All stochastic simulations intrinsically use a 6D model. Figure 6B show the “generalized separatrix” surface that delineates the various basins of attractions of the system. Specialized software to scan the phase volume and discover the location of this surface was developed (see Methods). Discovery of this surface requires sampling a discretized grid. The algorithm is highly computer time intensive. We note that the separatrix surface is comprised of three planar structures, each one made out of three leaves. The shape of this surface changes with asymmetry. More will be discussed below.

A new computational approach to study parameter asymmetry is needed

To begin a study of parameter asymmetry, the bifurcation diagram of the 3D version of DS2 vs. “b₂₁” was computed (data not shown for brevity). Parameter “b₂₁” is one of the several similar modifiers of repression parameters in the circuit that are specifically designed to study parameter-asymmetry. These parameters will be referred to as “asymmetry provider” parameters. The DS2 bifurcation behavior vs. “b₂₁” (or any other asymmetry provider parameter) is extremely intricate; the dynamics exhibits multiple saddle-node bifurcations. In order to better understand the DS2 circuit, it was quickly realized that new efficient numerical algorithms needed to be developed (see Methods). These algorithms were employed to discover the intersections of the three null-surfaces. At these intersection points, the global dynamics of the circuit is at rest: these are the fixed points of the dynamics. Since the surfaces are very intricate and the fixed point locations are often rendered difficult to discover due to very complicated and fast changing (with X₁, X₂, X₃ coordinates) null-surface shapes, a brand new and efficient method of fixed point finding was also devised (see Methods). As shown on Figure S1 (for which b₂₁=0.5), the location of the intersection of any two null-surfaces is first discovered numerically and is graphically denoted by elongated clouds of locator points. These intersections are referred to as “bi-null-surface intersections”. Any global fixed point of the dynamics must occur at the intersection of any two different such bi-null-surface intersections e.g., the “12”—intersection and the “13”—intersection together locate “123”—global fixed points. Thus the three possible bi-null-surface intersections (“12”, “13” and “23”) will redundantly locate the global fixed points of the dynamics. Stated differently, the global fixed points arise where two different-color clouds of locator points intersect, and the three possible different-color clouds will all intersect together at the fixed points of the dynamics.

Discrete parameter asymmetry in DS2

The computational approach described above is particularly useful to elucidate saddle-node bifurcations in the variation of “b₂₁”. Figure S1B shows that one of the bi-null-surface intersection clouds (black color) diverges up and away from the other two, thus causing the loss of a stable fixed point. Hence, the approach imparts new understanding of the effect of parameter asymmetry in the circuit. On Figure S1 (b₂₁=0.5) there are only two stable fixed points, rather than three on Figure 6 (b₂₁=1).

At a different higher value of asymmetry provider parameter “b₂₁”, b₂₁=0.85, as shown on Figure S2, the unique saddle of Figure S1 splits into two; so now two saddles co-exist in the system. The bi-null-surface intersections are particularly useful to better understand the changing dynamical topology of the system, thus yielding additional insight to the dynamics of the DS2 circuit in this new regime.

A common asymmetry-provider: “bcommon”

Because of the high complexity of studying the behavior of DS2 vs. every asymmetry-provider parameter in the system, the three key asymmetry providers were regrouped under a single asymmetry provider called “bcommon”. Specifically, b₂₁=b₁₃=b₃₂=bcommon. Figure 1B shows that “bcommon” follows the continuous directed clockwise path in the circuit that involves all state variables in a row. This method of installing common asymmetry will be emulated in other circuits so as to permit meaningful comparisons.

Dynamical behavior of DS2 vs. “bcommon”

In order to further understand the dynamical behavior of DS2 vs. asymmetry-provider parameter “bcommon”, the bifurcation behaviors of DS2 6D and of DS2 3D were computed. The results are shown on Figure 7. The 3D calculations (Figure 7C) reveal three basic regimes: oscillatory to the left, oscillation free in the center and oscillatory on the right. There are two Hopf bifurcations with corresponding limit cycles, and multiple saddle-node bifurcations in the center. On the left, the oscillatory regime terminates at homoclinic orbit/saddle node collisions (α, β). On the right, with increasing “bcommon”, the amplitude of the oscillatory regime steadily decreases from the (β, γ) homoclinic orbit/saddle node collision points, to eventually vanish at the second Hopf. The behavior seen in the 6D circuit (Figure 7A and B) is more intricate. From the first Hopf on the left arises a large unstable limit cycle that does not extend in the positive range of “bcommon”. From the second Hopf arises a stable limit cycle that ends in the homoclinic orbit / saddle node collisions (α, β). Markedly different from the 3D behavior however, is the oscillatory regime on the right that arises from the (β, γ) collision points but with oscillation amplitude ever growing with increasing “bcommon”. Thus 6D differs from 3D principally in that oscillations persist as “bcommon” increases.

In the paragraphs below, the impact of stochasticity and time-scale separation is studied as a function of parameter asymmetry “bcommon”.

Stochasticity and time-scale separation break commitment in DS2

Figure S3 shows the dynamical behavior of DS2 at bcommon=0.85. In this figure, “2” refers to a complex dynamical region, to be discussed later. Figure S3C shows details of stochastic and deterministic simulations. Depending on different initial seeds, different stable steady states (#3 and #4) are reachable from the same initial location in the phase volume. Stochastic run #3 and #4 reach steady state #3 and #4, respectively. This is completely consistent with finite stochasticity. Furthermore, depending on the amount of time-scale separation intrinsic to the 6D deterministic system, one particular stable steady state or another can be reached. Specifically, starting from the same location in the phase volume, a finite time-scale separation deterministic 6D track reaches stable fixed point #3, while the counterpart infinite time-scale separation 3D track reaches fixed point #1. The 3D deterministic track is, by construction, the infinite time-scale separation version of the 6D track. Effectively, the reachable steady state can be selected by the amount of time-scale separation in the system. Thus, time-scale separation in the DS2 circuit strongly impacts commitment. This is analogous to the situation in DS1.

On Figure 8, three stochastic tracks with differing seeds, but common initial location in the phase volume reach, in turn, any and all of the three stable fixed points. Here again, stochasticity impacts commitment. But in contrast to the situation above, because of sufficient intrinsic time-scale separation in the 6D deterministic system, the 3D and the 6D tracks are very close so both reach the same fixed point.

Parameter asymmetry breaks commitment in DS2

Figure S4 shows several time-series of 6D deterministic tracks with varying amount of asymmetry (different “bcommon”). All tracks start at the same location in the phase volume. Between bcommon=1.1 and bcommon=1.15, a switchover occurs between the upper and lower fixed points. Thus the amount of asymmetry in the system selects the reachable fixed point. The mechanism for the switchover to occur —for breaking commitment— is time-scale separation affecting the separatrix surface (shape and position), initial conditions remaining the same. Hence, this is the corresponding effect to the moving separatrix seen in DS1.

Parameter asymmetry and time-scale separation impact oscillations and break commitment in DS2

Figure S4 also shows that between bcommon=1.175 and bcommon=1.2, there is dynamical regime switchover from “escape to stable fixed points” to “oscillations” i.e. the attractor switches type, from static to dynamic. This is consistent with the 6D and 3D bifurcation analyses shown on Figure 7 (“bcommon” crossing points γ and δ). Crossing the onset of oscillations banishes commitment.

As mentioned earlier, the dynamical region in the center of the DS2 phase volume is complex. Here, this region is studied in more details. At bcommon=1.2, as shown on Figure 9, the system clearly presents a 3D stable limit cycle. The stability of the central fixed point #7 is subtle. The first eigenvalue is real and negative, and the first eigenvector direction is along the unit vector centered on the fixed point. Thus along that direction the dynamics is attractive to the fixed point. The other two eigenvalues are complex conjugate of each other and thus rotation about this first eigenvector arises everywhere else. Direct integration verifies that the dynamics is drawn to the surrounding stable limit cycle. However, since the real parts of the second and third eigenvalues are negative, very close to the first eigenvector, the dynamics must actually be attractive. Consequently, a small diameter unstable spiral must separate the two dynamical regions, but its size is below numerical accuracy for determination. For completeness, Figure S5 shows the agreement of 6D and 3D tracks; specifically showing this 6D limit stable cycle (Figure S5C). The green stochastic track on Panel D has 10^1/2 less intrinsic biochemical noise than the blue track on Figure S5A or S5B, but otherwise is similar. The computation shows that it does not, however, escape to fixed point #2. Instead, it reaches the surrounding limit cycle. Because it is less noisy, it already behaves more like the 6D deterministic track (yellow) that spirals out to the dynamic attractor. In the limit of zero intrinsic noise (thermodynamic limit), the agreement between stochastic and deterministic computations is expected to be exact. Similar to the situation described earlier, the switchover between static and dynamic attractor again breaks commitment.

The DS2 dynamics at bcommon=1.6 is shown in Figure 10. As is expected from the 3D bifurcation analysis shown in Figure 7C, there is no 3D stable limit cycle. This is because this value of “bcommon” lies above the second Hopf in the system. The stable limit cycle has already disappeared. However, as further shown in Figure 10C and 10D, the system does present a 6D deterministic limit cycle. Why is this important? In the absence of a limit cycle, the dynamics will be drawn to a stable fixed point, but in the presence of the limit cycle, starting from within its basin of attraction, the dynamics will instead be drawn to the dynamic attractor. Since the controlling parameter is time-scale separation, time-scale separation therefore breaks commitment.

At bcommon=2.0, as shown in Figure S6, there is only one single global (stable) fixed point of the dynamics, consistent with the single intersection of all the elongated clouds of locator points along the bi-null-surface intersections. Interestingly, the flat remnants of the three saddle regions are still visible near locations α, β, and γ (Figure S6B). Figure S6D shows three 6D deterministic tracks. The first one, (i) is seen to be falling, from outside, into a stable limit cycle. The second one, (ii) is seen to be falling into the same limit cycle, but from inside. The third one, (iii) is seen falling directly into the unique central stable fixed point. Clues to the 6D dynamics can be obtained from the 3D stability analysis. Because the first eigenvalue is real and negative, the system dynamics is therefore attractive to the fixed point anywhere on the first eigenvector. The first eigenvector lies in the unit vector direction, centered on the fixed point. The other two eigenvalues are complex conjugates of each other so, everywhere else, the dynamics exhibits rotation about the first eigenvector. Direct 6D integration shows that the dynamics is drawn to the surrounding limit cycle. As on Figure 9, real parts of second and third eigenvalues are actually negative, thus the dynamics very close to the axis of the first eigenvector must be attractive. Therefore, the existence of a small diameter unstable spiral around the first eigenvector can be inferred, but it is beyond numerical accuracy to determine. At bcommon=2.0, the rotational aspects of the 6D and 3D systems differ markedly. Figure 7 shows that in the low time-scale separation case, the dynamics presents a limit cycle; in the infinite time-scale separation case however, it does not. There, the dynamics presents a single stable fixed point. So again, time-scale separation breaks commitment.

In other regions of DS2, time-scale separation does not break commitment

For completeness, on the other side of the dynamics, at bcommon=0.5, both 3D and 6D bifurcation analyses (Figure 7) present stable limit cycles. The integration of the 3D and 6D systems shown on Figure S7 are in complete agreement with this. Contrasted to the regions above bcommon=1, here, no matter what time-scale separation is, the system presents a limit cycle, so this mechanism will not yield to breaking commitment.

The regulation of DS3 and DS4

The regulation schemes of DS3 and DS4 circuits (Figure 1C and 1D) are progressive generalizations of the two- and three-gene commitment circuits discussed above. Each additional gene introduced in the system is mutually repressive to all the other genes in the circuit, and is positively regulated to itself. All mRNA regulation is additive. In the DS3 circuit, asymmetry was installed in the only possible scheme that maintains its continuous end-to-end application, in a strict clockwise manner, starting from X₁ to any other gene, and beyond, except for one gene (here, chosen to be X₂) that must act to return the clockwise path back onto X₁. This can be seen to precisely mirror the asymmetry pattern applied in DS2. In the DS4 circuit, the same asymmetry scheme is used, except that in this five-gene network, there are two genes that return the path onto X₁: here X₂ and X₅. Thus, in all gene networks, parameter asymmetry was installed in a consistent manner allowing meaningful comparison.

Figure 11 shows the phase volume of DS3 in unit asymmetry (bcommon=1). The colored clouds of dots map the location of the four three-variable (“123”, “124”, “134” and “234”) null-surfaces in the 4D infinite time-space separation system (protein only). Their common intersections (redundantly) define a total of nine fixed points; 4 peripheral stable fixed points, 4 intervening unstable saddles and one central stable fixed point. The black dots are located at the global four-variable (“1234”) fixed points of the dynamics. Local linear analysis is provided on the figure. The generalization of previous circuits is clear. Computational limitations restricting the grid size already begin to affect the practicality of a full fine-grid map of the dynamics. For DS4, the additional dimension (DS4 infinite time-scale separation is 5-dimensional) prevents obtaining a full comprehensive map of the system; only local areas around fixed points are possible.

DISCUSSION

The focus of the present work is the impact of time-scale separation and parameter asymmetry in the mutually repressive strengths between all pairs of genes of fully connected core regulation networks delivering the commitment phenotype. It is interesting to notice that, for example, in the case of the SOS pathway in E. coli, the core genetic regulatory circuit involves as many as nine genes arranged in a fully-connected pattern [ 33]. In the case of the core transcriptional gene regulatory network of human and mouse embryonic stem-cells, the OCT4, SOX2 and NANOG genes are also presenting a fully connected network [ 34]. We wondered why Nature seems to have preferentially evolved more topologically complex core gene regulatory networks over simpler ones. By carefully studying progressively more complex idealized–hence tractable– core gene regulatory networks, we hoped to first, show how commitment is broken by parameter asymmetry and time-scale separation and, ultimately, to shed light on how necessary robustness may be a property of topologically more complex networks.

DS1 commitment circuit

The DS1 circuit is a simple commitment circuit between two fates. Impact of parameter asymmetry is profound. Figure 2 (C and D) clearly demonstrates that the circuit’s ability to deliver the phonotype is not robust to asymmetry. Here the phenotype is understood to be the switch over from mono-stability to bi-stability. Interestingly, when the system is symmetric (Figure 5, B and D), time-scale separation is immaterial; the commitment switchover location between fixed points, is independent of time-scale separation. But more interestingly, in the biologically significant situation of even modest parameter asymmetry (10%), A and C of Figure 5 show that time-scale separation has a profound effect on the switchover commitment location. So the DS1 system is seen to be both fragile to time-scale separation and parameter asymmetry in the sense that it is non-robust to phenotype delivery.

DS2 commitment circuit

Time-scale separation induces large (small) effects on the rightmost (leftmost) limit cycles

The effect of increasing the intrinsic time-scale separation between transcription and translation manifolds was studied via bifurcation analysis on the 6D system. Bifurcation curves were computed for reference, ×2, ×10 and ×1000 time-scale separation. The results are combined on Figure 12. For clarity, an unstable limit cycle that lies entirely in the unphysical negative “bcommon” range is omitted. The graph presents the behavior of the system for all biologically meaningful values of the asymmetry parameter “bcommon”. The dynamical region for which bcommon<1 is referred to as the “leftmost region”; the dynamical region for which bcommon>1 is referred to as the “rightmost region”. Limit cycles terminate at homoclinic/saddle collision points (α, β) on the left, and (γ, δ) on the right.

In the leftmost region, for all values of time-scale separation, there is only one Hopf. In the rightmost region, the number of on-scale Hopf bifurcations depends on the time-scale separation. For reference, ×2, ×10 and ×1000 time-scale separation, there are 0, 2, 2 and 1 rightmost Hopf bifurcations, respectively. In the time-scale separation regimes that present two Hopf bifurcations on scale (×2 and ×10), one of the Hopfs is at lower value of asymmetry than the other. As the time-scale separation is increased, there is little effect on the leftmost limit cycle except that its associated Hopf moves to slightly higher values of asymmetry (towards bcommon =1). However, concurrently, the lowest of the rightmost two Hopf bifurcations and the highest one move apart from each other. The lower bifurcation moves towards lower asymmetry. The higher bifurcation moves towards higher values of asymmetry. Thus, in contrast to the situation on the left, time-scale separation induces a dramatic effect on their limit cycles. In the reference time-scale separation, a limit cycle exists throughout much of the rightmost dynamics. That limit cycle’s asymmetry range of existence is terminated at the lower end by homoclinic orbit / saddle collisions (γ and δ points). At the upper end, the asymmetry range of existence extends all the way to the highest values of asymmetry. As the time-scale separation increases, the limit cycle range of existence splits into two separate regions. The end of one region attaches to the upper Hopf. The end of the other region attaches to the lower Hopf. This creates an asymmetry gap devoid of oscillations in between the two rightmost Hopfs (×2 and ×10). As time-scale separation is further increased, the upper rightmost Hopf bifurcation (along with its attached limit cycle) proceeds to move to higher values of asymmetry (×10) and, eventually, it disappears altogether from range (×1000). Thus, as timescale separation is increased, the large stable limit cycle on the right is effectively removed from the dynamics altogether. However, ×2, ×10 and ×1000 results also show that, upon increasing time-scale separation in the system, the lower Hopf bifurcation on the right of the dynamics slides to lower values of asymmetry, and, as it does, so does the end of its attached limit cycle. This limit cycle ranges from homoclinic/saddle collision points (γ and δ) all the way to the lower of the rightmost Hopf bifurcation location, where it terminates. For ×1000 time-scale separation, the behavior is the same as the behavior exhibited by the 3D system (Figure 7C) which has (by construction) infinite intrinsic time-scale separation.

Time-scale separation leaves the location of fixed points and their stability invariant

Whereas changing intrinsic time-scale separation induces effects on the limit cycles- slight effects on the limit cycle in the bcommon<1 range, and major effects on the limit cycle(s) in the bcommon>1 range- most of the dynamical features of the system remain invariant. In particular, the locations of the four homoclinic orbit/saddle collision points (α, β, γ, and δ) remain invariant. Even more strikingly, the locations of all the fixed points in the system (and their stability) remain invariant. Specifically, in Figure 12, for all time-scale separation, fixed point curves do not change. Thus, time-scale separation only affects the rotational aspects of the dynamics (Hopf bifurcations and their limit cycles are shifted), but it leaves other features invariant.

Time-scale separation impacts commitment by removing oscillations from the rightmost region

As time-scale separation increases, the region above the γ, δ collision points (Figure 12) is of particular interest. As stated above, while the locations of the fixed points and the collision points remain invariant, the rotational aspects of the dynamics markedly change. The large amplitude stable limit cycle asymmetry range of existence splits into two. As time-scale separation is increased, one end of the range moves to higher “bcommon” values and, as time-scale separation is further increased, the associated limit cycle moves completely out of range. Concurrently, the other end of the limit cycle asymmetry range of existence shifts to lower values of “bcommon”. This particular orbit becomes a permanent feature of the high time-scale separation dynamics. Naturally, this limit cycle is also found in the 3D infinite time-scale separation version of the system.

Overall, the most significant effect of increasing time-scale separation on the DS2 circuit dynamics is the effective removal of oscillations over an ever increasing range of parameter asymmetry (controlled by “bcommon”). This happens for “bcommon” above ~1.8 (saddle/node annihilation point). As oscillations are removed, it leaves the system in this region subject to only a stable fixed-point attractor. Thus, in this way, time-scale separation directly impacts commitment. In contrast, in the range of parameter asymmetry where there is co-existence of oscillations and fixed-point attractors (between the γ, δ points and bcommon~1.8) the presence of a dynamic attractor is expected to act as a stochastic dynamic mixer in system. Hence, as time-scale separation increases, the gradual effective removal of oscillations is expected to also remove stochastic mixing effects.

In sections below, the two flanking dynamical regions where oscillations and fixed points coexist are studied in more details.

Parameter asymmetry and time-scale separation impact commitment in DS2

The general separatrix surface delimits the various attractors in the system. Below, several studies reveal how the general separatrix surface changes with parameter asymmetry. Studies of its changing nature shed light on the impact of parameter asymmetry on commitment. Earlier, it was shown that time-scale separation impacts only the dynamic attractors. Therefore these studies also permit delineation of the concomitant impact of time-scale separation on commitment via limit cycles.

The flanking dynamical regions present coexisting oscillations with fixed point attractors

Because time-scale separation particularly influences rotational dynamics in the region of “bcommon” above the γ, δ collision points and below the saddle/node annihilation points (bcommon~1.8), detailed computations of the general separatrix surfaces in two contrasting parameter asymmetry regimes were performed. The first regime investigated is at bcommon=1.2 where both 3D (at infinite time-scale separation by definition) and 6D (set at reference time-scale separation for greater contrast) present limit cycles. The second regime is at bcommon=1.5 where 6D still presents a limit cycle, but 3D does not.

Parameter asymmetry and time-scale separation break commitment where oscillations coexist with fixed point attractors

Figure S8A shows an overlay of the 6D and 3D general separatrix surfaces at bcommon=1.2 (configuration #24; 6D: reference time-scale separation, 3D: infinite time-scale separation). Here, an attractive limit cycle and attractive fixed points of the dynamics coexist both in 6D and in 3D. Because the limit cycles in both finite and infinite time-scale separation are (in this region) of similar size, the two general separatrix surfaces are very close in shape and size. Figure S8B shows the basin of attraction of the dynamic attractor: the location in phase space of all X₁, X₂, and X₃ initial conditions for which the 6D dynamics is attracted to the limit cycle. The initial conditions for falling into the rotational attractor are confined to smaller values of (X₁, X₂) for X₃ small. As the X₃ level increases, the extent of allowed values of (X₁, X₂) grow lesser, matching the triangular shape of the limit cycle. The equivalent 3D figure is omitted because it is very similar: at each increased X₃ level, the extent of allowed (X₁, X₂) values is slightly smaller than that in 6D. The 3D dynamics is the infinite time-scale separation version of the 6D dynamics.

Figure S8 (C and D) show the corresponding situation at bcommon=1.5 (parameterization #24.1; 6D: reference time-scale separation, 3D: infinite time-scale separation). Figure S8C shows the overlay of the 6D and 3D general separatrix surfaces. In 3D, there is no limit cycle because bcommon=1.5 is located beyond the Hopf bifurcation. Irrespective of the lack of limit cycle in 3D, the two surfaces are nevertheless close in shape and size. Thus in 3D, the dynamics switches from the three attractive fixed points (1, 2 and 3) to the center fixed point (7). In 6D however, the switch is to the limit cycle (a stable but dynamic attractor). The stability of fixed point #7 in 6D has been discussed earlier. Except for initial locations along the eigenvector directions for which dynamics falls into #7, elsewhere, dynamics falls into the surrounding limit cycle. Panel D shows the phase volume locations for which the 6D dynamics falls into the limit cycle.

Figure S8 illustrates clearly how commitment to a fixed point changes by way of parameter asymmetry and time-scale separation. We see that the general separatrix surfaces are similar between 3D and 6D in either system. Yet by changing the parameter asymmetry from 1.5 to 1.2, a system committed to a stable state in infinite time-scale separation (3D) may switch to a dynamical state (limit cycle). This switch cannot occur in 6D (reference time scale separation) since both regimes present limit cycles.

The oscillations-free central dynamical region remains unaffected by time-scale separation

This region is of particular biological interest from the point of view of the commitment phenotype. Figure S9 shows an overlay of the generalized separatrix surfaces over the central dynamical region of DS2, computed for the reference time-scale separation case. However, since there is no oscillations in this region, time-scale separation does not affect it (see previous sections). Therefore, the computation is valid for all time-scale separations. In the symmetric case (bcommon=1), the separatrix surface is a striking three-leaved multi-planar surface. Each of its three leaves is, in turn, comprised of three layers. The two outer layers delimit the basin of attraction of the nearest surrounding stable fixed points. The middle layer is the attractive manifold of the central fixed point of the dynamics, in this case a saddle (in the symmetric case only). On Figure S9, one of these central layers is indicated by an arrow. The other two are omitted because they are hidden from this viewpoint. In the two asymmetric cases however (bcommon =0.79 and 1.15), the three leaves of the general separatrix surfaces are not planar and each one is single-layered. This is due to the fact that the central fixed point, in these cases, is an unstable spiral, so it has no attractive manifold.

In the oscillations-free central region, parameter asymmetry switches commitment among fixed point attractors by changing the shape of the general separatrix surface

The central dynamical region is devoid of oscillations. The shape deformation of the general separatrix surface with changing parameter asymmetry (“bcommon”) in DS2 is reminiscent to the situation seen in the two-gene version of the circuit (DS1). But here in DS2, we see that parameter asymmetry induces both a rotation of the surfaces about the eigenvector of the central fixed point, and a deformation of separatrix planes into curved surfaces. The changing surfaces with parameter asymmetry directly impact commitment.

In the flanking dynamical regions, noisy oscillations and parameter asymmetry break commitment

The flanking dynamical regions where oscillations coexist with stable fixed points are of particular interest because they reveal the manner in which the commitment phenotype is broken by parameter asymmetry. Figure 13 shows the behavior of ten independent stochastic tracks, all started at a common location in phase space: precisely on the axis of the eigenvector of the central fixed point (#7). The expectation is that, because of noise, rather than falling into the central fixed point, the tracks should pick up rotation and (noisily) coalesce into the limit cycle, thus becoming trapped in the dynamic attractor. However, the same noise that permits the dynamics to escape the attractive fixed point along the axis of the eigenvector also lets dynamics escape from the dynamic attractor. One might expect that all three fixed points should have equal probability of being selected, thus noise not selecting commitment. Something more subtle is observed however. In fact, at low bcommon=.7250 and high bcommon=1.2, noisy dynamics selects different fixed points, #3 and #1, respectively.

The reason for this unintuitive behavior is that, at finite time-scale separation, the axis of rotation of the 6D limit cycle and the 3D eigenvector counterpart do not quite line up, such that the approach to one or another saddle is favored. Stated differently, the protein subspace of the 6D eigenvector is not perfectly aligned with the (protein-only) 3D eigenvector. This, combined with the pitch of the outward spiral leads to preferential selection of one deflecting saddle over another, and consequently of its associated stable fixed point. It is by way of the intervening saddle that the stochastic track is deflected into an associated attractive fixed point. In Figure 13, the majority (six, in both cases) of the ten tracks fluctuates out of the dynamic attractor and heads for a preferred, but asymmetry-dependent, stable fixed point. The four tracks that do remain on the dynamic attractor (limit cycle) have been omitted from the figure, for clarity. This behavior has been verified for higher time-scale separation as well (data not shown for brevity). Thus, these simulations show that, in a noisy dynamical regime, parameter-asymmetry does break the commitment phenotype but, it does so in an unintuitive way. The onset of oscillations could be expected to spread the commitment equally between all three fixed points, but in fact, in this low-noise regime, because of the intricate noisy dynamics involved, commitment becomes tied to the amount of parameter asymmetry in the system. Note that in the noise-free system, tracks initiated precisely along the eigenvector belonging to fixed point #7 of Figure 13, as expected, head for the fixed point (data not shown for brevity).

More noise, more surprises: oscillations entrain noisy dynamics to break commitment

Figure 14 presents the dynamical behavior of the system for 10x the noise level, again at bcommon=.7250 (Figure 14A) and also at bcommon=.9 (Figure 14B). This noisy regime is achieved by scaling down the number of molecules and the volume of the system both by the same factor of 100, thus leaving biochemical concentrations unaltered. From Figure 12A, at bcommon=.7250 there is coexistence of a limit cycle with three attractive fixed points. But at bcommon=.9, the dynamics is absent of the limit cycle, presenting only the fixed points. The combination of a noisy regime with the driving effect of the limit cycle is striking. As seen on Panel A, one single stochastic trajectory visits all three fixed points in the system, systematically, in a counter-clockwise manner matching the direction of rotation of the limit cycle. This behavior was verified by following the stochastic track in time (data not shown for brevity). The stochastic dynamics is not composed of random hopping from one fixed point to any of the other two located on either side of it. Instead it is observed to be driven in a counter-clockwise cyclic dynamics. Specifically, the stochastic dynamics settles in one basin of attraction for some time, until it hops to the next one over, entrained by the nearby limit cycle, in a strict counter-clockwise direction. This process is repeated over and over again. Thus, presence of the parameter-asymmetry induced limit cycle breaks the commitment phenotype. This is further verified on Panel B which shows a close-by dynamical regime that lacks any limit cycle. There, 10 statistically independent tracks each reach one, and only one, of the fixed points. Once the dynamics has entered the basin of attraction of a particular fixed point, it permanently remains in the basin of attraction of that fixed point. Hence, in this contrasting situation lacking any limit cycle driving effect, the commitment phenotype is preserved.

Note that from an inspection of Panel A, one should expect that whether the commitment phenotype would be broken by the presence of the asymmetry-induced limit cycle and biochemical noise should depend on the relative sizes of the limit cycle and of the noise clouds over the basins of attraction of the fixed points. Further simulations (not shown for brevity) confirm this.

Figure 15 shows the dynamical behavior of the DS2 system for 10× the noise level, but at bcommon=1.6, hence on the other –higher- parameter-asymmetry side of the dynamics. This is also a regime for which there is coexistence of the three stable fixed points with a 6D limit cycle, but here in contrast, there is no corresponding 3D cycle. As in the previously studied case, the limit cycle not only causes the dynamics to break the commitment phenotype, but it also induces significant dwelling near the central region. In effect, the limit cycle induces the dynamical maintenance of an additional “quasi-stable” central fixed point. The stochastic trajectory first accesses the central “quasi-stable” fixed point region, dwells there for some time, then it hops over to the neighborhood of stable fixed point #3. After some while, the limit cycle influences the dynamics to hop, in clockwise order (matching the rotational direction of the limit cycle), first to stable fixed point #2, then on to stable fixed point #1, and then back on to stable fixed point #3 again. The hop order is: 3→2→1→3.

In summary, the presence of a parameter-asymmetry induced limit cycle and sufficient biochemical noise together break, in a time-ordered systematic manner, the simple commitment phenotype otherwise imparted by the topology of the circuit. Further, because time-scale separation controls the onset of oscillations, it also breaks the commitment phonotype by initiating the time-ordered visitation sequence. Biological conditions are neither expected to be parameter-symmetric, be noiseless or to present infinite time-scale separation.

The dynamics of the DS3 and DS4 commitment circuits

DS3 and DS4 are the four- and five-gene generalizations of the commitment circuits studied above. It is appropriate to discuss them together in the light of how oscillations might affect the commitment phenotype. Except for “bcommon” adjusted as discussed below, all other parameters remain the same as in DS2, configuration #12.1 (finite, reference time-scale separation).

Figure 16 shows a bifurcation diagram of DS3 8D, in reference finite time-scale separation. The dynamics presents two separate (non-overlapping) dynamical regimes. At lower asymmetry (below bcommon~1.44), saddle/node bifurcations arise. At much higher asymmetry (bcommon~6.24), a Hopf bifurcation gives rise to a stable limit cycle that quickly grows to large amplitude, as asymmetry increases. However, there is no overlap of limit cycle rotation with the dynamical regime of fixed points isolated by saddles. Hence, the commitment phenotype is un-impaired by oscillations. As time-scale separation increases to infinity, the limit cycle disappears from range (data not included for brevity). Figure S10 shows the situation at asymmetry bcommon=10. The reference finite time-scale separation 8D stable limit cycle is clearly visible. But the corresponding infinite time-scale separation 4D system is devoid of limit cycle; so in this situation, the dynamic attractor is replaced by a stable attractor and the 4D track spirals down to it.

Figure 17 shows a bifurcation diagram of the DS4 10D system, in reference time-scale separation. As is the case for the DS3 commitment circuit, the DS4 system presents a limit cycle that occurs at relatively high asymmetry (the Hopf bifurcation is located at bcommon~8.39, higher than for DS3). Again, this rotation does not overlap with a lower-asymmetry dynamical regime (located in the area from bcommon~1.3 to bcommon~ 2.1) that is dominated by saddle/node bifurcations. Thus again, in DS4 as is the case in DS3, rotation is incapable of spoiling the commitment phenotype of the circuit.

Interestingly however, DS4 presents a particularly un-attractive dynamical feature from the point of view of commitment. As shown in Figure 17B, a dynamical fold affords the commitment phenotype an opportunity to jump from a single stable fixed point (bcommon~1.79) to two stable fixed points, separated by a saddle (bcommon~1.76). This occurs because of the presence of a saddle/node bifurcation. Deterministic tracks were integrated for a system located just before the bifurcation (bcommon ~1.76) and another, just after the bifurcation (bcommon ~1.79). The data is omitted for brevity. In the former, tracks are attracted to the two distinct fixed points located on either sides of the saddle; in the latter, passed the saddle/node bifurcation, all tracks are attracted to the one residual single fixed point. Commitment is thus found fragile to parameter asymmetry because passing the bifurcation leads to very large jumps (here, 30% and 50%) in the expressed value of two out of three genes studied. A similar situation occurs when asymmetry drops from ~1.8 to ~1.76 (data also omitted for brevity), crossing another saddle/node bifurcation. The dynamical region where the fragility occurs is relatively small.

CONCLUSION

Parameter asymmetry breaks the commitment phenotype in DS1 by causing bimodality to become topologically disconnected from unimodality. Parameter asymmetry causes the appearance of a dynamical gap in the otherwise bimodal behavior of DS1. Figure 2 (C and D) shows that from a repression of b ~.3 to b ~.38, the circuit that exhibited bimodal behavior in the parameter symmetric case, is induced to remain high/low unimodal even in the mild 10% asymmetric case studied. Hence, from the biological standpoint, delivery of the commitment phenotype by the DS1 circuit is fragile with respect to parameter asymmetry. The same is found to be true due to time-scale separation (Figure 3 and Figure 5). Noise also breaks commitment (Figure 4).

Parameter asymmetry breaks the commitment phenotype in DS2 as well, but it only does so through the appearance of oscillations. Oscillations only matter to phenotype delivery in dynamical regions where it coexists with stable steady states. This happens from bcommon ~.7 to bcommon~.76 for all time-scale separations, and from bcommon~1.18 to bcommon~1.3 in the infinite time-scale separation case, or to bcommon~1.82 in the finite time-scale separation case. In the coexistence range for asymmetry less than unity, time-scale separation does not affect oscillations markedly. Only for asymmetry greater than unity does it do so markedly. But, remarkably, the central dynamical region (from bcommon ~.76 to bcommon ~1.18) is unaffected by oscillations. Therefore, in the central region, the commitment phenotype is immune to parameter asymmetry altogether. Since time-scale separation only affects oscillations, the commitment phenotype in the central region is also immune to time-scale separation. Hence the DS2 circuit is more robust to both parameter asymmetry and time-scale separation than DS1.

The DS3 and DS4 circuits both present oscillations at high asymmetry, but this rotation never co-exists with the dynamical regime of the commitment phenotype located at much lower asymmetry. Hence, contrary to DS2, oscillations cannot spoil the commitment phenotype. Furthermore, the infinite time-scale separation version of these circuits does not present any oscillations (on the scale of asymmetry studied). Thus, the DS3 and DS4 commitment phenotype delivery ability is immune to time-scale separation. So, both DS3 and DS4 can be thought of as more robust to parameter asymmetry than DS1 and DS2. In the case of DS4 however, this finding is somewhat tempered by the appearance of the fold effect depicted on panel B of Figure 17. The dynamical region where fragility appears is however small.

In this work, prompted by the observation that Nature seems to have preferentially evolved complex fully-connected gene regulatory networks [ 33, 34], we focused on the following question: From the point of view of the commitment phenotype delivery, could robustness to parameter asymmetry and time-scale separation be improving systematically with increasing network topology complexity? Stated differently, is fragility a feature of simpler networks? The current work presents evidence this is so. However, further studies on more complex topologies are needed to settle the question.

It is also interesting to reflect on the suggestive finding that oscillations sourced in parameter asymmetry combined with biochemical noise, together lead to systematic breaking of the commitment phenotype. In other words, natural conditions i.e. non-equality of gene suppression strengths combined with unavoidable biochemical noise together lead to systematics to be observed in the commitment time history of the system, thus effectively amounting to a possible mechanism for progressive differentiation.

The work presented here focused on the strength asymmetry between opposing repression forces in all the possible pairs of genes within relatively simple canonical fully-connected core regulatory gene networks delivering the commitment phenotype. As shown above, this already provides a richness of dynamical and stochastic behaviors deserving investigation because of their relevance to biology. Absent from the current concern however, is the role of auto-activation; all circuits are endowed with the same self-regulation strength on all their respective constituent genes. Not considered here also, is any asymmetry in the degradation strengths. These topics, and others, require much further study and will be the subject of future work.

MATERIAL AND METHODS

Null-surfaces, and bi-null-surface intersections

For any system studied herein, each of the coupled nonlinear differential equation expresses the time rate of change of one the state variables in terms of all the state-variables in the system. When set to zero, each differential equation therefore defines one of the steady-state manifolds of the system. In this work, each steady-state manifold is referred to as a “generalized null-surface” in analogy to the common two-dimensional system’s “nullcline”. Each generalized null-surface equation is however much too complex to solve analytically. Therefore, in order to find the generalized null-surface locations within a phase volume, or to find the locations of their mutual intersections, a special numerical algorithm was devised. For example, in the X₁, X₂, X₃ system, there are three generalized null-surfaces: f₁(X₁, X₂, X₃)=0, f₂(X₁, X₂, X₃)=0 and f₃(X₁, X₂, X₃)=0. The three surfaces are discovered numerically by finding the sets of X₁, X₂ and X₃ locations for which, in turn, f₁<ϵ, f₂<ϵ, and f₃<ϵ.

In this work, it was realized that much subtle and intricate information about dynamics was better rendered graphically by plotting the bi-null-surface intersections, in addition to the generalized null-surfaces. For example, here in the X₁, X₂ and X₃ system, the “12”, “13” and “23” intersections denote the 1 and 2, 1 and 3 and 2 and 3 null-surface intersections, respectively. The intersections are found numerically in a similar manner as the actual surfaces. They are plotted in the phase volume using different colored clouds of dots creating distinctive swirling patterns. At the bi-null-surface intersections (thus, within a cloud of a single color), two of the three state variables are at rest. Thus, any co-located pair of two-variable intersections (two different colored clouds intersecting) defines a global fixed point of the dynamics for which all three variables are at rest (e.g., “12” and “23” clouds intersecting define a “123” global fixed point). The redundancy offered by using three mutually intersecting clouds is found to be very useful in order to handle the high complexity of the null-surfaces’ changing shapes. Higher dimensionality systems (e.g., four, five) are handled in a straightforward generalization of the numerical approach described above.

Generalized separatrix

In the two-dimensional system, the two stable steady-state basins of attraction are separated by a curve called the separatrix. Starting a trajectory on either side of the separatrix will result in the trajectories reaching alternate fixed points; at the separatrix there is switchover of reached stable fixed point. Since the location on the phase plane where the switchover occurs is not known a priori, it is discovered by systematically starting multiple tracks on a grid of initial conditions covering the phase plane. This concept straightforwardly generalizes to systems of higher dimensionality. There, the separatrix becomes a surface embedded in the phase volume. To discover the location of such a generalized separatrix surface, multiple tracks are systematically started on a grid of initial conditions filling the phase volume. The reached attractor is recorded for each trajectory. The generalized separatrix surface location is therefore approximated by the switchover location on the grid.

Local linear analysis

At the numerically located fixed points of the dynamics, the analytically-determined (using Mathematica: Wolfram, Champaign, Illinois) Jacobian matrix of the system of coupled nonlinear ordinary differential equations is numerically evaluated. Eigenvalues and eigenvectors of the Jacobian matrix are computed using Matlab (The Mathworks, Natick, Massachusetts) and the local stability of each fixed point is determined based on the eigenvalues.

Stochastic simulations

The stochastic simulations were achieved by decomposing the dynamical systems into the set of their elemental discrete production and degradation events, and by randomly actuating them using the Gillespie algorithm [ 35– 38]. The behavior of multiple statistically independent realizations of the system is therefore governed by the underlying governing Master Equation of the dynamics [ 37, 38]. All simulations were developed, performed and analyzed within the Matlab (The Mathworks, Natick, Massachusetts) framework.

DS1 residency diagrams

Two-dimensional histograms of state variables are computed using extended duration stochastic trajectories throughout phase space. Color is used to indicate the base-10 logarithm of the contents that measures the probability of the dynamical system to visit the defined region of phase space. More details can be found in [ 11, 12].

DS1 4D-nullclines

The 4D-nullclines are computed by removing one ODE from the four-dimensional system, and solving the associated homogenous sub-problem. More details can be found in [ 11, 12].

Bifurcation analysis

Bifurcation analyses were performed using the Oscill8 continuation analysis software (http://oscill8.sourceforge.net/ ) and the MATCONT continuation toolbox in Matlab (The Mathworks, Natick, Massachussetts). Continuation of limit cycles was also performed using standalone software written in Matlab.

Model equations; finite time-scale separation version

DS1

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{1}}\times X_{\mathrm{1}}+t_{\mathrm{1}}\times mRNA{X}_{\mathrm{1}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{2}}\times X_{\mathrm{2}}+t_{\mathrm{2}}\times mRNA{X}_{\mathrm{2}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{1}}+a\times {\displaystyle \frac{{X_{\mathrm{1}}}^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{21}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{2}}+a\times {\displaystyle \frac{{X_{\mathrm{2}}}^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}},

DS2

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{1}}\times X_{\mathrm{1}}+t_{\mathrm{1}}\times mRNA{X}_{\mathrm{1}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{2}}\times X_{\mathrm{2}}+t_{\mathrm{2}}\times mRNA{X}_{\mathrm{2}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{3}}\times X_{\mathrm{3}}+t_{\mathrm{3}}\times mRNA{X}_{\mathrm{3}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{1}}+a\times {\displaystyle \frac{{X_{\mathrm{1}}}^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{21}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{31}}

\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{2}}+a\times {\displaystyle \frac{{X_{\mathrm{2}}}^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{32}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{3}}+a\times {\displaystyle \frac{{X_{\mathrm{3}}}^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{13}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{23}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}.

DS3

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{1}}\times X_{\mathrm{1}}+t_{\mathrm{1}}\times mRNA{X}_{\mathrm{1}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{2}}\times X_{\mathrm{2}}+t_{\mathrm{2}}\times mRNA{X}_{\mathrm{2}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{3}}\times X_{\mathrm{3}}+t_{\mathrm{3}}\times mRNA{X}_{\mathrm{3}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{4}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{4}}\times X_{\mathrm{4}}+t_{\mathrm{4}}\times mRNA{X}_{\mathrm{4}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{1}}+a\times {\displaystyle \frac{{X_{\mathrm{1}}}^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{21}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{31}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{41}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{2}}+a\times {\displaystyle \frac{{X_{\mathrm{2}}}^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{32}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{42}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{3}}+a\times {\displaystyle \frac{{X_{\mathrm{3}}}^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{13}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}

+b\times b_{\mathrm{23}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{43}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{4}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{4}}+a\times {\displaystyle \frac{{X_{\mathrm{4}}}^n}{S^n+{X_{\mathrm{4}}}^n}}+b\times b_{\mathrm{14}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{24}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{34}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}.

DS4

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{1}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{1}}\times X_{\mathrm{1}}+t_{\mathrm{1}}\times mRNA{X}_{\mathrm{1}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{2}}\times X_{\mathrm{2}}+t_{\mathrm{2}}\times mRNA{X}_{\mathrm{2}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{3}}\times X_{\mathrm{3}}+t_{\mathrm{3}}\times mRNA{X}_{\mathrm{3}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{4}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{4}}\times X_{\mathrm{4}}+t_{\mathrm{4}}\times mRNA{X}_{\mathrm{4}},

{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{5}}}{\mathrm{d}t}}=-kdeg{X}_{\mathrm{5}}\times X_{\mathrm{5}}+t_{\mathrm{5}}\times mRNA{X}_{\mathrm{5}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{2}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{2}}+a\times {\displaystyle \frac{{X_{\mathrm{2}}}^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}

+b\times b_{\mathrm{32}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{42}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}}+b\times b_{\mathrm{52}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{5}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{3}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{3}}+a\times {\displaystyle \frac{{X_{\mathrm{3}}}^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{13}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{23}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{43}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}}+b\times b_{\mathrm{43}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{5}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{4}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{4}}+a\times {\displaystyle \frac{{X_{\mathrm{4}}}^n}{S^n+{X_{\mathrm{4}}}^n}}+b\times b_{\mathrm{14}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{24}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{34}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{54}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{5}}}^n}},

{\displaystyle \frac{\mathrm{d}mRNA{X}_{\mathrm{5}}}{\mathrm{d}t}}=-kdegmRNA\times mRNA{X}_{\mathrm{5}}+a\times {\displaystyle \frac{{X_{\mathrm{5}}}^n}{S^n+{X_{\mathrm{5}}}^n}}+b\times b_{\mathrm{15}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{1}}}^n}}+b\times b_{\mathrm{25}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{2}}}^n}}+b\times b_{\mathrm{35}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{3}}}^n}}+b\times b_{\mathrm{45}}\times {\displaystyle \frac{S^n}{S^n+{X_{\mathrm{4}}}^n}}.

Model equations; infinite time-scale separation version

In all cases, to obtain the infinite time-scale separation equations, set the time rate of change of mRNA to zero to algebraically solve for the steady-state mRNA. Then, substitute this steady-state mRNA back into the differential equation for each state variable.

Parameter Configurations

The parameter configurations of the various models are detailed in Table 1, Table 2 and in Supplementary Tables.

SUPPLEMENTARY MATERIALS

The supplementary materials can be found online with this article at DOI 10.1007/s40484-015-0042-1

ACKNOWLEDGEMENTS

Thanks to Jianfeng Feng for hosting MT on a visit to Fudan University’s SCMS where some early ideas underlying this work were formed. Also thanks to his student Wenbo Sheng for assistance double checking one early calculation. Thanks to Michael Zhang for his kind hospitality to MT at UT Dallas. Thanks to Zhiyuan Li at UCSF for early discussions. MT designed the studies, wrote the code, performed the simulation, did the analysis, and wrote the paper. HX contributed to the consistency of the work, assisted with performing simulations and performed the migration of the figures to postscript. Many thanks to the anonymous reviewers for their useful suggestions. This work was supported by NSF grant #924296 to MT.

COMPLIANCE WITH ETHICS GUIDELINES

Marc Turcotte declares that he has no conflict of interest. Hongguang Xi declares that he has no conflict of interest. ƒThis article does not contain any studies with human or animal subjects performed by any of the authors.

References

[1] Waddington, C. H. (1957) The Strategy of the Genes. London: Routledge

[2] Ferrell, J. E. Jr. (2012) Bistability, bifurcations, and Waddington’s epigenetic landscape. Curr. Biol., 22, R458–R466

[3] Strogatz, S. H. (1994) Nonlinear Dynamics and Chaos. Cambridge: Perseus Books Publishing

[4] Jaeger, J., Monk, N. (2014) Bioattractors: Dynamical systems theory and the evolution of regulatory processes. J. Physiol., 592, 2267–2281

[5] Çağatay, T., Turcotte, M., Elowitz, M. B., Garcia-Ojalvo, J. and Süel, G. M. (2009) Architecture-dependent noise discriminates functionally analogous differentiation circuits. Cell, 139, 512–522

[6] Elowitz, M. B., Levine, A. J., Siggia, E. D. and Swain, P. S. (2002) Stochastic gene expression in a single cell. Science, 297, 1183–1186

[7] Süel, G. M., Garcia-Ojalvo, J., Liberman, L. M. and Elowitz, M. B. (2006) An excitable gene regulatory circuit induces transient cellular differentiation. Nature, 440, 545–550

[8] Süel, G. M., Kulkarni, R. P., Dworkin, J., Garcia-Ojalvo, J. and Elowitz, M. B. (2007) Tunability and noise dependence in differentiation dynamics. Science, 315, 1716–1719

[9] Thattai, M. and van Oudenaarden, A. (2004) Stochastic gene expression in fluctuating environments. Genetics, 167, 523–530

[10] Turcotte, M., Garcia-Ojalvo, J. and Süel, G. M. (2008) A genetic timer through noise-induced stabilization of an unstable state. Proc. Natl. Acad. Sci. USA, 105, 15732–15737

[11] Xi, H., Duan, L. and Turcotte, M. (2013) Point-cycle bistability and stochasticity in a regulatory circuit for Bacillus subtilis competence. Math. Biosci., 244, 135–147

[12] Xi, H., Yang, Z. and Turcotte, M. (2013) Subtle interplay of stochasticity and deterministic dynamics pervades an evolutionary plausible genetic circuit for Bacillus subtilis competence. Math. Biosci., 246, 148–163

[13] Li, C., Wang, E. and Wang, J. (2011) Landscape and flux decomposition for exploring global natures of non-equilibrium dynamical systems under intrinsic statistical fluctuations. Chem. Phys. Lett., 505, 75–80

[14] Li, C., Wang, E. and Wang, J. (2011) Landscape, flux, correlation, resonance, coherence, stability, and key network wirings of stochastic circadian oscillation. Biophys. J., 101, 1335–1344

[15] Li, C., Wang, E. and Wang, J. (2012) Landscape topography determines global stability and robustness of a metabolic network. ACS Synth Biol, 1, 229–239

[16] Li, C. and Wang, J. (2013) Quantifying Waddington landscapes and paths of non-adiabatic cell fate decisions for differentiation, reprogramming and transdifferentiation. J. R. Soc. Interface, 10, 20130787

[17] Li, C. and Wang, J. (2014) Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle. Proc. Natl. Acad. Sci. USA, 111, 14130–14135

[18] Li, C. and Wang, J. (2014) Quantifying the underlying landscape and paths of cancer. J. R. Soc. Interface, 11, 20140774

[19] Wang, J., Zhang, K., Xu, L. and Wang, E. (2011) Quantifying the Waddington landscape and biological paths for development and differentiation. Proc. Natl. Acad. Sci. USA, 108, 8257–8262

[20] Wu, W. and Wang, J. (2013) Landscape framework and global stability for stochastic reaction diffusion and general spatially extended systems with intrinsic fluctuations. J. Phys. Chem. B, 117, 12908–12934

[21] Wu, W. and Wang, J. (2013) Potential and flux field landscape theory. I. Global stability and dynamics of spatially dependent non-equilibrium systems. J. Chem. Phys., 139, 121920

[22] Xu, L., Zhang, F., Zhang, K., Wang, E. and Wang, J. (2014) The potential and flux landscape theory of ecology. PLoS One, 9, e86746

[23] Zhang F., Xu L., Zhang K., Wang E., Wang J., (2012) The potential and flux landscape theory of evolution. J. Chem. Phys., 137, 065102

[24] Beard, D. A. D., Babson, E., Curtis, E. and Qian, H. (2004) Thermodynamic constraints for biochemical networks. J. Theor. Biol., 228, 327–333

[25] Beard, D. A. and Qian H. (2008) Chemical Biophysics, Cambridge: Cambridge University Press

[26] Qian, H. and Cooper, J. A. (2008) Temporal cooperativity and sensitivity amplification in biological signal transduction. Biochemistry, 47, 2211–2220

[27] Qian, H. (2007) Phosphorylation energy hypothesis: open chemical systems and their biological functions. Annu. Rev. Phys. Chem., 58, 113–142

[28] Qian, H. and Beard, D. A. (2005) Thermodynamics of stoichiometric biochemical networks in living systems far from equilibrium. Biophys. Chem., 114, 213–220

[29] Qian, H., Beard, D. A. and Liang, S. D. (2003) Stoichiometric network theory for nonequilibrium biochemical systems. Eur. J. Biochem., 270, 415–421

[30] Ma, W., Trusina, A., El-Samad, H., Lim, W. A. and Tang, C. (2009) Defining network topologies that can achieve biochemical adaptation. Cell, 138, 760–773

[31] Zhang, J., Yuan, Z., Li, H. X. and Zhou, T. (2010) Architecture-dependent robustness and bistability in a class of genetic circuits. Biophys. J., 99, 1034–1042

[32] Snoussi, E. H. (1998) Necessary Conditions for Multistationarity and Stable Periodicity. J. Biol. Syst., 06, 3–9

[33] Gardner, T. S. and Faith, J. J. (2005) Reverse-engineering transcription control networks. Phys. Life Rev., 2, 65–88

[34] Chickarmane, V., Troein, C., Nuber, U. A., Sauro, H. M. and Peterson, C. (2006) Transcriptional dynamics of the embryonic stem cell switch. PLoS Comput. Biol., 2, e123

[35] Gillespie, D. T. (1976) A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. J. Comput. Phys., 22, 403–434

[36] Gillespie, D. T. (1977) Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem., 81, 2340–2361

[37] Gillespie Markov Processes, D. T.An Introduction for Physical Scientists, Academic Press, 1991

[38] Gillespie, D. T. (2007) Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem., 58, 35–55

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg