New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Oscillation patterns in negative feedback loops

Edited by Leo P. Kadanoff, University of Chicago, Chicago, IL, and approved February 6, 2007 (received for review December 5, 2006)
Abstract
Organisms are equipped with regulatory systems that display a variety of dynamical behavior ranging from simple stable steady states, to switching and multistability, to oscillations. Earlier work has shown that oscillations in protein concentrations or gene expression levels are related to the presence of at least one negative feedback loop in the regulatory network. Here, we study the dynamics of a very general class of negative feedback loops. Our main result is that, when a single negative feedback loop dominates the dynamical behavior, the sequence of maxima and minima of the concentrations exhibit a pattern that uniquely identifies the interactions of the loop. This allows us to devise an algorithm to (i) test whether observed oscillating time series are consistent with a single underlying negative feedback loop, and if so, (ii) reconstruct the precise structure of the loop, i.e., the activating/repressing nature of each interaction. This method applies even when some variables are missing from the data set, or if the time series shows transients, like damped oscillations. We illustrate the relevance and the limits of validity of our method with three examples: p53Mdm2 oscillations, circadian gene expression in cyanobacteria, and cyclic binding of cofactors at the estrogensensitive pS2 promoter.
Physiological processes in living cells exhibit a wide range of dynamical behavior, ranging from stable steady states [e.g., iron regulation (1) and response to DNA damage (2) in bacteria], multistability [e.g., lysis–lysogeny decision in temperate phage (3)], to oscillations. The most obvious examples of the latter are cell division and circadian (24 h) rhythms. Cellular processes are often coupled to the circadian clock, e.g., respiration and carbohydrate synthesis in cyanobacteria (4), which makes them periodic. Recently, faster “ultradian” oscillations with time periods of the order of hours have been found in several systems of interacting proteins, which influence the immune system (NFκB; refs. 5 and 6), apoptosis (p53; ref. 7), and development (Hes; ref. 8).
Theoretical studies of these oscillatory systems (9–12) describe the dynamics of the relevant variables (usually protein concentrations or gene expression levels) using differential equations. Very often, these models share two properties: the first is that the interactions between variables are monotone: proteins that activate a particular process will not change to repress that process at a different concentration, and vice versa. This is related to the properties of the biochemical processes acting on these slow timescales, most notably (but not only) transcription regulation. The second feature is that the system contains at least one negative feedback loop (i.e., a loop with an odd number of repressors). Indeed, a conjecture by Thomas (13), rigorously proven in refs. 14 and 15, states that, in a monotone system, at least one positive feedback loop is needed to have multistability (i.e., existence of multiple steady states), and at least one negative feedback loop is needed to have periodic behavior.
Feedback loops may thus be seen as the building blocks of the nontrivial dynamical behaviors of these systems; a network without loops can only reach a unique fixed point, regardless of the initial conditions. However, in general, oscillations in a particular system could depend on several factors, e.g., multiple feedback loops, time delays (11, 12), noise (16), or spatial effects (17). Oscillating time series observed in experiments may not have the precision or length required to discriminate between different models. For example, the p53Mdmd2 oscillations observed in singlecell fluorescence experiments, shown in Fig. 1 a, are consistent with several different models (7, 11). Nevertheless, we are going to show that oscillations like those of Fig. 1 a contain precious information about the real system.
Here, we propose a simple algorithm that can be used to deduce whether an oscillating time series is consistent with a single underlying negative feedback loop with no crosslinks, and if it is, it can also deduce whether each interaction is activating or repressing. This algorithm, described below, is based on the fact that the time order of the maxima and minima of the concentrations has a periodic pattern that uniquely identifies the interactions between variables, provided that the oscillations are produced by just one underlying negative feedback loop. The mathematical basis for the algorithm is laid out in the subsequent two sections. In the final section, we apply our algorithm to reconstruct the feedback loop from oscillating time series of two more biological systems; the examples also clarify the scope and limitations of the algorithm.
Extracting the Feedback Loop
As mentioned above, the time series in Fig. 1 a shows a specific order of maxima and minima of the concentrations. We can investigate this pattern by dividing the data set into intervals whose ends are determined by the occurrence of an extremal value (maximum or minimum) of a variable. In all of these intervals, marked by vertical dotted lines in Fig. 1 a, each variable shows an unchanging trend, either increasing or decreasing with time. We can therefore uniquely associate to each interval a “symbol” of the form (+, −,+, …), containing one sign for each variable, with a “+” meaning that that variable is increasing and a “−” meaning it is decreasing. In Fig. 1 a, each box corresponds to one such symbol (for convenience the signs are arranged vertically in all figures, but horizontally in the text). Thus, the continuous time series is converted into a discrete sequence of symbols, which we term the “symbolic dynamics” (18). The algorithm listed below can then determine whether the sequence is consistent with the dynamics of a single negative feedback loop, and if so, the precise order of activators and repressors in that loop. We emphasize that this algorithm is not at all specific to the p53Mdm2 example of Fig. 1; it can be applied to any oscillating time series, with any number of dynamical variables.
The algorithm:

