A common approach to the modeling of gene regulatory networks is to represent activating or repressing interactions using ordinary differential equations for target gene concentrations that include Hill function dependences on regulator gene concentrations. An alternative formulation represents the same interactions using Boolean logic with time delays associated with each network link. We consider the attractors that emerge from the two types of models in the case of a simple but nontrivial network: a figure-8 network with one positive and one negative feedback loop. We show that the different modeling approaches give rise to the same qualitative set of attractors with the exception of a possible fixed point in the ordinary differential equation model in which concentrations sit at intermediate values. The properties of the attractors are most easily understood from the Boolean perspective, suggesting that time-delay Boolean modeling is a useful tool for understanding the logic of regulatory networks.

The changing concentrations of mRNA and protein species in a living cell are determined by complex sequences of chemical interactions that cannot be modeled in all their biomolecular detail. In order to understand the fundamental function of a given network, one would like to construct a mechanistic model of the dynamics. Some form of direct simulation at the level of Gillespie algorithms for the production and degradation of key molecular components may be possible, but the basic logic of the network is best revealed if one can construct a deterministic model, supplemented by stochastic perturbations if necessary, that reproduces the observed behavior. Our aim here is to compare and contrast two rather different types of deterministic models: ordinary differential equation (ODE) networks based on continuous variables and autonomous Boolean networks (ABNs) based on binary variables and continuous time delays. The two types of model may be thought of as parsimonious choices within a large array of possibilities that includes time-delay models with multiple-state discrete variables and differential equation models with explicit time delays. For one particularly instructive network architecture, we find that when the models are constructed with the minimal features required to avoid nongeneric artifacts, the ODE and ABN attractors are qualitatively similar and that the ABN models are particularly amenable to analysis. In this case, the ABN analysis clearly distinguishes between two types of oscillators, frustration oscillators that rely on negative feedback to generate the oscillations and transmission oscillators in which a negative feedback loop serves only to regulate the duration of a pulse of activity that travels around a positive feedback loop.1 

We study the dynamics of the network shown in Figure 1, which consists of a single node A that participates in a positive feedback loop containing n nodes B1,,Bn and also in a negative feedback loop containing m nodes C1,Cm. Each variable may be thought of as representing the level of expression of a different mRNA species. Note that A has one activating input and one repressing input. We study here the case where activation of A occurs if and only if Bn is high and Cm is low; i.e., the case where the activator is required for expression, but the repressor is dominant.

FIG. 1.

The figure-8 network. Arrows indicate activating links. The bar indicates a repressing link.

FIG. 1.

The figure-8 network. Arrows indicate activating links. The bar indicates a repressing link.

Close modal

The figure-8 class of networks is chosen because of its simplicity and because it supports several distinct types of oscillations that exhibit interesting features. The class contains networks of the type suggested for modeling the transcriptional oscillator associated with the yeast cell cycle1 as well as the typical negative feedback oscillator, such as the repressilator.2 As we will see below, the negative feedback loop is necessary either as the primary source of oscillations or as a stabilizing factor when the oscillations are generated by the positive loop. The positive loop is not strictly necessary. It may be replaced by a single node Bn that is assumed to be on at all times. Parameters can be chosen for the figure-8 network, however, such that Bn does remain on at all times, which allows the attractor dynamics of the simple negative feedback oscillator to be viewed as a special case of attractor dynamics in a figure-8 model.

Our emphasis here is on qualitative features of the dynamics. More realistic models might include interpolating nodes representing protein levels and other post-transcriptional or post-translational steps, which can have the effect of reducing the values of n necessary for obtaining the attractors of interest in the ODE models. Our goal here is to clarify the connection between a certain class of ODEs that has been used in the literature2–4 and our ABN models, which represent a different branch of the literature.5–7 

ABN models, in general, are designed to implement directly the three core features of a regulatory system: (1) the architecture of causal links between regulatory elements; (2) the logic of the response of each element to its regulators; and (3) the times required for targets to respond to changes in the states of their regulators. The first two features are common to all Boolean models of regulatory networks. The third allows for dynamical behaviors that may be missed in Boolean models with synchronous or random asynchronous dynamics, and for the study of their stability against small timing fluctuations. We have included two features that are essential for avoiding artifacts that may arise in some types of time-delay Boolean models: a distinction between time delays for on and off signals removes marginally stable attractors that would not be observed in real systems, and the filtering of very short pulses avoids the unlimited buildup of sharp spikes of activity. The ABN model studied here is of the type described by Cheng et al. in Ref. 7. For present purposes, the simplicity of the network here (and lack of direct application to a particular biological system) allows a simpler presentation.

