Causal structure of oscillations in gene regulatory networks: Boolean analysis of ordinary differential equation attractors

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 ﬁgure-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 ﬁxed 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. V C 2013 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/1.4807733]

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

I. THE FIGURE-8 NETWORK
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 B 1 ; …; B n and also in a negative feedback loop containing m nodes C 1 ; …C m . 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 B n is high and C m is low; i.e., the case where the activator is required for expression, but the repressor is dominant.
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 cycle 1 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 B n that is assumed to be on at all times. Parameters can be chosen for the figure-8 network, however, such that B n 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 a) Author to whom correspondence should be addressed. Electronic mail: socolar@phy.duke.edu. URL: http://fds.duke.edu/db/aas/Physics/faculty/ socolar.
1054-1500/2013/23(2)/025104/9 V C Author(s) 2013 23, 025104-1 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 literature 2-4 and our ABN models, which represent a different branch of the literature. [5][6][7] A. ABN model 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 2 f0; 1g. A Boolean function F Y determines the value of y based on the values of its inputs. For the figure-8 network, F B i and F C j are single-input Boolean "copy" functions, while F A is a two-input NIF function: F A ðB n ; C m Þ ¼ 1 if and only if B n ¼ 1 and C m ¼ 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 s XY;1 ; s XY;2 ! 0, with s XY;1 ! r X;1 and s XY;2 ! r X;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. s XY;1 and s XY;2 may be thought of as the times required for signals to traverse the link. s XY;1 is the delay between the instant x updates from 0 to 1 and the instant this update affects y; s 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 where t x is a time infinitesimally later than the latest origination time of the signals transmitted from node B n or C m that have arrived at A before time t. These can be written as t b ¼ t À s B n A;r b and t c ¼ t À s C m A;r c , where r b 2 f1; 2g and r c 2 f1; 2g 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 t À s < r A;1 (or r A;2 ), then Aðt 0 Þ is set to 0 (or 1) for the entire interval s < t 0 < t and we say that the short pulse (or dip) is rejected. 1,7,8 Note that the restrictions s AY;1 ! r A;1 and s AY;2 ! r A;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 eventdriven 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][9][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.

B. Attractors of the ABN model
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 B n and C m . These input values represent the signals reaching the target A at time t, which are equal to b n ðt À s B n A;1 Þ, and c m ðt À s C m A;1 Þ for signals arriving from sources turning on, and b n ðt À s B n A;2 Þ, and c m ðt À s C m A;2 Þ for signals arriving from sources turning off (see Eq. (1)).
Let b 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 b 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 c 1 and c 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 b 2 > b 1 , which would cause the width of a pulse to grow, but pulses are annihilated while traversing the loop because s XY;2 < s 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 s XY;1 > s 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 C m turning off; deactivation of A is caused by C m turning on. In this case, B n is always on and system behaves like a simple negative feedback loop, which generates oscillations because the F A and all of the F C j 's cannot be simultaneously satisfied for any state of the system. For this reason, we call this a frustration oscillation. This case requires 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: II: Activation of A is caused by B n turning on; deactivation of A is caused by B n turning off. This case requires fine tuning to make 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 B n turning on; deactivation of A is caused by C m turning on. For the case shown, the condition for the existence of this attractor is Here, we have a pulse whose width is determined by the negative loop delay time c 1 . The period of the oscillation is T ¼ b 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 C m turning off; deactivation of A is caused by B n turning off. For the case shown, the condition for the existence of this attractor is Here, we have a dip whose width is determined by the negative loop delay time c 2 . The period of the oscillation is T ¼ b 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. (4) and (5) are mutually exclusive. For r A;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 r A;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 b 1 < b 2 . That is, for the given choice of F A , 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 ¼ b 1 . The same pattern could be obtained by adding T to c 1 and c 2 . (We cannot add only to c 1 because causality requires that the switch in C m 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 b 1 and b 2 are extremely long compared to c 1 þ c 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.

C. ODE model
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][3][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 f a and f r be activating and repressing Hill functions f a ðX; KÞ X K þ X and f r ðX; KÞ where is a cooperativity parameter that we take to be the same for all nodes for simplicity. The equations governing the concentrations are _ _ C jþ1 ¼ f a ðC j ; K Cj Þ À C jþ1 for 1 j m À 1 ; where each K x is a constant parameter between 0 and 1. For the numerical studies below, we assume for simplicity that all K Bi for i < n are equal to K AB , all of which control the activations of B nodes, and that K Bn , which controls the activation of A may be different. Similarly, we take K Cj for j < m to be equal to K AC , while K Cm 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.
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 s 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 t 1 be the time at which X ¼ K 1 , where K 1 is the threshold appearing in the Hill function governing Y, and let t 2 be the time at which Y ¼ K 2 , the threshold appearing in the Hill function governing the target of Y. Then our estimate for the activation delay s XY;1 is t 2 À t 1 . To determine the de-activation delay s XY;2 , we use the same procedure but initialize with W ¼ 0; For the two-input node A, we use the natural generalization of the single input procedure. To calculate s B n A;1 or s B n A;2 , we hold C m ¼ 0 and treat the B n -to-A link exactly as for a single input case. To calculate s C m A;1 , which is to the delay between C m turning on and its repressive effect on A, we hold B n ¼ 1, set C mÀ1 ¼ 1; C m ¼ 0, and A ¼ 1, then solve for the threshold crossing times for C m rising and A falling. To calculate s C m A;2 , which is to the delay between C m turning off and A no longer being repressed, we hold B n ¼ 1, set C mÀ1 ¼ 0; C m ¼ 1, and A ¼ 0, then solve for the threshold crossing times for C m 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 b 1 ; b 2 ; c 1 , and c 2 that correspond to a given set of ODE parameters. For example, b 1 is the sum of the activation delays computed for each of the links with targets B n , plus the activation delay for the link from B n to A. Table I shows the results. The ODE parameters are given in the figure captions.
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 r A;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 b 1 þ c 2 < b 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 ðK AC ; Þ parameter plane. Taking K AC ¼ 0:5 and analyzing the bifurcation at ¼ 4, one finds a period of 2p= ffiffi ffi 3 p , 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 2 12 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. 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; K Bn ¼ 0:45, and K Cm ¼ 0:5, as in Figure 5. For a range of values of K AB and K AC , 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: B i ¼ 1 for i 1 i i 2 , 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 i 1 and i 2 satisfying The thick lines in Figure 7 show boundaries obtained from the conditions of Eqs. (2) and (5) using b and c values determined from the ODEs as described above. The value of r A;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 b 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 b 2 becomes smaller than b 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 f A and f R , 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 f A and f R , 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

_
A ¼ ð1 À f r ðB n ; K Bn Þf a ðC m ; K Cm ÞÞ À A; (12) and parameter values are given in Figure 8.

II. CONCLUSIONS
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, U i , and down-regulation events, D i . The causes of U n and D n are traced back to previous rises or falls in the pulse train. In Figure 2, take U 2 to be the second (later time) rise in the time series of A, and U 1 to be the first rise, and similarly for D 2 and D 1 . There are four possibilities: I: U n is caused by D nÀ1 ; D n is caused by U n . In this case, both events are triggered by an inverted event, indicating that both are the result of propagation around a negative feedback loop. II: U n is caused by U nÀ1 ; D n is caused by D nÀ1 . This is the case where the negative feedback loop plays no role and fine tuning is required. III: U n is caused by the U nÀ1 ; D n is caused by U n . IV: U n is caused by the D nÀ1 ; D n is caused by D nÀ1 .
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 U n caused by U nÀ2 and D n caused by U n . 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.