Homogeneous Time Constants Promote Oscillations in Negative Feedback Loops.

Biological oscillators are present in nearly all self-regulating systems, from individual cells to entire organisms. In any oscillator structure, a negative feedback loop is necessary, but not sufficient to guarantee the emergence of periodic behaviors. The likelihood of oscillations can be improved by careful tuning of the system time constants and by increasing the loop gain, yet it is unclear whether there is any general relationship between optimal time constants and loop gain. This issue is particularly relevant in genetic oscillators resulting from a chain of different subsequent biochemical events, each with distinct (and uncertain) kinetics. Using two families of genetic oscillators as model examples, we show that the loop gain required for oscillations is minimum when all elements in the loop have the same time constant. On the contrary, we show that homeostasis is ensured if a single element is considerably slower than the others.

Timekeeping elements coordinate and synchronize most processes required to sustain life, from the physiology of individual cells to the daily rhythms of entire organisms.[1] The design principles underlying the operation of biomolecular clocks have been investigated by dissecting natural systems,[2−5] as well as by building de novo molecular oscillators in an effort to identify minimal requirements for periodic behaviors.[6−10] Experiments[6,2,7,9,11] and modeling[12−21] have established that a necessary requirement for a system to exhibit oscillations is the presence of a negative loop. Conversely, any negative feedback loop potentially leads to oscillations provided that the loop includes destabilizing features, for instance delaying elements associated with a large feedback gain.[7,22,23] Delay can be introduced by a variety of phenomena. In addition to transcription and translation steps, delay is increased by transport, mRNA splicing and stabilization, phosphorylation, and protein maturation.[24,25,43] While an explicit delay term is useful to aggregate these phenomena in a single, black-box parameter, in other instances it is possible to model their individual kinetics as a chain of interconnected subsystems.[17,25,26] In a linearized closed-loop system, a chain of slow subsystems has the same effects as an explicit delay in inducing oscillations.[22] The overall delay of the negative feedback chain depends on all the time constants of the dynamic elements. Given a sufficient delay, the larger the feedback gain, the likelier is the onset of persistent oscillations.

Despite the simplicity of these design principles, it is unclear how the individual time constants and the overall gain affect the emergence of oscillations. Are we equally likely to observe oscillations in a system with homogeneous versus nonhomogeneous time scales? How large should the loop gain be in either case? These questions are particularly relevant in the context of genetic oscillators, where time constants are determined by mRNA and protein degradation, transport, and processing rates, which may widely vary among the oscillator components, and a large loop gain is energetically expensive, because it depends on the production rate of the components.

In this paper, we address the following question in mathematical terms: Given a negative loop of first order elements, each associated with its own time constant, which is the choice of the time constants that requires the smallest gain to allow for persistent oscillations? We demonstrate that homogeneous time constants are the most favorable choice when a small loop gain is desired. In particular we prove that (1) the smallest negative feedback gain required to trigger oscillations is achieved when all time constant are equal; (2) the smallest gain is invariant under a homogeneous scaling of all time constants and the period of the oscillations is proportional to the scaling factor; (3) as a converse result, in a negative feedback loop, the best strategy to avoid oscillations is to have a single element of the chain that is much slower than all the others, and this fact explains why, in several pathways with negative feedback, the presence of a slow element ensures a robust nonoscillatory behavior.

We apply our results to well-known genetic oscillators, the Goodwin oscillator and a two-node (inhibitor-activator) oscillator, and we derive exact (necessary and sufficient) conditions for the existence of parameters that ensure oscillations.

Architecture of Candidate Negative Feedback Oscillators

As candidate oscillator architectures, we consider negative feedback loops of the following form:representing the series connections of n ordered subsystems in which any element has a positive influence, quantified by parameters b > 0 (with i = 2, ..., n), on the next one, while u is a perturbing input that triggers oscillations. Each ordinary differential equation (ODE) in the model in eq is suited to model phenomena such as production, conversion, processing, and degradation of molecular components (mRNA and proteins) interconnected in a regulatory chain.[27] The model can also capture enzymatic processes that operate at low substrate concentration relative to the binding affinity of the enzyme and substrate; in this regime, Michaelian or Hill-type reaction rates become approximately linear (first-order rates). A negative feedback loop is generated via the inhibitory effect of the last element in the chain on the first one, quantified by parameter b1 > 0. The parameter τ represents the time constant of process i (which can be seen as the speed of the reaction of the species x due to the regulatory effect of x).

