Article Outline

Article Information

PubMed

Related Articles

  • …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 10, 3407-3424, 15 May 2007

doi:10.1529/biophysj.106.096891

Biophysical Theory and Modeling

Dynamics of a Minimal Model of Interlocked Positive and Negative Feedback Loops of Transcriptional Regulation by cAMP-Response Element Binding Proteins

Hao Song, Paul Smolen, Evyatar Av-Ron, Douglas A. Baxter and John H. ByrneGo To Corresponding Author,  

Department of Neurobiology and Anatomy, W. M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas Medical School at Houston, Houston, Texas 77225

Address reprint requests to John H. Byrne, Dept. of Neurobiology and Anatomy, W. M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas Medical School at Houston, PO Box 20708, Houston, TX 77225. Tel.: 713-500-5602; Fax: 713-500-0623.

Abstract

cAMP-response element binding (CREB) proteins are involved in transcriptional regulation in a number of cellular processes (e.g., neural plasticity and circadian rhythms). The CREB family contains activators and repressors that may interact through positive and negative feedback loops. These loops can be generated by auto- and cross-regulation of expression of CREB proteins, via CRE elements in or near their genes. Experiments suggest that such feedback loops may operate in several systems (e.g., Aplysia and rat). To understand the functional implications of such feedback loops, which are interlocked via cross-regulation of transcription, a minimal model with a positive and negative loop was developed and investigated using bifurcation analysis. Bifurcation analysis revealed diverse nonlinear dynamics (e.g., bistability and oscillations). The stability of steady states or oscillations could be changed by time delays in the synthesis of the activator (CREB1) or the repressor (CREB2). Investigation of stochastic fluctuations due to small numbers of molecules of CREB1 and CREB2 revealed a bimodal distribution of CREB molecules in the bistability region. The robustness of the stable HIGH and LOW states of CREB expression to stochastic noise differs, and a critical number of molecules was required to sustain the HIGH state for days or longer. Increasing positive feedback or decreasing negative feedback also increased the lifetime of the HIGH state, and persistence of this state may correlate with long-term memory formation. A critical number of molecules was also required to sustain robust oscillations of CREB expression. If a steady state was near a deterministic Hopf bifurcation point, stochastic resonance could induce oscillations. This comparative analysis of deterministic and stochastic dynamics not only provides insights into the possible dynamics of CREB regulatory motifs, but also demonstrates a framework for understanding other regulatory processes with similar network architecture.

Introduction

The cyclic AMP (cAMP)-response element-binding protein (CREB) family of transcription factors is crucial for a variety of cellular processes 1, including induction of neuronal plasticity and formation of long-term memory (LTM), in both invertebrates and vertebrates 2,3,4,5, circadian rhythm 6,7,8, neuronal differentiation 9, addiction and depression 10, hormonal control of metabolic processes 11, and spermatogenesis 12. Thus, understanding the regulatory properties and dynamics of CREB is important for understanding many key biological functions.

The involvement and regulation of CREBs in the induction of LTM have been extensively studied 4,13. For example, in the mollusc Aplysia, two creb genes (creb1 and creb2) and their corresponding proteins have been characterized 14,15,16. CREB1 is a transcriptional activator necessary for induction of long-term synaptic facilitation, a mechanism for LTM. Upon exposure of sensory neurons to the neurotransmitter serotonin (5-HT), creb1 is activated 13,15 via activation of the protein kinase A (PKA) intracellular signaling pathways 4. CREB2 functions as a transcriptional repressor that poses inhibitory constraints on the induction and formation of LTM 16,17. Upon exposure to 5-HT, repression by CREB2 is relieved, possibly via phosphorylation of CREB2 by mitogen-activated protein kinase (MAPK) 16. CREB proteins regulate gene expression by binding to enhancer sequences termed cAMP-response elements (CREs). creb genes themselves may be regulated by CREs. For example, Mohamed et al. 14 recently characterized CREs in the promoter regions of Aplysia creb1 and creb2, and the mammalian CREB activator gene has CREs in its promoter region 18. The existence of these CREs suggests that a regulatory motif with interlocked positive and negative feedback loops may characterize transcriptional regulation by CREBs. In this motif, an activator such as CREB1 would further activate expression of its own gene via binding to CREs, forming a positive feedback loop. The activator would also activate expression of a repressor such as creb2. The repressor would, in turn, repress expression of the activator, closing a negative feedback loop. The repressor would also repress its own gene. In primary rat Sertoli cells 12, a transactivator form of CREB (similar to the activator CREB1) and a repressor isoform of CREB (similar to CREB2) are suggested to interact via such interlocked positive and negative loops, generating oscillations in gene expressions. This motif is illustrated in Fig. 1, where for definiteness, we term the activator CREB1 and the repressor CREB2.

Display large version of this figure
Figure 1
Schematic of the minimal model with interlocked positive and negative feedback loops of activator (CREB1) and repressor (CREB2).

The potential regulatory dynamics of CREB feedback have not been characterized. To explore these dynamics, we developed a minimal mathematical model to represent the feedback loops. Model dynamics can suggest experimentally testable hypotheses. For example, a prediction of the effect of increasing (or decreasing) the strength of positive (or negative) feedback might assist in designing protocols to regulate the induction and expression of CREB activators and repressors. Induction or activation of CREBs correlates with long-term memory formation not only in Aplysia, but also in mammals 19,20,21,22. Therefore, modulating CREB expression is likely to regulate the formation and persistence of LTM and other cellular processes mentioned above 1,2,3,4,5,6,7,8,9,10,11,12. In the minimal model, CREB1 and CREB2, respectively, denote activator and repressor CREB proteins.

Given that the actual regulatory dynamics of signaling pathways and gene network of CREBs have not yet been fully characterized, a minimal model representing the essential nonlinear network architecture of CREBs has several advantages. 1), The minimal model has relatively few parameters, which allows a systematic analysis of model dynamics and their robustness to parameter variation. 2), The minimal model is a core mathematical structure that can be extended to more detailed and biologically faithful models by incorporating other regulatory pathways once their kinetics are characterized. 3), Because the minimal model captures the essential nonlinearity embedded in the network structure, the dynamics may not change greatly upon the addition of other upstream and downstream regulatory pathways, and 4), Similar dynamics may characterize cellular processes based on similar regulatory motifs.

A number of minimal (or reduced) models have been developed to describe regulation of gene expression associated with the cell cycle, circadian rhythms, the tumor suppressor protein p53, and synthetic gene networks 23,24,25,26,27,28. Among these, only the model developed by Smolen et al. 25 has described positive and negative feedback loops of two transcription factors similar to Fig. 1. However, that model did not represent dimerization of the repressor, as occurs in the CREB family. That study also did not include bifurcation analysis to examine how dynamics depend on system parameters.

We focused on four essential questions about the dynamic properties of the interlocked feedback loops of CREBs: 1), Are bistability and oscillations possible? 2), How does the strength of positive and negative feedback affect the system dynamics? 3), How could time delays inherent in the regulation of expression of activator and repressor CREB proteins impact the dynamics? 4), How do stochastic fluctuations due to small numbers of molecules of CREBs impact the dynamics? Bifurcation analysis of the minimal model illustrates that both bistability and oscillations are possible given proper strengths of positive and negative feedback. Bistability and oscillations are fairly robust against variations of parameters. Time delays in regulation of CREB expression can alter a stable steady state to a stable limit cycle, or to chaotic dynamics. Stochastic simulations illustrate that a bimodal distribution of the numbers of CREB molecules characterizes the dynamics in the bistable region. The relative robustness of the HIGH and LOW stable states to molecular noise depends on the location of the unstable steady state between the HIGH and LOW states. A critical number of molecules is needed to sustain HIGH CREB expression for a day or longer, which may correlate with consolidation of LTM. A critical number of CREB molecules is also required to sustain robust oscillations. Stochastic resonance is also found. An optimal noise level induces a stochastic oscillation with the most reliable periodicity.