In an ABN, each expression level of the species Y is represented as a binary variable y{0,1}. A Boolean function FY determines the value of y based on the values of its inputs. For the figure-8 network, FBi and FCj are single-input Boolean “copy” functions, while FA is a two-input NIF function: FA(Bn,Cm)=1 if and only if Bn=1 and Cm=0.

For a link from input X to output Y, we refer to nodes X as the “source” and Y as the “target.” Associated with the link are two time delay parameters τXY,1,τXY,20, with τXY,1rX,1 and τXY,2rX,2, where the r's are positive constants representing minimum durations for generating a response in X. The link from X to Y represents a chain of events through which x may influence y. τXY,1 and τXY,2 may be thought of as the times required for signals to traverse the link. τXY,1 is the delay between the instant x updates from 0 to 1 and the instant this update affects y; τXY,2 is the corresponding delay when x updates from 1 to 0.

For node A, the time series a(t) is developed in chronological order as follows. At time t, we first compute

(1)

where tx is a time infinitesimally later than the latest origination time of the signals transmitted from node Bn or Cm that have arrived at A before time t. These can be written as tb=tτBnA,σb and tc=tτCmA,σc, where σb{1,2} and σc{1,2} are determined by whether the latest signal to arrive at A was a 0-to-1 or 1-to-0 update in each case.

Note that it is the origination time from the source node that enters the equation rather than the arrival time at node A. After a given source switches, it may switch back quickly enough that the trailing signal would reach the target before the leading one, in which case the first signal would never affect A. Physically this corresponds to a case in which the trailing signal annihilates the leading one before it reaches the target.

Next, we adjust the recent history of A to account for the fact that a target cannot respond to counteracting signals received in sufficiently rapid succession. Let s be the latest time before t that A switched from 0 to 1 (or 1 to 0). If a signal arrives at time t that causes A to switch back from 1 to 0 (or 0 to 1) and ts<rA,1 (or rA,2), then A(t) is set to 0 (or 1) for the entire interval s<t<t and we say that the short pulse (or dip) is rejected.1,7,8 Note that the restrictions τAY,1rA,1 and τAY,2rA,2 ensure that the adjustment to a(t) can be made in time to remove any signals generated by the switch at time s before they reach their targets.

The ABN dynamics can be simulated by an event-driven computer code. A time-ordered queue is maintained containing signals and arrival times at all targets, and the target receiving the next earliest signal is updated according to its logic function. Stochastic fluctuations may be modeled by adding a random increment to each event time as it is added to the queue.1,5

For present purposes, we neglect the dependence of the time delay on the recent history of the source and target; i.e., we do not attempt to model the fact that a node that has just recently turned off may turn on more quickly than one that has been off for a long time. In some physical systems, such as unclocked digital electronic circuits, history dependent effects can be important and the dynamics of the network may be chaotic.8–10 While it is possible to modify the ABN model to account for such effects, we do not consider them here. Thus the scope of the present paper is limited to periodic attractors.

As for any model with explicit time delays, the initial condition for a simulation must specify the state of the system over a time interval equal to the largest time delay. The present study, however, addresses only the attractor structure of the networks, so the details of the implementation of initial conditions are not relevant.

Regardless of the values of the time-delay parameters, our network has a trivial fixed point attractor in which all of the variables are equal to 0. In that state, all of the logic functions are satisfied. To characterize the other possibilities, we begin by outlining four distinct options for attractors in which A turns on and off exactly once per cycle. The options are depicted in the panels of Figure 2, each having a different causal structure. The figure shows time series for three quantities at node A: the output value a(t), and the input values coming from nodes Bn and Cm. These input values represent the signals reaching the target A at time t, which are equal to bn(tτBnA,1), and cm(tτCmA,1) for signals arriving from sources turning on, and bn(tτBnA,2), and cm(tτCmA,2) for signals arriving from sources turning off (see Eq. (1)).

FIG. 2.

Four types of oscillation in the ABN model. See text for explanations. Black dots indicate the events that cause the switch in A at a given time. (I) Frustration oscillation. (II) Marginally stable (unphysical) case. (III) Pulse transmission oscillation. (IV) Dip transmission oscillation.

FIG. 2.

Four types of oscillation in the ABN model. See text for explanations. Black dots indicate the events that cause the switch in A at a given time. (I) Frustration oscillation. (II) Marginally stable (unphysical) case. (III) Pulse transmission oscillation. (IV) Dip transmission oscillation.

Close modal