A similar negative feedback structure can be found in many oscillatory systems.[17] We now take the Laplace transform of these ODEs: we formally replace x(t) with X(s) and the derivative d/dt with the complex variable s, ẋ(t) → sX(s). After Laplace-transformation, the model in eq can be rewritten as a block-interconnection of elements:as shown in Figure .

Figure 1

Loop of n first order systems: block diagram.

The quantityis called the loop gain and has a fundamental role. It is the product of all the interaction strengths and thus represents the cumulative strength of the loop. It turns out that the characteristic polynomial depends on the product κ only, and not on the individual parameters b; hence, even if the individual rates b are changed, the system behavior remains the same as long as their product is unchanged (see the Supporting Information for the detailed derivation). The onset of oscillations in this negative feedback loop is therefore associated with two fundamental ingredients:

The time constants τ, which introduce an overall delay in the loop;

A sufficiently large feedback gain κ > 0.

In the next section we ask ourselves whether there is any ideal relationship between the loop gain and the time constants to achieve or avoid oscillatory behavior.

Influence of Time Constants on the Oscillatory Regime

We next investigate how the time constants τ influence the onset of persistent oscillations. We define τ = [τ1 τ2 ... τ], the vector of time constants, and consider the characteristic polynomial associated with the (linearized) system of ODEs:For κ = 0, the roots of p(s, τ) are λ = −1/τ, real and negative, hence the system response has an exponentially decreasing pattern. For large values of κ, p(s, τ) has complex roots, associated with oscillations.

The oscillations are damped if the roots have a negative real part. To have persistent oscillations, the roots of p(s, τ) must reach and cross the imaginary axis in the complex plane. This can happen only if(as discussed at the end of this section; see also ref (28)). Henceforth, we assume that the necessary condition n ≥ 3 is verified.

For n ≥ 3, let us increase κ. Then, there exists a critical gain κ* such that, for all κ > κ*, p(s, τ) has complex roots with positive real part (namely, the system becomes unstable). For κ = κ*, p(s, τ) has two purely imaginary roots ± jω*, while the other roots have negative real part. The limit value κ* is associated with the onset of an oscillation with frequency ω*/2π; we call ω* critical pulsation. Note that ω* ≠ 0: p(s, τ) cannot have 0 as a root for κ > 0, since p(0, τ) = κ + 1 ≠ 0.

We can formally define the critical gain κ* as the smallest value of κ for which p(s, τ) has a pair of purely imaginary roots (corresponding to the stability limit). The value κ* depends on the time constants τ and we can write

Which are the most favorable values of τ to promote oscillations? We address this question in terms of the minimum critical gain, by seeking a value τ* = [τ1* τ2* ... τ*] that minimizes the critical κ*(τ) enabling the onset of oscillations.

Problem: Find a value τ* that minimizes κ*(τ) in eq .

Main result: The problem is solved by a value τ* withThis result is proved in the Supporting Information (Theorem 1): our proofs are based on frequency analysis tools, linear algebraic tools and principles of convex optimization.

Therefore, an essential factor to promote oscillations in a negative feedback loop is the homogeneity of the time constants of the subsystems involved in the loop.

Further, we find that scaling the time constants influences exclusively the critical pulsation, without affecting the critical gain: when the time constants are scaled as τ → στ, for arbitrary σ > 0, the critical gain κ* is invariant, κ*(στ*) = κ*(τ*), while the critical pulsation scales proportionally to σ: ω* → σω* (cf. Corollary 1, Supporting Information).

Also, the critical gain κ* is a decreasing function of the number of elements in the loop (see the Supporting Information, Proposition 2, for details).

Our result (the critical gain that allows for oscillations is minimized when all the time constants are equal) indirectly suggests how to prevent a system from oscillating. This aspect is relevant in the context of biological and biochemical feedback loops, in all the situations where it is important to preserve homeostasis and oscillatory behaviors must be avoided. Being κ* a decreasing function of n, long feedback chains are more prone to instability, which can be of an oscillatory type. Hence, a natural question is which is the best strategy to avoid oscillations in the loop. Our result suggests that incongruous time constants lead to a robustly nonoscillatory behavior. Let us now consider the complementary question: assuming (without restriction) that the time constants are normalized aswhere Ttot is the overall loop delay, which is the best distribution of time constants to prevent oscillatory behaviors? We find that, roughly speaking, it is better to have the delay concentrated in a single subsystem (see the Supporting Information, Proposition 3, for details). Then, a robust strategy to prevent oscillations is, for instance, including in the loop a single subsystem that is much slower than the others, so that their time constant is negligible with respect to the slow part. This result also explains the previous statement about the necessity of condition n ≥ 3 to have persistent oscillations. Indeed, setting τ = 0 is mathematically equivalent to neglecting the ith process, since then 1/(1 + τs) = 1.