List the order in which the maxima and minima of the variables occur. For example, in Fig. 1 a, the order is p53 max, Mdm2 max, p53 min, Mdmd2 min, p53 max, Mdm2 max, …

If the variables occur in this list in an unchanging cyclic order, then this fixes the order of species in the loop, i.e., a variable activates or represses the one immediately following it in the list. Otherwise, a single negative feedback loop is inconsistent with the time series. Furthermore, if any variable changes sign more than twice in one period then, too, a single loop is inconsistent with the time series.

Construct the symbolic dynamics for the time series, with +/− symbols listed in the order obtained in step 2.

If the symbolic dynamics is not periodic, a single negative feedback loop is inconsistent with the time series. Otherwise, start with any variable and the one pointing to it (say, variables B and A) and note the steps where B changes sign. If at the previous step (before B changed), A had the same sign as B then A represses B, otherwise it activates B (see Fig. 1 b).

This procedure is repeated for each variable to obtain the effect of the preceding variable. For example, in Fig. 1, we conclude that p53 activates Mdm2, and Mdm2 represses p53. If the various sign changes of any variable give inconsistent conclusions, then a single negative feedback loop is inconsistent with the time series.

If the number of repressors in the loop is even, then a single negative feedback loop is inconsistent with the time series.
For p53Mdm2 oscillations, there are only two possibilities for a single negative feedback loop involving only these two. Either p53 activates Mdm2, which represses p53, or p53 represses Mdm2, which activates p53. The above algorithm applied to Fig. 1 a picks the former, for which experimental evidence already existed (19) (discussed further in Examples and Discussion). This is a particularly simple case because it involves only two variables. Our algorithm aids model selection much more when there are more than two variables involved, as shown in the two other examples discussed in the final section of the paper. First, we establish the mathematical basis for this algorithm.
A General Class of Negative Feedback Loops
Consider a feedback loop composed of an arbitrary number, N, of nodes, where each node can be either an activator or a repressor (see Fig. 2 c). Nodes could be genes, proteins, metabolites, or any other chemical species that could, directly or indirectly, activate or repress other nodes in the system. The equations we study are of the form x _{i} is a dynamical variable associated to node i, e.g., the concentration of the chemical species i, or the expression level of gene i. Henceforth, we will call i a chemical species and x _{i} its concentration. The vector field, whose components are g _{i}(x, y), contains all of the basal production, degradation, and possibly selfcatalytic terms, and specifies the interaction between variables. We denote explicitly with the superscript that each species i is either activated (A) or repressed (R) by the species i − 1 immediately preceding it in the loop. We allow for heterogeneity, meaning that the different species can be characterized by different production and degradation rates, and different interaction strengths, i.e., different functions g _{i} for each i. The functions can also be nonlinear. For instance, one way of implementing the threerepressor loop in Fig. 2 b would be for i = 1, 2, 3 (with i = 0 the same as i = 3). These equations model the basal production of each protein (c), the uniform dilution of each proteins by cell growth (γ) and the production of each protein (α) that is repressed by the one behind it in the circuit. The repression we have chosen to be of a standard Michaelis–Menten form, with halfmaximum at K. The Hill coefficient, h, models the cooperative prevention of transcription by h molecules of i − 1 binding at the promoter of i. This is a simplified version of the repressilator (20).
This example provides a single illustration of possible g _{i} functions. In fact, we do not constrain the g _{i} to be like Eq. 1 . The only restrictions on the functions g _{i} are the following two conditions:

All trajectories should be bounded and persistent, meaning that all of the concentrations should stay positive and finite in the time evolution.