Let β1 be the sum of the time delays associated with propagating a switch from 0 to 1 one full time around the positive loop of the figure-8 network and let β2 be the sum of the time delays associated with propagating a switch from 1 to 0 one full time around the positive loop. Similarly, let γ1 and γ2 be the times required to propagate 0-to-1 and 1-to-0 updates around the negative loop.

We will assume throughout that the effect of a full loop can be summarized by the sums of the time delays associated with the separate links in the loop. We do not analyze here all of the possibilities that may arise if, for example, we have β2>β1, which would cause the width of a pulse to grow, but pulses are annihilated while traversing the loop because τXY,2<τXY,1 for some links in the loop.

We have the following four cases. We assume here that variations in the time delays for links in the same loop are similar enough that either τXY,1>τXY,2 for all links or vice versa, so that we do not have a situation where a pulse that would grow in width when traveling around the whole loop actually shrinks (and possibly disappears) at some intermediate stage.

I: Activation of A is caused by Cm turning off; deactivation of A is caused by Cm turning on. In this case, Bn is always on and system behaves like a simple negative feedback loop, which generates oscillations because the FA and all of the FCj's cannot be simultaneously satisfied for any state of the system. For this reason, we call this a frustration oscillation. This case requires

(2)

This is also a sufficient condition under the assumption above that time delay variations on the links in a loop are sufficiently weak. The period T of the oscillations is determined by the negative loop delay times: T=γ1+γ2.

II: Activation of A is caused by Bn turning on; deactivation of A is caused by Bn turning off. This case requires fine tuning to make

(3)

and hence is only marginally stable. Deviation from this equality will cause the pulse width to grow or shrink from one cycle to the next and eventually lead either to the all off fixed point or a different type of oscillation.

III: Activation of A is caused by Bn turning on; deactivation of A is caused by Cm turning on. For the case shown, the condition for the existence of this attractor is

(4)

Here, we have a pulse whose width is determined by the negative loop delay time γ1. The period of the oscillation is T=β1, the time required for the leading edge of the pulse to propagation around the positive loop. We call this a pulse transmission oscillation.

IV: Activation of A is caused by Cm turning off; deactivation of A is caused by Bn turning off. For the case shown, the condition for the existence of this attractor is

(5)

Here, we have a dip whose width is determined by the negative loop delay time γ2. The period of the oscillation is T=β2, the time required for the leading edge of the dip to propagation around the positive loop. We call this a dip transmission oscillation.

We note that the conditions expressed in Eqs. (5) and (4) are mutually exclusive. For rA,2=0, Eq. (2) is also incompatible with either Eq. (4) or Eq. (5). Thus, for sufficiently small rejection time parameters, a given network cannot simultaneously support two of these types of stable oscillations. In cases where rA,2 is relatively large, however, we might expect to see coexistence of attractors corresponding to frustration oscillation and pulse or dip transmission oscillation.

We note also that each of the three stable attractors requires β1<β2. That is, for the given choice of FA, oscillations can be maintained only if the width of a pulse grows as it traverses the positive loop.

Completeness demands that we consider several other possible attractor structures. Consider the pulse transmission case (III) for example, which has period T=β1. The same pattern could be obtained by adding T to γ1 and γ2. (We cannot add only to γ1 because causality requires that the switch in Cm induced by turning A on must precede the switch induced by subsequently turning A off.) In this scenario, the basic logic of the oscillations remains the same: a growing pulse on the positive loop is cut back to its original size by a signal that has traversed the negative loop. We might expect the negative loop to be short enough so that the cutoff signal is the one generated by the leading edge of the pulse in question, but it is possible for it to have come from the previous pulse or even an earlier one.

Next consider a variation of case III in which β1 and β2 are extremely long compared to γ1+γ2. It would then be possible to have multiple pulses circulating on the positive loop periodically and, within our ABN formalism, any spacing between the pulses would be possible, resulting in a family of marginally stable attractors with multiple pulses in A per period. In some continuous models, these pulses might repel each other so that the system eventually settles into a stable pattern of pulses with a preferred spacing between them. An example can be found in the delay-differential equation model of Norrell et al.11 Similar effect may be expected for the case of dip transmission oscillations (IV).

It is interesting to note that our classification of attractors according to their causal structure is complementary to the identification of equivalence classes with different state transition diagrams. (See, for example, Glass.12) The latter approach focuses on the sequence in which different nodes switch states. What matters for our classification of the oscillation type, on the other hand, is only the time that signals arrive back at node A. By independently adjusting the delays on different links in an ABN, one can arrange to have identical symbolic sequences that correspond to different types of oscillations. That is, in our figure-8 example, the effects on A of its two inputs need not occur in the same order that the input nodes are updated. The state of the time delay system cannot be fully specified by a list of the node values at one instant, and therefore the symbolic dynamics cannot contain the information about which node is the cause of a given change at A.

