MODELING IN PHYSIOLOGY
Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory systems

Paul Smolen, Douglas A. Baxter, and John H. Byrne

Department of Neurobiology and Anatomy, University of Texas Medical School, Houston, Texas 77030

    ABSTRACT
Top
Abstract
Introduction
Discussion
Appendix
References

To examine the capability of genetic regulatory systems for complex dynamic activity, we developed simple kinetic models that incorporate known features of these systems. These include autoregulation and stimulus-dependent phosphorylation of transcription factors (TFs), dimerization of TFs, crosstalk, and feedback. The simplest model manifested multiple stable steady states, and brief perturbations could switch the model between these states. Such transitions might explain, for example, how a brief pulse of hormone or neurotransmitter could elicit a long-lasting cellular response. In slightly more complex models, oscillatory regimes were identified. The addition of competition between activating and repressing TFs provided a plausible explanation for optimal stimulus frequencies that give maximal transcription. Such optimal frequencies are suggested by recent experiments comparing training paradigms for long-term memory formation and examining changes in mRNA levels in repetitively stimulated cultured cells. In general, the computational approach illustrated here, combined with appropriate experiments, provides a conceptual framework for investigating the function of genetic regulatory systems.

genetic models; transcriptional regulation; optimal stimulus frequencies; control theory

    INTRODUCTION
Top
Abstract
Introduction
Discussion
Appendix
References

REGULATION OF GENE expression by signals from outside and within the cell plays important roles in many biological processes, including development (38), the cell cycle (33), hormone action (6), and neural plasticity (8). As the basic principles of genetic regulation have been characterized, it has become increasingly evident that nonlinear interactions, positive and negative feedback within signaling pathways, time delays, protein oligomerization, and crosstalk between different pathways need to be considered to fully understand genetic regulation. A conceptual problem arises, however, of how to predict the functional properties of these complex, nonlinear biochemical systems.

Genetic regulatory systems have often been modeled simply as networks of Boolean logical elements (43, 44). However, the range of nonlinear behaviors exhibited by such complex biochemical systems can be more thoroughly understood from an explicitly mathematical, dynamic systems approach (10, 13, 48). Such an approach allows the array of tools and concepts developed for analysis of systems of ordinary differential equations to be brought to bear. For example, bifurcation analysis (13) is used to determine the steady-state solutions, oscillatory solutions, and other attracting solutions of a system of ordinary differential equations as a function of parameters and thereby determine parameter values at which qualitative transitions in the dynamic behavior of these systems occur. This is done numerically by means of specialized software such as the package AUTO (5) or in simple cases algebraically. Also, by tracing the solution structure of a system as a function of parameter values, either through repetitive numerical integration or more rigorously via bifurcation analysis, one can quantify the sensitivity of these systems to parameter changes.

Consequences of applying the dynamic systems approach to modeling of genetic regulatory systems might include identification of multiple steady states of gene product concentrations. If brief perturbations can induce transitions between these states, a system with this structure could act as a switch. Transient exposure to a stimulus, like a growth factor, might elicit a long-lasting response, such as a switch from cellular quiescence to growth. In addition, the dynamic systems approach can identify key physiological control parameters to which the behavior of specific genetic regulatory systems is particularly sensitive. Such parameters might provide targets for pharmacological intervention. Oscillatory regimes and regimes in which transcriptional behavior is particularly sensitive to initial conditions (i.e., apparently chaotic) can also be identified, and the stability of such regimes can be examined. Mechanistic hypotheses can be tested by first using simulations to predict behavioral characteristics of a particular mechanism and then testing for these experimentally.

Modeling biochemical systems, with nonlinear elements and feedback pathways, as dynamic systems has an extensive history. For example, circadian rhythms (9) and glycolysis (10) have been modeled as biochemical oscillators. However, aside from work of Keller (19) modeling simple generic genetic systems, relatively little application of dynamic systems techniques has been made to modeling genetic regulatory systems specifically. Studies modeling a specific genetically regulated process (e.g., Ref. 36) focus on reproduction of experimental data pertaining to that process and not on the elucidation of dynamic properties generic to genetic regulatory systems.

To elucidate dynamic properties of genetic regulation, we have modeled genetic regulatory systems with different levels of complexity. First, in a confirmation and extension of Keller's earlier research (19), we consider a relatively simple model of transcription factors (TFs) subject to positive and negative autoregulation of their own transcription. Principles of operation that give rise to multistable and oscillatory dynamic behavior are discussed. Second, this model is extended to allow time-dependent phosphorylation of interacting transcriptional activator and repressor proteins. Simulations with this extended model identify possible mechanisms for generating maximal transcription at an optimal stimulus frequency. A similar system with interacting transcriptional activators and repressors is believed to mediate early stages of long-term memory (LTM) formation, and the formation of LTM can be greatest for an optimal frequency of stimulus presentation. Finally, we point out that the principles identified in these models are likely to apply to a variety of genetic regulatory systems.

Our models do not include stochastic thermal fluctuations in molecule numbers. Such fluctuations could, however, be important when small numbers of important molecular species, e.g., TFs bound to a limited number of DNA sites, are present. For example, McAdams and Arkin (26) have suggested that different phenotypes of prokaryotic systems could be selected by stochastic switching between alternative dynamic states. One can say intuitively that fluctuations could destabilize or switch among multiple steady states and would tend to randomize the amplitude and period of oscillatory solutions. In specific systems in which the numbers of individual molecular species could be estimated, further simulations would be necessary to quantify the importance of such effects.

    PRINCIPLES OF DYNAMIC GENE REGULATION ILLUSTRATED BY SPECIFIC MODELS

Responsive elements affecting regulatory genes can provide crosstalk between genetic regulatory systems and feedback within systems. We consider signal-transduction pathways in which stimuli lead to second messenger generation and phosphorylation of TFs, which in turn bind to DNA sequences known as responsive elements and thereby regulate the transcription of specific genes (18). The regulatory activity of TFs is often modulated by phosphorylation and by intermolecular interactions. For example, TFs often bind to DNA as homodimers or as heterodimers of different TF family members.

A well-known family of TFs is the Ca2+/adenosine 3',5'-cyclic monophosphate (cAMP)-responsive element binding protein/activating TF (CREB/ATF) family of homo- and heterodimers, which bind to Ca2+/cAMP-responsive elements (CREs). A second well-known group of TFs is the activating protein-1 (AP-1) family of heterodimers of Fos and Jun proteins, which bind to phorbol ester-responsive elements. Responsive elements that bind TFs have been found to affect the transcription of genes for diverse TFs such as Jun, Fos, and CREB (18, 29, 39). As a result, some TFs, such as Jun, autoregulate their own transcription (28). Moreover, responsive elements can mediate crosstalk between pathways. For example, members of the CREB/ATF family of TFs activate transcription of fos (30), and at least one ATF-Fos heterodimer binds to CREs, potentially activating CREB/ATF transcription in turn (14). Thus responsive elements affecting TF gene transcription can provide crosstalk and positive feedback. Responsive elements have also been found that regulate genes for potent transcriptional inhibitors. An example is the inducible Ca2+/cAMP-responsive early repressor (ICER) protein, whose transcription is increased on binding of phosphorylated dimers of CREB to a nearby CRE (31). Such responsive elements provide negative feedback loops. Given the above interactions, there is the possibility for rich dynamic activity.

To explore this possibility, we constructed a model (Fig. 1A) that captures the salient features of TF dimerization, binding, and phosphorylation-dependent regulation of transcription. Simplifications were made to obtain a model that can be appreciated intuitively. A single transcriptional activator, which we term TF-A, is considered as part of a pathway mediating a cellular response to a stimulus. The TF forms a homodimer that can bind to responsive elements (TF-REs). The tf-a gene incorporates one of these responsive elements, and when homodimers bind to this element TF-A transcription is increased. Binding to the TF-REs is independent of dimer phosphorylation. Only phosphorylated dimers, however, can activate transcription. The fraction of dimers phosphorylated is dependent on the activity of kinases and phosphatases whose activity can be regulated by external signals. Thus this model incorporates both signal-activated transcription and positive feedback on the rate of TF synthesis.