This comparative study of the deterministic bifurcation analysis and stochastic simulations also illustrates a framework to understand other cellular regulatory processes with similar network architectures. Examples include circadian rhythms, the cell cycle, and artificial gene networks.


Model development and computational methods

Minimal model development

The minimal model represents regulation of CREB1 and CREB2 expression by binding of these proteins to CREs near their own genes, degradation of CREB1 and CREB2, and basal synthesis of CREB1 and CREB2 (Fig. 1). A number of simplifications are made. The minimal model does not explicitly describe the translation of creb mRNAs, the phosphorylation of CREBs by kinases, the signaling pathways regulating phosphorylation, or the translocation of CREBs between nucleus and cytoplasm. A homogeneous, single-compartment cellular structure is assumed. The omitted steps are lumped into a single effective synthesis process of CREBs, as illustrated in Fig. 1.

These simplifications lead to a minimal model with two deterministic ordinary differential equations (ODEs):

(1)
(2)

In Eq. (1), the first term of the right-hand side represents the autoactivation of creb1 expression by CREB1 homodimer binding to a CRE. The dependence on CREB2 in the denominator represents repression of creb1 expression by competitive binding of CREB2 homodimer to the same CRE (for further details, see Supplementary Material ). The production rate of CREB1 saturates hyperbolically with CREB1 homodimer concentration, and decreases hyperbolically with CREB2 homodimer concentration. The Hill coefficients of 2 for [CREB1] and [CREB2] represent, in the simplest plausible manner, the requirement for two CREB1 or CREB2 monomers to form homodimers. By mass action, homodimer concentration should be more or less proportional to the square of monomer concentration. Equation (2), analogously, represents competitive binding of CREB1 and CREB2 homodimers to the CRE regulating CREB2 expression. These equations are identical to a previously published minimal model for competitive binding of two transcription factors 25, except that the Hill coefficient of CREB2 is 2 here (Eqs. (1)) and a basal synthesis rate is incorporated for CREB2.

In Eqs. (1), Vx and Vy are the maximum induced synthesis rates of CREB1 and CREB2. Kx and Ky are, respectively, the dissociation constants of CREB1 and CREB2 homodimers from the CRE. The second term in the right-hand side of Eq. (1) represents a first-order degradation of CREB1, where kdx is the degradation rate constant. A small basal synthesis rate of CREB1, rbas,x, is also assumed. In Eq. (2), kdy and rbas,y, have analogous meanings.

The purpose of the minimal model is to investigate the qualitative nonlinear dynamics of the interlocked feedback loops. Currently, there is little data to quantitatively determine the temporal dynamics of CREBs. The parameters in the model (Eqs. (1)) are assigned a set of basal values that have physiological relevance to the characteristic timescales of LTM and that lie within a physiologically feasible range. The basal values of the parameters in Eqs. (1) are Vx=0.4min−1, Vy=0.01min−1, Kx=5, Ky=10, kdx=0.04min−1, kdy=0.01min−1, rbas,x=0.003min−1, and rbas,y=0.002min−1. These basal parameter values were chosen to allow qualitative simulation of CREB1 phosphorylation dynamics observed in 5-HT stimulus protocols (see Fig. 6). Concentrations of CREB1 and CREB2 are expressed in arbitrary units, because the actual protein concentrations of CREBs in neurons are not well defined.

To investigate the robustness of the model dynamics to a change in the binding mechanisms of CREBs to CREs, a model variant was developed. This variant assumes that expression of CREB1 is regulated by noncompetitive binding of CREB1 and CREB2 to different CREs, and similarly for expression of CREB2. This mechanism is represented by another pair of ODEs (Eqs. (3)) (see Supplementary Material ):

(3)
(4)

Parameters in Eqs. (3) have the same meanings and basal values as given for Eqs. (1).


Delayed differential equations (DDEs)

Gene regulatory systems have ubiquitous time delays associated with translocation of mRNA and transcription factors between cytoplasm and nucleus and, thereby, transcription of genes 29,30,31,32. To explore the effects of time delays on the model of Eqs. (1), discrete delays for the synthesis of CREB1, τCREB1, and CREB2, τCREB2, were incorporated into terms that describe the synthesis of CREBs. The resulting DDEs are:

(5)
(6)

These DDEs were simulated as follows, the values of [CREB1] and [CREB2] at any given time t are stored and used after delays τCREB1 and τCREB2 to calculate the derivatives of [CREB1] and [CREB2], respectively.


Stochastic format of ODE model

The numbers of CREB1 and CREB2 molecules are likely to be in the range of tens to hundreds in a single neuron 33. Consequently, we investigated the stochastic effects of the finite number of molecules on the system dynamics. The Gillespie stochastic modeling algorithm 30,34 was implemented to simulate a stochastic version of the minimal model (Eqs. (1)). The representative reactions and their respective propensity functions are listed in Table 1.

Table 1 Stochastic version of the minimal model
ReactionPropensity function
(7)
(8)
(9)
(10)
(11)
(12)
NCREB1 and NCREB2 in Eqs. (7) represent the number of molecules of CREB1 and CREB2. In the second column, the propensities of occurrence of each reaction step are derived from kinetic terms in Eqs. (1), and are scaled by multiplication of the system volumetric parameter Ω. Eqs. (7) describe the production and degradation propensity function, respectively, of CREB1. Eq. (9) describes the basal synthesis propensity of CREB1. Eqs. (10), for CREB2, are analogous to Eqs. (7). Parameters in Eqs. (7) are the same as for the deterministic ODEs (Eqs. (1)), and have the same basal values.

Reaction probabilites, or propensities for the synthesis of CREB1 and CREB2, are modeled by use of the Hill-type rate expression in the Gillespie algorithm. An alternative approach would be to decompose the synthesis of CREB1/2 into a series of elementary reactions (unimolecular or bimolecular) and use only rates for such elementary reactions in the Gillespie algorithm. We chose the former approach, for two reasons. First, recent comparisons between both approaches in simulations of gene regulation and protein phosphorylation underlying circadian rhythms found that composite Michaelis-Menten type rate expressions gave very similar results to elementary reaction rates 35,36,37. Second, retaining the basic form of the minimal model in both cases allows for a direct comparison of the minimal model dynamics predicated by deterministic and stochastic simulations. This approach also avoided substantially increasing the number of reactions and differential equations.


Computational methods

For deterministic simulation, ODEs were numerically integrated using the fourth-order Runge-Kutta algorithm in the numerical analysis package XPPAUT 38. Numerical bifurcation analyses of the ODEs were performed with XPPAUT. DDEs were numerically integrated by three software tools: DDE solver in XPP, dde23 in Matlab 39, and a program implementing forward Euler algorithm for DDEs. Differences of computational results with these DDE integrators were not significant.

For stochastic simulation, the Gillespie algorithm was programmed in FORTRAN. The random number generator used was from the IMSL statistics library. Power spectrum density (PSD) analysis to characterize the spectrum (or period) features of stochastic time series was performed in Matlab.



Results

Varying the strength of negative feedback generates four different bifurcation diagrams, in which bistability and oscillations are observed