We now turn to the question of whether the different types of oscillation identified by the ABN analysis account qualitatively for the behavior observed in ODE models. The class of ODE models that we might want to consider is not well-defined, so we do not have rigorous theorems showing that all ODE attractors must be of one of the types identified above. We can, however, ask whether a generic (and popular) class of ODEs that is often used to model biological networks does in fact exhibit each of the three types of oscillation — frustration, pulse transmission, and dip transmission — in some parameter regimes. We show here that the answer is “Yes” by exhibiting examples of each and explaining the relationship between the ODE and ABN parameters.

We consider a class of models employing Hill functions of the type typically used in models of gene regulation,2–4 which are known to be capable of supporting nontrivial oscillatory attractors.13 Our choice of Hill functions is motivated primarily by the fact that they have been used in the literature on gene regulatory networks. We take the degradation terms to be linear for convenience. The effects illustrated here are not sensitive to the values we choose for the cooperativity exponents, though smaller exponents typically require longer loops to sustain oscillations. With no cooperativity (i.e., the case of Michaelis-Menten dynamics, where all exponents are equal to unity), we find no oscillations. (Oscillations based on Michaelis-Menten dynamics are known to be possible when degradation terms are nonlinear.14 For a review of mechanisms for generating oscillations in biochemical systems, see Ref. 15.)

Let fa and fr be activating and repressing Hill functions

(6)

where ν is a cooperativity parameter that we take to be the same for all nodes for simplicity. The equations governing the concentrations are

(7)
(8)
(9)
(10)
(11)

where each Kx is a constant parameter between 0 and 1. For the numerical studies below, we assume for simplicity that all KBi for i<n are equal to KAB, all of which control the activations of B nodes, and that KBn, which controls the activation of A may be different. Similarly, we take KCj for j<m to be equal to KAC, while KCm may be different.

The attractors shown below are obtained by setting the parameters as indicated and integrating the equations numerically (using the ndsolve function of mathematica): Figure 3 shows a frustration oscillation; Figure 4 shows a pulse transmission oscillation; Figure 5 shows a dip transmission oscillation.

FIG. 3.

An attractor showing frustration oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.4,KCm=0.5,KAB=0.3,KAC=0.35. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the pulse that results in B9 remaining on at all times.

FIG. 3.

An attractor showing frustration oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.4,KCm=0.5,KAB=0.3,KAC=0.35. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the pulse that results in B9 remaining on at all times.

Close modal
FIG. 4.

An attractor showing pulse transmission oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.45,KCm=0.5,KAB=0.43,KAC=0.38. One initial condition in the basin of attraction has Bi=1 for i=5,6,7,8,9 and all other variables set to 0. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. In this case, A is turned on by a rise in Bn and turned off by a rise in Cm. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the pulse.

FIG. 4.

An attractor showing pulse transmission oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.45,KCm=0.5,KAB=0.43,KAC=0.38. One initial condition in the basin of attraction has Bi=1 for i=5,6,7,8,9 and all other variables set to 0. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. In this case, A is turned on by a rise in Bn and turned off by a rise in Cm. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the pulse.

Close modal
FIG. 5.

An attractor showing dip transmission oscillations. The parameter values are n=9,m=6,ν=5,KBn=0.3,KCm=0.7,KAB=0.44,KAC=0.45. Top: The thick line shows A; the dashed line Bn; and the thin line Cm. The thick line shows A; the thin line Bn, and the dashed line Cm. In this case, A is turned on by a drop in Cm and turned off by a drop in Bn. Bottom: The time series for each of the nodes Bi on the positive loop showing that the pulse width (rather than the dip width) still grows.

FIG. 5.

An attractor showing dip transmission oscillations. The parameter values are n=9,m=6,ν=5,KBn=0.3,KCm=0.7,KAB=0.44,KAC=0.45. Top: The thick line shows A; the dashed line Bn; and the thin line Cm. The thick line shows A; the thin line Bn, and the dashed line Cm. In this case, A is turned on by a drop in Cm and turned off by a drop in Bn. Bottom: The time series for each of the nodes Bi on the positive loop showing that the pulse width (rather than the dip width) still grows.

Close modal

