Effects of macromolecular transport and stochastic fluctuations on dynamics of genetic regulatory systems

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

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


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

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


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

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 beta -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 Ca2+/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
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

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 Ca2+/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 sigma 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. 1A). 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 kf, which is proportional to TF-A phosphorylation. At negligible dimer concentration, the synthesis rate is Rbas. TF-A is eliminated with a rate constant kd. 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 kf and kd, which set the time scale for [TF-A] equilibration, are fairly rapid (e.g., kd = 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 kf and kd 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 tau  tau 1 + tau 2, with tau 1 the time taken for the transcription of tf-a mRNA and its movement to translation and tau 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.


View larger version (20K):
[in this window]
[in a new window]
 
Fig. 1.   A: phosphorylated dimers of TF-A activate tf-a transcription when bound to specific DNA sequences (TF-REs). Degradation is also indicated (kd). B: bistability in model in A. For 0.6 min-1 < kf < 2.5 min-1, 2 stable steady-state solutions of [TF-A] exist (top and bottom of d[TF-A]/dt = 0 curve) with an unstable solution between them (middle, dashed portion of curve). Outside this region there is a single steady-state solution. Other parameters are as follows: Rbas = 0.01 min-1, kd = 0.1 min-1, and Kd = 10 nM2. See Glossary and RESULTS for definition of abbreviations.



View larger version (15K):
[in this window]
[in a new window]
 
Fig. 2.   Responses of model with discrete delays. A: model of Fig. 1A exhibits a "staircase" transition between states. Starting from a lower steady state of Eq. 2 with Rbas = 0.02 min-1, kd = 0.2 min-1, Kd = 10 nM2, tau  = 120 min, and delta  = 20 min; kf is increased at t = 200 min from 2 to 20 min-1. After approximately the delay tau , a small step in [TF-A] occurs. This is followed by successive steps to a new steady state. B: state-dependent responses to perturbations in model of Fig. 1A. Equation 2 is initially in upper steady state of Fig. 1B, with kf = 1 min-1, tau  = 120 min, and delta  = 20 min. At t = (210 - tau ) min (arrow a), kf is increased to 10 min-1 for 2 min, generating a large increase in [TF-A] at t = 210 min. At t = (400 - tau ) min, kf is decreased to a new baseline of 0.1 min-1 (bar). Thus, at t = 400 min, [TF-A] decreases most of the way toward a new, low steady state. At t = (500 - tau ) min (arrow b), kf is increased to 10 min-1 for 2 min. A large increase in [TF-A] results at t = 500 min. At t = (800 - tau ) min (arrow c), kf is again increased to 10 min-1 for 2 min. An imperceptible increase in [TF-A] occurs at t = 800 min. * Small increases in [TF-A] that are "echoes" of larger excursions (see results).

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
<FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR>
= <FENCE><FR><NU><IT>k</IT><SUB>f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>d</SUB></DE></FR></FENCE> (<IT>t</IT> − &tgr;) − <IT>k</IT><SUB>d</SUB>[TF-A] + R<SUB>bas</SUB> (1)
with Kd the dissociation constant of dimer from TF-REs. The first term (in braces) on the right-hand side is evaluated at a time tau  previous to the time when d[TF-A]/dt 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]/dt 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
<FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>k</IT><SUB>f</SUB>⟨[TF-A]⟩<SUP>2</SUP></NU><DE>⟨[TF-A]⟩<SUP>2</SUP> + <IT>K</IT><SUB>d</SUB></DE></FR> − <IT>k</IT><SUB>d</SUB>[TF-A] + R<SUB>bas</SUB> (2)
where the angle brackets average [TF-A] over a window with a width of delta  min centered a time tau  before the time at which d[TF-A]/dt 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]/dt = 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. 1B illustrates, for kf 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 kf by which [TF-A] activates its own transcription. A change in kf 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, Rbas, of the tf-a promoter could be similarly interpreted and also can give a state transition. The model of Fig. 1A without delay would switch between states on a brief (~10-min) perturbation (42). However, with a discrete delay tau , or with a narrowly distributed delay with a mean of tau , the perturbation must be applied for considerably longer than tau . 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 2A illustrates the nature of the transition between steady states of [TF-A] on an increase in kf. 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 tau , as reflected in the parameter values chosen for Fig. 2A (equilibration time constant = 1/kd = 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. 2A, 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. 2A, if the width delta  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 delta  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. 2B a perturbation of kf at arrow a gives a large excursion of [TF-A] above its original value (the upper steady state of Fig. 1B) at t = 210 min. This same perturbation would give only a minute excursion of [TF-A] above the lower steady state of Fig. 1B (too small to see on the scale of Fig. 2B). 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 kf), the response to a subsequent change in kf depends strongly on whether the interval between the two changes is less than or greater than the delay tau . Figure 2B also illustrates this point. An abrupt decrease in kf 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 kf was applied within one delay time tau  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 kf applied after a time greater than tau  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 kf at t = 380 min, at which time [TF-A] was still high. The high [TF-A] combined with the increased kf 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 kf, 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 tau , 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 delta  of the distributed delay. Typically, significant echoes remain evident, even with delta  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
<FR><NU>d[<IT>tf-a</IT> mRNA]</NU><DE>d<IT>t</IT></DE></FR>
= <FR><NU><IT>k</IT><SUB>1,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>d</SUB></DE></FR> − <IT>k</IT><SUB>1,d</SUB>[<IT>tf-a</IT> mRNA] + R<SUB>bas</SUB> (3)
<FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR>
= <IT>k</IT><SUB>2,f</SUB>[<IT>tf-a</IT> mRNA](<IT>t</IT> − &tgr;) − <IT>k</IT><SUB>2,d</SUB>[TF-A] (4)
The 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 k1,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. 3A). 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: kRNA mRNAN and kp PN, 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
<FR><NU>dmRNA<SUB>1</SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>dmRNA<SUB>1</SUB></NU><DE>d<IT>t</IT><SUB>diffusion</SUB></DE></FR> − <IT>k</IT><SUB>RNA</SUB> mRNA<SUB>1</SUB> + R<SUB>transcription</SUB> (5)
<FR><NU>dP<SUB><IT>N</IT></SUB></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU>dP<SUB>N</SUB></NU><DE>d<IT>t</IT><SUB>diffusion</SUB></DE></FR> − <IT>k</IT><SUB>P</SUB>P<SUB><IT>N</IT></SUB> + R<SUB>translation</SUB> (6)
Integration of models with diffusion was done by the forward-Euler method.


