## Abstract

To predict the dynamics of genetic regulation, it may be necessary to consider macromolecular transport and stochastic fluctuations in macromolecule numbers. Transport can be diffusive or active, and in some cases a time delay might suffice to model active transport. We characterize major differences in the dynamics of model genetic systems when diffusive transport of mRNA and protein was compared with transport modeled as a time delay. Delays allow for history-dependent, non-Markovian responses to stimuli (i.e., “molecular memory”). Diffusion suppresses oscillations, whereas delays tend to create oscillations. When simulating essential elements of circadian oscillators, we found the delay between transcription and translation necessary for oscillations. Stochastic fluctuations tend to destabilize and thereby mask steady states with few molecules. This computational approach, combined with experiments, should provide a fruitful conceptual framework for investigating the function and dynamic properties of genetic regulatory systems.

- transcriptional regulation
- transport delays
- circadian oscillations
- multistability
- genetic modeling

computational modeling of biochemical processes often assumes a homogenous cytoplasmic medium. However, to predict actual dynamics, it may be necessary to consider that transport of macromolecules between the nucleus and the cytoplasm can take a time on the order of hours, and the mechanism of transport will therefore help determine dynamics on this time scale. There are two basic types of macromolecular transport: passive diffusion and active transport along cytoskeletal elements mediated by motor proteins such as kinesins or dyneins (17, 26). A first approximation for modeling active transport might be obtained by assuming a discrete time delay for movement of macromolecules from their place of synthesis to the location where they exert an effect. Consideration of qualitative differences in the behavior of models incorporating diffusion vs. a time delay can be expected to yield insights into the dynamics of cellular processes that incorporate diffusional vs. active transport of macromolecules. If a discrete delay is too drastic a simplification, another approach is to assume a distributed delay, with the derivative of a concentration dependent on an integral over a specified range of previous time (27).

Genetic regulation is a process in which the mechanism of transport can be expected to have a profound influence. Regulation of gene expression by signals from outside and within the cell plays important roles in many biological processes, including development (38), hormone action (9), and neural plasticity (3, 16, 37). Recently, the formation and movement of single β-actin mRNA transcripts have been visualized by fluorescent in situ hybridization (10). In about one-half of the cases examined, the movement of transcripts away from the transcription site appeared to follow specific tracks. This finding suggests that an active transport mechanism might be operating to direct mRNAs along cytoskeletal elements. In the remaining cases, however, the transcripts appeared to simply diffuse away from their site of formation. These results motivated us to compare the dynamics of model genetic regulatory systems incorporating diffusional vs. active modes of transport for mRNA and protein.

Previous work has examined the dynamics of models of genetic regulatory systems that use time delays to model the translocation of molecules within the cell (2, 29-31). These studies focus largely on determining conditions for steady states to lose stability and oscillatory solutions to be concurrently formed. Relatively little work has compared the effect of time delays with the effect of diffusion. Exceptions are the studies of Busenberg and Mahaffy (5) and Mahaffy and Pao (31), which considered the stability of steady states in models of genetic control by repression that incorporate diffusion and delays. As discrete delays were increased, steady states could lose stability to oscillatory solutions. Slowing diffusion, by contrast, damped oscillations.

The present study focuses on comparing the dynamic characteristics of transitions between coexisting steady states in models in which macromolecular transport is dominated either by diffusion or by active transport. Active transport was modeled as a discrete or narrowly distributed time delay. We identified several clear qualitative differences that could in principle be detected experimentally. We found that active transport can produce “staircase” transitions of a series of steps in transcription rate or in nuclear concentrations of macromolecules. Moreover, a system dominated by active transport can possess a type of short-term memory, such that the amplitude of a response to a brief perturbation depends on the length of time since a prior change in transcription rate. Also, in a system dominated by active transport, a single increase in transcription rate due to a brief stimulus can give a series of distinct subsequent “echoing” increases in macromolecular concentrations. We also compared the effects of diffusion with those of active transport on the capability for oscillations in systems with positive and negative feedback. A delay tended to create an oscillatory attractor with a period similar to the delay, whereas diffusion tended to damp oscillations.

To illustrate specific applications of such models, we considered two genetic systems, the dynamics of which are strongly dependent on macromolecular transport and time delays: circadian oscillators and genes responsible for long-term synaptic facilitation (LTF). For circadian oscillators, we found that time delays are of value for encapsulating complex biochemical processes but need to be supplemented by more detailed biochemical models to address specific issues. For LTF, it was recently suggested that transport of a messenger protein from active synapses to the nucleus may provide a signal for induction of essential genes (8). We found that including active synapse-to-nucleus transport of a messenger protein greatly sharpened this signal. However, this conclusion was sensitive to parameter variations within the physiological range.

Finally, we considered the effects of stochastic fluctuations in molecule numbers. We extended the previous results of McAdams and Arkin (34) by finding that stochastic fluctuations can preferentially destabilize and thereby mask the existence of steady states characterized by low concentrations.

### Glossary

- CRE
- Ca
^{2+}/cAMP response element - CREB
- CRE binding protein
*D*- Diffusion coefficient
- LTF
- Long-term facilitation
- LTM
- Long-term memory
- mRNA
- messenger RNA
- TF
- Transcription factor
- TF-A
- Transcriptional activator
- TF-R
- Transcriptional repressor
- TF-RE
- TF responsive element

## RESULTS

### Active Transport (i.e., Models With Delays) Allows Qualitatively Novel Responses to Stimuli

Genetic regulatory systems often are activated by signal transduction pathways, in which stimuli (e.g., hormones or neurotransmitters) lead to phosphorylation of transcription factors (TFs), which in turn bind to DNA sequences known as responsive elements and thereby regulate the transcription of specific nearby genes (21). Some TFs, such as Jun and possibly Ca^{2+}/cAMP response element (CRE) binding protein (CREB), autoregulate their own transcription (35, 44). Ubiquity of genetic autoregulation in even relatively simple organisms is suggested by an inventory of*Escherichia coli*ς^{70} promoter regulation, identifying 21 regulatory proteins that repress their own synthesis and 4 that activate their own synthesis (7).

#### Stimulus-response properties of models with transport represented by a time delay.