Fig. 1 represents the minimal model. In the positive feedback loop, binding of the CREB1 homodimer to CRE further activates the production of CREB1. In the negative feedback loop, CREB2 repressor homodimers bind to CRE, inhibiting the production of CREB2. In addition, CREB1 activates CREB2 expression; whereas CREB2 represses CREB1 expression. Therefore, the positive feedback loop of CREB1 and negative feedback loop of CREB2 are interlocked.

We proceeded to use Eqs. (1) to investigate how dynamics vary with the strength of positive feedback by CREB1 for different strengths of negative feedback by CREB2. Varying Vx, the strength of autoactivation by CREB1 of CREB1 synthesis, alters the strength of positive feedback. Therefore, we depicted the dynamics with one-parameter bifurcation diagrams of [CREB1] versus Vx. Figure 2AD, illustrate four different bifurcation diagrams, each at a different fixed value of Vy, the production rate of repressor CREB2. Each diagram therefore corresponds to a different strength of the negative feedback in which CREB2 represses synthesis of both CREBs.

Display large version of this figure
Figure 2
Bifurcation diagrams of [CREB1] versus Vx at four values of Vy and graphs depicting period of oscillation as a function of Vx. Solid lines depict stable steady states, dashed lines depict unstable steady states; solid (open) circles denote maximum and minimum values of [CREB1] on stable (unstable) limit-cycles. Codimension-1 singular points are marked as LP, limit point (or saddle-node point); sub-H, subcritical Hopf bifurcation point; sup-H, supercritical Hopf point; HC, homoclinic point; SNIC, saddle-node-invariant-circle bifurcation; and LP-PO, limit point on periodic orbits. (A) Bifurcation diagram at Vy=0.01min−1. (B) Bifurcation diagram at Vy=0.4min−1. (C1) Bifurcation diagram at Vy=0.8min−1. (C2) Period of oscillations at Vy=0.8min−1. (D1) Bifurcation diagram at Vy=2min−1. (D2) Period of oscillations at Vy=2min−1. Other parameters have basal values: Kx=5; Ky=10; kdx=0.04min−1; kdy=0.01min−1; rbas, x=0.003min−1; and rbas, y=0.002min−1.

In Figure 2A, for a very weak negative feedback strength (Vy=0.01min−1), bistable steady states of [CREB1] arise in the range of Vx bounded by two limit points (LP1 and LP2). Outside of this range only a single steady state of [CREB1] exists. Upon increasing Vy to 0.4min−1 (Figure 2B), oscillations bifurcate from the upper steady state of the bistability region at a subcritical Hopf (sub-H) point (Figure 2B, oscillation amplitudes marked by open circles). These unstable oscillations coalesce with the unstable middle branch of steady states connected by LP1 and LP2 (Figure 2B, dashed line), generating a homoclinic (HC) bifurcation point. The upper steady-state solution is unstable between LP1 and sub-H, but becomes stable beyond sub-H for larger values of Vx. Upon further increasing Vy to 0.8min−1 (Figure 2C), the sub-H point moves beyond the bistability range of Vx bounded by LP1 and LP2 (Figure 2C1). Now, as sub-H is crossed from lower to higher values of Vx, a stable limit-cycle oscillation (solid circles) is replaced, in a narrow region, by the coexistence of a stable limit-cycle oscillation and a stable steady state. At the limit point of periodic orbits (LP-PO), stable and unstable limit-cycle oscillations coalesce. Decreasing Vx from LP-PO, the amplitude of oscillation decreases, but the period lengthens (Figure 2C2). When Vx decreases to LP2, the limit-cycle oscillation disappears at a saddle-node-invariant-circle (SNIC) bifurcation point, and is replaced by a stable steady state. The period of oscillation at SNIC tends to infinity (Figure 2C).

As Vy is increased to 2.0min−1 (Figure 2D), the two limit points (LP1 and LP2) both move together and coalesce at a SNIC. Consequently, bistability disappears. Oscillations are preserved, but now arise on the left from a supercritical Hopf bifurcation point (sup-H). The sub-H in Figure 2D1 is similar to that in Figure 2C1. The main difference between sub-H and sup-H is that upon crossing sup-H, an oscillation with small amplitude arises from a stable steady state, whereas upon crossing sub-H, an oscillation with large amplitude abruptly arises from a stable steady state. Figure 2D2 illustrates that the period of limit-cycle oscillation decreases from sup-H to sub-H. Thus, the interlocked feedback loops of CREBs might support limit-cycle oscillation with a broad range of periods, e.g., from hours to days (Figure 2C2D2). Alternatively, by scaling all the parameters containing units of time, any biologically reasonable oscillation period can be obtained.

How many qualitatively different one-parameter bifurcation diagrams of [CREB1] versus Vx may be attained at different Vy? This question can be answered by inspection of the two-parameter bifurcation diagram in the (Vx, Vy)-parameter plane (Figure 3A), which is constructed by continuation of the loci of five different types of codimension-1 singular points identified in Figure 2AD, namely, LP (LP1, LP2), Hopf (sub-H, sup-H), SNIC, HC, and LP-PO. XPPAUT was used to numerically determine and continue the loci of these singular points. Four distinct bifurcation diagrams of [CREB1] versus Vx are observed. The typical value of Vy for each case is marked A–D in Figure 3A. These values correspond to the one-parameter bifurcation diagrams of Figure 2AD, respectively. The loci of SNIC and LP2 are overlapping in Figure 3A because SNIC and LP2 occur at identical Vx (Figure 2C1). The loci of HC and LP-PO are close to the sub-H in Figure 3A. Two codimension-2 singular points, a cusp point and a Bogdanov-Takens (BT) bifurcation point, are identified in Figure 3A. At the cusp point (Vx, Vy)=(0.736min−1, 1.627min−1), the LP1 and LP2 limit points coalesce. Therefore, upon crossing the cusp point, three steady states appear from a single steady state. The BT bifurcation, at (Vx, Vy)=(0.218min−1, 0.057min−1), is caused by coalescence of a limit point (here LP1) and a Hopf bifurcation point (sub-H). Upon crossing BT, a Hopf bifurcation (sub-H) appears, corresponding to the appearance of an oscillating solution.

Display large version of this figure
Figure 3
(Vx, Vy) bifurcation diagram with corresponding distinct phase diagrams of [CREB1] versus [CREB2]. (A) Bifurcation diagram with loci of six codimension-1 singular points: LP, sub-H, sup-H, SNIC, HC, and LP-PO. These loci delineate seven regions, denoted by I, II, III, IVa, IVb, Va, and Vb. Two codimension-2 singular points exist: a cusp point and a Bogdanov-Takens (BT) bifurcation point. Values of Vy denoted by A–D correspond to the Vy values in Figure 2AD, respectively. Basal parameter values are otherwise used. (B) Phase diagrams corresponding to the labeled regions in A. Stable (unstable) steady states are denoted by solid (open) circles, stable (unstable) limit-cycle oscillations by solid (dashed) lines.