View larger version (18K):
[in this window]
[in a new window]
 
Fig. 3.   Dynamics of model based on transport modeled as diffusion. A: schematic of model. Transcription of mRNA takes place in innermost spherical shell, and translation takes place in outermost shell (shell boundaries indicated by dashed circles). B: at t = 100 min, increasing Rbas from 0.08 to 3.75 nM/min for 140 min (bar) gives a smooth state transition to a new stable state. Other parameters are as follows: ktranscription = 3.75 nM/min, kRNA = 0.004 min-1, kp = 0.08 min-1, Kd = 10 nM2, ktranslation = 2.0 min-1, N = 10, and DRNA = DP = 5 µm2/s. C: oscillations of TF-A protein and mRNA in model incorporating diffusion as well as positive and negative feedback via Eqs. 11 and 12. Concentrations in nuclear region are shown. Parameter values are as follows: N = 10, ksynA = 10.5 nM/min, Rbas = 0.4 nM/min, kARdeg = 0.01 min-1, ksynR = 3.0 nM/min, kRRdeg = 0.002 min-1, ksynP = 2.22 min-1, kdegP = 0.06 min-1, Kd = 10 nM2, KR,d = 2.0 nM, and DRNA = DP = 5 µm2/s.

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 µm2/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
R<SUB>transcription</SUB> = <FR><NU><IT>k</IT><SUB>transcription</SUB> [TF-A]<SUP>2</SUP><SUB>1</SUB></NU><DE>[TF-A]<SUP>2</SUP><SUB>1</SUB> + <IT>K</IT><SUB>d</SUB></DE></FR> + R<SUB>bas</SUB> (7)
R<SUB>translation</SUB> = <IT>k</IT><SUB>translation</SUB> [<IT>tf-a</IT> mRNA]<SUB><IT>N</IT></SUB> (8)
where ktranscription represents a maximal rate of transcription, ktranslation is a proportionality constant, Rbas represents a basal transcription rate in the absence of TF-A, and Kd represents the dissociation constant for TF-A dimer binding to the tf-a promoter.