We added a discrete or distributed delay for macromolecular transport to a model utilizing a single TF that activates its own transcription (i.e., positive feedback; Fig.1
*A*). The behavior of this model without delay was considered previously (42). The factor, TF-A, forms a homodimer that can bind to responsive elements (TF-REs). The *tf-a* gene incorporates a TF-RE, and when homodimers bind to this element,*tf-a* transcription is increased. Responses to stimuli are modeled by varying the degree of TF-A phosphorylation. The transcription rate saturates with TF-A dimer concentration to a maximal rate*k*
_{f}, which is proportional to TF-A phosphorylation. At negligible dimer concentration, the synthesis rate is R_{bas}. TF-A is eliminated with a rate constant *k*
_{d}. Binding processes are considered comparatively rapid, so the concentration of dimer is proportional to [TF-A]^{2}. The variable [TF-A] is the nuclear concentration of TF-A. The rate constants *k*
_{f} and*k*
_{d}, which set the time scale for [TF-A] equilibration, are fairly rapid (e.g., *k*
_{d} = 0.1 min^{−1}) in simulations below (Fig. 2). This assumption of fairly rapid equilibration of [TF-A] would probably not be reasonable for overall cellular [TF-A], because the equilibration time would be on the order of the degradation time for TF-A protein. However, a short time scale for equilibration is more likely for nuclear [TF-A]. This is because the rate constants *k*
_{f} and*k*
_{d} include implicitly entrance and exit of TF-A protein from the relatively small nuclear volume and are thus larger than those governing the dynamics of overall cellular [TF-A]. The model incorporated a time delay τ = τ_{1} + τ_{2}, with τ_{1} the time taken for the transcription of *tf-a* mRNA and its movement to translation and τ_{2}the time required for movement of TF-A protein to the nucleus. This delay appears between any change in the level of nuclear TF-A and the appearance in the nucleus of TF-A synthesized in response to that change.

These considerations yield a model consisting of a single delay differential equation for the concentration of TF-A monomer in the nucleus. For discrete delay
Equation 1with *K*
_{d} the dissociation constant of dimer from TF-REs. The first term (in braces) on the right-hand side is evaluated at a time τ previous to the time when d[TF-A]/d*t* is computed.

The simplest reasonable distributed delay assumes that [TF-A] is averaged over a time interval, and this average is used to calculate d[TF-A]/d*t* at a later time. This corresponds to assuming that the times required for individual *tf-a* mRNAs to be transcribed and translated are equally likely to lie anywhere within a specific interval and never lie outside it. The same assumption is made for TF-A protein movement. This results in
Equation 2where the angle brackets average [TF-A] over a window with a width of δ min centered a time τ before the time at which d[TF-A]/d*t* is computed. For integration of *Eqs. 1
* and *
2
*, the program XPP (B. Ermentrout, University of Pittsburgh) was used, with the Gear integration method selected. For other simulations, the simple forward-Euler method was implemented in FORTRAN, with storage of delayed variables for use in subsequent calculations.

With zero delay, we have shown (42) that the genetic regulatory system described by *Eq. 1
* is bistable. There are two steady-state solutions, or fixed points, for [TF-A], which are asymptotically stable, in that if [TF-A] is slightly perturbed from these solutions, it will relax back. At these solutions, d[TF-A]/d*t* = 0. Therefore, a fixed point in the absence of delay will remain a fixed point when finite delay is assumed, and no new fixed points can be created by adding a delay. As Fig. 1
*B*illustrates, for*k*
_{f} between 0.6 and 2.5 min^{−1}, there are three values of [TF-A] that are steady states of the system. For the middle (unstable) steady state, a small perturbation of [TF-A] will grow until [TF-A] moves to either of the other two steady states. Dimerization of TF-A to activate its own transcription is essential for bistability in this system. In the variant of *Eq. 1
* where only first powers of [TF-A] are present, corresponding to activation of transcription by TF-A monomers, there is only a single nonzero fixed point.

Large perturbations of [TF-A] can switch this model between stable steady states. Such perturbations can be induced by brief changes in the strength*k*
_{f} by which [TF-A] activates its own transcription. A change in*k*
_{f} could correspond to a change in the fraction of phosphorylated TF-A. The resulting state transitions could correspond physiologically to short-lived stimuli, such as exposure to a neurotransmitter or hormone, leading to long-lasting changes in the levels of proteins. A brief change in the basal transcription rate, R_{bas}, of the*tf-a* promoter could be similarly interpreted and also can give a state transition. The model of Fig.1
*A* without delay would switch between states on a brief (∼10-min) perturbation (42). However, with a discrete delay τ, or with a narrowly distributed delay with a mean of τ, the perturbation must be applied for considerably longer than τ. This difference comes about because, with no delay, positive feedback of increased [TF-A] on TF-A synthesis can occur immediately and accentuate the effect of even a brief perturbation. With a delay, however, positive feedback cannot begin to act until the delay has passed. Thus inclusion of a delay does not change the property of multiple coexisting steady states, but the nature of the perturbations required to induce transitions and the nature of the transitions themselves are strongly affected.