Further inspection of Figure 3A illustrates that the (Vx, Vy)-parameter plane is divided into seven regions by the loci of different bifurcation points, as depicted by Regions I–III, IVa, IVb, Va, and Vb. The seven corresponding phase diagrams of the dynamic solutions for [CREB2] and [CREB1] are illustrated in Figure 3B. In Region I of Figure 3A, only a single steady state with low concentrations of CREB1 and CREB2 exists (Figure 3B,I). In Region II of Figure 3A, similar to Region I, only a single steady state may exist, but with higher concentrations of CREB1 and CREB2. This is due to the stronger positive feedback (Vx) in Region II than in region I as well as the weaker negative feedback (Vy) in Region II. In Region-III, bounded by LP1, LP2, and sub-H, three steady states arise, but only the lower steady state is stable to small perturbations in [CREB1] or [CREB2] (Figure 3B, III). In Region IV (IVa and IVb) of Figure 3A, bounded by LP1, LP2, sub-H, and horizontal axis of Vx, three steady states exist (Figure 3B, IVa and IVb). Bistability exists, in which the upper and lower steady states are stable to small perturbations, but the middle steady state is unstable. The only difference between the phase diagrams depicted in Figure 3B, IVa and IVb, is that the upper state in the phase diagram of IVb is surrounded by an unstable limit cycle. In the phase diagram of Figure 3B, Va, a stable limit cycle surrounds an unstable steady state. In Figure 3B, Vb, a stable limit cycle surrounds an unstable limit cycle, which further surrounds a stable steady state.

In summary, by manipulating the strength of the positive and negative feedbacks, three classes of system dynamics may be attained: a single stable steady state (Figure 3B, I–III), bistable steady states (Figure 3B, IVa and IVb), and limit-cycle oscillations (Figure 3B, Va and Vb). The relatively large regions of bistability (Region IV) and oscillatory dynamics (Region V) in Fig. 3 suggest that parameters with values in these regions may have some physiological significance in regulating CREB dynamics.


Bistability and oscillations in the minimal model are robust to variations in parameters and binding mechanisms

One way to investigate the robustness (or sensitivity) of system dynamics is to vary system parameters and observe changes in the size and location of the parameter regions where particular dynamics exist. The small number of parameters in the minimal model allows a systematic analysis. Fig. 4 displays the two-parameter bifurcation diagram in the (Vx, Vy) plane after doubling the values of the six parameters in Eqs. (1). Since homoclinic bifurcations (HC) and limit points on periodic orbits (LP-PO) are very close to the subcritical Hopf (sub-H) point in these parameter ranges, only the Hopf points (sub-H and sup-H) and the limit points (LP1 and LP2) are plotted.

Display large version of this figure
Figure 4
Impact of doubling the values of parameters on the two-parameter bifurcation diagram in the (Vx, Vy) plane. Solid lines depict loci of limit points (LP1, LP2), dashed lines depict loci of Hopf bifurcation points (sup-H, sub-H). Thin lines correspond to parameters with basal values, thick lines correspond to doubling a parameter value. (A) Kx is doubled. (B) Ky is doubled. (C) kdx. (D) kdy. (E) rbas,x. (F) rbas,y.

Figure 4AB, depicts the result of doubling the dissociation constants of CREB1 and CREB2 from CREs (Kx and Ky, respectively). Doubling Kx versus doubling Ky gives opposite effects on the location of the loci of limit points (LP1 and LP2) and Hopf bifurcation points (sub-H and sup-H), with doubling of Kx moving LPs and Hopf points to larger values of Vx and Vy. This shift is produced because an increase in Kx increases the dissociation of CREB1 from CREs, thereby decreasing the production of CREB1 and CREB2, so that larger values of Vx and Vy are needed to compensate. In contrast, an increase in Ky increases dissociation of CREB2 from CREs, thereby increasing the production of CREB1 and CREB2. Figure 4CD, illustrates the result of doubling the degradation rate constants of CREB1 and CREB2 (kdx and kdy). Increases in these parameters also have opposite effects on the location of LPs and Hopf points. One significant effect of increasing kdy is an elimination of oscillatory dynamics (Figure 4D) because at large kdy, rapid degradation of CREB2 significantly weakens the strength of the negative feedback loop that is required for oscillations. Figure 4EF, display the location of LPs and Hopf points when doubling the basal synthesis rates of CREB1 and CREB2 (rbas,x and rbax,y). Doubling rbas,x significantly shrinks both regions of bistability and oscillation, whereas doubling rbas,y has little effect on these regions. Taken together, the bistability and oscillatory dynamics are robust to variations of system parameters.

We next examined whether bistability and oscillations are robust to a change in the mechanism of CREBs binding to CREs. The model of Eqs. (1) describes competitive binding of CREB1 and CREB2. The model variant of Eqs. (3) describes noncompetitive binding. Fig. 5 compares the regions of bistability and oscillation for the two binding mechanisms. Changing the mechanism has negligible impact on the limit points (LP1 and LP2), but significantly affects the location of the Hopf bifurcations (sub-H and sup-H). For noncompetitive binding, only sup-H bifurcations exist (thick sup-H curve), in contrast with competitive binding (thin sup-H and sub-H curves). Limit-cycle oscillations for noncompetitive binding occur only in the region bounded by the thicker sup-H curve and LP2, which is much smaller in size than the corresponding region for competitive binding (bounded by the thin curves sub-H, sup-H, and LP2). For noncompetitive binding, the sup-H bifurcation points coalesce and disappear at a degenerate Hopf point (DH) of codimension 2. Fig. 5 illustrates that for both binding mechanisms, the size of the region corresponding to bistability and oscillations is fairly large, indicating that the dynamics are relatively robust to a change in binding mechanism.

Display large version of this figure
Figure 5
Bifurcation diagrams in the (Vx, Vy) plane. Thin curves correspond to competitive binding (Eqs. (1)), thick curves correspond to noncompetitive binding (Eqs. (3)). Solid curves depict loci of LP, dashed lines depict loci of Hopf bifurcation points (sub-H, sup-H). Three codimension-2 singular points are identified as cusp point, Bogdanov-Takens (BT) bifurcation, and degenerate Hopf bifurcation (DH). The basal set of parameter values was used.

The minimal model can qualitatively simulate several experimental findings

Bistability of creb expression might help to explain several experimental findings concerning the dynamics of CREB1 and CREB2 in protocols that can lead to LTM formation in Aplysia. Empirically, five pulses of treatment with 5-HT induce long-term strengthening or facilitation of synapses between Aplysia sensory neurons and motor neurons, whereas three pulses do not 40. This long-term synaptic facilitation (LTF) correlates both with increased CREB expression 14 and with long-term memory formation. Empirically, injection of CREB1 into sensory neurons lowers the threshold for LTF induction from five 5-HT pulses to one pulse 15. In contrast, injection of CREB2 prevents LTF from being elicited by five pulses 16.

To phenomenologically simulate pulsed 5-HT application, the synthesis rate constant of CREB1 (Vx) in Eqs. (1) was briefly increased. Vx was increased from its basal value (0.4min−1) to 3.7min−1 to correspond to a [5-HT] elevation from 0 to 10μM. A bistable bifurcation diagram of [CREB1] and [CREB2] versus Vx was constructed using Eqs. (1) (Figure 6A). This diagram exhibits two stable high and low branches (solid lines) and an unstable middle branch (dashed line). An elevation in Vx can cause [CREB1] to transit from the LOW to the HIGH state. The HIGH state is stable when [5-HT]=0. After 5-HT is removed, [CREB1] will not return to a basal level, but will remain in the HIGH state. Such a bifurcation diagram is called an “irreversible switch”. Two thresholds exist in this bifurcation diagram. One threshold concerns the bifurcation parameter (Vx)—the limit point at Vx≈0.67min−1. The second threshold concerns the variable [CREB1]—the unstable middle steady state of [CREB1] (∼0.43) at [5-HT]=0.