Bistability could be obtained (Fig. 3B) as for the version without diffusion (i.e., Eqs. 3 and 4). For parameters as in Fig. 3B and for small Rbas (0.08 min-1), the model is bistable for 3.5 < ktranscription < 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 ktranscription or Rbas (Fig. 3B) 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. 3B). 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 Rbas and ktranscription. 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 tau  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 APPENDIX, 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 APPENDIX) to rule out changes in stability of fixed points as tau  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 APPENDIX). 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. 1A 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
<FR><NU>d[TF-A]</NU><DE>d<IT>t</IT></DE></FR>
= <FENCE><FR><NU><IT>k</IT><SUB>1,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>1,d</SUB> (1 + [TF-R]/<IT>K</IT><SUB>R,d</SUB>)</DE></FR> </FENCE> (<IT>t</IT> − &tgr;)
− <IT>k</IT><SUB>1,d</SUB>[TF-A] + R<SUB>bas</SUB> (9)
<FR><NU>d[TF-R]</NU><DE>d<IT>t</IT></DE></FR>
= <FENCE><FR><NU><IT>k</IT><SUB>2,f</SUB>[TF-A]<SUP>2</SUP></NU><DE>[TF-A]<SUP>2</SUP> + <IT>K</IT><SUB>2,d</SUB>(1 + [TF-R]/<IT>K</IT><SUB>R,d</SUB>)</DE></FR> </FENCE> (<IT>t</IT> − &tgr;)
− <IT>k</IT><SUB>2,d</SUB>[TF-R] (10)
Parameters in Eqs. 9 and 10 are analogous to those in Eq. 1. KR,d is the dissociation constant of TF-R monomers from TF-REs. Robust oscillations are readily generated by this model when tau  = 0 (42). Discrete delays of an order reasonable for macromolecular transport (tau  ~ 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 tau  = 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 delta  >=  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
R<SUB>synA</SUB> = <FR><NU><IT>k</IT><SUB>synA</SUB>[TF-A]<SUP>2</SUP><SUB>1</SUB></NU><DE>[TF-A]<SUP>2</SUP><SUB>1</SUB> + <IT>K</IT><SUB>d</SUB> (1 + [TF-R]<SUB>1</SUB>/<IT>K</IT><SUB>R,d</SUB>)</DE></FR> + R<SUB>bas</SUB> (11)
R<SUB>synR</SUB> = <FR><NU><IT>k</IT><SUB>synR</SUB>[TF-A]<SUP>2</SUP><SUB>1</SUB></NU><DE>[TF-A]<SUP>2</SUP><SUB>1</SUB> + <IT>K</IT><SUB>d</SUB>(1 + [TF-R]<SUB>1</SUB>/<IT>K</IT><SUB>R,d</SUB>)</DE></FR> (12)
Degradation of all macromolecular species is assumed to take place in every shell. Degradation of tf-a and tf-r mRNA is given by kARdeg[tf-a mRNA]N and kRRdeg[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 ksynP and kdegP. With D of ~5 µm2/s, oscillatory behavior was readily obtained (Fig. 3C). 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. 3C, oscillations were not suppressed until D for all species was lowered from 5 to ~0.2 µm2/s.

We repeated the simulation of Fig. 3C, 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, x2 = 4Dt/pi , 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 µm2/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. 4A, 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 Atot. The rate of transcription is a sigmoidal function of free activator with a Hill coefficient of 3, maximal velocity vRNA, and half-maximal transcription occurring at a free activator concentration Kact. The rate of translation is proportional to mRNA concentration with proportionality constant ktrans. No delay is assumed between synthesis of protein and its repression of transcription. However, a delay tau  of 90 min is assumed between transcription and synthesis of protein. The states of increasing protein phosphorylation are denoted P0,...,Pm. The concentration of each phosphorylation state requires a separate differential equation; however, the kinetics of each phosphorylation are assumed identical (Fig. 4B). The model simulated circadian oscillations as illustrated in Fig. 4B, 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. 4B, 10 phosphorylations were assumed.


View larger version (19K):
[in this window]
[in a new window]
 
Fig. 4.   Circadian oscillations generated by a model incorporating sequential protein phosphorylations. A: schematic of model. Transcription of a circadian gene is inhibited by its protein product P. Parameter values not given in text are as follows. All phosphorylations are irreversible and first order with identical rate constant kphos = 3.0 h-1, and degradation of fully phosphorylated protein is irreversible and Michaelis-Menten with maximal reaction velocity (Vmax) = 1.7 µM/h and Michaelis-Menten constant (Km) = 0.02 µM. Degradation of mRNA is Michaelis-Menten with Vmax = 12 µM/h and Km = 0.8 µM. Binding of P to TF-A is bimolecular with forward and reverse rate constants 20 µM-1 · h-1 and 14 h-1, respectively. Other parameters are as follows: vRNA = 40 µM/h, ktrans = 0.3 µM/h, Atot = 0.21 µM, and Kact = 0.1 µM. B: circadian oscillations with period close to 24 h are simulated with parameter values in A.

Not only were many phosphorylations essential for obtaining circadian oscillations, but also the delay tau  was required. Only relaxation to a steady state was seen if tau  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 tau  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 µm3 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 µm2/s was assumed (39, 50); 2 µm2/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 dt; i.e., the mean step length xavg = <RAD><RCD>(4<IT>D</IT>/&pgr;)d<IT>t</IT></RCD></RAD> (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 5B compares simulations incorporating slow diffusion (D = 2 µm2/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 µm2/s), inclusion of active transport produced only a modest improvement in the peak CREB level (a factor of ~2; Fig. 5C). 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 µm2/s), virtually no phosphorylated CREB arrived at the nucleus if diffusive transport alone was assumed.


View larger version (21K):
[in this window]
[in a new window]
 
Fig. 5.   A: schematic of model neuron used to simulate transport of protein from an active synaptic zone, assumed to be located at end of a dendrite, to nucleus. B: time courses for number of protein molecules present in nucleus after placement of 50 molecules at end of a dendrite. At t = 20 h, molecules were placed and allowed to move. Movement was by diffusion combined with active transport (top curve) or by diffusion alone (bottom curve). Relatively slow diffusion (D = 2 µm2/s) was assumed. Population of molecules decayed with a time constant of 10 h. Simulation time step was fixed at 1 s. At each time step and for each molecule, a random number was chosen between 0 and 1, and if that number was less than deterministic probability of molecule decaying, as determined from time constant and time step, that molecule was eliminated from simulation. C: same as B, except D is raised to 8 µm2/s.

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
R<SUB>transcription</SUB> = <FR><NU><IT>k</IT><SUB>1,f</SUB>[TF-A]<SUB>dimer</SUB></NU><DE>[TF-A]<SUB>dimer</SUB> + <IT>K</IT><SUB>d</SUB></DE></FR> + R<SUB>bas</SUB> (13)
R<SUB>degRNA</SUB> = <IT>k</IT><SUB>1,d</SUB>[<IT>tf-a</IT> mRNA] (14)
R<SUB>translation</SUB> = <IT>k</IT><SUB>2,f</SUB>⟨[<IT>tf-a</IT> mRNA]⟩ (15)
R<SUB>degP</SUB> = <IT>k</IT><SUB>2,d</SUB>[TF-A]<SUB>monomer</SUB> (16)
The expression for Rtranslation 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 Rmonomer-dimer = kf([TF-A]monomer)2 and Rdimer-monomer = kb[TF-A]dimer, respectively. The reciprocals of the deterministic rates are the average time intervals between reactions. Denote such a time interval by Tavg. 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 Delta t centered at a later time t is (13)
<IT>P</IT>(<IT>t</IT>) = <FR><NU>1</NU><DE><IT>T</IT><SUB>avg</SUB></DE></FR> exp <FENCE><FR><NU>− <IT>t</IT></NU><DE><IT>T</IT><SUB>avg</SUB></DE></FR> </FENCE> &Dgr;<IT>t</IT> (17)

At each time step Delta 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 Delta t and the average rate (Eqs. 13-16), we assumed that the process occurred once. If the products of Delta 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 Rbas, and the other has [tf-a mRNA] high and its synthesis rate close to k1,f. Figure 6A illustrates an example in which the deterministic model is stable in the lower steady state until a temporary increase in Rbas 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. 6B). This behavior occurred with physiologically reasonable parameter values, similar to those used by McAdams and Arkin (34). For example, in the simulation of Fig. 6B, the maximal mean transcription rate k1,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 delta  was varied over the range 5-100 min. However, assuming a discrete delay (delta  = 0) resulted in the appearance of much larger and spurious fluctuations due to amplification by fluctuations exactly one delay time earlier (not shown).


View larger version (25K):
[in this window]
[in a new window]
 
Fig. 6.   Stochastic fluctuations in molecule numbers can destabilize steady states of genetic regulatory systems. A: without stochastic fluctuations, deterministic model (Eqs. 13-16) is bistable. Initial levels of RNA and protein are low (<1 nM) and steady. Parameters are as follows: k1,f = 0.2 nM/min, k1,d = 0.02 min-1, k2,f = 0.2 nM/min, k2,d = 0.03 min-1, Rbas = 0.0015 nM/min, kf = 2 nM/min, kb = 20 min-1, tau  = 90 min, and Kd = 10 nM. At t = 10 h (arrow), Rbas is increased to 0.1 nM/min for 500 min. Increase causes a transition to a stable state of high RNA and protein levels. B: initial state of A is spontaneously destabilized when stochastic fluctuations are incorporated. Parameters as in A, except 1) delta  = 5 min and 2) to convert nanomolar concentrations to molecule numbers within a nuclear region assumed to have a 2-µm radius, k1,f, Rbas, and Kd were multiplied by a numerical factor of 19.2 nM-1 and kf was divided by this factor.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

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. 2B). 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 kf and kd 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. 2A) 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 congruent  1/2 (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 congruent  1/2 the variance in the distance moved by individual molecules will be 0.08 µm2. A variance as small as 0.08 µm2 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.


    APPENDIX
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

Here we demonstrate the mathematical points stated in RESULTS. 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 Kd and kf to 1. Rbas 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
<FR><NU>d<IT>X</IT></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>X</IT><SUB>del</SUB></NU><DE><IT>X</IT><SUB>del</SUB> + 1</DE></FR> − <IT>aX</IT> (18)
Here we have introduced a notation we use throughout the APPENDIX. For a variable X, a discrete delay tau , and present time t, Xdel triple-bond  X(t - tau ), whereas X triple-bond  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 tau  >=  0, then in general, if g(X,Xdel) 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(lambda t) into the linearized equation (see Ref. 27 for details)
<FR><NU>∂<IT>g</IT></NU><DE>∂<IT>X</IT><SUB>del</SUB></DE></FR> exp (− &lgr;&tgr;) + <FR><NU>∂<IT>g</IT></NU><DE>∂<IT>X</IT></DE></FR> − &lgr; = 0 (19)
The derivatives in Eq. 19 are to be evaluated at the fixed point. At the nonzero fixed point of Eq. 18, Eq. 19 gives a2exp(-lambda tau - a - lambda  = 0. If the real part of the eigenvalue lambda  is negative (respectively, positive), the fixed point is stable (respectively, unstable). A switch in stability is only possible when the real part of lambda  is 0 for some tau  > 0.

Substituting a pure imaginary lambda  = ki into Eq. 19 and evaluating at the fixed point of Eq. 18 give two simultaneous equations
cos (<IT>k</IT>&tgr;) = <FR><NU>1</NU><DE><IT>a</IT></DE></FR>
sin (<IT>k</IT>&tgr;) = − <FR><NU><IT>k</IT></NU><DE><IT>a</IT><SUP>2</SUP></DE></FR>
Using the Pythagorean identity gives
<FR><NU><IT>k</IT><SUP>2</SUP></NU><DE><IT>a</IT><SUP>4</SUP></DE></FR> = 1 − <FR><NU>1</NU><DE><IT>a</IT><SUP>2</SUP></DE></FR>
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 tau .

For Eq. 1, nondimensionalization of [TF-A] and t and neglect of Rbas gives
<FR><NU>d<IT>X</IT></NU><DE>d<IT>t</IT></DE></FR> = <FR><NU><IT>X</IT><SUP>2</SUP><SUB>del</SUB></NU><DE><IT>X</IT><SUP>2</SUP><SUB>del</SUB> + 1</DE></FR> − <IT>aX</IT> (20)
X = 0 is a fixed point of Eq. 20, and in addition, there are two positive fixed points if 0 < a < 1/2. If tau  = 0 the upper positive fixed point is stable, the lower is unstable. If tau  > 0, we again look for a purely imaginary lambda  that solves Eq. 19. For lambda  = ki, we have
cos (<IT>k</IT>&tgr;) = − <FR><NU>∂<IT>g</IT></NU><DE>∂<IT>X</IT></DE></FR><FENCE><FR><NU>∂<IT>g</IT></NU><DE>∂<IT>X</IT><SUB>del</SUB></DE></FR></FENCE>
sin (<IT>k</IT>&tgr;) = − <IT>k</IT><FENCE><FR><NU>∂<IT>g</IT></NU><DE>∂<IT>X</IT><SUB>del</SUB></DE></FR></FENCE>
partial g/partial Xdel must be greater than -(partial g/partial X) to admit a real solution (k,tau ). For Eq. 20, -(partial g/partial X) = a, and for a < 1/2, numerical calculation shows that partial g/partial Xdel > -(partial g/partial X) for only the lower positive fixed point. Thus only this fixed point could switch stability as tau  varies. The other fixed points must remain stable for all tau .

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

This average equals the difference of the average of X2/(X2 + 1) from time 0 to time (t - tau ) and the average of aX from time tau  to time t. For sufficiently large t, one can neglect the relatively small intervals from 0 to tau  and from (t - tau ) to t when taking these averages. Thus the average of [X2/(X2 + 1)] - aX from time tau  to time t must approach zero for large t. However, X lies within (X1,X2) from time tau  to time t. We recall that 1) g(X,Xdel) triple-bond  [X2/(X2 + 1)] - aX = 0 at the lower fixed point and 2) partial g/partial Xdel > -(partial g/partial X) for all X between this fixed point and X2. These conditions imply that [X2/(X2 + 1)] - aX is positive and monotonically increasing from X1 to X2, which in turn implies that the average of [X2/(X2 + 1)] - aX from time tau  to time t cannot approach zero for large t. This contradiction implies that, rather than staying within (X1,X2) for all t, X must eventually increase above X2 and above its original upward perturbation from the lower nonzero fixed point. Thus this fixed point is unstable, regardless of the value of tau . We conclude that none of the fixed points of Eq. 20 can switch stability as tau  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 dX/dt triple-bond  g(X,Xdel) is the sum of two terms, a synthesis term dependent only on Xdel 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) partial g/partial Xdel > -(partial g/partial X) [respectively, partial g/partial Xdel < -(partial g/partial X)] at those fixed points that are unstable (respectively, asymptotically stable) in the absence of delay.


    ACKNOWLEDGEMENTS

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


    FOOTNOTES

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.

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