A fundamental consequence of our results is the following: in a negative feedback loop, the presence of a single slow element is an effective strategy to preserve stability and prevent oscillations. Indeed, several negative feedback loops in nature are practically always stable, and this can be explained by noting that the involved time constants are very different. For example, biologically, degradation rates of mRNA and proteins may vary in a broad range and, as foreseen by our results, this variability could contribute to stabilizing negative loops, making it difficult to achieve oscillations. Yet, we point out that, when the kinetics of a molecular species are much slower or much faster than the rest of the system, they can be simply eliminated from the model via time scale separation methods;[27] this type of drastic time scale variability may affect the capacity for oscillations when the system dimension collapses below 3.

Examples

Our results allow us to derive analytical bounds in the parameter space for the oscillatory regions of general families of genetic oscillators. These bounds are also numerically verified in the following sections via random parameter sampling.[29]

Goodwin Oscillator

The well-known Goodwin oscillator[30] is associated with the following equations:

The model is characterized by the number n of stages, by the cooperativity (Hill) coefficient N, by the apparent dissociation constant K, and by the rate constants a and b. Rates a and b can model protein translation and degradation, mRNA processing phenomena, or protein phosphorylation/dephosphorylation.[31] All these parameters are positive. As shown in the Supporting Information (Section 2.1), the system admits a single positive equilibrium. The emergence of sustained oscillations depends on the choice of the parameters, which influences the values of both the variables at steady-state and of the entries of the Jacobian matrix. By applying our main result to the Goodwin oscillator model, we discover that there exists at least one choice of the parameter values that leads to oscillations (namely, for which the linearized system admits complex eigenvalues with nonnegative real part) if and only if(the full derivation is in the Supporting Information, Section 2.3). The characteristic equation, equating the characteristic polynomial to zero, isHere, x̅ is the steady state value of x derived by solving the equilibrium equationfrom which existence and uniqueness of the equilibrium are proven in the Supporting Information (Section 3.1). If we could arbitrarily choose the rates a and b, then we could achieve any positive steady-state x̅. Hence, the gain in eq ,could take any value in the interval: 0 < κ(x̅) < N. The value of the ith time constant for this system is τ = 1/b. Hence, our result tells us that the (minimum) critical gain is obtained by setting b1 = b2 = ... = b = b (which leads to equal time constants), henceThe condition in eq is necessary and sufficient for eq to admit imaginary solutions. As exemplified in Figure , if the condition in eq is satisfied, oscillations are possible if the rates a are large enough with respect to the rates b; this guarantees a large equilibrium value x̅, hence κ(x̅) is larger than the critical gain.

Figure 2

Goodwin oscillator: equilibrium conditions for different ratios a/b. The orange line represents the first expression in eq , while the blue lines represent the second expression in eq for different values of the ratios a/b. Their intersections give, on the vertical axis, the equilibrium values of x for various choices of a/b.

Note that, if the oscillator includes three stages (n = 3), the condition in eq becomes and the minimum value of the Hill coefficient N to have oscillations is N = 8, consistently with the results in ref (32).

In Figure we numerically compute the oscillatory region in the N–n space when a = a and b = b, for all i. In each panel, the black line represents the condition in eq converted to an equality, and delimits the region where oscillations can occur. As a further numerical experiment, we fixed a = a for all i and randomly generated different values of the parameters b, taken from different distributions (normal and uniform distribution) with the same expected value E[b] = 1 and variance ϵ. With this sampling method, the total delay is not necessarily constant for all samples: the degradation rates are randomly generated and are drawn from different distributions, which all have the same average, but which can be more or less spread depending on the value of the variance ϵ. The lower the variance, the more homogeneous are the degradation rates b. Figure shows the fraction of oscillating samples (parameter choices for which characteristic polynomial has positive-real-part roots) as a function of ϵ. As predicted by our analytical results, decreasing the variance increases the likelihood of oscillations. Choosing homogeneous time constants favors oscillations, but also other design decisions are important: for instance, the overall loop gain needs to be high enough.

Figure 3

Oscillatory domain of the Goodwin oscillator with K = 1 in the (n, N) plane, for various choices of homogeneous rates, a = a and b = b for all i. The black line represents ; hence, eq is satisfied in the whole region above. We compute the solutions of eq : red dots indicate an oscillatory behavior, blue dots indicate no oscillations, while gray dots indicate that no equilibrium can be found computationally due to numerical problems.