View larger version (18K):
[in this window]
[in a new window]
 
Fig. 1.   Models of genetic regulation with positive and negative autoregulatory feedback loops. Genes are indicated by italic type, proteins by uppercase lettering. A: the transcription factor (TF) TF-A activates transcription when phosphorylated (P) and bound as a dimer to specific responsive-element DNA sequences (TF-REs). Degradation of TFs is also indicated. TF-A can stimulate transcription of later genes. For simplicity, phosphorylation and dimerization are not shown as separate kinetic steps. B: schematic resulting from addition of a 2nd TF, TF-R, which represses transcription by competing with TF-A dimer for binding to TF-REs; k1,d and k2,d, rate constants for degradation of TF-A and TF-R, respectively.

It is further assumed that the transcription rate of tf-a saturates hyperbolically with TF-A dimer concentration, to a maximal rate k1,f, which is itself proportional to the degree of TF-A phosphorylation (Eq. 1). This proportionality implies that a singly phosphorylated TF-A dimer is one-half as effective at activating transcription as is a doubly phosphorylated dimer. A basal rate of synthesis of activator (r1,bas) is present at negligible dimer concentration. TF-A is degraded with first-order kinetics, with a rate constant k1,d. Binding processes are considered comparatively rapid and close to equilibrium, so the concentration of homodimer is proportional to the square of TF-A monomer concentration ([TF-A]2). Here and subsequently, translation of mRNA into protein is not explicitly modeled. These simplifications give a model with a single ordinary differential equation for the concentration of TF-A
 <FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>k</IT><SUB>1,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>1,d</SUB></DE></FR> − <IT>k</IT><SUB>1,d</SUB> [TF-A] + r<SUB>1,bas</SUB> (1)
Here K1,d is the dissociation constant of TF-A dimer from TF-REs. Concentrations and concentration-based parameters such as K1,d are dimensionless because effective concentrations of TFs and of DNA elements in the nucleus are generally not known. Numerical integration of this and other models was carried out with Gear's adaptive step size method (7).

Transient phosphorylation of TFs can regulate the dynamics of genetic systems by switching among multiple stable states. One intriguing aspect of the simple gene regulatory system described by Eq. 1 is the emergence of bistability, i.e., two steady-state solutions for [TF-A] (Fig. 2A, left). At these solutions, d[TF-A]/dt = 0. These solutions are termed stable because if [TF-A] is slightly perturbed from these solutions, it will relax back. For a relatively wide range of parameters, there is one such solution with [TF-A] low and its synthesis rate close to r1,bas and another with [TF-A] high and its synthesis rate close to k1,f. Moreover, appropriate, larger perturbations can switch the model between these states. It is assumed that only properly phosphorylated dimers can activate transcription, as is the case for the TFs AP-1 and CREB (12, 42). Signal-dependent activation of protein kinases or phosphatases induces changes in the fraction of dimers phosphorylated. These changes are modeled as transient increases or decreases in k1,f. These perturbations induce upward or downward transitions, respectively, between the steady states. Such transitions could correspond physiologically to brief stimuli, such as exposure to a hormone, leading to long-lasting increases or decreases in the levels of particular proteins.


View larger version (17K):
[in this window]
[in a new window]
 
Fig. 2.   Complex dynamics of a genetic regulatory system, Eq. 1. A, left: perturbations switch system between steady states. TF-A concentration ([TF-A]) is initially in a low steady state. Parameter values are maximal rate of TF-A synthesis (k1,f) = 10 min-1, basal rate of TF-A synthesis (r1,bas) = 0.1 min-1, k1,d = 1 min-1, and dissociation constant of TF-A dimer from TF-REs (K1,d) = 10. Brief changes in k1,f (50 min-1 for t between 30 and 31 min, 1 min-1 for t between 60 and 63 min) give small excursions of [TF-A] from low steady state (dotted circles mark imperceptible transients in [TF-A]). At t = 90 min, k1,f is increased to 100 min-1 for 1 min, switching [TF-A] to an upper steady state. At t = 120 min, k1,f is increased to 50 min-1 for 1 min. At t = 150 min, k1,f is decreased to 1 min-1 for 3 min. [TF-A] then returns to lower steady state. A, right: k1,f is kept at 5 min-1 except for brief perturbations (here taken as identical to those at left). Only transient excursions of [TF-A] occur. Note change in scale for [TF-A] between left and right. B: bifurcation analysis of multistability in this model as a function of k1,f; other parameters are as in A, left. For 6.1 min-1 < k1,f < 25.1 min-1, there exist 2 stable steady-state solutions of [TF-A] (solid portions of d[TF-A]/dt curve) with an unstable solution between them (dashed curve). Outside this region there is only a single steady-state solution.

The transient responses of this system are also state dependent (Fig. 2A, left). With the system in the upper steady state, a perturbation of k1,f which gave only a minute excursion of [TF-A] above the lower steady state (dotted circle at left) now gives a very large excursion above the upper steady state. This difference in response magnitude following a state transition is a model for "priming" of a system to respond more vigorously to subsequent stimuli.

The range of parameters permitting bistability can be more precisely delineated by bifurcation analysis. Figure 2B, derived algebraically, demonstrates this technique for varying k1,f in Eq. 1. For each value of k1,f between 6.1 min-1 and 25.1 min-1, there are three values of [TF-A] that are steady states of Eq. 1. The upper and lower states are those illustrated in Fig. 2A. For the middle steady state, a small perturbation of [TF-A] will grow until [TF-A] moves to either of the other two steady states. Thus the middle state is unstable and would not be expected to be observed experimentally. For k1,f < 6.1 min-1 or > 25.1 min-1, Eq. 1 has only one stable steady-state solution for [TF-A].

Figure 2A, right, demonstrates that, for parameters outside the permissible range supporting multiple stable steady states, brief perturbations in k1,f only yield transient excursions in [TF-A], which return to the single steady state.

Combined positive and negative feedback, or time delays, can generate oscillations and complex transients in genetic systems. The model of Fig. 1A is similar to a model analyzed by Keller (case B of Ref. 19), which also exhibits bistability. Indeed, in that work, five other models based on autoactivation or autorepression of transcription by binding of gene product to DNA were also analyzed. All these models exhibit bistability over significant parameter ranges. Thus biochemical architectures capable of supporting bistability may be common in genetic regulatory systems. However, models similar to that of Fig. 1A have only one type of feedback: either positive or negative. Without both types of feedback, such a model cannot be expected to support oscillations. For example, positive feedback can act to drive [TF-A] to high levels, but then there is no process to bring [TF-A] back down. To investigate the effect of negative feedback, we introduced a protein, TF-R, that represses transcription by binding to TF-REs. Its rate of synthesis is increased by binding of the TF-A dimer to a TF-RE (Fig. 1B). TF-R competitively inhibits binding of TF-A dimers to TF-REs; thus it inhibits the transcription of the genes tf-a and tf-r. The equations for this model are
<FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>k</IT><SUB>1,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>1,d</SUB>(1 + [TF-R]/<IT>K</IT><SUB>R,d</SUB>)</DE></FR>
− <IT>k</IT><SUB>1,d</SUB> [TF-A] + r<SUB>1,bas</SUB> (2)
<FR><NU>d[TF-R]</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>k</IT><SUB>2,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>2,d</SUB>(1 + [TF-R]/<IT>K</IT><SUB>R,d</SUB>)</DE></FR>
− <IT>k</IT><SUB>2,d</SUB> [TF-R] (3)