Received 4 December 1998; accepted in final form 20 May 1999.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
RESULTS
DISCUSSION
REFERENCES
APPENDIX

1.   An der Heiden, U. Delays in physiological systems. J. Math. Biol. 8: 345-364, 1979[Medline].

2.   Banks, H., and J. Mahaffy. Global asymptotic stability of certain models for protein synthesis and repression. Q. Appl. Math. 36: 209-221, 1978.

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

4.   Blumenfeld, H., L. Zablow, and B. Sabatini. Evaluation of cellular mechanisms for modulation of Ca2+ transients using a mathematical model of fura-2 Ca2+ imaging in Aplysia sensory neurons. Biophys. J. 63: 1146-1164, 1992[Abstract].

5.   Busenberg, S., and J. Mahaffy. Interaction of spatial diffusion and delays in models of genetic control by repression. J. Math. Biol. 22: 313-333, 1985[Medline].

6.   Cantor, C., and P. R. Schimmel. Biophysical Chemistry. San Francisco, CA: Freeman, 1980, vol. 2.

7.   Collado-Vides, J., B. Magasanik, and J. D. Gralla. Control site location and transcriptional regulation in Escherichia coli. Microbiol. Rev. 55: 371-394, 1991.

8.   Crino, P., K. Khodakhah, K. Becker, S. Ginsberg, S. Hemby, and J. Eberwine. Presence and phosphorylation of transcription factors in developing dendrites. Proc. Natl. Acad. Sci. USA 95: 2313-2318, 1998[Abstract/Free Full Text].

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