Display large version of this figure
Figure 6
Bifurcation diagram of [CREB1] and [CREB2] versus Vx and time courses of [CREB1]. (A) Diagram depicting an irreversible switch of [CREB1] (thick lines) and [CREB2] (thin lines) versus Vx. (B) Five pulses of 5-HT induce persistent [CREB1] elevation; whereas three pulses of 5-HT only induce transient [CREB1] elevation. (C) An imposed elevation of [CREB1] is simulated by choosing initial conditions [CREB1]t=0=0.3 and [CREB2]t=0=0.2 (its basal steady state). [CREB1] elevation is paired with a single pulse of 5-HT. Persistent [CREB1] elevation results (solid line). Only a transient [CREB1] elevation results if five pulses of 5-HT are paired with elevation of CREB2 to 2.0 ([CREB1]t=0 is at basal steady state) (dashed line). Parameters are at basal values except when [5-HT]=10μM, Vx=3.7min−1.

If Vx is elevated above the first threshold for a sufficient duration, [CREB1] will increase to the HIGH state and will remain there when Vx returns to basal levels. Figure 6B illustrates that five simulated pulses of 10μM [5-HT] (Vx increased to 3.7min−1 for 5min, interpulse interval=15min) induce persistent [CREB1] elevation. However, three pulses cannot, in qualitative agreement with experiment 40.

The physiological significance of the [CREB1] threshold at [5-HT]=0 correlates with the CREB1 elevation required to induce the persistent HIGH state of [CREB1]. The simulation displayed as the solid line in Figure 6C pairs a single pulse of 5-HT (5min, [5-HT]=10μM) with an imposed elevation of [CREB1] (to 0.3) such as might be produced by injecting CREB1 into a neuron. This pairing induces long-lasting CREB1 elevation. Because increased CREB1 expression correlates with LTF, this simulation is consistent with the observation that injection of CREB1 allows LTF following a single 5-HT pulse 15. In the simulation, the imposed elevation of [CREB1] is near but below the threshold of CREB1 (middle, unstable steady state of [CREB1], Figure 6A). Thus, a single pulse of 5-HT suffices to move [CREB1] up past the threshold to elicit persistent CREB1 elevation. In contrast, the simulation displayed as the dashed line in Figure 6C pairs an imposed elevation of [CREB2] (to 2.0) with five pulses of 5-HT. With [CREB2]=2.0, only a lower stable steady state of [CREB1] exists; therefore, only a transient elevation of CREB1 results. This simulation is consistent with the empirical observation that injection of CREB2 blocks 5-HT-induced LTF 16. In the model, the elevated CREB2 inhibits CREB1 expression sufficiently to overwhelm positive feedback, eliminating the stable HIGH state of [CREB1].


Time delays in the expression of CREB1 or CREB2 can destabilize steady states and oscillations

Gene regulatory systems are characterized by time delays between regulation of transcription and appearance of functional gene product. Delays are created by processes such as translocation of mRNA and transcription factors between the cytoplasm and nucleus. Differential equations containing discrete time delays have often been used to provide approximate, qualitative models of biological systems with delay 29,30,31,32, and we adopt this technique here. The addition of discrete delays to ordinary differential equations such as Eqs. (1) allows for qualitative predictions of how delays affect dynamics without greatly altering the form or complexity of the model. However, we note that biological delays are not actually discrete. To model more precisely the processes underlying delays, such as macromolecular diffusion or active transport, partial differential equations would be required.

To qualitatively explore the effects of time delays on the model of Eqs. (1), discrete delays for the synthesis of CREB1, τCREB1, and of CREB2, τCREB2, were incorporated into the terms that describe the production processes of CREBs in Eqs. (1). The resulting DDEs are given as Eqs. (5). With low levels of negative feedback (Vy=0.01min−1, other parameters as in Figure 2A), neither τCREB1 or τCREB2 can induce a limit-cycle oscillation. For small Vy, the stability of the steady states in Figure 2A (lower and upper stable, middle unstable) is not affected by introducing τCREB1 or τCREB2. Negative feedback is required to sustain oscillations irrespective of time delays. With more negative feedback (Vy≥0.4min−1, Figure 2BD) either τCREB1 or τCREB2 can destabilize a steady state and generate a limit-cycle oscillation. Fig. 7 illustrates how introducing τCREB1 and τCREB2 impacts the stability of the steady-state and oscillatory solutions in Figure 2D1. As shown in Figure 7A1, starting from the stable limit cycle of Figure 2D1 (Vy=2min−1, Vx=1min−1), the amplitude of oscillations is rapidly reduced by increasing τCREB1. The stable oscillation disappears at a Hopf bifurcation, τCREB1=9min, and is replaced by a stable steady state for 9min≤τCREB1≤100min. This state is destabilized by another Hopf bifurcation at τCREB1=100min. The amplitude of the resulting stable oscillations varies in a complicated way for τCREB1>100min, as does the period (Figure 7A1, inset). Figure 7A2 illustrates the impact of τCREB2 on the limit cycle of Figure 2D1. For 0≤τCREB2<100min, τCREB2 rapidly increases the amplitude of oscillations; however, when τCREB2>100min, the amplitude reaches a plateau. The inset illustrates that the period varies linearly with τCREB2.

Display large version of this figure
Figure 7
Bifurcation diagrams illustrating effects of time delays. Diagrams were each constructed by a large number of individual simulations. Each point corresponds to a single simulation that was run sufficiently long to determine the attractor dynamics, whether a steady state or a limit cycle. For limit cycles, maxima and minima are plotted. Parameter values are as in Figure 2D1 (Vy=2min−1). (A) At Vx=1min−1, oscillations exist with no time delays. (A1) Bifurcation parameter τCREB1 (at τCREB2=0); destabilization of oscillation or stable steady state is through Hopf bifurcations. (Inset) Dependence of period on τCREB1. (A2) Bifurcation parameter τCREB2 (at τCREB1=0). (Inset) Dependence of period on τCREB2. (B) Vx=0.5min−1, a steady state of low [CREB1] exists with τCREB1=0 and τCREB2=60min. Increasing τCREB1 destabilizes the steady state through a Hopf bifurcation. (Inset) Dependence of period on τCREB1. (C) Vx=2.8min−1, a steady state of high [CREB1] exists with τCREB1=τCREB2=0. (C1) Bifurcation parameter τCREB1 (at τCREB2=0). (C2) Bifurcation parameter τCREB2 (at τCREB1=0). (Inset) Expansion of the τCREB2 axis.

For a low stable steady state of [CREB1] in Figure 2D1 (Vy=2min−1, Vx=0.6min−1) neither τCREB1 nor τCREB2 separately can destabilize this state. However, a proper combination of τCREB1 and τCREB2 can destabilize it (Figure 7B). At τCREB2=60min, increasing τCREB1 to 180min destabilizes the lower steady state. The amplitude of oscillations significantly increases with increasing τCREB1 and reaches a plateau for τCREB1>600min. The period varies linearly with τCREB1 (Figure 7B, inset).

Finally, starting from the stable steady state with high [CREB1] in Figure 2D1 (Vy=2min−1, Vx=2.8min−1), introduction of either τCREB1 or τCREB2 can destabilize the state. As illustrated in Figure 7C1 with τCREB2=0, an increase in τCREB1 to 180min induces an oscillation via Hopf bifurcation from the steady state. A further increase in τCREB1 to ∼850min or above generates a time course of [CREB1] that, although oscillatory, exhibits substantial random variation in oscillation amplitudes (broad dark area in Figure 7C1). Such a time course suggests chaotic dynamics. For τCREB1=0 and τCREB2 increasing, when τCREB2 reaches ∼3min (Figure 7C2, inset), an oscillatory solution bifurcates from the steady state. Increasing τCREB2 increases the period linearly (not shown).

