|Home | About | Journals | Submit | Contact Us | Français|
Robust oscillations have recently been observed in a synthetic gene network composed of coupled positive and negative feedback loops. Here we use deterministic and stochastic modeling to investigate how a small time delay in such regulatory networks can lead to strongly nonlinear oscillations that can be characterized by “degrade and fire” dynamics. We show that the period of the oscillations can be significantly greater than the delay time, provided the circuit components possess strong activation and tight repression. The variability of the period is strongly influenced by fluctuations near the oscillatory minima, when the number of regulatory molecules is small. We explore the effect of the positive feedback loop on the robustness of these oscillations.
Recently, we constructed a synthetic two-gene oscillator which exhibited highly robust and tunable oscillations . The design of the oscillator was motivated by the original theoretical proposal  in which oscillations were due to the interplay of positive and negative feedback loops operating at different time scales. However, experimental results suggested that the core of the oscillator was a single multi-stage NFB circuit (Fig. 1), while the PFB loop helped to increase the robustness and tunability of oscillations. In that study  we also developed a detailed biochemical model of this circuit, which explicitly incorporated multiple steps leading from gene transcription to formation of mature protein multimers and protein-DNA interactions. Such models are useful for making testable predictions, however they are not transparent enough to permit analytical investigation and thus gaining qualitative insights into the mechanisms controlling the circuit dynamics. In particular, it was not clear how relatively small (3-5 min) delays caused by the intermediate steps in protein synthesis can lead to rather long-period (20-60 min) oscillations. In this Letter, we replace the feedback cascade (Fig. 1) by a single reaction with a fixed deterministic delay time. It is well known that delayed auto-repression can lead to oscillatory gene expression [3-8]. Here we demonstrate analytically that in a strongly nonlinear regime this system exhibits what we term degrade-and-fire (DF) oscillations, in analogy with integrate-and-fire models . An essential feature of this regime is that the period of oscillations can be arbitrary long compared to the delay time. In this regime, we are able to derive analytic estimates for the mean period and its variance.
While the results are mostly presented for the NFB-only case, we also address the role of an additional PFB loop, which was present in the experimental realization of the oscillator . Details of these calculations can be found in the Supplementary Information (SI) .
We consider a system with only two reactions: production of a repressor protein and its degradation. These reactions are characterized by the rates and to produce or degrade repressor, respectively. The production rate is negatively regulated by the repressor itself, and this reaction is assumed to take a finite time τ from the start to completion of a mature repressor protein, , where r(t) is the number of repressor molecules at time t, rτ(t) r(t − τ), is a Hill function governing NFB . Repressor is degraded enzymatically at a maximum rate γr, and we also assume a small first-order degradation due to dilution with rate β, with R0 characterizing the concentration of the enzyme.
We simulated these two stochastic reactions using the extension of the Gillespie algorithm  proposed in  to treat delayed reactions. For sufficiently large α γr C0/τ > R0/τ  the system exhibits strongly nonlinear degrade-and-fire oscillations in which the number of protein molecules rises rapidly from zero, and then slowly decays back to zero, after which the process repeats. Figure 2a shows a typical trajectory r(t).
DF oscillations consist of two main phases. Each oscillation begins with a “production” phase. During this phase the existing repressor concentration is too small to preclude transcription and nascent repressor is produced (to become active a time τ later). The second phase, which we call the “degradation” phase, primarily consists of the degradation of functional repressor. For the majority of the degradation phase, production is tightly repressed, and transcription begins only when the level of the repressor falls near the threshold value C0. The onset of the production phase can be associated with time ton defined by condition F(ron) = γr (Fig. 2b). However, due to production delay, new repressor molecules do not appear en masse until time ton + τ. This causes a rapid rise of the repressor level. Once the repressor level significantly exceeds C0 again, production of new repressor ceases, and the concentration of repressor reaches a maximum within a time τ later, after which a new decay phase commences, and so on. Because of the stochasticity of the biochemical reactions, the amplitudes of individual pulses vary, but their shapes are similar. The properties of DF oscillations can be analyzed in terms of these alternating degradation and production phases.
In the deterministic limit, the mass-action kinetics of delayed production and degradation can be expressed by a delay-differential equation
This equation describes a sequence of nearly triangular pulses of repressor similar to that of Fig. 2a, however with a constant amplitude and duration. The key to analytically solving Eq. (1) is to recognize that, in this strongly nonlinear regime of large-amplitude pulses, the production of repressor is negligibly small during most of the degradation phase. For simplicity of analysis, we further assume R0 → 0 and β = 0, in which case the decay becomes a zeroth-order reaction. While this assumption is not essential, it greatly simplifies the subsequent analysis. Note that in this approximation we have to augment the corresponding simplified model by the natural condition r(t) ≥ 0. In this case the solution during the degradation phase is linear in time,
Here t0 is the time at which r first reaches zero (see Fig. 2b). This law for linear decay is essentially independent of the amplitude of the previous pulse. Thus, individual pulses are decoupled from each other.
The structure of the pulse at t > t0 can be reconstructed iteratively based on the knowledge of the decay phase solution (2) at t < t0. In the integral form
for some times t, ti. Eq. (3) advances the repressor trajectory in intervals of time duration τ, except when r(t) = 0, which should be treated separately. With this expression, we can first determine the trajectory for ton + τ < t < ton + 2τ using ti = ton + τ and then in a similar manner compute the trajectory at ton + 2τ < t < ton+3τ. During this second iteration, the solution reaches a maximum at which r(t) is sufficiently large to neglect the production term, and the new degradation phase commences. Thus, only two iterations of Eq. (3) from time t = ton + τ are sufficient to determine the magnitude P for a pulse of repressor.
The degradation of a burst of repressor takes an approximate time TD = P/γr and is typically the dominant contribution to the oscillation period. A simple, but useful, estimate for P is the production P0 that occurs during t0 < t < ton + τ when r(t) = 0, i.e. when the production rate is maximal. P0 is given by
As seen from this expression, the period of oscillations can be arbitrarily large compared with the delay time τ, depending on the production and decay rates of the protein. This estimate for the period has a maximum at a certain γr and becomes zero when γr → α, ron/τ. A more accurate period estimate can be obtained using the integral expression (3), see SI. Fig. 3a compares both these estimates to the results of stochastic simulations and direct integration of the deterministic Eq. (1).
The case of nonzero but small β (βτ 1) can also be analyzed. Production bursts of repressor are only slightly perturbed in this case. However, if the number of produced proteins P γr/β, the initial degradation of the burst is exponential rather than linear, and so the time TD(β) to degrade repressor to zero can be significantly less than for β = 0. It is easy to show that βTD(β) ≈ ln (1 + βTD(0)), where TD(0) ≈ P/γr. Thus, the period of oscillations in growing cells cannot be much larger than the dilution time β−1. On the other hand, the β-correction is small for βTD(0) 1.
Stochasticity in gene networks may arise due to randomness of discrete, molecular reactions (intrinsic noise) and randomness associated with fluctuating environments (extrinsic noise) [14-16]; here we only address the former. Enhanced variability, e.g. due to transcriptional or translational bursting , is also not treated here but can be included readily.
Stochastic DF oscillations can be approximated by a rapid burst of repressor with random magnitude P that is subsequently degraded through a zeroth-order reaction with rate γr (again, we assume β = 0 and R0 → 0). The random time for degradation to zero of P repressor molecules satisfies a gamma distribution. Allowing the time to complete a burst of repressor (magnitude P) to be sharply-defined (e.g. 3τ + ton − t0), the period variance can be shown to consist of two parts (see SI)
where P and ΔP2 denote mean and variance of the burst amplitude. The first term in Eq. (5) is due to production variability and the second term is due to stochasticity of the degradation process.
Analogous to the deterministic case, the statistics of the burst amplitude can be found by integrating birth and death events for r(t) from time ton + τ to some later time, e.g. ton + 3τ, assuming that the statistics of trajectories for earlier times are known. Assuming that the probability of r = 0 for t > ton + τ is small, the equations for the evolution of the mean repressor concentration r(t) and the correlation k (r(t2), r(t1)) read
with ti + τ ≥ t2 ≥ t1 ≥ ti, υ(r) = F(r) − γr, B(r) = F(r) + γr, and k(x1,x2) x1x2 − x1x2 (see SI for details). The RHS of Eq. 7 comprises two parts. The integral over B(rτ) is a diffusive contribution to correlations due to uncorrelated production and degradation events [18, 19]. The remaining terms reflect variability due to the amplification of the randomness in past repressor trajectories in the course of new repressor synthesis.
For very small ron, the production rate is effectively binary: it is zero when rτ > 0 and α when rτ = 0. In this case the amplified variability is negligible, and the diffusive fluctuations are dominant. It is easy to see in this case that P obeys Poisson statistics (assuming α + γr ≈ α), i.e. ΔP2 ≈ P, such that . This Poisson-like period variability becomes relatively small for large P and is independent of the details of how the pulse of repressor rises. However, for finite ron, the amplified history fluctuations can lead to a significant period variability even when P is very large. In fact, simulations show that in many cases the variability is much larger than the Poisson contribution (see Fig. 3b,d).
Equations (6), (7) can be solved perturbatively by writing r(t) = r(t) + δr(t), with δr(t) being a small stochastic correction with zero mean. This approximation is used to investigate oscillations in the limit γrτ ron . After substituting Taylor expansions of v(r) and B(r) in Eqs. (6), (7), we can compute correlation functions and variances sequentially based on the knowledge of correlation functions and variances at earlier times (see SI). These analytical results are compared to numerical results in Fig. 3b,d. As seen from the figure, at sufficiently large γr, the analytical theory agrees well with numerics. Furthermore at large γr, the main contribution to the period variability is provided by the Poisson-like fluctuations. However, for moderate and small γr (or for large α), the main contribution to period variability is provided by amplified history fluctuations. Note that the period variance has a maximum near the value of γr where repressor first fails to consistently reach zero.
The genetic oscillator built in  featured two coupled feedback loops, one negative and another positive. The two proteins are encoded by two genes which are under the control of identical hybrid promoters. These promoters are activated by one of the proteins (activator) and repressed by the other (repressor). It was demonstrated experimentally that the dual feedback loop design provided enhanced robustness and tunability of oscillations as compared with a NFB-only design. To address this issue theoretically, we extend our delayed model system to include two additional reactions of production and degradation of activator a. Following the experiment, we assume that transcription of r and a genes is controlled by identical promoters, but the delay in producing mature activator proteins τa is different from τ, so repressor and activator production rates are, respectively, with F(r, a) α(f−1 + a/C1)/[(1 + a/C1)(1 + r/C0)2]. Here k is the ratio of the production rates for activator and repressor, and f is the relative activation strength (f = 1 corresponds to the NFB-only circuit). Repressor and activator degradation occur independently with zero-order rates γr and γa, respectively. In the deterministic limit, the dual feedback system is described by a pair of delay-differential equations (see SI).
In the DF limit, the inclusion of the PFB does not drastically alter the basic staging of phases in NFB oscillations in Fig. 2b. PFB with high f and low τa can produce robust and large amplitude oscillations that maintain a relatively low coefficient of variation (C.V.), i.e. with C.V. similar to the NFB-only case with basal rate α/f (Fig. 4). Oscillation amplitude and period tend to decrease as the activator delay becomes comparable to repressor delay, since late-arriving activator allows repressor to be produced at the low basal rate α/f during most of the production phase, but the C.V. can be remarkably insensitive to τa. However, for other choices of parameters, the presence of PFB can actually lead to the increased C.V. as compared to the NFB-only case (see SI).
Synthetic gene circuits designed around a core delayed NFB circuit offer a promising strategy for generating robust genetic oscillations. Here we have presented an analysis of a simple model for such oscillators. A NFB-only system with strong tightly repressible promoter and slow degradation can produce degrade-and-fire oscillations with a period much larger than the delay time and with relatively small period variability. The main source of oscillation variability is the short protein production phase when the number of repressor molecules is low. An additional fast PFB loop in this model was found to be able to increase the mean period and maintain low C.V. compared to an analogous NFB-only system. This result is similar to recent findings regarding the role of positive feedback in biological oscillations . However, it should be noted that while other coupled positive and negative feedback biological oscillator models [2, 22] rely on a separation of time scales between the two components to create relaxation oscillations, here the mechanism is based on the presence of delay in the negative feedback loop and can function even in the absence of PFB. Note that a different mechanism of long-period oscillations in cell dynamics have been recently investigated in .
We believe that the results presented here are also relevant for many naturally occurring gene oscillators, such as circadian clocks and Hes1 system [23-25]. While they are generally made up of many interdependent genes and proteins, one can often identify a core NFB loop similar to Fig. 1.
We thank M. Mackey for illuminating discussions. This work was supported by NIH grants GM69811-01 and GM082168-02.