10.   Femino, A. M., F. Fay, K. Fogarty, and R. H. Singer. Visualization of single RNA transcripts in situ. Science 280: 585-590, 1998[Abstract/Free Full Text].

11.   Gekakis, N., L. Saez, A. Delahaye-Brown, M. Myers, A. Sehgal, M. W. Young, and C. J. Weitz. Isolation of timeless by PER protein interaction: defective interaction between timeless protein and long-period mutant PERL. Science 270: 811-815, 1995[Abstract].

12.   Gillespie, D. T. The multivariate Langevin and Fokker-Planck equations. Am. J. Phys. 64: 1246-1257, 1996.

13.   Gillespie, D. T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 61: 2340-2361, 1977.

14.   Goldbeter, A. A model for circadian oscillations in the Drosophila period (PER) protein. Proc. R. Soc. Lond. B Biol. Sci. 261: 319-324, 1995[Medline].

15.   Guo, S., and K. J. Kemphues. Par-1, a gene required for establishing polarity in C. elegans embryos, encodes a putative Ser/Thr kinase that is asymmetrically distributed. Cell 81: 611-620, 1995[Medline].

16.   Hevroni, D., A. Rattner, M. Bundman, D. Lederfein, A. Gabarah, M. Mangelus, M. A. Silverman, H. Kedar, C. Naor, M. Kornuc, T. Hanoch, R. Seger, L. E. Theill, E. Nedivi, G. Richter-Levin, and Y. Citri. Hippocampal plasticity involves extensive gene induction and multiple cellular mechanisms. J. Mol. Neurosci. 10: 75-98, 1998[Medline].