Parameters in Eq. 2 have the same meaning as in the model of Eq. 1. KR,d is the dissociation constant of TF-R monomers from TF-REs. Parameters in Eq. 3 are analogous to those in Eq. 2 (k2,f, k2,d, and K2,d are maximal synthesis rate, degradation rate constant, and dissociation constant from TF-REs of TF-A dimers). Robust oscillations are readily generated by this model (Fig. 3A). Both TFs oscillate approximately in phase, with a period on the order of 1 h. Moreover, physiologically plausible variations in parameters can give disproportionate changes in oscillation amplitude and mean, thus providing for a signal-dependent modulation of the oscillations. In Fig. 3A, at t = 160 min the maximal rate of TF-A synthesis k1,f is reduced from 10.5 min-1 to 10 min-1. This modest parameter change gives a large change in the pattern of oscillation.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 3.   A: addition of repressor protein dynamics (Fig. 1B) to simple model of Fig. 1A allows oscillations to occur. Parameter values are k1,f = 10.5 min-1, r1,bas = 0.4 min-1, k1,d = 1 min-1, K1,d = 10, maximal rate of TF-R synthesis (k2,f) = 0.3 min-1, k2,d = 0.2 min-1, K2,d = 10, and KR,d = 0.2. At t = 160 min, k1,f is decreased to 10 min-1. This change induces a large damping of oscillations. B: bifurcation analysis showing oscillatory solution in [TF-A] as a function of k1,f; other parameters are as in A. For k1,f < 9.9 min-1 or k1,f > 10.8 min-1 there is only the stable steady state shown; for k1,f between these values there is a stable oscillatory state. Solid curves illustrate maxima and minima of oscillations, which surround an unstable steady state (dashed curve).

Bifurcation analysis can precisely determine the parameter range over which oscillations exist. We used the numerical software package AUTO (5) to trace the amplitude of oscillations as a function of parameters. Figure 3B gives an example. When k1,f is varied, with all other parameter values as in Fig. 3A, oscillations exist only in the interval 9.9 min-1 < k1,f < 10.8 min-1.

Biological oscillations with periods on the order of hours might be explained by such genetic models. Examples could be oscillatory secretion of hormones (3), such as human growth hormone (35) or leutinizing hormone-releasing hormone (21). In secretory cells, the genetic regulatory system utilizing CREB and related TFs has been postulated to yield oscillations of transcription rate (47). A negative feedback loop on gene expression, dependent on TF phosphorylation, has also been postulated within a molecular model for circadian rhythm generation (9).

In various models of biological phenomena, e.g., population dynamics (32), time delays serve as another way to generate oscillations or complex transients. In genetic regulatory systems, there are ubiquitous time delays associated with, for example, the transport of mRNA and proteins between the cytoplasm and the nucleus. To briefly explore the effect of introducing time delays into a simple model, a delay of several minutes was introduced in the model of Fig. 1A. The delay was between changes in TF-A concentration and the resultant changes in the rate of formation of new TF-A due to tf-a transcription. This converted Eq. 1 into a delay differential equation, which was integrated by the fourth-order Runge-Kutta method. As a result of incorporating the delay, a brief increase in k1,f leads to 10 or more large oscillations in the TF concentrations, with a period comparable to the delay, before the system stabilizes in the upper state (not shown).

Protein oligomerization can underlie complexity, such as multistability and oscillations, in the dynamics of genetic regulatory systems. Dimerization of TFs is essential for bistability and oscillations in these models. Figure 4 illustrates a graph of d[TF-A]/dt vs. [TF-A] for the simple model of Eq. 1. The synthesis rate of TF-A is a nonlinear function of [TF-A]. This allows the graph to have a complex shape, with three steady states of [TF-A] with d[TF-A]/dt = 0. The nonlinearity is due to the square of [TF-A] in Eq. 1, which represents dimer concentration. The middle state with d[TF-A]/dt = 0 is unstable to small perturbations in [TF-A], whereas the lower and upper states are stable and represent the two steady states illustrated in Fig. 2A, left.


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 4.   TF dimerization is important for bistability. Solid curve, d[TF-A]/dt vs. [TF-A] for model of Eq. 1, with parameters as in Fig. 2A. The 3 crossings of 0 on ordinate correspond to lower, middle, and upper steady states (SS1-SS3, respectively). Dashed curve, d[TF-A]/dt vs. [TF-A] for model obtained by replacing 2nd powers of [TF-A] in Eq. 1 with 1st powers; parameter values are as in Fig. 2A. In this model, SS1 and SS3 have been eliminated; single steady state remaining is SS2.

Qualitatively different results are obtained with a variant of our model in which just monomers of TF-A bind to the TF-REs. When Eq. 1 is changed so that the activation of tf-a transcription is proportional simply to [TF-A], the graph of d[TF-A]/dt vs. [TF-A] is qualitatively altered (Fig. 4) such that only one stable state at which d[TF-A]/dt = 0 remains.

We also find that in the model with repressor TF-R, Eqs. 2 and 3, dimerization of TF-A is essential for oscillations in TF concentration. If transcription of TFs depends only on the concentration of monomeric TF-A, the system always settles to a steady state. Thus it can be inferred that protein oligomerization provides an important mechanism for achieving complex dynamics in transcription.

Stimulus frequency can be decoded in complex ways into variations in response strength. The above models do not consider further dynamic behaviors of genetic regulatory systems that could result from concurrent modification of transcriptional activator and repressor efficacy by stimulus-dependent phosphorylation. We therefore added an additional level of complexity to the model of Eqs. 2 and 3. The fractions of repressor and activator proteins that are phosphorylated are made dynamic variables. Only when phosphorylated can these proteins affect transcription. Different kinetics of phosphorylation and dephosphorylation for the two TFs can allow for different sensitivities to stimuli, such that stimuli sufficient to saturate the phosphorylation of one TF may not do so for the other.

Such a model might be able to simulate experiments that have demonstrated optimal stimulus frequencies for activation, or repression, of transcription. For example, transcription of the cell adhesion molecule L1 in cultured neurons is strongly repressed by imposed, continual 0.1-Hz electrical stimulation but not repressed significantly by 0.3-Hz stimulation (17). Also, c-fos transcription in cultured neurons is enhanced almost 200% by bursts of 6 electrical stimuli at 10 Hz with an interburst interval of 1 min but not significantly enhanced by bursts of 12 stimuli with an interburst interval of 2 min (41). Continuous 0.1-Hz stimuli gave a 70% enhancement. A model in which an intermediate intensity or frequency of stimulation phosphorylated and activated one TF only, whereas a higher intensity of stimulation activated also a second TF that counteracted the effect of the first, might explain these phenomena. Moreover, a similar model might apply to initial steps in LTM formation. It has recently been proposed that changes in the relative phosphorylation of transcriptional activators and repressors may be important for induction of transcription required for the formation of LTM (49). The relationship between stimulus frequency and amount of LTM formation, and by inference amount of transcription, appears sometimes to be nonmonotonic. In Drosophila, spaced olfactory stimulus presentations, with a relatively long interstimulus interval (ISI), yield much more LTM than do massed presentations (a short ISI) even given the same total training time (thus more massed presentations) (46). In this system, there is some optimal stimulus frequency for formation of LTM, and by inference for activation of transcription. Optimal stimulus frequencies appear also to exist for types of task learning by humans (20, 24).

We developed a model with two antagonistic TFs for the purpose of testing qualitatively to what extent a single model could provide a unified explanation of these results concerning optimal stimulus frequencies and groupings. In the model, the dimeric character of TF-A and TF-R allows different efficacies of activation, or repression, by singly vs. doubly phosphorylated dimers. As before, monomer-dimer equilibria and the existence of heterodimers are neglected. Only homodimers of TF-A and TF-R, which compete for binding to responsive elements (TF-REs), are considered. Phosphorylation of dimers of TF-R is considered necessary for binding to TF-REs. Phosphorylation does not affect binding of TF-A to TF-REs, in accordance with data for the specific transcriptional activator CREB implicated in the formation of LTM (37). The amount of gene transcription during a given time interval is assumed to be proportional to the concentration of phosphorylated, TF-RE-bound TF-A integrated over that interval. The overall flow of our model from stimulus to the transcription of genes that mediate a cellular response is diagrammed in Fig. 5A, and the binding and phosphorylation processes are schematized in more detail in Fig. 5B. The corresponding equations as given in the APPENDIX are more complex than Eqs. 2 and 3 because phosphorylation of TF-R affects binding to DNA so that separate differential equations are needed for each TF-R species.