Figure 2
*A* illustrates the nature of the transition between steady states of [TF-A] on an increase in *k*
_{f}. There is a series of five steps, each separated by the delay, caused by successive increases in positive feedback due to previous jumps in [TF-A]. The appearance of staircase steps does require that the time scale of [TF-A] equilibration be shorter than the transport delay τ, as reflected in the parameter values chosen for Fig. 2
*A*(equilibration time constant = 1/*k*
_{d} = 10 min). Otherwise, equilibration cannot occur rapidly enough to allow a sharp relaxation to a plateau of concentration, which accounts for the shape of the step. Indeed, in Fig. 2
*A*, one could as well plot the actual rate of*tf-a* transcription, which would have to itself be undergoing well-defined staircase steps, as newly synthesized TF-A protein is actively transported to the vicinity of the*tf-a* gene, where it can rapidly activate transcription. Such steps in transcription rate might be experimentally visualized by fluorescent in situ hybridization (10).

For *Eq. 2
* with parameter values as in Fig. 2
*A*, if the width δ of the distributed delay is less than ∼30 min, steps in [TF-A] are evident during the transition from the lower to the upper state, whereas if δ is greater than ∼45 min, only a smooth transition is seen (result not shown). Thus, if the distribution of times required for individual mRNA molecules to be transported and translated and for the corresponding TF-A protein molecules to move to the nucleus was not too broad, distinct steps in *tf-a*transcription rate might be expected.

The transient responses of the model of *Eqs.1
* and *
2
* are also state dependent. In Fig. 2
*B* a perturbation of *k*
_{f} at*arrow a* gives a large excursion of [TF-A] above its original value (the upper steady state of Fig. 1
*B*) at*t* = 210 min. This same perturbation would give only a minute excursion of [TF-A] above the lower steady state of Fig. 1
*B* (too small to see on the scale of Fig.2
*B*). This difference in response magnitude after a state transition is a model for “priming” a system to respond more vigorously to subsequent stimuli (42). In addition, the presence of a delay gives the system a type of “memory.” The memory can be described as follows: after a sustained change in the degree of TF-A phosphorylation (i.e., a change in *k*
_{f}), the response to a subsequent change in*k*
_{f} depends strongly on whether the interval between the two changes is less than or greater than the delay τ. Figure2
*B* also illustrates this point. An abrupt decrease in*k*
_{f} from 1 to 0.1 min^{−1} at*t* = 280 min gave a rapid transition of [TF-A] to a lower steady state at*t* = 400 min. If a brief increase in*k*
_{f} was applied within one delay time τ after the decrease, at *arrow b*, then a large excursion in [TF-A] resulted at *t* = 500 min, even though just before this excursion [TF-A] had decreased most of the way to the lower steady state. The same brief increase in*k*
_{f} applied after a time greater than τ had elapsed, at *arrow c*, gave an excursion in [TF-A] too small to see at *t* = 800 min. The memory is due to the delay for a change in mRNA synthesis to propagate to a corresponding change in nuclear [TF-A]. The large [TF-A] excursion at *t* = 500 min was due to a brief increase in*k*
_{f} at*t* = 380 min, at which time [TF-A] was still high. The high [TF-A] combined with the increased*k*
_{f} to strongly activate *tf-a* transcription, with the resulting large peak in [TF-A] not occurring until*t* = 500 min.

The response to a single perturbation can “echo” many times when a discrete delay characterizes transport. A brief increase in*k*
_{f}, insufficient to give a state transition, will elicit a brief perturbation in [TF-A]. However, this response, in turn, briefly increases the synthesis of *tf-a* mRNA, thus yielding a second, although diminished, increase in [TF-A]. Within a significant range of parameters, many successive echoes can be obtained. A larger brief initial perturbation can cause successive echoes that grow until the system undergoes a permanent state transition. Distinct echoes require that the time scale for [TF-A] equilibration be shorter than the delay τ, so that each perturbation relaxes close to baseline before the following perturbation occurs. As discussed above, this is more likely for nuclear than for overall cellular [TF-A], and also for the actual rate of *tf-a* transcription, which is here thought of as undergoing repeated distinct increases, as newly synthesized TF-A protein is actively transported to the vicinity of the *tf-a* gene. With*Eq. 2
* the echo perturbations are fairly robust to the width δ of the distributed delay. Typically, significant echoes remain evident, even with δ on the order of 45 min.

To examine whether the qualitative dynamic properties found above apply to a somewhat more detailed and realistic model and to allow for comparison with subsequent simulations that assume diffusive transport of mRNA and protein (see below), the model of*Eq. 1
* was extended to explicitly include *tf-a* mRNA. Translation was modeled as a simple first-order process, and a delay was included between mRNA transcription and changes in nuclear [TF-A]. The differential equations are
Equation 3
Equation 4The dynamics of this model are qualitatively similar to those of *Eq. 1
*. In particular, if the time scale characterizing the variation in*tf-a* mRNA concentration ([*tf-a* mRNA]) is rapid, the model reduces to that of *Eq. 1
*. For a relatively wide range of parameters, there is bistability. Perturbations in*k*
_{1,f} have to be of significant length to cause state transitions. Repeated echoes in mRNA and protein levels can be elicited by a single perturbation.

#### Stimulus-response properties of models with diffusion.

To compare the above dynamics with those produced by an analogous model with diffusive transport, we considered a cell divided into*N* spherical shells, with the innermost region, *shell 1*, a small sphere at the center (Fig.3
*A*). To obtain insights into the dynamics, it sufficed to take*N* = 10, because larger values of*N* did not alter the qualitative conclusions discussed below. Equations for diffusion with spherical symmetry (4) were used. In each shell, there are first-order degradation terms for mRNA and protein:*k*
_{RNA}mRNA_{N} and*k*
_{p}P_{N}, respectively. Transcription of mRNA is assumed to occur only in *shell 1* and translation of protein only in the outermost*shell N*, leading to the equations
Equation 5
Equation 6Integration of models with diffusion was done by the forward-Euler method.

A total cell radius of 10 μm was assumed for most of the simulations; thus a shell thickness was 1 μm. To determine plausible values for diffusion coefficients (*D*), we considered recent experimental observations of diffusion of dextran (39) and green fluorescent protein (50). Unless otherwise noted,*D* = 5 μm^{2}/s was used. Transcription depends on TF-A dimer, as in *Eq. 1
*; translation is simply proportional to the concentration of mRNA. Thus in *Eqs. 5
* and *
6
* we take
Equation 7
Equation 8where *k*
_{transcription} represents a maximal rate of transcription,*k*
_{translation} is a proportionality constant, R_{bas}represents a basal transcription rate in the absence of TF-A, and*K*
_{d} represents the dissociation constant for TF-A dimer binding to the*tf-a* promoter.

Bistability could be obtained (Fig.3
*B*) as for the version without diffusion (i.e., *Eqs. 3
* and *
4
*). For parameters as in Fig.3
*B* and for small R_{bas} (0.08 min^{−1}), the model is bistable for 3.5 <*k*
_{transcription}< 10. If concentrations are in the nanomolar range, these parameter values correspond to a maximal gene transcription rate of ∼10 mRNAs/min, which is reasonable for a strongly expressed gene (20, 47). Perturbations of appreciable length in*k*
_{transcription} or R_{bas} (Fig.3
*B*) were needed to bring about state transitions, because newly synthesized mRNA must be translated and the protein must reenter the nucleus before it can bind to DNA and initiate positive feedback. When state transitions did occur, the time course of protein and mRNA concentrations always remained smooth (Fig.3
*B*). There was no staircase pattern of steps to increasing levels. As with delay, bistability was observed only when dimerization of TF-A is assumed necessary for its effect on transcription.

We also examined responses to brief perturbations in R_{bas} and*k*
_{transcription}. With physiologically reasonable parameter values, a brief increase in either parameter gave a relatively abrupt increase in mRNA and protein levels followed by a slow decline (not shown). The time scale of the decline is set by the slow mRNA degradation rate. There was a significant lag of ∼30 min between peak mRNA and peak protein levels, which is due to the time scale of protein equilibration set by the slow protein degradation rate constant.

With diffusive transport, as opposed to a simple delay, the spreading out of any excess of mRNA or protein in one region disrupts the distinct peaks that are necessary for propagation of echo perturbations (i.e., the occurrence of more than one significant perturbation of protein or mRNA levels after a brief perturbation in a parameter). No echoes were observed with physiologically reasonable parameter values.

### Time Delays and Diffusion Suppress Oscillations Present When Macromolecular Transport Is Neglected; Time Delays, but Not Diffusion, Can Create New Oscillatory Modes

We compared the effect of diffusion with the effect of delay on the capability for stable oscillations in models of genetic regulation. These models incorporate positive and negative feedback via TF-As and TF-Rs.

#### Oscillations with only positive or only negative feedback.

In our earlier study (42), we noted that models analogous to*Eq. 1
* with only positive or only negative feedback do not have oscillatory solutions when the time taken for macromolecular transport is neglected. However, with a discrete delay, it has been proven that models with only negative feedback often possess oscillatory solutions (1). What if only positive feedback is present with a delay: can there be stable oscillations?

We are not aware of a general theorem for this case. However, we considered whether stable oscillations can be supported by*Eq. 1
* or by the simpler analog in which transcription is activated by monomeric TF-A. If oscillatory solutions came into existence as τ was increased from zero, they would be expected to appear via Hopf bifurcations from fixed points (46). Such a bifurcation requires a concomitant change in stability of the fixed point: stable to unstable or vice versa. As demonstrated in the
, however, the fixed points of *Eq. 1
* and the analog have the same stability properties with and without delay. In addition, our arguments appear to generalize (see) to rule out changes in stability of fixed points as τ is varied and creation of stable oscillatory solutions via Hopf bifurcations for all models in a class defined by three conditions. First, the models are represented by a single first-order differential equation comprised of a degradation term and a synthesis term for a single variable. Second, a discrete delay affects only the synthesis term. Third, the fixed points satisfy two not very restrictive constraints (see). Furthermore, if addition of a discrete delay cannot destabilize a fixed point, neither can it be destabilized by a biologically reasonable distributed delay (27). Thus*Eq. 2
* and any model in the above class are not expected to exhibit stable oscillations for distributed delay.

#### Oscillations with positive and negative feedback.

To investigate the effect of concurrent negative and positive feedback, we introduced into the model of Fig.1
*A* an additional gene,*tf-r*, the transcription rate of which is increased by binding of the TF-A dimer to a TF-RE. TF-R monomer represses transcription by competitively inhibiting binding of TF-A dimers to TF-REs. The equations for this model are
Equation 9
Equation 10Parameters in *Eqs. 9
* and *
10
* are analogous to those in*Eq. 1
*.*K*
_{R,d} is the dissociation constant of TF-R monomers from TF-REs. Robust oscillations are readily generated by this model when τ = 0 (42). Discrete delays of an order reasonable for macromolecular transport (τ ∼ 120 min) suppress this oscillatory pattern. However, a new oscillatory attractor with a period on the order of the delay is created. For example, if τ = 120 min, a lengthy transient of complex oscillations is seen (not shown). After ∼150 h, the transient evolves to a stable limit cycle. For distributed delay, if δ ≥ 60 min, these oscillations were abolished.

To contrast the above results with the dynamics obtained by assuming diffusional transport of TF-A and TF-R protein and mRNA, we modeled diffusion of these species within a spherically symmetrical cell. Transcription rates for *tf-a* and*tf-r* mRNA were given by
Equation 11
Equation 12Degradation of all macromolecular species is assumed to take place in every shell. Degradation of *tf-a* and*tf-r* mRNA is given by*k*
_{ARdeg}[*tf-a*mRNA]_{N} and*k*
_{RRdeg}[*tf-r*mRNA]_{N}, respectively. The rate of protein synthesis was directly proportional to the mRNA concentration in the outermost shell, and the rate of degradation of protein was directly proportional to protein concentration. The proportionality constants are*k*
_{synP} and*k*
_{degP}. With*D* of ∼5 μm^{2}/s, oscillatory behavior was readily obtained (Fig. 3
*C*). The period of the oscillations was ∼30 h. Further simulations found that slowing diffusion by lowering *D* always acts to suppress oscillations. For suppression, diffusion had to be slowed considerably, because the time scale for movement of macromolecules between nucleus and cytoplasm had to become comparable to that for other chemical processes (i.e., macromolecular degradation and synthesis). For the simulation of Fig.3
*C*, oscillations were not suppressed until *D* for all species was lowered from 5 to ∼0.2 μm^{2}/s.

We repeated the simulation of Fig. 3
*C*, reducing *D* between the second and third spherical shells to simulate the effect of a nuclear membrane. The modified coefficient was determined by using the equation that describes diffusion in the radial direction,*x*
^{2} = 4*Dt*/π, where*x* is the mean distance traveled in time *t* given*D* (6). We assumed that a molecule moves one shell width (1 μm), crossing the membrane, in a time on the order of 1 min. Therefore, we reduced*D* between the second and third spherical shells to 0.04 μm^{2}/s for protein and mRNA. This reduction sufficed to abolish oscillations.

### Dynamics of Two Specific Genetic Systems Are Strongly Influenced by Time Delays

#### Molecular processes underlying circadian oscillations.

Models incorporating negative-feedback loops have been used to simulate circadian rhythms (14, 25). Circadian “clocks” involve negative feedback of one or a few core genes on their own expression. In*Drosophila* and mammals, accumulation of PER and TIM proteins represses *per* and *tim*transcription. Later, the repression is relieved by PER and TIM degradation subsequent to multiple slow phosphorylation steps, so that the next cycle can begin. In the fungus*Neurospora* the FRQ protein exhibits analogous dynamics (36). Recent data should allow for the development of more comprehensive models for these circadian oscillators. In particular, prior models have assumed that phosphorylation of the core gene product must precede the onset of transcriptional repression (14,25). However, autorepression by FRQ protein of its own transcription occurs rapidly after *frq* induction and may therefore be independent of slow FRQ phosphorylation and degradation (36). Autorepression by PER may also be independent of slow PER phosphorylation (11). We have carried out preliminary simulations with a generic model of a circadian rhythm generator to examine whether independence of transcriptional repression from phosphorylation implies that any specific biochemical features are necessary for the production of robust circadian oscillations.

Multiple protein phosphorylations could proceed sequentially (e.g.,*phosphorylation 1* must precede*phosphorylation 2*) or independently. We simulated both possibilities. In the model schematized in Fig.4
*A*, a protein undergoes an obligatory series of sequential phosphorylations before degradation. A key aspect of the model is that all species of protein are equally capable of repressing their own transcription. Repression occurs by binding and sequestration of a putative transcriptional activator, the total concentration of which is a fixed parameter A_{tot}. The rate of transcription is a sigmoidal function of free activator with a Hill coefficient of 3, maximal velocity*v*
_{RNA}, and half-maximal transcription occurring at a free activator concentration*K*
_{act}. The rate of translation is proportional to mRNA concentration with proportionality constant *k*
_{trans}. No delay is assumed between synthesis of protein and its repression of transcription. However, a delay τ of 90 min is assumed between transcription and synthesis of protein. The states of increasing protein phosphorylation are denoted P_{0},...,P_{m}. The concentration of each phosphorylation state requires a separate differential equation; however, the kinetics of each phosphorylation are assumed identical (Fig. 4
*B*). The model simulated circadian oscillations as illustrated in Fig.4
*B*, which displays time courses of mRNA and of total protein concentration. Thus circadian rhythms can be generated by a model in which all species of protein are equally capable of repressing their own transcription. Circadian oscillations could only be simulated if many phosphorylation states (>8) were included. For the simulation of Fig.4
*B*, 10 phosphorylations were assumed.

Not only were many phosphorylations essential for obtaining circadian oscillations, but also the delay τ was required. Only relaxation to a steady state was seen if τ was removed. This result is not surprising, because the slow phosphorylations are not within the negative-feedback loop of protein synthesis followed by transcriptional repression. Oscillations within a pure negative-feedback loop require a delay within the loop or a very high Hill coefficient of feedback along with at least three variables in the loop (1).

The above simulations assuming sequential phosphorylations were compared with analogous simulations in which the individual phosphorylations were allowed to proceed independently of each other (not shown). Oscillations were also obtained in the latter case. However, in both cases, oscillations were always of a very short period unless many (∼10) slow phosphorylations were assumed necessary for degradation. Thus not only the delay τ but also the slow phosphorylations are required for simulation of long-period (∼24-h) circadian oscillations. Therefore, these simulations yield the specific prediction that many obligatory phosphorylations must precede protein degradation in the circadian clocks of*Neurospora* and*Drosophila*.

#### Molecular processes underlying LTF.

CREB protein can bind to CREs to induce the transcription of immediate-early genes crucial for LTF and the formation of long-term memory (LTM) (32, 48). Active transport of CREB protein appears to occur in neurites. Fluorescently tagged CREB protein perfused into dendrites of hippocampal neurons is rapidly and unidirectionally transported to the nucleus, whereas CREB protein perfused into the nucleus remains there (8). Transport of newly phosphorylated CREB protein from active neuronal synapses to the nucleus is hypothesized to provide a signal for the transcription of genes necessary for LTF (8). Active transport of CREB protein might be preferred over passive diffusion for signal transmission, because it could preserve concentration peaks. Suppose synaptic activity causes a brief episode of local CREB protein translation and phosphorylation. If CREB protein spread by passive diffusion throughout the neuronal volume, only a very small and slow change in its concentration might occur at the nucleus. However, we have seen that active transport, which might be modeled as a time delay, could preserve concentration peaks, so that a sharp peak of phosphorylated CREB protein might arrive at the nucleus in response to a brief episode of synaptic CREB protein translation. This peak of CREB protein could then initiate the transcription of genes important for synaptic modification.

To place these considerations on a more quantitative footing, transport of phosphorylated CREB was simulated within a simple neuronal geometry. A cubic “soma” with a volume of 1,000 μm^{3} was considered, with four attached “dendrites,” each 1 μm thick and 350 μm long. Hippocampal pyramidal cells can have dendrites of similar length and mean thickness (19). However, this simple morphology only suffices for a preliminary estimation of whether adding active transport to diffusive transport, both with rates similar to experimental data, is likely to significantly sharpen a pulse of messenger protein arriving in the soma after a synaptic event. More detailed modeling would be required for any specific cell type.

All dendrite lengths were measured along the*x* coordinate, dendrites were modeled as long boxes with *y* and*z* dimensions of 1 μm, and diffusive movement was allowed along *x*,*y*, and*z* coordinates.*D* in the range of 2–8 μm^{2}/s was assumed (39, 50); 2 μm^{2}/s is rather slow diffusion; however, diffusion in fine dendrites might be expected to be slowed by organelles and other obstructions. For each coordinate, spatial step lengths were chosen from a Gaussian distribution with a mean determined by *D* and the simulation time step d*t*; i.e., the mean step length*x*
_{avg} =
(6). Reflection of particles at boundaries was implemented. At each boundary, reflection only requires retrograde movement along the coordinate perpendicular to the boundary; the steps in the other two coordinates are unaltered. The retrograde movement was chosen, such that the sum of absolute values of the anterograde and retrograde movements gave the total step length that would have occurred in the absence of reflection. In some simulations, active transport was superimposed on diffusion. Active transport occurred along the*x* coordinate as a constant movement toward the soma of 300 μm/h, similar to some reported rates of fast mRNA transport within dendrites (45). CREB was assumed to be dephosphorylated with a time constant on the order of 10 h. CREB molecules were placed at the end of one dendrite, corresponding to the abrupt phosphorylation of CREB that might be engendered by a synaptic event. The arrival of phosphorylated CREB in the soma and its accumulation and dephosphorylation were simulated.

Figure 5
*B*compares simulations incorporating slow diffusion (*D* = 2 μm^{2}/s, *bottom trace*) with simulations incorporating active transport superimposed on slow diffusion (*top trace*). For generation of each of these traces, 50 CREB molecules were placed at the end of a dendrite at*t* = 20 h. It is evident that inclusion of active transport increased the peak level of somatic CREB by a factor of 2–5. The time required to reach the peak was decreased considerably, from ∼4 h to 2 h. However, for faster diffusion (*D* = 8 μm^{2}/s), inclusion of active transport produced only a modest improvement in the peak CREB level (a factor of ∼2; Fig. 5
*C*). The time to peak was decreased only slightly. These simulations indicate that the presence of active synapse-to-nucleus transport of a putative messenger protein, such as phosphorylated CREB, could considerably sharpen the nuclear signal for synaptic activity compared with slow diffusive transport. We note that this conclusion is sensitive to parameter variations within a physiologically reasonable range. Sometimes, if diffusion is rapid, the presence of active transport would give little advantage. However, the simple geometry used in these simulations underestimates the dilution of messenger protein that would be expected to occur by diffusion in neurons with complex morphologies in the absence of active transport. Also, if the dephosphorylation of CREB proceeded more rapidly, e.g., with a time constant of 3 h, the advantage of including active transport became greater (not shown). For slower diffusion (*D* = 2 μm^{2}/s), virtually no phosphorylated CREB arrived at the nucleus if diffusive transport alone was assumed.

### Stochastic Fluctuations in Molecule Numbers Can Mask the Existence of Stable Steady States

Stochastic fluctuations are generally significant in reacting systems when relatively small numbers of molecules are present (22, 34). A thorough analysis is difficult, because transport should be treated stochastically. Here, as a preliminary step, stochasticity was added to biochemical reactions in the model of *Eqs.3
* and *
4
* for synthesis and degradation of a TF that activates its own transcription, but transport was still modeled as a distributed delay. Given physiologically reasonable parameters that allow multiple stable steady states, do fluctuations due to stochasticity cause state transitions?

The variables [*tf-a* mRNA] and [TF-A] were reinterpreted as molecule numbers rather than as concentrations by rescaling parameter values. Deterministic rates for synthesis and degradation processes are given by
Equation 13
Equation 14
Equation 15
Equation 16The expression for R_{translation} incorporates a distributed delay for macromolecular transport. To include fluctuations due to association and dissociation of TF-A monomers, we dropped the assumption that TF-A dimer concentration is proportional to [TF-A]^{2}. Deterministic rates for association and dissociation were assigned R_{monomer-dimer} =*k*
_{f}([TF-A]_{monomer})^{2}and R_{dimer-monomer} =*k*
_{b}[TF-A]_{dimer}, respectively. The reciprocals of the deterministic rates are the average time intervals between reactions. Denote such a time interval by *T*
_{avg}. If a particular biochemical reaction occurs at*t* = 0, the probability*P*(*t*) that the next reaction of that type will occur within a specific short time interval Δ*t* centered at a later time *t* is (13)
Equation 17

At each time step Δ*t*, a separate random number was chosen for each elementary process of synthesis or degradation of mRNA or protein and TF-A association and dissociation. Each random number was drawn from a uniform distribution on {0,1}. For any elementary process, if the random number was less than the product of Δ*t* and the average rate (*Eqs. 13-16
*), we assumed that the process occurred once. If the products of Δ*t* and the average rates are kept small (<0.1), then this scheme approximates *Eq.17
*.

*Equations 13-16
* were also used to formulate ordinary differential equations for [*tf-a* mRNA], [TF-A]_{monomer}, and [TF-A]_{dimer}. Simulations with this “deterministic” model were compared with simulations that included fluctuations. The deterministic model exhibits bistability for a significant range of parameters. One stable solution has [*tf-a* mRNA] low and its synthesis rate close to R_{bas}, and the other has [*tf-a* mRNA] high and its synthesis rate close to*k*
_{1,f}. Figure6
*A*illustrates an example in which the deterministic model is stable in the lower steady state until a temporary increase in R_{bas} switches it to a new stable state. The stochastic variant of the model was then initialized near the lower steady state. Within ∼40 h, random fluctuations accumulated and enabled a transition to the upper state (Fig.6
*B*). This behavior occurred with physiologically reasonable parameter values, similar to those used by McAdams and Arkin (34). For example, in the simulation of Fig.6
*B*, the maximal mean transcription rate *k*
_{1,f} of 3.8 mRNA molecules per minute is reasonable for a strongly activated promoter. A reverse transition, from the upper to the lower state, was not seen even after 400 h of simulated time. Upward transitions within 10–100 h of simulated time, with no downward transitions after 400 h, were seen for five different random number sequences. Thus the lower stable state can be “masked” by fluctuations, in that the system would not ordinarily be observed near it. These results were essentially unchanged when δ was varied over the range 5–100 min. However, assuming a discrete delay (δ = 0) resulted in the appearance of much larger and spurious fluctuations due to amplification by fluctuations exactly one delay time earlier (not shown).

## DISCUSSION

We have considered the dynamics of model genetic regulatory systems in which TFs activate or repress their own transcription. Differences have been illustrated that depend on whether transport of mRNA and protein is described by diffusion or by active transport modeled as a time delay. Active transport was modeled by a narrowly distributed time delay. For example, a time delay leads to a type of memory. After a sustained change in transcription rate, the amplitude of the excursion in protein concentration due to a brief change in transcription rate depends strongly on whether the interval between the two changes is less than or greater than the delay (Fig.2
*B*). This memory is due to the delay for a change in mRNA synthesis to propagate through translation to a corresponding change in nuclear protein level. Also, given appropriate parameters (large *k*
_{f} and *k*
_{d} in*Eq. 1
*), models with time delays can exhibit staircase-like transitions from one steady state of macromolecular concentrations to another in response to a sustained stimulus (Fig.2
*A*) or repeated perturbations of concentrations in response to a single brief stimulus. If newly synthesized TF was actively transported to the vicinity of its own gene, then staircase transitions in transcription rate or repeated perturbations of this rate after a single stimulus might be visualized by fluorescent in situ hybridization (10). Observation of such phenomena would constitute strong evidence for active transport of mRNA and protein. It does not seem that a system reliant on diffusional transport would exhibit either of these phenomena.

The dynamics exemplified by Fig. 2 might help determine the relative efficacies of different training paradigms in forming LTM (33, 42). If transcription of CREB or of another TF that activates its own transcription is essential for plasticity underlying the formation of LTM, then stimuli separated by short time intervals might be expected to be relatively ineffective at forming LTM. A longer time interval after the first stimulus might be required to allow for TF transcription, translation, and translocation to the nucleus. Thus a larger amount of nuclear TF would be available at the arrival of the second stimulus, and, when activated, that TF could strongly promote the transcription of itself and of other genes essential for LTM formation. Such considerations could help explain the usual superiority of temporally spaced, as opposed to massed, training sessions in forming LTM.

Two issues arise concerning the validity and importance of models incorporating different transport mechanisms. First, what conditions might justify modeling macromolecular transport with a time delay? Second, how useful is this type of modeling, in conjunction with experiments, for understanding the dynamics of genetic systems?

### Can Active or Diffusive Transport Be Plausibly Modeled as a Narrowly Distributed Time Delay?

Numerous examples of active mRNA and protein transport in the cytoplasm are now known, including microtubule-dependent movement of*Vg1* mRNA in*Xenopus* oocytes (49), the myosin-dependent segregation of the lineage-determining PAR-1, PAR-2, and PAR-3 proteins in *Caenorhabditis elegans* (15), microtubule-dependent localization of several maternal mRNAs for proteins that direct embryonic development (such as *bicoid, bicuadal-D*, and*oskar*) in*Drosophila* oocytes (24), and retrograde transport of proteins with a nuclear localization signal in*Aplysia* neurons (40).

Recent experimental studies of active transport by motor proteins suggest that the distribution of times taken by individual macromolecules to be transported a given intracellular distance could be quite narrow. The kinetics of motion along microtubules have recently been examined by use of latex beads with single kinesin molecules attached. The motion is comprised of single steps ∼8 nm long (18, 41). The distribution of distances traveled in a given time can be characterized by a randomness parameter*r*, defined as the variance divided by the product of the mean and the single-step distance (41). At >1 mM ATP, similar to that in vivo, *r* ≅ ½ (41). The variance can be estimated for typical intracellular distances of transport. With use of the definition of*r*, if a population of macromolecules moves a mean distance of 20 μm and if the movement is driven by kinesin with a single step length of 8 nm, then with*r* ≅ ½ the variance in the distance moved by individual molecules will be 0.08 μm^{2}. A variance as small as 0.08 μm^{2} could be modeled as a discrete time delay, because for small variances the ratio of the variance in arrival times to the mean travel time is the same as the ratio of the variance in position to the mean distance moved. In contrast, analogous calculations demonstrate that diffusive transport could not be modeled as a time delay, because the variance in the distance moved by individual macromolecules is much greater.

However, it is not generally known whether successive mRNAs from a given gene are transported to the same translation site or closely grouped sites. If successive mRNAs were directed to a variety of translation sites at different distances from the gene, there would be a distribution of times to reach the sites. Also, it can be expected that some diffusional movement would be required even if active transport predominated, because each mRNA would probably need to move along more than one cytoskeletal element to reach its destination, and “jumps” between elements might be via diffusion. A model formulation using a distributed delay to describe transport might be an appropriate approximation for such complications. In our simulations, oscillations were abolished if the width of the distribution was 20–30% of the average delay. Repeated perturbations and transitions via concentration steps were not abolished until somewhat larger distribution widths. Whether distributions this narrow are reasonable characterizations of intracellular transport is an issue requiring experimental investigation.

### Specific Issues for Further Investigation by Modeling and Experiment

Our modeling of stochastic fluctuations in genetic regulatory systems, which illustrates that such fluctuations could preferentially destabilize and mask the existence of steady states with low concentrations of macromolecules, could certainly be further developed. One issue of interest is whether fluctuations tend to be significantly smaller or to have a different frequency spectrum when diffusional transport dominates over active transport. It appears that if simulations such as those of McAdams and Arkin (34) are to predict the degrees of variability in the behavior of actual genetic systems, transport and perhaps its stochastic nature will often have to be included. However, the inclusion of transport may be less necessary in prokaryotic systems such as those modeled by McAdams and Arkin. If transport is diffusional, the small dimensions of many prokaryotic cells and the lack of a nuclear membrane imply that, on the time scale of minutes, concentrations of species not rapidly degraded can be regarded as homogenous.

Transport that is partly diffusive and partly active could be described in a more quantitative, spatially resolved manner by use of a stochastic ordinary differential (Langevin) equation in combination with a partial differential (Fokker-Planck) equation (12). A Langevin equation can be formulated to describe the average and fluctuating components of motion of an actively transported macromolecule. The corresponding Fokker-Planck equation describes the evolution of the spatial distribution of an ensemble of macromolecules. It has a convective term for the mean flow due to active transport and a diffusive term for thermal fluctuations. An additional term could be added to describe ordinary diffusion. Given estimates of kinetic parameters, numerical simulations using this equation could predict the variance in distances traveled by members of an ensemble of macromolecules. Related statistical quantities, such as a correlation function for overlap of concentration perturbations due to separate genetic induction events, could also be predicted.

On a more general level, the diversity of transcription factors and their interactions suggest that behaviors such as those considered here (e.g., long-lasting state transitions in response to perturbations or oscillations dependent on time delays) will be identified. Thus the dynamic principles illustrated here are likely to be important in phenomena where regulation of transcription has an essential role, such as development or the formation of LTM. It has recently been proposed (28) that epigenetic, heritable changes in gene expression after exposure to chemicals might play a role in carcinogenesis. Such changes correspond dynamically to perturbations of genetic regulatory systems from one steady state to another, such as we have modeled. An outstanding issue will be to determine whether the parameters of particular genetic systems in vivo are permissive for specific types of dynamic behavior. To this end, we suggest that the rate of movement and the variance in distance moved in a given time should be examined for a few TFs of particular importance and for their mRNA transcripts. This can help determine whether diffusional or active transport dominates in particular systems. We believe that as the dynamic behaviors of gene networks are explored empirically the present work can provide a conceptual framework for the interpretation of such experiments.

## Acknowledgments

We thank C. Canavier and R. Butera for comments on the manuscript.

## Appendix

Here we demonstrate the mathematical points stated inresults. First, we show that the fixed points of *Eq. 1
* and of its analog where first powers of [TF-A] replace second powers have the same stability properties with discrete delay as without. In the analog, nondimensionalization of [TF-A] and time successively sets*K*
_{d} and*k*
_{f} to 1. R_{bas} is also assumed to be small and has little effect on the dynamics about the nonzero fixed point (simulations have verified this for several combinations of parameter values) and can be neglected. The analog then reduces to
Equation 18Here we have introduced a notation we use throughout the. For a variable*X*, a discrete delay τ, and present time *t*,*X*
_{del} ≡*X*(*t*− τ), whereas *X* ≡*X*(*t*).*Equation 18
* has a single nonzero asymptotically stable fixed point at*X* = (1 −*a*)/*a*.

Biological relevance requires *X* > 0 at this fixed point, so *a* < 1. If τ ≥ 0, then in general, if*g*(*X,X*
_{del}) denotes the right-hand side of a first-order delay differential equation, a characteristic equation governing stability of the fixed point is derived by linearizing the differential equation about the fixed point and then substituting*X*(*t*) = exp(λ*t*) into the linearized equation (see Ref. 27 for details)
Equation 19The derivatives in *Eq. 19
*are to be evaluated at the fixed point. At the nonzero fixed point of*Eq. 18, Eq. 19
* gives*a*
^{2}exp(−λτ) − *a* − λ = 0. If the real part of the eigenvalue λ is negative (respectively, positive), the fixed point is stable (respectively, unstable). A switch in stability is only possible when the real part of λ is 0 for some τ > 0.

Substituting a pure imaginary λ = *ki*into *Eq. 19
* and evaluating at the fixed point of *Eq. 18
* give two simultaneous equations
Using the Pythagorean identity gives
which cannot be satisfied if*a* < 1, as is required for*X* > 0 at the fixed point. Thus the fixed point remains stable for all τ.

For *Eq. 1
*, nondimensionalization of [TF-A] and *t* and neglect of R_{bas} gives
Equation 20
*X* = 0 is a fixed point of*Eq. 20
*, and in addition, there are two positive fixed points if 0 < *a* < ½. If τ = 0 the upper positive fixed point is stable, the lower is unstable. If τ > 0, we again look for a purely imaginary λ that solves *Eq. 19
*. For λ =*ki*, we have
∂*g*/∂*X*
_{del}must be greater than −(∂*g*/∂*X*) to admit a real solution (*k*,τ). For*Eq. 20
*, −(∂*g*/∂*X*) = *a*, and for*a* < ½, numerical calculation shows that ∂*g*/∂*X*
_{del}> −(∂*g*/∂*X*) for only the lower positive fixed point. Thus only this fixed point could switch stability as τ varies. The other fixed points must remain stable for all τ.

We consider further whether the lower nonzero fixed point could become stable for any τ. Suppose at *time 0*an upward perturbation in *X* is made from this point, small enough to remain within the neighborhood where ∂*g*/∂*X*
_{del}> −(∂*g*/∂*X*) (a similar argument applies for downward perturbations). The perturbation is free until time τ to relax toward the fixed point. Denote the value of *X* at τ by*X*
_{1}. At time τ,*X* will begin to increase. This occurs because the positive contribution to d*X*/d*t*,
/(
+ 1), abruptly becomes greater than the negative contribution*aX*
_{1}. Later, if*X* “tries” to decrease below*X*
_{1}, it will be unable to, because
/(
+ 1) will remain greater than*aX*
_{1} (since*X*
_{del} >*X*
_{1}). We now choose a value *X* =*X*
_{2}, above*X*
_{1} and the peak of the small perturbation but still within the region where ∂*g*/∂*X*
_{del}> −(∂*g*/∂*X*). Suppose *X* never increases above*X*
_{2} for subsequent*t*. If so, one can average d*X*/d*t*from time τ to time *t* and, since*X* must always remain in the interval (*X*
_{1},*X*
_{2}), this average must approach zero for large*t*.

This average equals the difference of the average of*X*
^{2}/(*X*
^{2}+ 1) from *time 0* to time (*t* − τ) and the average of*aX* from time τ to time*t*. For sufficiently large*t*, one can neglect the relatively small intervals from 0 to τ and from (*t* − τ) to*t* when taking these averages. Thus the average of [*X*
^{2}/(*X*
^{2}+ 1)] − *aX* from time τ to time *t* must approach zero for large *t*. However,*X* lies within (*X*
_{1},*X*
_{2}) from time τ to time *t*. We recall that *1*)*g*(*X,X*
_{del}) ≡ [*X*
^{2}/(*X*
^{2}+ 1)] − *aX* = 0 at the lower fixed point and *2*) ∂*g*/∂*X*
_{del}> −(∂*g*/∂*X*) for all *X* between this fixed point and*X*
_{2}. These conditions imply that [*X*
^{2}/(*X*
^{2}+ 1)] − *aX* is positive and monotonically increasing from*X*
_{1} to*X*
_{2}, which in turn implies that the average of [*X*
^{2}/(*X*
^{2}+ 1)] − *aX* from time τ to time *t* cannot approach zero for large *t*. This contradiction implies that, rather than staying within (*X*
_{1},*X*
_{2}) for all *t*,*X* must eventually increase above*X*
_{2} and above its original upward perturbation from the lower nonzero fixed point. Thus this fixed point is unstable, regardless of the value of τ. We conclude that none of the fixed points of *Eq.20
* can switch stability as τ varies.

The arguments of the preceding three paragraphs appear, for a class of models, to rule out switches of stability of fixed points when the length of the discrete delay is varied. This class satisfies the following conditions: *1*) there is a single dependent variable *X*, and d*X*/d*t*≡*g*(*X,X*
_{del}) is the sum of two terms, a synthesis term dependent only on*X*
_{del} and a degradation term dependent only on *X*;*2*) all fixed points are either asymptotically stable or unstable in the absence of delay; and*3*) ∂*g*/∂*X*
_{del}> −(∂*g*/∂*X*) [respectively, ∂*g*/∂*X*
_{del}< −(∂*g*/∂*X*)] at those fixed points that are unstable (respectively, asymptotically stable) in the absence of delay.

## Footnotes

Address for reprint requests and other correspondence: J. H. Byrne, Dept. of Neurobiology and Anatomy, W. M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas-Houston Medical School, PO Box 20708, Houston, TX 77225 (E-mail:jbyrne{at}nba19.med.uth.tmc.edu).

This work was supported by National Institutes of Health Grants T32 NS-07373 and R01 RR-11626 and by Texas Higher Education Coordination Board Grant 011618-048.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

- Copyright © 1999 the American Physiological Society