To understand the behavior of the ODE system using the insights gained from the ABN models, we need to extract the relevant ABN time delay parameters from the ODEs. To make a plausible correspondence between the two, we use the following procedure to determine τXY,1. If X is the only input to Y, let W be the activator of X. Initialize the ODE system with W=1,X=0, and Y = 0, then, holding W constant, solve for X(t) and Y(t). Let t1 be the time at which X=K1, where K1 is the threshold appearing in the Hill function governing Y, and let t2 be the time at which Y=K2, the threshold appearing in the Hill function governing the target of Y. Then our estimate for the activation delay τXY,1 is t2t1. To determine the de-activation delay τXY,2, we use the same procedure but initialize with W=0,X=1,Y=1.

For the two-input node A, we use the natural generalization of the single input procedure. To calculate τBnA,1 or τBnA,2, we hold Cm=0 and treat the Bn-to- A link exactly as for a single input case. To calculate τCmA,1, which is to the delay between Cm turning on and its repressive effect on A, we hold Bn=1, set Cm1=1,Cm=0, and A = 1, then solve for the threshold crossing times for Cm rising and A falling. To calculate τCmA,2, which is to the delay between Cm turning off and A no longer being repressed, we hold Bn=1, set Cm1=0,Cm=1, and A = 0, then solve for the threshold crossing times for Cm rising and A falling.

Note that there is one element of ambiguity associated with the fact that A has two outputs. In principle, the two activating functions could have substantially different thresholds and it would not be clear which one should count for defining A as “on.” One natural choice is to use the higher of the two thresholds for defining when A turns on, but the lower of the two for defining when A turns off; i.e., we call A “on” only if it is high enough to activate both of its targets, and call it “off” only when it drops low enough to de-activate both of them.

Using these rules, we can estimate the values of the ABN parameters β1,β2,γ1, and γ2 that correspond to a given set of ODE parameters. For example, β1 is the sum of the activation delays computed for each of the links with targets Bn, plus the activation delay for the link from Bn to A. Table I shows the results. The ODE parameters are given in the figure captions.

Table I.

ABN parameters extracted from ODEs. The period listed is the actual period measured from numerical integration of the ODEs. The error indicated is the percent difference between the measured period and the period expected from the ABN analysis.