View larger version (23K):
[in this window]
[in a new window]
 
Fig. 5.   Schematic of TF-A and TF-R regulatory model. A: overall flow of model from stimulus to gene transcription. A stimulus induces phosphorylation of both TF-A and TF-R, which compete for binding to TF-RE. In some simulations (dashed box), product of gene regulated by TF-A and TF-R itself represses transcription of a further gene product (e.g., L1). B: binding and phosphorylation. Dimers of TF-A and of TF-R are phosphorylated in response to a stimulus. TF-A binds to TF-REs irrespective of phosphorylation but can only activate transcription when phosphorylated. TF-R binds competitively to TF-REs but only when phosphorylated.

TF phosphorylation and dephosphorylation were simulated with Michaelis-Menten kinetics. For example, TF-A phosphorylation has a time-dependent maximal velocity kA,f (t) and a Michaelis constant KA,ph; TF-R phosphorylation has a maximal velocity kR,f (t). Relatively low values were chosen for the Michaelis constants of phosphorylation. This choice introduces zero-order ultrasensitivity, defined generally as an amplification of the response of a reversible covalent modification system to perturbations in enzymatic activity when the opposing enzymes operate in a zero-order kinetic regime (11). Specifically, a modest increase in the maximal velocity of phosphorylation of either TF, from somewhat less than to somewhat greater than the maximal velocity of dephosphorylation of that for TF, results in a large increase in the fraction of the TF phosphorylated. Further details of equations and parameters are given in the APPENDIX. In general, the maximal velocities of TF phosphorylation would be dynamic variables whose values depend on the activation of kinases, located in the nucleus, by stimulus applications. These stimuli are transduced in a possibly complex, time-dependent manner from the cell membrane (where stimuli are sensed) to these kinases. However, for simulation of ISIs on the order of seconds, as are utilized in Ref. 17, we assume that temporal averaging occurs during this transduction so that different stimulus frequencies correspond to different constant values of the maximal velocities of TF phosphorylation.

As illustrated in Fig. 6A, when either kA,f (t) or kR,f (t) is varied alone, the transcription rate varies monotonically. However, let us now specialize to the case in which the first-order dephosphorylation rate constant for TF-A is less than that for TF-R and increase together both phosphorylation maximal velocities (keeping them equal). One then sees a sharp rise in transcription rate to a peak as the maximal velocity for TF-A phosphorylation passes that for TF-A dephosphorylation and most of the TF-A becomes phosphorylated, followed by a sharp decline in transcription rate as the maximal velocity of TF-R phosphorylation exceeds that for dephosphorylation and most of the TF-R becomes phosphorylated. Figure 6B illustrates that the sharpness of this tuning curve depends strongly on the Michaelis constants for phosphorylation and dephosphorylation being small enough to allow rapid shifts in TF phosphorylation states. Larger Michaelis constants give a flatter curve. Also, we investigated whether TF dimerization was important for the qualitative dynamics of this model as had been found for the models of Fig. 1. Figure 6B illustrates that, indeed, if TF-R is assumed to be monomeric, the sharpness of the tuning curve is greatly reduced.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 6.   Dependence of transcription rate (in units of min-1) in model of Fig. 5 on rate constants for TF phosphorylation. The 1st set of parameter values in APPENDIX is used, except as noted. A: increasing only forward phosphorylation rate constant for TF-A (kA,f) gives a monotonically increasing transcription rate [bullet ; forward rate constant for TF-R phosphorylation (kR,f) fixed to equal backward rate constant for TF-R phosphorylation (kR,b)], and increasing only kR,f for TF-R gives a monotonically decreasing transcription rate (open circle ; kA,f fixed at 0.1 min-1), but increasing both rate constants while keeping their ratio at 1 yields an optimum (black-square). B: peak of transcription rate (black-square, same as corresponding curve in A) is greatly reduced by increasing Michaelis constants for TF-A and TF-R phosphorylation and dephosphorylation to 10 (open circle ). Also, sharpness is greatly reduced by assuming monomeric TF-R (bullet ), with all other parameters unchanged, i.e., phosphorylated TF-R monomers bind to DNA with on and off rate constants k2,f and k2,b, and Rtot is a TF-R monomer concentration. C: with transcription of cell adhesion molecule L1 added, L1 transcription and [L1] as a function of stimulus frequency have pronounced minima. kR,f and kA,f exhibit small oscillations about mean values. These values, in units of min-1, are both assumed to be equal to value of stimulus frequency in Hz.

We have qualitatively simulated the data of Itoh et al. (17) that demonstrate suppression of transcription of the cell adhesion molecule L1 by continuous 0.1-Hz but not 0.3-Hz stimuli. As indicated in Fig. 5A, an extra kinetic step is necessary, with TF-A and TF-R regulating the transcription of a gene whose protein product represses the transcription of the L1 gene. This converts the maximum of transcription of the gene regulated by TF-A and TF-R into a minimum of L1 transcription. We assume that, in response to stimuli applied at the membrane, maximal velocities of TF phosphorylation vary on a time scale considerably longer than the intervals between stimuli in the protocols of Itoh et al. (17). Therefore, after an initial transient, these velocities exhibit a steady-state behavior of small oscillations about mean values. Figure 6C illustrates that, with these mean values taken as directly proportional to stimulus frequency, the model predicts that 0.1-Hz stimulation produces an ~90% reduction in L1 transcription from the basal rate. In contrast, 0.3-Hz stimulation gives only a 7% reduction. Stimulation at 0.1 Hz suffices to phosphorylate TF-A only, which activates the transcription of the repressor of L1 transcription. Stimulation at 0.3 Hz phosphorylates TF-R as well, so the repressor of L1 transcription has its own transcription repressed. Small values for the Michaelis constants of TF-R phosphorylation and dephosphorylation are necessary. Only then can the model give a virtually complete cancellation of strong transcriptional repression in response to a threefold increase in stimulus frequency (from 0.1 to 0.3 Hz).

As previously mentioned, we had also hoped to simulate experiments concerning c-fos transcription (41). Given a fixed average stimulus frequency, bursts of 6 stimuli at 10 Hz repeated every minute were reported to yield much more transcription than evenly spaced 0.1-Hz stimuli or bursts of 12 stimuli at 10 Hz repeated every 2 min. However, our current model cannot simulate these results. If it is assumed, as above, that velocities of TF phosphorylation are approximately proportional to stimulus frequency, the stimulus paradigms would be expected to yield, on average, the same velocities of phosphorylation of the TFs, and the same transcription rates, because the number of stimuli averaged over time is the same in all paradigms. A more detailed model of stimulus coupling to nuclear events, considering nonlinear kinetics of particular second messenger systems, may be required.

To explain why massed stimulus presentations are less effective than spaced presentations in producing LTM in Drosophila, Yin et al. (49) proposed the same generic mechanism that is considered here. An intermediate intensity or frequency of stimulation phosphorylates and activates a TF that activates transcription of genes essential for LTM formation, while a higher frequency of stimulation activates also a second TF that counteracts the effect of the first. However, rather than assuming fixed average values for phosphorylation velocities, Yin et al. (49) assumed that the net dephosphorylation rate for the repressor TF-R is faster than that of the activator TF-A during ISIs. Then, during spaced stimuli, the net difference (phosphorylated activator - phosphorylated repressor) becomes large during the long ISIs, but during massed stimuli TF-A phosphorylation is always approximately canceled out by TF-R phosphorylation. The kinetic scheme of Fig. 5 is again used to test this hypothesis. Because the ISIs are now on the order of minutes rather than seconds, we assume that each stimulus abruptly resets the phosphorylation rate constants to maximal values that decay exponentially.