Monotonicity: all of the g _{i}(x, y) have to be monotonically decreasing functions of the first argument. Moreover, the g ^{R}(x, y) are decreasing functions of the second argument, whereas the g ^{A}(x, y) are increasing functions of the second argument.
Condition 1 is imposed to ensure that concentrations of the species are well defined and cannot grow infinitely, a biologically plausible constraint: typically for regulatory networks, g _{i} is bounded from above (i.e., there is a maximum possible rate of production) and is dominated by the degradation terms for sufficiently large concentrations, thus ensuring condition 1. Condition 2 implies that activators of a specific process are activators at all concentrations (and similarly for repressors). In other words, we exclude regulation like CI in lambda phage which can activate the P _{RM} promoter at low concentrations, but repress it at high concentrations (21). Another example is the galactose regulator GalR, which at high concentrations of galactose is an activator of the promoter galP2 but in the absence of galactose forms a DNA loop, in which state galP2 is completely repressed (22). However, such examples are relatively rare, and we exclude them from the class of networks we study.
The monotonicity condition can be used to prove that, if the number of repressors in the loop is even, there can be multiple fixed points. This is necessary for bistable or multistable systems, as previously analyzed in ref. 23. In the following, we focus on the case of an odd number of repressors, i.e., negative feedback. Then, as shown in supplementary material, there is one and only one fixed point. Furthermore, a linear stability analysis shows that the transition to instability is necessarily a Hopf bifurcation, which implies that near the transition point there exists a periodic orbit [see supporting information (SI) Appendix ]. However, this stability analysis allows one to study the dynamics only locally, both in the phase space, i.e., close to the fixed point, and in parameter space, i.e., close to the bifurcation value. In the next section, we show that in general, whether the fixed point is stable or unstable, there are restrictions on the trajectory of the system.
Symbolic Dynamics
Our argument is the direct consequence of how the nullclines, i.e., the N manifolds defined by the equations g _{i}(x _{i}, x _{i−1}) = 0, partition the phase space (the positive orthant x _{i} > 0 ∀i). Each nullcline separates two fully connected regions, one in which the ith component of the field, g _{i}, is positive and one in which it is negative. Furthermore, all these manifolds have exactly one point in common, the unique fixed point x*. The phase space is thus the union of 2^{N} sectors, each characterized by the signs of the components of the field. These sectors correspond precisely to the symbols (+, −, −,+, …, −) defined previously.
Note that there cannot be an attractor entirely contained in the interior of a sector because the field does not change sign, the trajectory is bounded and there is only one fixed point. Thus, either the trajectory of the system spirals in toward a stable fixed point, or, if the fixed point is unstable, it will keep crossing from one sector to another indefinitely. In the first case, the symbolic dynamics is a finite sequence of symbols, which ends when the trajectory stops crossing sector boundaries. In the latter case, it will be an infinite sequence of symbols. In either case, adjacent symbols in the sequence will differ by only one sign change.
The key point of our argument is that any of the boundaries can be crossed in just one specific direction, due to the monotonicity of the functions g _{i}(x _{i}, x _{i−1}). This means that not all possible sign changes are allowed. In Fig. 3, we illustrate this using the example of the twospecies negative feedback loop of Fig. 2 a, one repressor and one activator. Fig. 3 b shows the only four transitions possible in this system. Thus, starting from any initial condition, the symbolic dynamics will follow the order shown in Fig. 3 b until the trajectory converges to the stable fixed point and there are no further sign changes. This example gives a good pictorial idea of the fact that the nullclines behave like oneway doors.
In the general Nspecies case, the same phenomenon occurs, with the symbolic dynamics obeying the following rules:

If the variable (i − 1) represses i, the nullcline i can be crossed if g _{i} and g _{i−1} have the same sign.