Oscillation typeβ1β2γ1γ2PeriodError (%)
Frustration (Fig. 34.20 12.40 2.74 2.18 3.62 −16 
Pulse transmission (Fig. 46.73 9.01 2.44 2.41 6.62 −1.6 
Dip transmission (Fig. 56.60 9.21 6.30 5.31 8.95 −2.9 
Oscillation typeβ1β2γ1γ2PeriodError (%)
Frustration (Fig. 34.20 12.40 2.74 2.18 3.62 −16 
Pulse transmission (Fig. 46.73 9.01 2.44 2.41 6.62 −1.6 
Dip transmission (Fig. 56.60 9.21 6.30 5.31 8.95 −2.9 

Two points are immediately apparent. First, the inequalities relevant to the different cases (Eqs. (2), (4), and (5)) are respected in every case, even if we neglect rA,2. Because these are mutually exclusive conditions, we have a strong indication that none of the ODE parameter sets shown here could support one of the other types of oscillations.

Second, the periods of both transmission oscillations are accurately accounted for by the ABN analysis, but the period of the frustration oscillation is not. The accuracy of our method for estimating absolute time delays (and hence the oscillation period) depends on the variables ranging well above and below their threshold values. Visual inspection of the figures shows that this is a reasonable assumption for the transmission oscillators, but not for the frustration oscillation. Lengthening the negative feedback loop by including more nodes (and lengthening the positive loop to keep β1+γ2<β2) increases the accuracy of the correspondence between the ABN description and the ODE behavior. Setting m = 6 and n = 19 for the frustration oscillator, while leaving all other parameters the same as for Figure 3, produces well developed oscillations and an accuracy of 3% in the ABN estimate of the period.

For cases in which the frustration oscillation has relatively low amplitude, a better estimate of the period may be obtained from analysis of the Hopf bifurcation that gives rise to it. In the present case, there is a Hopf bifurcation occurring on a curve in the (KAC,ν) parameter plane. Taking KAC=0.5 and analyzing the bifurcation at ν=4, one finds a period of 2π/3, which is quite close to the observed period.12 Identification of the relevant Hopf bifurcations in the full figure-8 system with generic parameter values is more difficult. Our point in the present context is that the ABN analysis, which can be applied to more complex networks, works reasonably well.

The possibility of multiple attractors in the ABN due to nonzero short pulse rejection times is reflected in the ODE system, as shown in Figure 6. While a rigorous study of the basins of these attractors is difficult due to the high dimensionality of the state space, we have made a rough survey by considering all possible initial conditions in which each node begins at either 0 or 1. Of the 212 choices of initial conditions, we find that 2525 lead to the frustration oscillation, 745 lead to the pulse transmission oscillation, and the remaining 826 lead to the trivial all-0 attractor. These numbers remained identical when the initial values were 0.2 or 0.8 rather than 0 or 1. This is a case where extracting time delays for the ABN do not yield sufficiently accurate values to predict the stability of pulse transmission oscillations because the pulse does not reach high values sufficiently close to saturation. Nevertheless, the behavior is consistent with the general ABN analysis that predicts the possibility of multistability due to short pulse rejection.

FIG. 6.

An attractor showing frustration oscillations (top) and pulse transmission oscillations (bottom) for the same parameter values. The parameter values are n=9,m=2,ν=5,KBn=0.45,KCm=0.5,KAB=0.42,KAC=0.37. The initial condition for the frustration oscillations was A=1,Bi=1 for all i, and C1=C2=0. The initial condition for the pulse transmission oscillations was A=B1=1,Bi=0 for all i2, and C1=C2=0. The thick line shows A; the thin lines show the Bi's.

FIG. 6.

An attractor showing frustration oscillations (top) and pulse transmission oscillations (bottom) for the same parameter values. The parameter values are n=9,m=2,ν=5,KBn=0.45,KCm=0.5,KAB=0.42,KAC=0.37. The initial condition for the frustration oscillations was A=1,Bi=1 for all i, and C1=C2=0. The initial condition for the pulse transmission oscillations was A=B1=1,Bi=0 for all i2, and C1=C2=0. The thick line shows A; the thin lines show the Bi's.

Close modal

It is instructive to study the regions in a 2-dimensional slice of parameter space where different types of oscillations may occur. Here, we consider the figure-8 network with m=6,n=9,ν=5,KBn=0.45, and KCm=0.5, as in Figure 5. For a range of values of KAB and KAC, we run the ODE system to determine whether it supports dip transmission or frustration oscillations and plot the results in Figure 7. The light gray region supports stable dip transmission oscillations (type IV in the classification above). The darker gray region supports frustration oscillations. The transition between the two is a smooth one. As the boundary between them is crossed the attractor itself does not undergo a discrete shift, but the relative timing of the signals from the B and C loops that turn A on changes sign. The overlap of the two regions (darkest gray) corresponds to oscillations where one cannot clearly define which signal comes first. The figure was made by running the system at an array of different parameter values for a single choice of the initial condition: Bi=1 for i1ii2, and all other variables set to 0. The boundary between the two types of oscillation was found to be nearly identical for all choices of i1 and i2 satisfying i2i1>2.

FIG. 7.

A phase diagram for the ODE system. The parameter values are n=9,m=6,ν=5,KBn=0.45,KCm=0.5. Shaded regions indicate different types of oscillation in the ODE. The thick dashed line outlines the region where time delay values extracted from the ODE system satisfy the condition for frustration oscillations (type I). The thick solid line outlines the region where conditions for dip transmission oscillations (type IV) are satisfied. The white regions in the figure indicate that the system always goes to the all off state. The black disk marks a small region where period-2 oscillations were observed.

FIG. 7.

A phase diagram for the ODE system. The parameter values are n=9,m=6,ν=5,KBn=0.45,KCm=0.5. Shaded regions indicate different types of oscillation in the ODE. The thick dashed line outlines the region where time delay values extracted from the ODE system satisfy the condition for frustration oscillations (type I). The thick solid line outlines the region where conditions for dip transmission oscillations (type IV) are satisfied. The white regions in the figure indicate that the system always goes to the all off state. The black disk marks a small region where period-2 oscillations were observed.

Close modal

The thick lines in Figure 7 show boundaries obtained from the conditions of Eqs. (2) and (5) using β and γ values determined from the ODEs as described above. The value of rA,2 is taken to be 2.0. This is a rough estimate determined empirically by examining the width of a dip in A that just reaches the threshold value for parameters in the middle of the plot. Note that the region of overlap of the two types of oscillation corresponds roughly to the ambiguous region in the ODE system. We also found a small region, marked by the black disk, where oscillations occur in which the peaks in A alternate in shape. Analysis of these period-2 (or higher) oscillations is beyond our present scope.

In general, there is good agreement between the ODE behaviors and ABN expectations. Minor shifts discrepancies in the boundaries of the regions are to be expected due to short pulse rejection effects or small discrepancies between time delay parameters extracted as described above and the actual delays arising in the ODE dynamics when variables do not fully saturate their high or low values. The failure to observe dip transmission oscillations in the upper left corner of the figure is also due to variations in the time delays associated with pulses of different amplitude. The ABN analysis does not account for the fact that β2 depends on the height of the pulse in A. Because the putative attractor in this region has a narrow pulse, the amplitude of that pulse is low enough that β2 becomes smaller than β1 and the oscillations collapse to the all off state. Discrepancies of this type generally become more prevalent for shorter loop lengths, where the attractor periods are short enough that nodes do not have time to rise to their saturation values.

Previous studies have emphasized the need for long time delays in regulatory oscillators. In the Elowitz-Leibler model of the repressilator (which is a frustration oscillator), protein creation and degradation equations were added to the system in order to capture the oscillatory dynamics.2 From our present perspective, the protein dynamics simply serves to lengthen the delay time for propagation of a pulse around the loop enough to allow elements to vary with sufficient amplitude. The explicit representation of protein variables is not necessary if the loop is made longer.

Norrell et al. studied a different mechanism for lengthening the loop propagation times: inserting explicit delays into the differential equations.11 Using a slightly different form for fA and fR, they studied frustration oscillations and pulse transmission oscillations, but did not address the distinct possibility of dip transmission oscillations.

Finally, it is worth emphasizing that the distinction between pulse transmission and dip transmission is not simply a matter of symmetry; that is, the dip transmission oscillations are not just pulse transmission oscillations with the on and off states exchanged. If that were the case, we would have a dip that grows in width as it traverses the positive loop, but Figure 5 clearly shows that it is pulses (not dips) that grow in the dip transmission oscillator. The on-off symmetry is broken by the Hill function forms for fA and fR, but this is merely a quantitative effect that determines the parameter domains where oscillation is possible. The more important symmetry breaking in the figure-8 system is the logic function for the two-input element A. If the default state (with both inputs off) were taken to yield A = 1 and the activating input were dominant, we could obtain oscillations in cases where dips grow rather than pulses. The language becomes a bit cumbersome: it might be best to refer to these cases as “anti-pulse transmission” and “anti-dip transmission” oscillations. Figure 8 shows an anti-pulse transmission oscillator, where the ODE system is the same as above except that Eq. (7) is replaced by

(12)

and parameter values are given in Figure 8.

FIG. 8.

An attractor showing anti-pulse transmission oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.55,KCm=0.5,KAB=0.52,KAC=0.55. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. Note that the width of the dip for Bn is larger than for A. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the dip (the anti-pulse).

FIG. 8.

An attractor showing anti-pulse transmission oscillations. The parameter values are n=9,m=2,ν=5,KBn=0.55,KCm=0.5,KAB=0.52,KAC=0.55. Top: The thick line shows A; the thin line Bn; and the dashed line Cm. Note that the width of the dip for Bn is larger than for A. Bottom: The time series for each of the nodes Bi on the positive loop showing the growth of the dip (the anti-pulse).

Close modal

This study serves to illustrate a sense in which ABN modeling can be used to identify distinct classes of oscillatory solutions of ODE systems of a type often used to model activating and repressing regulatory interactions. The time series of a selected node on a periodic attractor is viewed as two interleaved doubly infinite sequences of up-regulation events, Ui, and down-regulation events, Di. The causes of Un and Dn are traced back to previous rises or falls in the pulse train. In Figure 2, take U2 to be the second (later time) rise in the time series of A, and U1 to be the first rise, and similarly for D2 and D1. There are four possibilities:

  • Un is caused by Dn1; Dn is caused by Un. In this case, both events are triggered by an inverted event, indicating that both are the result of propagation around a negative feedback loop.

  • Un is caused by Un1; Dn is caused by Dn1. This is the case where the negative feedback loop plays no role and fine tuning is required.

  • Un is caused by the Un1; Dn is caused by Un.

  • Un is caused by the Dn1; Dn is caused by Dn1.

From the ABN analysis, it is clear that stable oscillations require some causative role of the negative feedback loop, consistent with a conjecture by Thomas that oscillation requires negative feedback.16,17

The above classification could clearly be extended to cases where the cause of a rise or fall goes back to earlier events. For example, a variation on III could have Un caused by Un2 and Dn caused by Un. These cases correspond to having more than one pulse propagating on at least one of the loops at any given time. Such complex oscillations appear to be unlikely candidates for biological implementation, but may present interesting possibilities for electronic and opto-electronic systems.18 

Because negative feedback seems to be necessary for stability while positive feedback does not, it is natural to suppose that oscillators found in nature will be of the frustration type. It may be important, however, to be aware of other possibilities, and there is some evidence that a transcriptional oscillator underlying the cell cycle of yeast is a pulse transmission oscillator.1 To our knowledge, no examples of dip transmission oscillations have been reported, but it is not clear whether there is a fundamental reason for this, be it evolutionary or functional.

This work was supported by the National Science Foundation Grant No. DMS-10-68602 and the National Institutes of Health Grant No. P50-GM081883.

1.
V.
Sevim
,
X.
Gong
, and
J. E. S.
Socolar
, “
Reliability of transcriptional cycles and the yeast cell-cycle oscillator
,”
PLOS Comput. Biol.
6
,
e1000842
(
2010
).
2.
M. B.
Elowitz
and
S.
Leibler
, “
A synthetic oscillatory network of transcriptional regulators
,”
Nature
403
,
335
338
(
2000
).
3.
G.
von Dassow
,
E.
Meir
,
E. M.
Munro
, and
G. M.
Odell
, “
The segment polarity network is a robust developmental module
,”
Nature
406
,
188
192
(
2000
).
4.
N. T.
Ingolia
, “
Topology and robustness in the Drosophila segment polarity network
,”
PLOS Biol.
2
,
e123
(
2004
).
5.
R.
Thomas
and
R.
D’Ari
,
Biological Feedback
, 1st ed. (
CRC Press
,
Boca Raton, FL
,
1990
).
6.
M.
Ghil
,
I.
Zaliapin
, and
B.
Coluzzi
, “
Boolean delay equations: A simple way of looking at complex systems
,”
Physica D
237
,
2967
2986
(
2008
).
7.
X.
Cheng
,
M.
Sun
, and
J. E. S.
Socolar
, “
Autonomous Boolean modelling of developmental gene regulatory networks Autonomous Boolean modelling of developmental gene regulatory networks
,”
J. R. Soc., Interface
10
,
20120574
(
2013
).
8.
R.
Zhang
,
H. de S.
Cavalcante
,
Z.
Gao
,
D.
Gauthier
,
J.
Socolar
,
M.
Adams
, and
D.
Lathrop
, “
Boolean chaos
,”
Phys. Rev. E
80
,
045202
(
2009
).
9.
L.
Glass
,
T. J.
Perkins
,
J.
Mason
,
H. T.
Siegelmann
, and
R.
Edwards
, “
Chaotic dynamics in an electronic model of a genetic network
,”
J. Stat. Phys.
121
,
969
994
(
2005
).
10.
H. L. D. D. S.
Cavalcante
,
D. J.
Gauthier
,
J. E. S.
Socolar
, and
R.
Zhang
, “
On the origin of chaos in autonomous Boolean networks
,”
Philos. Trans. R. Soc. A
368
,
495
513
(
2010
).
11.
J.
Norrell
,
B.
Samuelsson
, and
J.
Socolar
, “
Attractors in continuous and Boolean networks
,”
Phys. Rev. E
76
,
046122
(
2007
).
12.
L.
Glass
, “
Combinatorial and topological methods in nonlinear chemical kinetics
,”
J. Chem. Phys.
63
,
1325
(
1975
).
13.
L.
Glass
and
J.
Pasternack
, “
Prediction of limit cycles in mathematical models of biological oscillations
,”
Bull. Math. Biol.
40
,
27
44
(
1978
).
14.
R. D.
Bliss
,
P. R.
Painter
, and
A. G.
Marr
, “
Role of feedback inhibition in stabilizing the classical operon
,”
J. Theor. Biol.
97
,
177
193
(
1982
).
15.
J.
Tyson
, “
Computational cell biology
” in
Biological Oscillations
, edited by
C. P.
Fall
,
E. S.
Marland
,
J. M.
Wagner
, and
J. J.
Tyson
(
Springer
,
New York
,
2002
), Chap. 9.
16.
R.
Thomas
, “
On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations
,” in
Numerical Methods in the Study of Critical Phenomena
, Vol. 9, Springer Series in Synergetics, edited by
J.
Dora
,
J.
Demongeot
, and
B.
Lacolle
(
Springer
,
Berlin, Heidelberg
), pp.
180
193
.
17.
M.
Kaufman
,
C.
Soulé
, and
R.
Thomas
, “
A new necessary condition on interaction graphs for multistationarity
,”
J. Theor. Biol.
248
,
675
85
(
2007
).
18.
I.
Kanter
,
Y.
Aviad
,
I.
Reidler
,
E.
Cohen
, and
M.
Rosenbluh
, “
An optical ultrafast random bit generator
,”
Nature Photon.
4
,
58
61
(
2010
).