As Fig. 7 demonstrates, our model predicts an optimal stimulus frequency for transcription, and by inference for LTM formation, when maximal velocities for TF dephosphorylation are chosen in accordance with the hypothesis of Yin et al. (49). We also found (not shown) that qualitatively similar results are obtained if both TFs are dephosphorylated at identical rates and it is assumed instead that the phosphorylation rate of TF-R is slower than that of TF-A during exposure to a stimulus. Then TF-R is again only able to become highly phosphorylated during massed stimuli. In addition, alternative kinetic schemes utilizing only one TF were also found to give an optimal stimulus frequency for transcription. One such model variant postulates both activating and inhibiting phosphorylation sites on TF-A, with the inhibiting site only becoming significantly phosphorylated by massed stimuli. Another model variant relies on competing kinase and phosphatase activities. In principle, these model variants could also explain the aspects of L1 transcriptional regulation simulated above (Fig. 6C). It may be inferred that the existence of two competing processes, such as activator and repressor phosphorylation, that have different sensitivity to stimuli and opposing effects on transcription of a specific gene could provide a general mechanism for tuning the response of a genetic system to an optimum stimulus frequency.


View larger version (16K):
[in this window]
[in a new window]
 
Fig. 7.   Nonmonotonic dependence of transcription rate on stimulus frequency in model of Fig. 5 with kinetic parameters consistent with hypothesis of Yin et al. (49) for explaining greater efficacy of spaced stimuli in formation of LTM (APPENDIX, 2nd set of parameter values). A: time course of transcription rate (in units of min-1) during spaced stimuli. B: massed stimuli are used, and in all other respects simulation and graph are as in A, including time scale. Comparing A and B demonstrates that, over 100 min, 100 massed stimuli [interstimulus interval (ISI) = 1 min] produce considerably less transcription, and by inference less LTM formation, than 8 spaced stimuli (ISI = 15 min). C: dependence of transcription rate on ISI for 2 cases. Top curve, dephosphorylation rate constant for TF-A (kA,b; 1.0 min-1) < kR,b (7.0 min-1), as in hypothesis of Yin et al. (49). Bottom curve, dephosphorylation rate constants are equal (1.0 min-1).

    DISCUSSION
Top
Abstract
Introduction
Discussion
Appendix
References

Biochemical nonlinearities such as dimerization, feedback loops, and time delays are common in genetic regulatory systems (10). Our results indicate that incorporating these features into models of relatively simple genetic regulatory systems can give rise to complex dynamic activity and nonmonotonic dependence of response strength on stimulus. Thus the dynamic principles illustrated are likely to be important in phenomena in which regulation of transcription has an essential role, such as development and the formation of LTM.

It is of value to compare the method of modeling genetic regulatory systems as sets of ordinary differential equations, used herein, with other approaches. Genetic regulatory systems have also been modeled as networks of Boolean logical elements. In this approach, genes are considered as either on or off and biochemical connectivities are reduced to logical update rules for determining the set of genes that will be on at a given time step as a function of the genes that were on at the preceding time step. Relatively large time steps are used in numerical simulations. By consideration of the mathematical properties of such networks, significant insights concerning the expected dynamics of genetic systems have been obtained, e.g., by Thomas and colleagues (44, 45). These authors point out, for example, the importance of negative feedback loops for maintaining homeostasis in levels of gene products and of positive feedback loops for allowing multiple stable states of these levels. More concretely, Somogyi and Sniegoski (43) have developed a computational method for efficiently modeling large genetic networks in a Boolean fashion, which is designed to determine connectivity relations from simultaneously measured experimental expression time courses of sets of genes.

However, although the dynamic systems approach is more computationally intensive than the Boolean approach, we prefer the former because it is a more physically correct representation. Not only are there quantitative differences, but, more importantly, it has been found that simple biochemical models, when expressed as a network of Boolean logical elements instead of as a system of ordinary differential equations, can exhibit qualitatively different, and spurious, attractors, i.e., different stable oscillatory solutions (1). Although it has not been specifically demonstrated that more complex biochemical and/or genetic Boolean network models also suffer from this flaw, it seems plausible that a significant number would. However, it should be noted that in investigating complex systems that would require more differential equations than could be computationally simulated in a reasonable amount of time, applying the Boolean network approach while keeping in mind the above caveat may be the best practical alternative.

Some of the phenomena illustrated in this work can also be dealt with by modeling within the framework of electrical circuit analysis, as McAdams and Shapiro (27) have done with a recent model of the lambda  phage-Escherichia coli genetic regulatory system. These authors constructed a model of this system in which many nonlinear biochemical processes were represented as Boolean logical switches (the product of a reaction such as transcription of a given gene is either present or absent, depending on the value of an effector variable). However, other biochemical reactions were still modeled as continuous input-output relations and numerically integrated as ordinary differential equations. Time delays were also incorporated in some kinetic steps between effector concentration changes and rate changes. This method is therefore a hybrid of the continuous approach based on ordinary differential equations and the approach based on modeling genes as Boolean logical elements.

One specific application of the techniques discussed in this work could be to model the genetic regulatory system responsible for initial steps in LTM formation in greater detail. Signals that raise levels of cAMP, such as pulses of neurotransmitter, induce phosphorylation of the TF CREB (39). CREB can then induce the transcription of immediate early genes crucial for neuronal plasticity and formation of LTM (2, 34, 50). Proteins related to CREB, such as CREB2, are transcriptional repressors that bind to the same CRE sequence as CREB (39). They are generally phosphorylated by the same signals that phosphorylate CREB. Although the functional relevance of repressor phosphorylation has not yet been established, the basic architecture of the model of Fig. 5 appears to be present.

Recent data, obtained in Aplysia, actually suggest that phosphorylation of CREB2 by mitogen-activated protein kinase reduces its repressing activity (25). This would contradict the mechanism for generation of an optimal stimulus frequency hypothesized by Yin et al. (49) and simulated in Fig. 7, which relies on phosphorylation enhancing repression. An optimal stimulus frequency for transcription in this system might still, however, be explained in terms of two competing processes that have different sensitivity to stimuli and opposing effects on transcription. For example, CREB has both activating and inhibiting phosphorylation sites (Ser-133 and Ser-142, respectively), and, if the inhibiting site only became significantly phosphorylated in response to massed stimuli, an optimum could result.

Additional biochemical elements of the CREB genetic system might allow qualitatively new dynamics. For example, positive feedback exists via binding of CREB to CREs affecting its own transcription. Negative feedback exists in the form of a repressor protein, ICER, whose transcription is induced by CREB. These feedback loops could create multistability or oscillatory behavior. There is some empirical indication for complex dynamics mediated by these elements. Oscillations in CREB mRNA have been reported in secretory cells, and the above feedback loops have been proposed as essential components of the oscillatory mechanism (47).

More generally, other biochemical architectures that have been observed in genetic regulatory systems but not captured in any of our models so far include convergence of different signaling pathways through distinct TFs onto a particular gene [e.g., TFs of the ETS family and pituitary TF Pit-1 onto the prolactin gene (15)], heterodimerization of TFs in two separate pathways with a common third TF [e.g., thyroid hormone receptor and peroxisome proliferator-activated receptor heterodimerize with the retinoic acid receptor (16)], and conditional regulation by a single TF [e.g., enhancement of basal transcription of the dopamine beta -hydroxylase gene along with suppression of cAMP-induced transcription by the TF YY1 (40)]. Such architectures may be expected to provide alternative mechanisms for generating dynamic phenomena such as multistability and oscillations.