17.   Hirokawa, N. The mechanisms of fast and slow transport in neurons: identification and characterization of the new kinesin superfamily motors. Curr. Opin. Neurobiol. 7: 605-614, 1997[Medline].

18.   Hua, W., E. Young, M. Fleming, and J. Gelles. Coupling of kinesin steps to ATP hydrolysis. Nature 388: 390-393, 1997[Medline].

19.   Ishizuka, N., W. M. Cowan, and D. G. Amaral. A quantitative analysis of the dendritic organization of pyramidal cells in the rat hippocampus. J. Comp. Neurol. 6: 17-45, 1995.

20.   Iyer, V., and K. Struhl. Absolute mRNA levels and transcriptional initiation rates in Saccharomyces cerevisiae. Proc. Natl. Acad. Sci. USA 93: 5208-5212, 1996[Abstract/Free Full Text].

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

22.   Keizer, J. Statistical Thermodynamics of Equilibrium Processes. New York: Springer-Verlag, 1987.

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

24.   King, M. L. Molecular basis for cytoplasmic localization. Dev. Genet. 19: 183-189, 1996[Medline].

25.   Leloup, J. C., and A. Goldbeter. A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J. Biol. Rhythms 13: 70-87, 1998[Abstract/Free Full Text].

26.   Li, M. G., M. McGrail, M. Serr, and T. S. Hays. Drosophila cytoplasmic dynein, a microtubule motor that is asymmetrically localized in the oocyte. J. Cell Biol. 126: 1475-1494, 1994[Abstract].