If the variable (i − 1) activates i, the nullcline i can be crossed if g _{i} and g _{i−1} have opposite signs.
This is the main mathematical result of this paper, and leads directly to the algorithm described earlier; its derivation is described in more detail in supplementary material. To emphasize the point further, the direction in which the nullcline g _{i} = 0 can be crossed at a given point depends on the position of that point relative to only one other nullcline, g _{i−1} = 0. It does not depend on any other nullcline. In simple words, the rules encode the fact that an increasing activator can make the affected concentration increase (but not decrease), whereas an increasing repressor can make the affected concentration decrease (but not increase). Note that the allowed transitions apply to any trajectory, even transients. Thus, if one is analyzing oscillatory time series it doesn't matter whether the oscillations are sustained, or the measurement is of transients or damped oscillations.
To determine the general scenario which is compatible with these rules, consider that when a nullcline is crossed, the symbolic dynamics makes a transition between two states differing by just a single sign. We say that there is a mismatch between two adjacent signs if the nullcline depending on these two variables can be crossed according to the rules defined above. The effect of crossing the nullcline i is to remove the mismatch between i and i − 1. If there was also a mismatch between i and i + 1, it too is removed, otherwise it is created. For a negative feedback loop, we can show that the number of mismatches has to be odd, and cannot increase.
Now consider what happens if the fixed point is unstable. If there is just one mismatch, it can only keep traveling around the loop, in the direction of the loop arrows. This means that the symbolic trajectory is periodic of period N. In the general case, we can visualize the symbolic dynamics as several mismatches traveling around the loop in the same direction. Whenever two mismatches “hit” each other, they annihilate. Eventually, the number of mismatches will reach some limit, where each mismatch stays safely distant from the others. The length of the loop limits how many mismatches can, in principle, coexist; for example, only one mismatch can survive if N < 4. In practice, even in long loops, it is likely that only one mismatch survives, and we will restrict to this case in the following.
An interesting consequence of this periodicity is that any of the N nullclines can be used to define a Poincaré map for the dynamical system. Periodic oscillations in the symbolic dynamics translate into a fixed point or a stable periodic orbit of the Poincaré map (ref. 24 proves that quasiperiodicity and chaos are impossible for such systems). In this general class of systems, when the fixed point is unstable, the dynamics is oscillatory with well defined properties: each of the concentrations has exactly one maximum and one minimum during a time period of the symbolic dynamics, and the fact that the mismatch travels in the direction of the feedback loop implies that the sequence of maxima and minima has to follow the order of the species in the loop. From the particular order observed, it is also possible to argue which species acts as an activator and which as a repressor. Furthermore, the observation of a time series which is incompatible with the symbolic dynamics rules allows one to exclude a dynamics of the form of Eq. 2 , generally suggesting a topology more complicated than a simple feedback loop or more subtle effects like time delays and nonmonotonic regulation. In other words, the algorithm described previously is a straightforward consequence of the rules followed by the symbolic dynamics
Note that our method works even if one does not measure the time series of all of the species belonging to the loop. The algorithm gives a coherent conclusion about the overall sign between the variables: for example, a variable A will appear as an activator of a variable B if there is an even number of “unobserved” repressing links between them (see SI Appendix ). The following examples will further clarify these points.
Examples and Discussion
We now apply the above ideas to extract information about the loop structure from experimentally observed time series of three systems: p53Mdm2 oscillations in mammalian cells, circadian expression of kai genes in Synechocystis cyanobacteria, and cyclic binding of protein cofactors with DNA at the estrogensensitive pS2 promoter in human breast cancer cells.
Our first example is the well known p53Mdm2 negative feedback loop, already discussed in the introduction. The tumor repressor protein, p53, activates transcription of the mdm2 gene (19). Mdm2, once produced, binds to p53 preventing it from acting as a transcription factor, and subsequently ubiquitinates it, which enhances its proteolytic breakdown (19). This negative feedback loop is precisely the structure our algorithm predicts from the oscillating concentrations of p53 and Mdm2 in Fig. 1. In a couple of cases, there is some ambiguity about the order of maxima and minima. However, the pattern is clarified by comparing two separate time intervals in which the symbolic sequence is unambiguous. Note that both regions exhibit the same periodic symbolic sequence. Thus, the observed oscillations are consistent with a dynamics of the form of Eq. 2 . However, there must be at least one other unobserved species taking part in the loop, because the fixed point is always stable for N = 2. Indeed, several threevariable models of p53–Mdmd2 oscillations have been examined, which assume the third variable to be either an Mdm2 precursor (e.g., Mdm2 mRNA) or a third protein that interacts with p53 or Mdm2 (7). Although time delays (11) and other feedback loops (25) present in this system could also be important, we can conclude that a model like Eq. 2 is a serious candidate for a zeroorder model, being simple and reproducing correctly the qualitative behavior of the components with the correct interaction signs.
The next example involves circadian oscillations of gene expression in cyanobacteria. Cyanobacteria are the only bacterial species with a circadian clock and several of their cellular functions appear to be under circadian control (4). In Synechococcus elongatus, a cluster of three genes, kaiA,B,C, were found to be essential: deletion of any of these genes eliminated the oscillations (26). Fig. 4 a shows circadian rhythms in the expression levels of genes coding for homologues of the KaiA,B,C proteins in one Synechocystis strain. The symbolic dynamics is consistent with a threevariable feedback loop (Fig. 4 a Inset), where kaiA activates kaiC1, which represses kaiB3, which, in turn, activates kaiA. The first two of these predicted interactions exist in Synechococcus (26), whereas the third is a new prediction for how kaiA is brought into the loop. Note that our analysis only provides the sign of this interaction. It does not reveal the molecular mechanism of the interaction, nor whether the interaction is direct or through intermediate steps.
Finally, we consider the cyclic binding of cofactors to the estrogensensitive pS2 promoter. A coordinated sequence of binding and unbinding events modifies the DNA packing and nucleosome structure to enable transcription to proceed (28). This is a case where no model exists and not all of the proteins involved have been identified. Our method is particularly suited for such a case because it does not matter whether the dynamics of only a subset of the proteins involved is available. Fig. 4 b shows oscillations in the binding of 4 proteins, after the addition of estradiol (28). Estradiol receptor (ER) binds estradiol and is required for initiating transcription. Pol II is the RNA polymerase that transcribes the gene. TRIP1 is a component of the APIS proteasome subunit, whereas HDAC1 is involved in deacetylation of histones (28). The symbolic dynamics is consistent with the model shown in Fig. 4 b Inset. In this case, each variable measures the amount (in arbitrary units) of bound protein at the pS2 promoter. The predicted links indicate how a bound protein affects the probability of binding (or of remaining bound) of another one in the sequence. For example, the link from ER to Pol II indicates that ER, when bound at the promoter, increases the recruitment probability of Pol II. Ref. 29 models the dynamics of Fig. 4 b using a positive feedback loop, requiring >200 intermediate steps, which has only activating links. Our analysis suggests, however, that a negative feedback loop is a plausible hypothesis as the cause of oscillations, and predicts the existence of a repressive link between HDAC1 and ER.
Acknowledgments
We thank Ian Dodd for critical reading of the manuscript and many useful suggestions. This work was supported by the Danish National Research Foundation.
Footnotes
 ^{§}To whom correspondence should be addressed. Email: pigo{at}nbi.dk