It is unlikely that all transcriptionally regulated genes will exhibit the behaviors illustrated by our models. However, the diversity of TFs and their interactions suggest that behaviors such as these will be identified. Indeed, our models are simplifications of the actual kinetic schemes characterizing genetic systems. MacLeod (23) has recently proposed that epigenetic, heritable changes in gene expression following exposure to chemicals might play a role in carcinogenesis. Such changes would correspond dynamically to perturbations of genetic regulatory systems from one steady state to another.

An outstanding major issue for future investigation will be to determine whether the parameters of specific genetic systems in vivo are permissive for specific types of dynamic behavior. Experiments to help determine this might include monitoring transcription of transfected reporter gene constructs, with defined promoters subject to regulation by TFs, in cultured cells during specific patterns of hormone or neurotransmitter application. A prolactin promoter-luciferase gene construct has been used to provide real-time quantification of promoter activity in cultured secretory cells (4). Another relevant technique is polymerase chain reaction amplification and quantitation of specific mRNAs from tissue samples (43); however, this technique does not resolve dynamics at the single cell level. We believe that, as the dynamic behaviors of gene networks are explored empirically, the present work can provide a conceptual framework for the analysis and interpretation of such experiments.

    APPENDIX
Top
Abstract
Introduction
Discussion
Appendix
References

Details of equations and parameters for simulations with the model of Fig. 5 follow. We make some assumptions consistent with experimental results concerning the dynamics of a specific system with competing TFs that is thought to mediate the initial steps in LTM formation. In particular, we have been guided by analyses of competitive interactions of the transcriptional activator CREB and related repressors for their target DNA sequences, i.e., CRE sites. It is assumed that 1) total amounts of TF-A and TF-R remain constant (2), 2) phosphorylation does not affect binding of TF-A to CRE sites (37), and 3) singly phosphorylated TF-A dimers have one-half the activity of doubly phosphorylated ones, with unphosphorylated dimers inactive (22).

We denote the concentration of free CRE sites by Gfree. For brevity, in equations the activator TF-A is denoted by A and the repressor TF-R by R. [AP] is used to denote the concentration of phosphorylated TF-A sites, and [AA] is used to denote the concentration of free TF-A dimers. [AAB] denotes the concentration of TF-A dimer bound to DNA. Atot denotes the total concentration of TF-A dimers. The total concentration of repressor dimers is Rtot. [RR] is the concentration of free, unphosphorylated TF-R dimers. RRP and RRPP denote single or double phosphorylation; RRPPB denotes bound, doubly phosphorylated R dimers. Because phosphorylation of TF-A is independent of binding to DNA, the rate of change of phosphorylated TF-A can be described by a single differential equation for the concentration of phosphorylated TF-A sites
<FR><NU>d[AP]</NU><DE>d<IT>t</IT></DE></FR> 
= <IT>k</IT><SUB>A,f</SUB> (<IT>t</IT>) <FR><NU>2A<SUB>tot</SUB> − [AP]</NU><DE>2A<SUB>tot</SUB> − [AP] + <IT>K</IT><SUB>A,ph</SUB></DE></FR> − <IT>k</IT><SUB>A,b</SUB> <FR><NU>[AP]</NU><DE>[AP] + <IT>K</IT><SUB>A,deph</SUB></DE></FR> (A1)
where KA,deph is the Michaelis constant for TF-A dephosphorylation and kA,b is the backward rate constant for TF-A phosphorylation.

A single differential equation describes the rates of change of free and bound TF-A dimers because of conservation of total dimers. The association (forward) and dissociation (backward) rate constants are k1,f and k1,b, respectively
<FR><NU>d[AA]</NU><DE>d<IT>t</IT></DE></FR> = −<IT>k</IT><SUB>1,f</SUB> [AA] G<SUB>free</SUB> + <IT>k</IT><SUB>1,b</SUB> [AAB] (A2)
where [AAB] = Atot - [AA].

Separate equations are needed for the rates of change of each species of TF-R because binding and phosphorylation are not independent. Total phosphorylation and dephosphorylation rates are first expressed in terms of site concentrations, and then, in the differential equations, fractions of these total rates appropriate to each molecular species are used
r<SUB>R,ph</SUB> = <IT>k</IT><SUB>R,f </SUB>(<IT>t</IT>) <FR><NU>2[RR] + [RRP]</NU><DE>2[RR] + [RRP] + <IT>K</IT><SUB>R,ph</SUB></DE></FR> (A3)
r<SUB>R,deph</SUB> = <IT>k</IT><SUB>R,b</SUB> <FR><NU>2[RRPP]+[RRP]</NU><DE>2[RRPP]+[RRP]+<IT>K</IT><SUB>R,deph</SUB></DE></FR> (A4)
<FR><NU>d[RR]</NU><DE>d<IT>t</IT></DE></FR> = − r<SUB>R,ph</SUB> <FR><NU>2[RR]</NU><DE>2[RR] + [RRP]</DE></FR>
+ r<SUB>R,deph</SUB> <FR><NU>[RRP]</NU><DE>2[RRPP] + [RRP]</DE></FR> (A5)
<FR><NU>d[RRP]</NU><DE>d<IT>t</IT></DE></FR> = r<SUB>R,ph</SUB> <FR><NU>2[RR]</NU><DE>2[RR] + [RRP]</DE></FR>
 − r<SUB>R,deph</SUB> <FR><NU>[RRP]</NU><DE>2[RRPP] + [RRP]</DE></FR>
− r<SUB>R,ph</SUB> <FR><NU>[RRP]</NU><DE>2[RR] + [RRP]</DE></FR>
+ r<SUB>R,deph</SUB> <FR><NU>2[RRPP]</NU><DE>2[RRPP] + [RRP]</DE></FR> (A6)
<FR><NU>d[RRPP]</NU><DE>d<IT>t</IT></DE></FR> = r<SUB>R,ph</SUB> <FR><NU>[RRP]</NU><DE>2[RR] + [RRP]</DE></FR>
 − r<SUB>R,deph</SUB> <FR><NU>2[RRPP]</NU><DE>2[RRPP] + [RRP]</DE></FR>
 − <IT>k</IT><SUB>2,f</SUB> [RRPP] G<SUB>free</SUB>
 + <IT>k</IT><SUB>2,b</SUB> [RRPPB] (A7)
<FR><NU>d[RRPPB]</NU><DE>d<IT>t</IT></DE></FR> = <IT>k</IT><SUB>2,f</SUB> [RRPP] G<SUB>free</SUB> − <IT>k</IT><SUB>2,b</SUB> [RRPPB] (A8)
where KR,ph and KR,deph are Michaelis constants for TF-R phosphorylation and dephosphorylation, respectively, and kR,b and k2,b are the backward rate constants for TF-R phosphorylation and TF-R binding to DNA, respectively.

The rate of transcription of the target gene for whose promoter region TF-A and TF-R compete (rRep) is taken as proportional to the concentration of bound TF-A dimers multiplied by the fraction of phosphorylated TF-A sites, with a rate constant kRep
r<SUB>Rep</SUB> = <IT>k</IT><SUB>Rep</SUB> <FR><NU>[AAB] [AP]</NU><DE>A<SUB>tot</SUB></DE></FR> (A9)
There is a conservation condition on the total number of DNA binding sequences
G<SUB>free</SUB> = G<SUB>tot</SUB> − [AAB] − [RRPPB] (A10)
where Gtot is the total number of CRE sites.