27.   MacDonald, N. Biological Delay Systems: Linear Stability Theory. Cambridge, UK: Cambridge University Press, 1989.

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

29.   Mahaffy, J. Cellular control models with linked positive and negative feedback and delays. I. The models. J. Theor. Biol. 106: 89-102, 1984[Medline].

30.   Mahaffy, J., D. Jorgensen, and R. van der Heyden. Oscillations in a model of repression with external control. J. Math. Biol. 30: 669-691, 1992[Medline].

31.   Mahaffy, J., and C. V. Pao. Models of genetic control by repression with time delays and spatial effects. J. Math. Biol. 20: 39-57, 1984[Medline].

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

33.   Mauelshagen, J., C. M. Sherff, and T. J. Carew. Differential induction of long-term synaptic facilitation by spaced and massed applications of serotonin at sensory neuron synapses of Aplysia californica. Learn. Mem. 5: 246-256, 1998[Abstract/Free Full Text].

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

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

36.   Merrow, M. W., N. Garceau, and J. C. Dunlap. Dissection of a circadian oscillation into discrete domains. Proc. Natl. Acad. Sci. USA 94: 3877-3882, 1997[Abstract/Free Full Text].

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

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

39.   Sabry, J., T. O'Connor, and M. W. Kirschner. Axonal transport of tubulin in Ti1 Pioneer neurons in situ. Neuron 14: 1247-1256, 1995[Medline].