Further simulations illustrated qualitatively similar effects of these time delays on the stability of the HIGH and LOW steady states in the bistable regions of Figure 2B, and of oscillatory solutions in Figure 2C 1. In summary, both τCREB1 and τCREB2 can destabilize either stable steady states or stable limit-cycle oscillations. Increasing τCREB1 may generate chaotic dynamics (Figure 7C1). In contrast, increasing τCREB2 appears only to generate regular oscillations.


For plausible average molecule numbers, bistability is preserved when stochastic fluctuations are considered

Stochastic molecular noise may have a significant impact on nonlinear system dynamics 41,42,43. To what extent can deterministic bifurcation diagrams predict the behavior when stochastic noise due to fluctuating copy numbers of molecular species is present? In this section, the relationship between deterministic bifurcation analysis and the dynamics with noise is examined. The Gillespie algorithm is used to simulate the dynamics with molecular noise. The propensity associated with each reaction step in Table 1 is a stochastic version of each kinetic term in Eqs. (1). At each time step, the algorithm randomly chooses which possible reaction event occurs, updates the system state accordingly, and also determines the time step to the next reaction event. In Table 1, the parameter Ω transforms the concentration units into molecule numbers, and the kinetic rate constants are expressed in terms of molecule numbers. Ω represents the system volume.

As illustrated in Fig. 8, deterministic bifurcation diagrams are able to predict the outcomes of stochastic simulations of the model in the bistability region. Figure 8A repeats the bifurcation diagram of Figure 2A with Vy=0.01min−1 (low negative feedback), in which there is a large bistable region. From this diagram, specific values of Vx were chosen for stochastic simulations. These were 1), Vx=0.1min−1, with a single low steady state of [CREB1] (Figure 8A, dash-dotted line B); 2), Vx=0.18min−1; 3), Vx=0.2min−1; 4), Vx=0.22min−1; 5), Vx=0.24min−1; and 6), Vx=0.7min−1. Values 2–5 all lie within the bistable region of Figure 8A (Figure 8A, inset, lines C–F); and value 6 corresponds to a single high steady state. For these values of Vx, the outcomes of stochastic simulations are illustrated in Figure 8BG, respectively.

Display large version of this figure
Figure 8
Relation between the deterministic bifurcation diagram and stochastic simulations in or near the region of bistability. (A) Bifurcation diagram of [CREB1] versus Vx at Vy=0.01min−1. When 0.177min−1≤Vx≤0.68min−1, there exist two stable steady states (solid line) and an unstable steady state (dashed line). (Inset) Expansion of Vx axis. (B–G) Stochastic simulation at six values of Vx corresponding to lines marked B–G in A. (B1–G1) Time courses of the number of CREB1 molecules. (B2–G2) Steady-state probability distributions. (B) Vx=0.1min−1. (C) Vx=0.18min−1. (D) Vx=0.2min−1. (E) Vx=0.22min−1. (F) Vx=0.24min−1. (G) Vx=0.7min−1. The volume factor Ω=10.

At Vx=0.1, the stochastic time course in Figure 8B1 exhibits a fluctuation between 0 and 10 in the number of CREB1 molecules. Figure 8B2 illustrates the corresponding steady-state probability distribution of the number of CREB1 molecules. For this case, corresponding to a single deterministic steady state of low [CREB1], the probability distribution has a unimodal shape. At Vx=0.18min−1, two stable steady states and one unstable steady state exist in Figure 8A. In the presence of molecular noise, CREB1 molecule numbers flip between the neighborhoods of the two stable steady states, as illustrated by the stochastic time series in Figure 8C1. Most commonly, CREB molecule numbers reside around the LOW state, and only occasionally near the HIGH state. Consequently, the steady-state probability distribution of the number of CREB1 molecules (Figure 8C2) exhibits a bimodal distribution 44,45 with a higher hump at a small number of CREB1 (∼1–2) and a lower hump at CREB1∼15.

When Vx is successively increased to 0.2min−1, 0.22min−1, and 0.24min−1, the corresponding CREB1 time series (Figure 8D1F1) illustrate that CREB1 molecule numbers reside progressively more commonly near the HIGH state. The steady-state bimodal probability distributions (Figure 8D2F2) shows a progressive shift toward the HIGH state. Upon further increasing Vx to 0.7min−1, which corresponds deterministically to a single steady state with high [CREB1], CREB1 molecule numbers are found to reside near the HIGH state (Figure 8G1), with a unimodal probability distribution (Figure 8G2).

Why do CREB1 molecule numbers sometimes preferentially reside around the LOW state (e.g., Figure 8C), and other times around the HIGH state (e.g., Figure 8F)? One important factor determining these preferences can be found by inspecting the bistability region of the deterministic bifurcation diagram. Inspection of the dash-dotted line C (Figure 8A, inset) illustrates that the middle unstable steady state is vertically closer to the upper steady state than to the lower steady state, whereas on line F (Figure 8A, inset), the middle unstable steady state is vertically closer to the lower steady state. The unstable steady state plays the role of threshold for the noise-induced transitions between the two stable steady states (i.e., the difference in the CREB1 molecule numbers between either stable steady state and the middle state is a barrier that has to be overcome for the occurrence of state transitions). The larger the barrier, the more robust the steady state is to molecular noise. On line F (Figure 8A, inset) one can predict that the LOW state is less robust to noise. This is the outcome of the stochastic simulation, with CREB1 molecule numbers usually residing near the HIGH state (Figure 8F2). From lines C–F (Figure 8A, inset), the distance from the HIGH state to the middle unstable state becomes larger, predicting that the residence time near the HIGH state should become larger. Stochastic simulations illustrate this, with the probability distribution peak shifting to the HIGH state (Figure 8CF). The deterministic bifurcation diagram helps one to predict the rough shape of the probability distribution and the relative residence time of CREB molecule numbers near stable states.

Next, the first transition time (FTT) 46,47 was employed to examine how the robustness of bistable steady states changes with the numbers of molecules of each species, which scale with the volumetric parameter Ω. FTT is defined as the time interval required for the system to transit between steady states (e.g., FTT (HIGH to LOW) represents the time interval at which a system initially at the HIGH state flips for the first time to the LOW state). A simulated FTT depends on the seed of the random number generator. Therefore, an average FTT over a number of stochastic simulation runs must be used to characterize the average lifetime of a stable steady state. Fig. 9 illustrates FTTs for HIGH to LOW and vice versa. In Figure 9A, for discrete values of Ω, FTTs (LOW to HIGH) are plotted as solid circles and FTTs (HIGH to LOW) as open circles. Each vertical set of small circles consists of 100 individual FTTs for each Ω, and within each vertical set, the large circle corresponds to the average FTT for that set. Figure 9AB, differ in that Vx=0.18min−1 for Figure 9A (Figure 8A, line C) and Vx=0.24min−1 for Figure 9B (Figure 8A, line F). Figure 9AB, illustrates that the average FTTs (HIGH to LOW) and (LOW to HIGH) both always increase with increasing Ω. Increasing Ω corresponds to increasing average molecule numbers, thereby approaching a deterministic system. Thus, random transitions between the HIGH and LOW states become rare. For each of the four sets of FTTs, the average FTTs (larger circles) are fit reasonably well by straight lines, as illustrated. These linear fits of the logarithmic plots indicate that the average FTTs are all exponentially dependent upon Ω, in agreement with Bialek 48.

Display large version of this figure
Figure 9
First transition time and the average of 100 simulations of the bistable steady states (HIGH and LOW) versus Ω. (A) FTT and the mean of the HIGH and LOW states versus Ω at Vx=0.18min−1, Vy=0.01min−1. Small solid circles denote individual FTTs (LOW to HIGH), and large solid circles are the averages of 100 vertically adjacent FTTs. Small open circles denote individual FTTs (HIGH to LOW), and large open circles are average FTTs. Solid line denotes the exponential fit of the average FTT (LOW to HIGH) versus Ω, and dashed line denotes the exponential fit of the average FTT (HIGH to LOW) versus Ω. (B) Similar to A, but at Vx=0.24min−1, Vy=0.01min−1. Other parameters are as in Fig. 8.