For modeling the data of Itoh et al. (17) indicating an optimal frequency for repression of transcription of the cell adhesion molecule L1, an additional kinetic step is needed. The target gene for TF regulation is assumed to express a protein Rep that represses transcription of the L1 gene. L1 transcription is assumed to proceed at a basal rate rL1 in the absence of Rep. L1 transcription only occurs if Rep is not bound to a promoter for the L1 gene. Rep dimers bind to this promoter with dissociation constant KRep
<FR><NU>d[Rep]</NU><DE>d<IT>t</IT></DE></FR> = r<SUB>Rep</SUB> − [Rep] (A11)
<FR><NU>d[L1]</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>r<SUB>L1</SUB></NU><DE>(1 + [Rep]<SUP>2</SUP>/<IT>K</IT><SUB>Rep</SUB><SUP>2</SUP>)</DE></FR> − [L1] (A12)
In simulations of the formation of LTM, we posited that each stimulus immediately set kA,f (t) and kR,f (t) to maximal values kA,max and kR,max, respectively. After a stimulus, kA,f (t) and kR,f (t) decayed to zero with time constants tau 2 and tau 1, respectively.

All simulations used parameter values from one of the two following sets. Concentrations are left dimensionless due to lack of sufficient experimental data. Parameters marked "varies" have values given in the text or in Fig. 6 or 7.

For the simulations of Fig. 6
<AR><R><C>    G<SUB>tot</SUB> = 0.1</C><C>    A<SUB>tot</SUB> = 1.0</C><C>    R<SUB>tot</SUB> = 3.0</C></R><R><C>    <IT>k</IT><SUB>A,f</SUB> and <IT>k</IT><SUB>R,f</SUB> vary</C><C>    <IT>K</IT><SUB>A,ph</SUB> = 1.0</C><C>    <IT>K</IT><SUB>A,deph</SUB> = 1.0</C></R><R><C>    <IT>k</IT><SUB>A,b</SUB> = 0.1 min<SUP>−1</SUP></C><C>    (Fig. 6, <IT>A</IT> and <IT>B</IT>),   0.05 min<SUP>−1</SUP></C><C>    (Fig. 6<IT>C</IT>)</C></R><R><C>    <IT>k</IT><SUB>R,b</SUB> = 0.18 min<SUP>−1</SUP></C><C>    <IT>k</IT><SUB>1,f</SUB> = 10.0 min<SUP>−1</SUP></C><C>    <IT>k</IT><SUB>1,b</SUB> = 10.0 min<SUP>−1</SUP></C></R><R><C>    <IT>K</IT><SUB>R,ph</SUB> = 0.5</C><C>    <IT>K</IT><SUB>R,deph</SUB> = 0.5</C><C>    <IT>k</IT><SUB>2,f</SUB> = 10.0 min<SUP>−1</SUP></C></R><R><C>    <IT>k</IT><SUB>2,b</SUB> = 1.0 min<SUP>−1</SUP></C><C>    r<SUB><IT>L1</IT></SUB> = 1.0</C><C>    <IT>K</IT><SUB>Rep</SUB> = 0.11</C></R><R><C>    <IT>k</IT><SUB>Rep</SUB> = 5.0 min<SUP>−1</SUP></C></R></AR>
For the simulations of Fig. 6C, kA,f and kR,f were assumed to execute small oscillations about the mean values in the figure legend. In the absence of data to construct a kinetic model for these oscillations, we merely assumed sinusoidal oscillations with a frequency equal to the stimulus frequency and an amplitude of 5% of the mean value.

For the simulations of Fig. 7
<AR><R><C>    G<SUB>tot</SUB>= 0.1</C><C>    A<SUB>tot</SUB> = 1.0</C><C>    R<SUB>tot</SUB> = 3.0</C></R><R><C>    <IT>K</IT><SUB>A,ph</SUB> = 20.0</C><C>    <IT>K</IT><SUB>A,deph</SUB> = 20.0</C><C>    <IT>k</IT><SUB>A,b</SUB> = 0.7 min<SUP>−1</SUP></C></R><R><C>    <IT>k</IT><SUB>R,b</SUB> = 7.0 min<SUP>−1</SUP></C><C>    <IT>k</IT><SUB>1,f</SUB> = 10.0 min<SUP>−1</SUP></C><C>    <IT>k</IT><SUB>1,b</SUB> = 10.0 min<SUP>−1</SUP></C></R><R><C>    <IT>K</IT><SUB>R,ph</SUB> = 10.0</C><C>    <IT>K</IT><SUB>R,deph</SUB> = 10.0</C><C>    <IT>k</IT><SUB>2,f</SUB> = 10.0 min<SUP>−1</SUP></C></R><R><C>    <IT>k</IT><SUB>2,b</SUB> = 1.0 min<SUP>−1</SUP></C><C>    ISI varies</C><C>    <IT>k</IT><SUB>A,max</SUB> = 30 min<SUP>−1</SUP></C></R><R><C>    <IT>k</IT><SUB>R,max</SUB> = 30 min<SUP>−1</SUP></C><C>    &tgr;<SUB>1</SUB> = 4.0 min</C><C>    &tgr;<SUB>2</SUB> = 4.0 min</C></R></AR>

    ACKNOWLEDGEMENTS

We thank Pramod Dash, Ron Dror, and Shogo Endo for comments on an earlier draft of the paper and B. Ermentrout for use of his XPP program for simulations.

    FOOTNOTES

This work was supported by Office of Naval Research Grant N0014-95-1-0579, by National Institutes of Health Grants K05-MH-00649, T32-NS-07373, and R01-RR-11626, and by Texas Higher Education Coordination Board Grant 011618-048.

Address for reprint requests: J. H. Byrne, Dept. of Neurobiology and Anatomy, University of Texas Medical School, 6431 Fannin St., Suite 7.046, Houston, TX 77030.

Received 30 April 1997; accepted in final form 27 October 1997.

    REFERENCES
Top
Abstract
Introduction
Discussion
Appendix
References

1.   Bagley, R. J., and L. Glass. Counting and classifying attractors in high dimensional dynamical systems. J. Theor. Biol. 183: 269-284, 1996[Medline].

2.   Bartsch, D., M. Ghirardi, P. Skehel, K. Karl, S. Herder, M. Chen, C. Bailey, and E. R. Kandel. Aplysia CREB2 represses long-term facilitation: relief of repression converts a transient facilitation into long-term functional and structural change. Cell 83: 979-992, 1995[Medline].

3.   Brabant, G., K. Prank, and C. Schofl. Pulsatile patterns in hormone secretion. Trends Endocrinol. Metab. 3: 183-190, 1992.

4.   Castano, J., R. Kineman, and L. S. Frawley. Dynamic monitoring and quantification of gene expression in single, living cells: a molecular basis for secretory cell heterogeneity. Mol. Endocrinol. 10: 599-605, 1996[Abstract].

5.   Doedel, E. AUTO: a program for the automatic bifurcation analysis of autonomous systems. Congr. Num. 30: 265-284, 1981.

6.   Edwards, D. R. Cell signaling and the control of gene transcription. Trends Pharmacol. Sci. 15: 239-244, 1994[Medline].

7.   Gear, C. The numerical integration of ordinary differential equations. Math. Comput. 21: 146-156, 1967.

8.   Ghosh, A., and M. Greenberg. Calcium signaling in neurons: molecular mechanisms and cellular consequences. Science 268: 239-247, 1995[Medline].

9.   Goldbeter, A. A model for circadian oscillations in the protein PER. Proc. R. Soc. Lond. B Biol. Sci. B261: 319-324, 1995[Medline].

10.   Goldbeter, A. Biochemical Oscillations and Cellular Rhythms. Cambridge, UK: Cambridge University Press, 1996.

11.   Goldbeter, A., and D. Koshland. An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. USA 78: 6840-6844, 1991.

12.   Gonzalez, G., P. Menzel, J. Leonard, W. Fischer, and R. Montminy. Characterization of motifs which are critical for activity of the cyclic AMP-responsive transcription factor CREB. Mol. Cell. Biol. 11: 1306-1312, 1991[Medline].

13.   Guckenheimer, J., and P. Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Heidelberg, Germany: Springer-Verlag, 1983.