Figure 4

The fraction of oscillating samples of the Goodwin oscillator is largest when randomly drawn degradation rates are homogeneous. We simulated the model with K = 1 and n = 5; in the top panels N = 8, while in the bottom panel N = 10. We took a = a and randomly generated rates b with expected value E[b] = 1 and variance ϵ. (In each simulation, the randomly generated parameters are kept constant during all the integration steps of the ODEs.) We show the fraction of oscillating samples as a function of the variance ϵ when the rates b are taken from a normal distribution (top left, and bottom) and a uniform distribution (top right). In all panels, 1000 parameter samples are drawn per data point.

A Two-Node Oscillator

Consider a two-node oscillator given by the feedback interconnection of an activated module and an inhibited module:[9,48−50]The ODEs of variables r1 and r2 represent mRNA dynamics, and p1, p2 represent protein translation. As earlier, N is a Hill coefficient, and K1 and K2 are apparent dissociation constants. Parameters α1 and α2 are maximal mRNA transcription rates, and β1, β2 are mRNA degradation rates. Finally γ1, γ2 and δ1, δ2 are protein translation and degradation rates. This model can serve as a coarse-grained representation of a variety of molecular clocks. Many genetic oscillators result from the interconnected dynamics of inhibitor-activator elements,[25] such as the p53-mdm-2[33] and the IκB-NF-κB[34] oscillators; this architecture has also been demonstrated in artificial in vitro transcriptional oscillators.[9,35,36] Here we assume that the mRNA dynamics evolve on a time scale that is comparable to that of the proteins; this means that the order of the model cannot be reduced via time scale separation arguments. On the basis of our main result, we can prove that there exists at least one choice of parameters for which the system exhibits sustained oscillations (namely, the linearized system has complex eigenvalues with a nonnegative real part) if and only ifThe linearized system has characteristic equationwhere p̅1 and p̅2 are the steady-state values of p1 and p2. The loop gainranges in the interval 0 < κ(p̅1, p̅2) < N2; it is a decreasing function of p̅1, and an increasing function of p̅2. In view of our main result, the minimum critical gain guaranteeing oscillations is achieved when β1 = β2 = δ1 = δ2 = τ. The corresponding critical equation istherefore no unstable complex eigenvalues, hence no oscillations, can be found unless N > 2 (the full derivation is in the Supporting Information, Section 3.3). Numerical simulations illustrating and confirming this analytical result are in Figure . We have also generated random instances of the oscillator to show how increasing the variance in the delay decreases the chances of oscillatory behavior (Figure ). We note that in this particular example, if the mRNA dynamics are much faster than the protein dynamics, the model can be reduced to include exclusively the protein kinetics; in that case, oscillations cannot be achieved for any choice of the parameters.

Figure 5

Oscillatory regime of the two-node oscillator. We compute the solutions of eq , with α1 = α2 = α, β1 = β2 = β, γ1 = γ2 = 1, δ1 = δ2 = 1 and K1 = K2 = 1, and indicate with red dots an oscillatory behavior, with blue dots no oscillations. Parameter sets that give rise to oscillations cannot be found for N = 2, while they are easy to find for N > 2.

Figure 6

The fraction of oscillating samples of the two-node oscillator is largest when randomly drawn degradation rates are homogeneous. We simulated the model with N = 3, γ1 = γ2 = 1, K1 = K2 = 1, α1 = α2 = α and randomly generated (β, δ) with expected value E[(β, δ)] = (1, 1) and variance ϵ. (In each simulation, the randomly generated parameters are kept constant during all the integration steps of the ODEs.) We show the fraction of oscillating samples as a function of the variance ϵ when (β, δ) are taken from a normal distribution (left) and a uniform distribution (right). 1000 parameter samples are drawn per data point.

Conclusion and Discussion

We have demonstrated that homogeneity of the time constants within a negative feedback loop can facilitate the emergence of oscillations, in that it minimizes the critical loop gain (minimum gain to achieve oscillations). We have also shown that scaling of the (uniform) time constants influences the critical frequency, but does not affect the critical gain. A converse result is that a candidate oscillator can be stabilized, i.e., oscillations cannot occur, by increasing a single (arbitrarily chosen) time constant of the loop with respect to the others. The negative feedback architecture we consider is general, and it can be specialized to model many biomolecular oscillators.[25]