The average FTTs (LOW to HIGH) are much longer than the average FTTs (HIGH to LOW) for the case represented in Figure 9A, in which the deterministic HIGH state is closer to the unstable middle state than is the deterministic LOW state (Figure 8A, inset, line C). At Ω=10, the average number of CREB1 molecules is ∼20, and the average FTT (LOW to HIGH) (or the average lifetime of the LOW state) is ∼850h, whereas the average FTT (HIGH to LOW) is only ∼1.5h. Even at Ω=50 (∼100 CREB1 molecules on average), the average lifetime of HIGH is only ∼4h. However, with increased strength of positive feedback (Vx=0.24min−1) (Figure 8A, line F), there is a dramatic shift in the FTTs. The average FTT (HIGH to LOW) now becomes much larger than the average FTT (LOW to HIGH), as illustrated in Figure 9B. At Ω=8, the average number of CREB1 molecules is only ∼15, but nevertheless the average FTT (HIGH to LOW) is ∼450h. This means that ∼15 CREB1 molecules can sustain a HIGH state with an average lifetime of tens of days in this case. Thus, at a higher strength of positive feedback, small average molecule numbers suffice to maintain a persistent state of elevated CREB1 expression. As discussed previously, such a state might correlate physiologically with prolonged gene induction by CREB1 and consequent establishment of long-term memory.

We further examined how negative feedback impacts the robustness of the bistable steady states. By increasing the strength of negative feedback (e.g., varying Vy from 0.01min−1 to 0.06min−1), the distance of the HIGH state from the threshold becomes smaller for a fixed Vx=0.22min−1 (data not shown). Correspondingly, the average FTT (HIGH to LOW), decreases from ∼60h to ∼5h. For a fixed but strong negative feedback (Vy=0.07min−1), increasing Vx increases the distance of HIGH from the threshold and increases the average FTT (HIGH to LOW), as in Fig. 9. Thus the robustness of the “memory” (HIGH) state against noise can be increased substantially by decreasing negative feedback strength as well as by increasing positive feedback strength.


Molecular noise degrades deterministic oscillations, but can also induce oscillations via stochastic resonance

What molecule numbers are required for oscillations predicted by deterministic ODEs to persist in the presence of molecular noise? To examine this issue, we used the Gillespie method to simulate the effect of noise on the limit cycle of Figure 2C1 (Vx=0.9min−1, Vy=0.8min−1). When a small number of CREB molecules (Ω=10) are present in the system, the stochastic time course (Figure 10A, left) of the number of CREB1 molecules fluctuates between 1 and 140. The corresponding phase diagram of CREB1 versus CREB2 (Figure 10A, middle) exhibits an attractor very smeared out by noise, and the PSD versus frequency (Figure 10A, right) is rather broad. Therefore, at these low average molecule numbers, reproducible oscillations have been eliminated by noise. When Ω is increased to 50 (Figure 10B), the number of CREB1 molecules and the corresponding phase diagram exhibit less noisy oscillatory dynamics. The PSD exhibits a much narrower frequency distribution with a peak near 0.1h−1. Thus, with average molecule numbers in the low hundreds, regular oscillations are preserved despite noise. This range of molecule numbers was found to be sufficient to preserve oscillations in previous reduced models of gene regulation 29,30. When Ω increases to 200 (Figure 10C), fluctuations are still significant. However, the regularity of periodicity is prominent, as illustrated by the sharp PSD peak at ∼0.1h−1 (Figure 10C, right). With Ω=500, the oscillation is very close to the deterministic limit cycle (Figure 10D).

Display large version of this figure
Figure 10
Characterization of stochastic oscillatory dynamics at four values of Ω. (Left column) Time series of the number of CREB1 molecules. (Middle column) Molecule-number trajectories of CREB1 versus CREB2. (Right column) Power spectrum densities (arbitrary unit) versus frequency (h−1) of the time courses. (A) Ω=10. (B) Ω=50. (C) Ω=200. (D) Ω=500. Parameter values correspond to deterministic oscillations: Vx=0.9min−1, Vy=0.8min−1 (Figure 2C1), other parameters at basal values.

The above simulations illustrate degradation of oscillations by noise. However, in a phenomenon termed stochastic resonance, molecular noise has also been shown, in a variety of systems, to induce oscillations when the system is at a deterministic steady state but close to a Hopf bifurcation point 49,50. We examined whether the minimal model exhibits this phenomenon. Figure 11A illustrates the time courses of stochastic oscillations in the number of CREB1 molecules for Ω=50, 100, and 1000, with model parameters corresponding to a deterministic steady state (Vx=0.75min−1) (Figure 2D1) close to a supercritical Hopf bifurcation point (Vx=0.78min−1) (Figure 2D1). Large fluctuations in CREB1 molecule numbers occur for these values of Ω. Inspection of these time courses does not suffice to determine whether an oscillatory component is significant. It is necessary to construct PSD distributions.

Display large version of this figure
Figure 11
Molecular noise can yield stochastic resonance when a steady state (Vx=0.75min−1, Vy=2min−1) is close to a supercritical Hopf bifurcation (Vx=0.78min−1, Vy =2min−1, Figure 2D1). (A) Time courses of CREB1 at Ω=50, 100, and 1000. (B) Smoothed power spectrum density (PSD) versus frequency of stochastic oscillations at the values of Ω in A. The points A–C, with Ω=100, illustrated in B, are used to compute the effective SNR (Eq. (13)). (C) Effective SNR versus Ω. The existence of a peak at Ω=100 demonstrates stochastic resonance.

Here, the three PSDs were first computed as in Fig. 10. However, it is not possible to use these noisy PSDs to quantitatively characterize the periodicity of the stochastic oscillations. Therefore, the procedure proposed by Hou and Xin 50 was used to characterize the periodicity of the stochastic time courses in Figure 11A. The PSDs were smoothed by replacing each original point with the value obtained by averaging over that point and the 50 closest original points. These smoothed PSDs are illustrated in Figure 11B. The PSD for Ω=50 shows a weak but significant peak at a frequency of ∼0.06h−1. The PSD for Ω=100 shows a much sharper peak at about the same frequency. This peak illustrates that a significant oscillatory component is present in the Ω=100 time course. In contrast, for very large Ω (1000), the PSD fails to show a peak, signifying that as deterministic dynamics are approached, the periodic oscillation component disappears, and only random fluctuations about the stable steady state persist. These PSDs illustrate stochastic resonance, with a strong oscillatory component for the dynamics of Ω=100 in particular. An effective signal/noise ratio (SNR) can be defined by the following equation 50:

(13)

In Eq. (13), PSD(·) denotes the PSD level at a given point, and f(·) denotes the frequency at a given point. Point B corresponds to the maximum PSD, point A corresponds to the minimum PSD at the low-frequency side of point B, and point C is on the high-frequency side of point B and is determined by the condition PSD(C)=PSD(B)/e. This effective SNR characterizes the height of the peak in the PSD normalized by its relative width. As demonstrated in Figure 11C, there is a maximum in the SNR versus Ω curve in the vicinity of Ω=100, which corresponds to an optimal system size for inducing internal stochastic resonance.



Discussion

The minimal model is useful to delineate the dynamics of a putative gene regulatory motif involving interlocked positive and negative feedbacks