Author contributions: S.P., S.K., and M.J. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/cgi/content/full/0610759104/DC1.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵
 ↵
 ↵
 ↵

↵
 Hoffmann A ,
 Levchenko A ,
 Scott ML ,
 Baltimore D

↵
 Nelson DE ,
 Ihekwaba AEC ,
 Elliott M ,
 Johnson JR ,
 Gibney CA ,
 Foreman BE ,
 Nelson G ,
 See V ,
 Horton CA ,
 Spiller DG ,
 et al.

↵
 GevaZatorsky N ,
 Rosenfeld N ,
 Itzkovitz S ,
 Milo R ,
 Sigal A ,
 Dekel E ,
 Yarnitzky T ,
 Liron Y ,
 Polak P ,
 Lahav G ,
 et al.

↵
 Hirata H ,
 Yoshiura S ,
 Ohtsuka T ,
 Bessho Y ,
 Harada T ,
 Yoshikawa K ,
 Kageyama R

↵
 Leloup JC ,
 Goldbeter A

↵
 Krishna S ,
 Jensen MH ,
 Sneppen K

↵
 Tiana G ,
 Jensen MH ,
 Sneppen K
 ↵

↵
 Thomas R
 Gardiner CW
 ↵
 ↵

↵
 Forger DB ,
 Peskin CS

↵
 Lauzeral J ,
 Halloy J ,
 Goldbeter A

↵
 Strogatz S

↵
 Wu X ,
 Bayle JH ,
 Oldon D ,
 Levine AJ
 ↵

↵
 Dodd IB ,
 Perkins AJ ,
 Tsemitsidis D ,
 Egan JB
 ↵

↵
 Angeli D ,
 Ferrell JE ,
 Sontag ED
 ↵
 ↵

↵
 Ishiura M ,
 Kutsuna S ,
 Aoki S ,
 Iwasaki H ,
 Andersson CR ,
 Tanabe A ,
 Golden SS ,
 Johnson CH ,
 Kondo T

↵
 Kucho K ,
 Okamoto K ,
 Tsuchiya Y ,
 Nomura S ,
 Nango M ,
 Kanehisa M ,
 Ishiura M
 ↵
 ↵
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Applied Mathematics
Biological Sciences
Related Content
 No related articles found.
Cited by...
 Coupled Positive and Negative Feedbacks Produce Diverse Gene Expression Patterns in Colonies
 Dramatic reduction of dimensionality in large biochemical networks owing to strong pair correlations
 Variation in Candida albicans EFG1 Expression Enables HostDependent Changes in Colonizing Fungal Populations