PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Rev Lett. Author manuscript; available in PMC 2010 August 20.
Published in final edited form as:
PMCID: PMC2924583
NIHMSID: NIHMS225830

Delay-induced Degrade-and-Fire Oscillations in Small Genetic Circuits

Abstract

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 [1]. The design of the oscillator was motivated by the original theoretical proposal [2] 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 [1] 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 [9]. 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.

FIG. 1
An auto-repressor genetic oscillator described in Ref. [1]. Transcriptional activity of the gene is regulated by the concentration of the repressor protein through several sequential kinetic steps. In this work, these steps are replaced by an explicitly ...

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 [1]. Details of these calculations can be found in the Supplementary Information (SI) [10].

NFB-only system

We consider a system with only two reactions: production of a repressor protein and its degradation. These reactions are characterized by the rates K+(r) and K(r) 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, K+(r)(t)=F(rτ(t)), where r(t) is the number of repressor molecules at time t, rτ(t) [equivalent] r(t − τ), F(r)αC02/(C0+r)2 is a Hill function governing NFB [11]. Repressor is degraded enzymatically at a maximum rate γr, and we also assume a small first-order degradation due to dilution with rate β, K(r)(t)=γrr(t)[R0+r(t)]1+βr(t) with R0 characterizing the concentration of the enzyme.

We simulated these two stochastic reactions using the extension of the Gillespie algorithm [12] proposed in [7] to treat delayed reactions. For sufficiently large α [greater, similar] γr [dbl greater-than sign] C0/τ > R0/τ [13] 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).

FIG. 2
(a) DF oscillations in the stochastic NFB-only system, with parameters α = 300, γr = 80, C0 = 10, τ = 1, β = 0.1, R0 = 1. (b) A magnified view of the deterministic trajectory near the production phase in the limit β ...

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.

Deterministic analysis

In the deterministic limit, the mass-action kinetics of delayed production and degradation can be expressed by a delay-differential equation

r˙(t)=αC02(C0+rτ)2γrr(t)R0+r(t)βr(t)
(1)

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,

r(t)γr(tt0),t<t0
(2)

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

r(t)r(ti)=tiτtτdtα(1+r(t)C0)2γr(tti)
(3)

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

P0(αγr)[τC0γr(αγr1)]
(4)

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).

FIG. 3
(a) Stochastic mean period (diamonds), deterministic period (dash-dot line) and two analytic approximations, T0 = 2τ + P0/γr (Eq. (4), dashed line) and an analytic estimate based on Eq. (2) (solid line), vs. γr for the NFB-only ...

The case of nonzero but small β (βτ [double less-than sign] 1) can also be analyzed. Production bursts of repressor are only slightly perturbed in this case. However, if the number of produced proteins P [dbl greater-than sign] γ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) [double less-than sign] 1.

Stochastic analysis

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 [17], 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τ + tont0), the period variance can be shown to consist of two parts (see SI)

ΔT2=ΔP2γr2+Pγr2
(5)

where left angle bracketPright angle bracket and left angle bracketΔP2right angle bracket 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 left angle bracketr(t)right angle bracket and the correlation k (r(t2), r(t1)) read

r(t)r(ti)=titdtυ(rτ(t))
(6)

k(r(t2),r(t1))k(r(ti),r(ti))=tit1dtB(rτ(t))+tit2dt2tit1dt1k(υ(rτ(t2)),υ(rτ(t1)))+tit1dtk(r(ti),υ(rτ(t)))+tit2dtk(r(ti),υ(rτ(t)))
(7)

with ti + τt2t1ti, υ(r) = F(r)γr, B(r) = F(r) + γr, and k(x1,x2) [equivalent] left angle bracketx1x2right angle bracketleft angle bracketx1right angle bracketleft angle bracketx2right angle bracket (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. left angle bracketΔP2right angle bracketleft angle bracketPright angle bracket, such that ΔT22P/γr2. 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) = left angle bracketr(t)right angle bracket + δr(t), with δr(t) being a small stochastic correction with zero mean. This approximation is used to investigate oscillations in the limit γrτ [dbl greater-than sign] ron [20]. 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.

Coupled positive-negative feedback oscillator

The genetic oscillator built in [1] 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, K+r=F(rτ,aτ),K+a=kF(rτa,aτa) with F(r, a) [equivalent] α(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).

FIG. 4
Mean period (a) and coefficient of variation (b) of the stochastic coupled positive-negative feedback system vs. 1/f, with different values of τa (labeled) and fixed α averaged over 5000 samples. These results are compared with the NFB-only ...

Closing remarks

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 [21]. 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 [8].

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.

Supplementary Material

Supp

Acknowledgments

We thank M. Mackey for illuminating discussions. This work was supported by NIH grants GM69811-01 and GM082168-02.

References

1. Stricker J, Cookson S, Bennett M, Mather W, Tsimring L, Hasty J. Nature. 2008;456:516. [PubMed]
2. Hasty J, Dolnik M, Rottschafer V, Collins JJ. Phys Rev Lett. 2002;88:148101. [PubMed]
3. MacDonald N, Theor J. Biol. 1977;67:549. [PubMed]
4. Bliss RD, Painter PR, Marr AG, Theor J. Biol. 1982;97:177. [PubMed]
5. Mackey MC, Nechaeva IG. Phys Rev E. 1995;52:3366. [PubMed]
6. Monk N. Curr Biol. 2003;13:1409. [PubMed]
7. Bratsun D, Volfson D, Tsimring LS, Hasty J. Proc Natl Acad Sci USA. 2005;102:14593. [PubMed]
8. Pujo-Menjouet L, Bernard S, Mackey M, Hasty J. SIAM J Appl Dyn Syst. 2005;4:312.
9. Knight BW. J Gen Phys. 1972;59:734. [PMC free article] [PubMed]
10. Supplementary Information supplied by EPAPS
11. While nonlinearity of the delayed negative feedback loop is important, the qualitative features of protein oscillations are weakly sensitive to the specific form of the Hill function, see SI for more details
12. Gillespie DT. J Phys Chem. 1977;81:2340.
13. For very large α [greater, similar] γr((γrτ/C0)+1)2 there exists another regime in which the repressor levels fail to reach zero, the amplitude of oscillations becomes small and the period is approximately 4τ
14. Elowitz MB, Levine AJ, Siggia ED, Swain PS. Science. 2002;297:1183. [PubMed]
15. Ozbudak E, Thattai M, Kurtser I, Grossman A, van Oudenaarden A. Nature Gen. 2002;31:69. [PubMed]
16. Raser J, O’Shea E. Science. 2005;309:2010. [PMC free article] [PubMed]
17. Golding I, Paulsson J, Zawilski SM, Cox EC. Cell. 2005;123:1025. [PubMed]
18. Gardiner C. Handbook of Stochastic Methods. Springer-Verlag; New York: 2004.
19. Kepler TB, Elston TC. Biophys J. 2001;81:3116. [PubMed]
20. This approximation is not valid at r = 0, however this region mostly contributes to the Poisson-like variability and can be treated separately
21. Tsai TY, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE. Science. 2008;321:126. [PMC free article] [PubMed]
22. Wang J, Xu L, Wang E. Proc Natl Acad Sci USA. 2008;105:12271. [PubMed]
23. Barkai N, Leibler S. Nature. 2000;403:267. [PubMed]
24. Goldbeter A. Nature. 2002;420:238. [PubMed]
25. Bernard S, Čajavec B, Pujo-Menjouet L, Mackey M, Herzel H. Phil Trans Roy Soc A. 2006;364:1155. [PubMed]