Different types of biochemical regulatory motifs (or minimal building blocks of biochemical circuitry), such as feedback and feedforward loops, have been suggested 51,52. The interlocked feedback loops of Fig. 1 represent a combination of minimal motifs, which can exhibit more complex behaviors than the individual components. For example, upon modulation of the relative strengths of positive and negative feedback (e.g., by manipulating the production rate of CREB1 activator or CREB2 repressor), the minimal model can exhibit complex behaviors, such as a bistable genetic toggle switch and oscillations of gene expression (Fig. 3).

The dynamic complexity arises from the nonlinear interplay of the different strengths of positive and negative feedback loops. Although one can construct a switch with just a single-gene positive feedback or double-negative feedbacks 52,53, and an oscillator with a single-gene negative feedback 54, it may be possible to use the CRE enhancer sequences and creb genes to construct a genetic circuit element that can transduce stimuli into multiple types of responses, manifesting both switch and oscillatory behaviors by modulating the relative strengths of the interlocked positive and negative feedbacks. Therefore, this interlocked positive and negative feedback circuit could act as a “multifunctional converter” in gene circuit design.

Interlocked positive and negative feedback loops play essential roles in determining system dynamics in various biological contexts. For example, positive and negative feedbacks of Cdc2-cyclin B and APC support oscillations in the cell cycle regulatory circuit 55; and positive and negative feedback associated with tumor suppressor protein p53 and Mdm2 may demonstrate bistability or oscillatory dynamics 56. Minimal models similar to that presented here, representing basic regulatory network structures, will be useful to systematically delineate complex behaviors of such biochemical architectures.


Deterministic bifurcation analysis and stochastic modeling provide complementary dynamic information

Simulations illustrated that a critical number of molecules is required to maintain a reproducible oscillation in the presence of molecular noise (e.g., Ω=200, Figure 10C, average number of CREB1∼600). The model also exhibits stochastic resonance. Molecular noise can move a system from a stable steady state into an oscillatory regime, if the steady state is near a deterministic Hopf bifurcation point (Fig. 11). An optimal noise level was found to yield stochastic resonance with a sharply peaked power spectrum (Ω=100) (Figure 11C).

In bistable regions of bifurcation diagrams (Figure 2AA and Figure 8AA), the system relaxes to either the HIGH or LOW state, and cannot transit between states unless forced by a large external perturbation. However, with molecular noise, molecules can flip between states, which leads to a bimodal probability density distribution of molecule numbers. Bimodal distributions have been experimentally validated in several biological systems 57,58,59.

Our investigation of the effect of stochastic noise on bistability illustrates how the deterministic bifurcation diagrams may be used to qualitatively determine the shape of bimodal distributions and the relative robustness of the stable states against molecular noise. By inspection of the bistable bifurcation diagram (e.g., Figure 8A), the relative distances from the threshold (middle unstable steady state) to the HIGH steady state versus the LOW steady state determine the shape of the bimodal probability density distribution. For example, if (e.g., Figure 8A, inset, line C), the threshold is relatively lower for the CREB1 molecule numbers to fluctuate down from the HIGH state to the LOW state, as compared with fluctuating up from the LOW state to the HIGH state. Therefore, the CREB1 molecule number more commonly resides near the LOW state, which has the larger peak in the probability distribution (Figure 8C2) and the average FTT (LOW to HIGH) is larger than the average FTT (HIGH to LOW) (Figure 9A). On the contrary, if the HIGH state is more robust and has a larger peak in the probability distribution (Figure 8F2).

This simple criterion, used alone, is only reliable if the distances of the two stable states from the unstable middle state are substantially different. For line D in Figure 8A, the HIGH state is only slightly, not substantially, farther from the threshold. Although for the probability distribution (Figure 8D), the peak for the HIGH state is lower, the area under the HIGH peak appears somewhat greater. In such a case, the average FTT (HIGH to LOW) and FTT (LOW to HIGH) at some Ω is further required to predict the shape of bimodal probability distribution, and to characterize the robustness of the HIGH and LOW states against molecular noise. For example, an approximately equal average FTT (HIGH to LOW) and FTT (LOW to HIGH) indicates an approximately identical area under the HIGH and LOW humps in the probability distribution.

Kummar et al. 60 recently proposed a measure to characterize the difference between deterministic and stochastic simulated behaviors. The divergence property (quantified by the sum of Liapunov exponents) of a stable or oscillating solution provides a measure of the robustness of the solution to molecular noise, with a more negative divergence corresponding to greater robustness. Applied to bistability, an equal divergence value of the HIGH and LOW states would indicate that the HIGH and LOW states have similar attractor properties, corresponding to similar areas below these peaks in the probability distribution. Therefore, this approach can provide a complementary method to predict the shape of a bimodal probability distribution in a region of bistability.

Further stochastic simulations illustrated that either increasing the strength of positive feedback (increasing Vx) or decreasing the strength of negative feedback (decreasing Vy) tends to increase the robustness and average lifetime of the HIGH state to molecular fluctuation. As discussed previously, such a state might correlate physiologically with prolonged gene induction by CREB1 and consequent establishment of long-term memory. For example, an average ∼1-h lifetime of the “memory” state with weak positive feedback (Vx=0.18min−1, Ω=8) (Figure 9A) increases to ∼450h given only a modest increase in positive feedback (Vx=0.24min−1, Ω=8) (Figure 9B). Similarly, decreasing negative feedback by decreasing Vy from 0.06min−1 to 0.01min−1 increases the average lifetime of the memory state from ∼5h to ∼40h (Vx=0.22min−1, Ω=10, data not shown). These results suggest that pharmacological interventions that vary the strength of positive and negative feedback may improve the formation of LTM, if long-lasting [CREB] elevation contributes to LTM. For example, negative feedback may be decreased by phosphorylation of CREB2 by MAPK 16. In addition, the putative positive feedback in which CREB1 activates its own gene requires phosphorylation of CREB1 14, and MAPK activation following 5-HT exposure correlates with phosphorylation of CREB1 61. These considerations suggest that pharmacological intervention enhancing and prolonging MAPK activation is likely to enhance CREB1 activation and the formation of LTM. A more detailed model of the dynamics of MAPK and CREB1 in Aplysia following 5-HT application has similarly predicted that prolongation of MAPK activation would enhance LTM 62.

In previous studies 62,63, we suggested an alternative positive feedback loop that might exhibit a bistable switch behavior contributing to the persistence of LTM. In this loop, PKA phosphorylates a transcription factor, thereby increasing the expression of an ubiquitin hydrolase (Ap-Uch). Ap-Uch then acts to increase PKA activity, closing the loop. Active PKA also phosphorylates CREB1. This phosphorylation is necessary for CREB1 to induce creb1 expression, thus sustaining the CREB1 positive feedback loop. Phosphorylated CREB1 may also induce creb2, thus sustaining the CREB2 negative feedback loop. Therefore, the feedback loops of CREBs examined in this study would be downstream effectors of the PKA/Ap-Uch feedback loop. In this way, all three feedback loops may contribute to the induction and persistence of CREB expression and LTM. Additional feedback loops may play a role in regulating CREB expression. Inducible cAMP early repressor (ICER) can bind to CREs and thereby repress the activity of its own expression, forming a negative feedback loop 18,64. In Aplysia, CREB1a (an isoform of CREB1) is the activator, whereas CREB1b (an alternatively spliced isoform of CREB1) is a repressor 15. CREB1b might form another negative feedback loop that interferes with the expression of CREBs. It would be of interest to incorporate these additional feedback loops into the minimal mo