## Abstract

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

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 Ca^{2+}/adenosine 3′,5′-cyclic monophosphate (cAMP)-responsive element binding protein/activating TF (CREB/ATF) family of homo- and heterodimers, which bind to Ca^{2+}/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 Ca^{2+}/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.1
*A*) 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.

It is further assumed that the transcription rate of*tf-a* saturates hyperbolically with TF-A dimer concentration, to a maximal rate*k*
_{1,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 (r_{1,bas}) is present at negligible dimer concentration. TF-A is degraded with first-order kinetics, with a rate constant*k*
_{1,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
Equation 1Here *K*
_{1,d} is the dissociation constant of TF-A dimer from TF-REs. Concentrations and concentration-based parameters such as*K*
_{1,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.2
*A*,*left*). At these solutions, d[TF-A]/d*t* = 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 r_{1,bas} and another with [TF-A] high and its synthesis rate close to*k*
_{1,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*k*
_{1,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.

The transient responses of this system are also state dependent (Fig.2
*A*,*left*). With the system in the upper steady state, a perturbation of*k*
_{1,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. Figure2
*B*, derived algebraically, demonstrates this technique for varying*k*
_{1,f} in*Eq*.*
1
*. For each value of*k*
_{1,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. 2
*A*. 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*k*
_{1,f} < 6.1 min^{−1} or > 25.1 min^{−1},*Eq*. *
1
*has only one stable steady-state solution for [TF-A].

Figure 2
*A*,*right*, demonstrates that, for parameters outside the permissible range supporting multiple stable steady states, brief perturbations in*k*
_{1,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. 1
*A* 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.1
*A* 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.1
*B*). 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
Equation 2
Equation 3

Parameters in *Eq*.*
2
* have the same meaning as in the model of *Eq*.*
1
*.*K*
_{R,d} is the dissociation constant of TF-R monomers from TF-REs. Parameters in*Eq*. *
3
*are analogous to those in *Eq*.*
2
*(*k*
_{2,f},*k*
_{2,d}, and*K*
_{2,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.3
*A*). 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. 3
*A*, at*t* = 160 min the maximal rate of TF-A synthesis *k*
_{1,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.

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 3
*B* gives an example. When*k*
_{1,f} is varied, with all other parameter values as in Fig.3
*A*, oscillations exist only in the interval 9.9 min^{−1} <*k*
_{1,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.1
*A*. 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*k*
_{1,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]/d*t* 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]/d*t* = 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]/d*t* = 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. 2
*A*,*left*.

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]/d*t* vs. [TF-A] is qualitatively altered (Fig. 4) such that only one stable state at which d[TF-A]/d*t* = 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.5
*A*, and the binding and phosphorylation processes are schematized in more detail in Fig. 5
*B*. The corresponding equations as given in the
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.

TF phosphorylation and dephosphorylation were simulated with Michaelis-Menten kinetics. For example, TF-A phosphorylation has a time-dependent maximal velocity*k*
_{A,f} (*t*) and a Michaelis constant *K*
_{A,ph}; TF-R phosphorylation has a maximal velocity *k*
_{R,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
. 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.6
*A*, when either *k*
_{A,f} (*t*) or*k*
_{R,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 6
*B* 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. Figure6
*B* illustrates that, indeed, if TF-R is assumed to be monomeric, the sharpness of the tuning curve is greatly reduced.

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.5
*A*, 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. Figure6
*C* 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. 6
*C*). 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.

## DISCUSSION

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 λ 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 β-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.

## Acknowledgments

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.

## Appendix

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 G_{free}. 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. A_{tot} denotes the total concentration of TF-A dimers. The total concentration of repressor dimers is R_{tot}. [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
Equation A1where *K*
_{A,deph}is the Michaelis constant for TF-A dephosphorylation and*k*
_{A,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*k*
_{1,f} and*k*
_{1,b}, respectively
Equation A2 where [AAB] = A_{tot} − [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
Equation A3
Equation A4
Equation A5
Equation A6
Equation A7
Equation A8where *K*
_{R,ph}and *K*
_{R,deph} are Michaelis constants for TF-R phosphorylation and dephosphorylation, respectively, and*k*
_{R,b} and*k*
_{2,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 (r_{Rep}) 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*k*
_{Rep}
Equation A9There is a conservation condition on the total number of DNA binding sequences
Equation A10where G_{tot} 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 r_{L1} 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*K*
_{Rep}
Equation A11
Equation A12In simulations of the formation of LTM, we posited that each stimulus immediately set *k*
_{A,f} (*t*) and*k*
_{R,f} (*t*) to maximal values*k*
_{A,max} and *k*
_{R,max}, respectively. After a stimulus,*k*
_{A,f} (*t*) and*k*
_{R,f} (*t*) decayed to zero with time constants τ_{2} and τ_{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
For the simulations of Fig. 6
*C*,*k*
_{A,f} and*k*
_{R,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

## Footnotes

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.

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.

- Copyright © 1998 the American Physiological Society