14.   Hai, T., and T. Curran. Cross-family dimerization of transcription factors Fos, Jun and ATF/CREB alters DNA binding specificity. Proc. Natl. Acad. Sci. USA 88: 3720-3724, 1991[Abstract].

15.   Howard, P., and R. Maurer. A composite Ets/Pit-1 binding site in the prolactin gene can mediate transcriptional responses to multiple signal transduction pathways. J. Biol. Chem. 270: 20930-20936, 1995[Abstract/Free Full Text].

16.   Hunter, J., A. Kassam, C. Winrow, R. Rachubinski, and J. Capone. Crosstalk between the thyroid hormone and peroxisome proliferator-activated receptors in regulating peroxisome proliferator-responsive genes. J. Mol. Cell. Endocr. 116: 213-221, 1996.

17.   Itoh, K., B. Stevens, M. Schachner, and R. D. Fields. Regulated expression of the neural cell adhesion molecule L1 by specific patterns of neural impulses. Science 270: 1369-1372, 1995[Abstract].

18.   Karin, M. Signal transduction from the cell surface to the nucleus through the phosphorylation of transcription factors. Curr. Opin. Cell Biol. 6: 415-424, 1994[Medline].

19.   Keller, A. Model genetic circuits encoding autoregulatory transcription factors. J. Theor. Biol. 172: 169-185, 1995[Medline].

20.   Kientzle, M. J. Properties of learning curves under varied distributions of practice. J. Exp. Psychol. 36: 187-211, 1946.

21.   Knobil, E. Remembrance: the discovery of the hypothalamic gonadotropin-releasing hormone pulse generator and of its physiological significance. Endocrinology 131: 1005-1006, 1992[Medline].

22.   Loriaux, M., R. Brennan, and R. Goodman. Modulatory function of CREB-CREMa heterodimers depends upon CREMa phosphorylation. J. Biol. Chem. 269: 28839-28843, 1994[Abstract/Free Full Text].

23.   MacLeod, M. A possible role in chemical carcinogenesis for epigenetic, heritable changes in gene expression. Mol. Carcinog. 15: 241-250, 1996[Medline].

24.   Madigan, S. Intraserial repetition and coding processes in free recall. J. Verb. Learn. 17: 493-507, 1969.

25.   Martin, K., D. Michael, J. Rose, M. Barad, A. Casadio, H. Zhu, and E. R. Kandel. MAP kinase translocates into the nucleus of the presynaptic cell and is required for long-term facilitation in Aplysia. Neuron 18: 899-912, 1997[Medline].

26.   McAdams, H., and A. Arkin. Stochastic mechanisms in gene expression. Proc. Natl. Acad. Sci. USA 94: 814-819, 1997[Abstract/Free Full Text].

27.   McAdams, H., and L. Shapiro. Circuit simulation of genetic networks. Science 269: 650-656, 1995[Medline].

28.   McKnight, S., and Y. Yamamoto. Transcriptional Regulation. Cold Spring Harbor, NY: Cold Spring Harbor Laboratory Press, 1992.

29.   Meyer, T., G. Waeber, J. Lin, W. Beckmann, and J. Habener. The promoter of the gene encoding cAMP response element binding protein contains cAMP response elements: evidence for positive autoregulation of gene transcription. Endocrinology 132: 770-780, 1993[Abstract].

30.   Moens, U., N. Subramaniam, B. Johansen, and J. Aarbakke. The c-fos cAMP-response element: regulation of gene expression by a beta -2-adrenergic agonist, serum and DNA methylation. Biochim. Biophys. Acta 1173: 63-70, 1993[Medline].

31.   Molina, C., N. Foulkes, E. Lalli, and P. Sassone-Corsi. Inducibility and negative autoregulation of CREM: an alternative promoter directs the expression of ICER, an early response repressor. Cell 75: 875-886, 1993[Medline].

32.   Murray, J. Mathematical Biology. Heidelberg, Germany: Springer-Verlag, 1989.

33.   Okayama, H., A. Nagata, S. Jinno, and H. Murakami. Cell cycle control in fission yeast and mammals: identification of new regulatory mechanisms. Adv. Cancer Res. 69: 17-62, 1996[Medline].

34.   O'Leary, F., J. Byrne, and L. Cleary. Long-term structural remodeling in Aplysia sensory neurons requires de novo protein synthesis during a critical time period. J. Neurosci. 15: 3519-3525, 1995[Abstract].

35.   Prank, K., M. Kloppstech, S. Nowlan, T. Sejnowski, and G. Brabant. Self-organized segmentation of time series: separating growth hormone secretion in acromegaly from normal controls. Biophys. J. 70: 2540-2547, 1996[Abstract].

36.   Reinitz, J., E. Mjolsness, and D. Sharp. Model for cooperative control of positional information in Drosophila by bicoid and maternal hunchback. J. Exp. Zool. 271: 47-56, 1995[Medline].

37.   Richards, J., H. Bachinger, R. Goodman, and R. Brennan. Analysis of the structural properties of cAMP-responsive element binding protein (CREB) and phosphorylated CREB. J. Biol. Chem. 271: 13716-13723, 1996[Abstract/Free Full Text].

38.   Rossant, J., and N. Hopkins. Of fin and fur: mutational analysis of vertebrate embryonic development. Genes Dev. 6: 1-13, 1992[Medline].

39.   Sassone-Corsi, P. Transcription factors responsive to cAMP. Ann. Rev. Cell Dev. Biol. 11: 355-377, 1995.[Medline]

40.   Seo, H., C. Yang, H. Kim, and K. Kim. Multiple protein factors interact with the cis-regulatory elements of the proximal promoter in a cell-specific manner and regulate transcription of the dopamine beta -hydroxylase gene. J. Neurosci. 16: 4102-4112, 1996[Free Full Text].

41.   Sheng, H., R. D. Fields, and P. Nelson. Specific regulation of immediate-early genes by patterned neuronal activity. J. Neurosci. Res. 35: 459-467, 1993[Medline].

42.   Smeal, T., B. Binetruy, D. Mercola, A. Grover-Bardwick, G. Heidecker, U. Rapp, and M. Karin. Oncoprotein-mediated signalling cascade stimulates c-Jun activity by phosphorylation of serines 63 and 73. Mol. Cell. Biol. 12: 3507-3513, 1992[Abstract].

43.  Somogyi, R., and C. Sniegoski. The gene expression matrix: towards the extraction of genetic network architectures. In: Proceedings of the Second World Congress of Nonlinear Analysis. Amsterdam: Elsevier Science. In press.

44.   Thomas, R., and R. d'Ari. Biological Feedback. Boca Raton, FL: CRC, 1990.

45.   Thomas, R., D. Thieffry, and M. Kaufman. Dynamical behaviour of biological regulatory networks. I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull. Math. Biol. 57: 247-276, 1995[Medline].

46.   Tully, T., T. Preat, S. Boynton, and M. Del Vecchio. Genetic dissection of consolidated memory in Drosophila melanogaster. Cell 79: 35-47, 1994[Medline].

47.   Walker, W., L. Fucci, and J. Habener. Expression of the gene encoding transcription factor CREB: regulation by follicle-stimulating hormone-induced cAMP signaling in primary rat sertoli cells. Endocrinology 136: 3534-3545, 1995[Abstract].

48.   Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Heidelberg, Germany: Springer-Verlag, 1990.

49.   Yin, J., M. Del Vecchio, H. Zhou, and T. Tully. CREB as a memory modulator: induced expression of a dCREB2 activator isoform enhances long-term memory in Drosophila. Cell 81: 107-115, 1995[Medline].

50.   Yin, J., and T. Tully. CREB and the formation of long-term memory. Curr. Opin. Neurobiol. 6: 264-267, 1996[Medline].


AJP Cell Physiol 274(2):C531-C542
0363-6143/98 $5.00 Copyright © 1998 the American Physiological Society