The gain of the biomolecular feedback loops we consider is proportional to the ratio of production and degradation rates of its components, and to the cooperativity coefficient of regulatory molecules (Hill coefficient N). Maintaining the lowest gain that can yield oscillations in the network is therefore tantamount to operating the circuit with minimum consumption of transcription and translation resources, minimum kinase activity, as well as with minimum copy number of regulators. This energy-efficient scenario can be achieved when the time constants in each process are similar (degradation, transport, and processing rates). This requirement may be easy to satisfy if these time constants are globally regulated for all components (for instance, mRNA and protein degradation).

It must be pointed out that, in our analysis, we have considered systems consisting of a single negative feedback loop: although this is a very common structure for biological oscillators,[17,25] it is not the only one. For this particular structure, we have argued that

Short negative loops have a stabilizing effect (which makes the onset of sustained oscillations less likely, because a higher loop gain is needed);

Long negative loops can favor the onset of sustained oscillations, and the most favorable case is that in which the time constants of the system in the loop are similar;

If one or two time constants are significantly larger than the others, then the long loop actually behaves as a short one and the probability of having sustained oscillations is smaller (because a higher loop gain is needed for the onset of oscillations).

However, if several feedback loops are concurrently present, our analysis does not apply, and the above statements are no longer true. In particular, our findings are valid in the absence of self-catalytic reactions (i.e., of positive self-loops). In the presence of a positive self-loop, also a single negative loop involving two nodes only can be easily destabilized. For instance, consider the systemwhere a, b, c and d are positive parameters. As shown by our results, this loop cannot be destabilized. However, this system can exhibit sustained oscillations if a can be negative (take for instance a = −d): in this case our results no longer hold, because now x1 is self-catalytic.

In practice, positive self-loops are not very common. Yet, a positive feedback loop can result from a chain of reactions. Indeed, different oscillator architectures are based on the coexistence of positive and negative loops (it is important to stress that the presence of a negative loop is necessary for the onset of oscillations[14,15]). In particularAn example is the genetic network present in the bread mold Neurospora crassa, which has been shown in ref (37) to be a successful biological oscillator; we investigate this case study in the Supporting Information, where we show that the result proposed in ref (37) are fully consistent with our analysis.

A (possibly short) negative loop can be destabilized by the concurrent presence of another loop that is positive.

Negative-feedback oscillators are very common in nature and appear also very robust. For instance, the Hes1 and Hes7 oscillators in mammalian embryos consist of a negative autoregulation loop where Hes protein represses its mRNA production.[38] These oscillators could be modeled taking n = 2 in eq ; however, in this case the ODE solution does not admit sustained oscillations even for very large values of N. Addition of an explicit delay, discrete or distributed, to Hes autoregulatory models yields oscillatory solutions for physically acceptable values of N(39,40) (Hes1 and Hes7 are dimers). Similar observations can be made for the p53-mdm-2 and the IκB-NF-κB oscillators.[41] The explicit delay, which captures mRNA processing and transport steps,[42] could be alternatively modeled as a chain of intermediate subsystems; while the number and kinetics of these steps are unknown, our results suggest that oscillations would be more likely to occur if they had similar time scales. Interestingly, the Hes1 oscillator requires nearly identical mRNA and protein half-lives (≈ 23 min) to operate.[43]

Our results are particularly relevant for the design of artificial negative feedback clocks. While the architecture of the Goodwin oscillator is attractive due to its simplicity, it has been difficult to build synthetic examples without including positive feedback, high Hill coefficients, or additional nonlinearities to destabilize the system.[6,7,16,44] The mathematical models developed to capture the dynamics of these artificial oscillators often assume similar degradation rates for all the mRNA and protein species. However, recent experiments on the famous repressilator circuit suggest that protein degradation rates in the original design might have been subject to temporal fluctuations caused by competition for shared proteolytic machinery, occurring due to the presence of protein degradation tags meant to reduce their half-life.[45] Removal of the degradation tags resulted in more regular (although slower) oscillations at the population level. It is possible that in the absence of degradation tags dilution (due to cells dividing) becomes the dominant time constant, which should be uniform for all the repressor proteins. In light of our results, a more homogeneous protein half-life could explain the improved robustness of the oscillations.

In conclusion, our work highlights that the variability of time constants within negative feedback oscillators could have under-appreciated effects on the dynamics; better estimation of these parameters could help explain the robustness of many natural oscillators. Conversely, we expect that the construction or improvement of artificial oscillators could be facilitated by ensuring that the modules being interconnected evolve with similar time scales.

Methods

The formal proofs of our results, which employ mathematical tools from dynamical systems and systems and control theory, as well as the detailed mathematical analysis of the proposed examples, are in the Supporting Information.