40.   Schmied, R., and R. T. Ambron. A nuclear localization signal targets proteins to the retrograde transport system, thereby evading uptake into organelles in Aplysia axons. J. Neurobiol. 33: 151-160, 1997[Medline].

41.   Schnitzer, M. J., and S. M. Block. Kinesin hydrolyses one ATP per 8-nm step. Nature 388: 386-389, 1997[Medline].

42.   Smolen, P., D. A. Baxter, and J. H. Byrne. Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory systems. Am. J. Physiol. 274 (Cell Physiol. 43): C531-C42, 1998[Abstract/Free Full Text].

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

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

45.   Wallace, C. S., G. L. Lyford, P. F. Worley, and O. Stewart. Differential intracellular sorting of immediate-early gene mRNAs depends on signals in the mRNA sequence. J. Neurosci. 18: 26-35, 1998[Abstract/Free Full Text].

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

47.   Wolffe, A. P., J. Glover, S. Martin, M. Tenniswood, J. Williams, and J. Tata. Deinduction of transcription of Xenopus albumin genes and destabilization of mRNA by estrogen in vivo and in hepatocyte cultures. Eur. J. Biochem. 146: 489-496, 1985[Abstract].

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

49.   Yisraeli, J. K., S. Sokol, and D. A. Melton. A two-step model for the localization of maternal mRNA in Xenopus oocytes: involvement of microtubules and microfilaments in the translocation and anchoring of Vg1 mRNA. Development 108: 289-298, 1990[Abstract].

50.   Yokoe, H., and T. Meyer. Spatial dynamics of GFP-tagged proteins investigated by local fluorescence enhancement. Nat. Biotechnol. 14: 1252-1256, 1996[Medline].


Am J Physiol Cell Physiol 277(4):C777-C790
0002-9513/99 $5.00 Copyright © 1999 the American Physiological Society