Stochastic Delay Accelerates Signaling in Gene Networks

The creation of protein from DNA is a dynamic process consisting of numerous reactions, such as transcription, translation and protein folding. Each of these reactions is further comprised of hundreds or thousands of sub-steps that must be completed before a protein is fully mature. Consequently, the time it takes to create a single protein depends on the number of steps in the reaction chain and the nature of each step. One way to account for these reactions in models of gene regulatory networks is to incorporate dynamical delay. However, the stochastic nature of the reactions necessary to produce protein leads to a waiting time that is randomly distributed. Here, we use queueing theory to examine the effects of such distributed delay on the propagation of information through transcriptionally regulated genetic networks. In an analytically tractable model we find that increasing the randomness in protein production delay can increase signaling speed in transcriptional networks. The effect is confirmed in stochastic simulations, and we demonstrate its impact in several common transcriptional motifs. In particular, we show that in feedforward loops signaling time and magnitude are significantly affected by distributed delay. In addition, delay has previously been shown to cause stable oscillations in circuits with negative feedback. We show that the period and the amplitude of the oscillations monotonically decrease as the variability of the delay time increases.


Introduction
Gene regulation forms a basis for cellular decision-making processes and transcriptional signaling is one way in which cells can modulate gene expression patterns [1]. The intricate networks of transcription factors and their targets are of intense interest to theorists because it is hoped that topological similarities between networks will reveal functional parallels [2]. Models of gene regulatory networks have taken many forms, ranging from simplified Boolean networks [3,4], to full-scale, stochastic descriptions simulated using Gillespie's algorithm [5].
The majority of models, however, are systems of nonlinear ordinary differential equations (ODEs). Yet, because of the complexity of protein production, ODE models of transcriptional networks are at best heuristic reductions of the true system, and often fail to capture many aspects of network dynamics. Many ignored reactions, like oligomerization of transcription factors or enzyme-substrate binding, occur at much faster timescales than reactions such as transcription and degradation of proteins. Reduced models are frequently obtained by eliminating these fast reactions [6][7][8][9]. Unfortunately, even when such reductions are done correctly, problems might still exist. For instance, if within the reaction network there exists a linear (or approximately linear) sequence of reactions, the resulting dynamics can appear to be delayed. This type of behavior has long been known to exist in gene regulatory networks [10].
Delay differential equations (DDEs) have been used as an alternative to ODE models to address this problem. In protein production, one can think of delay as resulting from the sequential assembly of first mRNA and then protein [10][11][12]. Delay can qualitatively alter the local stability of genetic regulatory network models [13] as well as their dynamics, especially in those containing feedback. For instance, delay can lead to oscillations in models of transcriptional negative feedback [11,[14][15][16][17][18], and experimental evidence suggests that robust oscillations in simple synthetic networks are due to transcriptional delay [19,20].
Protein production delay times are difficult to measure in live cells, though recent work has shown that the time it takes for transcription to occur in yeast can be on the order of minutes and is highly variable [21]. Still, transcriptional delay is thought to be important in a host of naturally occurring gene networks. For instance, mathematical models suggest that circadian oscillations are governed by delayed negative feedback systems [22,23], and this was experimentally shown to be true in mammalian cells [24]. Delay appears to play a role in cell cycle control [25,26], apoptosis induction by the p53 network [27], and the response of the NFkB network [15]. Delay can also affect the stochastic nature of gene expression, and the relation between the two can be subtle and complex [28][29][30][31].
In this study, we examine the consequences of randomly distributed delay on simple gene regulatory networks: We assume that the delay time for protein production, T, is not constant but instead a random variable. If f T denotes the probability density function (PDF) of T, this situation can be described deterministically by an integro-delay differential equation [32] of the form _ x x(t)~ð g x(t),x(t{s) where x is a positive definite state vector of protein concentrations, and g is a vector function representing the production and degradation rates of the proteins. Note that processes that do not require protein synthesis (like dilution and degradation) will depend on the instantaneous, rather than the delayed, state of the system. Therefore g is in general a function of both the present and past state of the system. Equation (1) only holds in the limit of large protein numbers [32]. As protein numbers approach zero, the stochasticity associated with chemical interactions becomes non-negligible. Here, we address this issue by expanding on Eq. (1) using an exact stochastic algorithm that takes into account variability within the delay time [32]. We further use a queueing theory approach to examine how this variability affects timing in signaling cascades. We find that when the mean of the delay time is fixed, increased delay variability accelerates downstream signaling. Noise can thus increase signaling speed in gene networks. In addition, we find that in simple transcriptional networks containing feed-forward or feedback loops the variability in the delay time nontrivially affects network dynamics.
Queueing theory has recently been used to understand the behavior of genetic networks [33][34][35]. Here we are mainly interested in dynamical phenomena to which the theory of queues in equilibrium used in previous studies cannot be applied. As we explain below, gene networks can be modeled as thresholded queueing systems: Proteins exiting one queue do not enter another queue, as would be the case in typical queueing networks. Rather, they modulate the rate at which transcription is initiated, and thus affect the rate at which proteins enter other queues.

Distributed delay in protein production
The transcription of genetic material into mRNA and its subsequent translation into protein involves potentially hundreds or thousands of biochemical reactions. Hence, detailed models of these processes are prohibitively complex. When simulating genetic circuits it is frequently assumed that gene expression instantaneously results in fully formed proteins. However, each step in the chain of reactions leading from transcription initiation to a folded protein takes time ( Figure 1). Models that do not incorporate the resulting delay may not accurately capture the dynamical behavior of genetic circuits [17]. While earlier models have included either fixed or distributed delay [32,36,37], here we examine specifically the effects of delay variability on transcriptional signaling.
In one recent study, Bel et al. studied completion time distributions associated with Markov chains modeling linear chemical reaction pathways [38]. Using rigorous analysis and numerical simulations they show that, if the number of reactions is large, completion time distributions for an idealized class of models exhibit a sharp transition in the coefficient of variation (CV, defined as the standard deviation divided by the mean of the distribution), going from near 0 (indicating a nearly deterministic completion time) to near 1 (indicating an exponentially distributed completion time) as system bias moves from forward to reverse.
However, it is possible, and perhaps likely, that the limiting distributions described by Bel et al. do not provide good approximations for protein production. For instance, when the number of rate limiting reactions is small, but greater than one, the distribution of delay times can be more complex. Moreover, linear reaction pathways only represent one possible and necessarily simplified reaction scheme. Protein production involves many reaction types that are nonlinear and/or reversible, each of which is influenced by intrinsic and extrinsic noise [39], and these Numerous reactions must occur between the time that transcription starts and when the resulting protein molecule is fully formed and mature. Though we call this phenomenon ''transcriptional'' delay, there are many reactions after transcription (such as translation) which contribute to the overall delay. (B) The creation of multiple proteins can be thought of as a queueing process. Nascent proteins enter the queue (an input event) and emerge fully matured (an output event) some time later depending on the distribution of delay times. Because the delay is random, it is possible that the order of proteins entering the queue is not preserved upon exit. (C) In a transcriptionally regulated signaling process the time it takes for changes in the expression of gene 1 to propagate to gene 2 depends on both the distribution of delay times, f T , and the number of transcription factors needed to overcome the threshold of gene 2, R. doi:10.1371/journal.pcbi.1002264.g001

Author Summary
Delay in gene regulatory networks often arises from the numerous sequential reactions necessary to create fully functional protein from DNA. While the molecular mechanisms behind protein production and maturation are known, it is still unknown to what extent the resulting delay affects signaling in transcriptional networks. In contrast to previous studies that have examined the consequences of fixed delay in gene networks, here we investigate how the variability of the delay time influences the resulting dynamics. The exact distribution of ''transcriptional delay'' is still unknown, and most likely greatly depends on both intrinsic and extrinsic factors. Nevertheless, we are able to deduce specific effects of distributed delay on transcriptional signaling that are independent of the underlying distribution. We find that the time it takes for a gene encoding a transcription factor to signal its downstream target decreases as the delay variability increases. We use queueing theory to derive a simple relationship describing this result, and use stochastic simulations to confirm it. The consequences of distributed delay for several common transcriptional motifs are also discussed.
reactions may impact the delay time distribution in complicated ways. Therefore, we do not try to derive the actual shape of f T , but examine the effects its statistical properties have on transcriptional signaling. To do this, we represent protein production as a delayed reaction of the form g ?
where g is the gene, and transcription is initiated at rate l(t,P), which can depend explicitly on both time and protein number, P. After initiation, it takes a random time, T, for a protein to be formed. Note that the presence of time delay implies that scheme (2) defines a non-Markovian process. Such processes can be simulated exactly using an extension of the Gillespie algorithm (See Methods and [28,32]). If the biochemical reaction pathway that leads to functional protein is known and relatively simple, direct stochastic simulation of every step in the network is preferable to simulation based on scheme (2). From the point of view of multi-scale modeling, however, paradigm (2) is useful when the biochemical reaction network is either extremely complex or poorly mapped, since one needs to know only the statistical properties of T.

Protein formation as a queueing system
In the setting of scheme (2), first assume that l(t,P) does not depend on P, and protein formation is initiated according to a memoryless process with rate l(t). A fully formed protein enters the population a random time T after the initiation of protein formation. We assume that the molecules do not interact while forming; that is, the formation of one protein does not affect that of another. Each protein therefore emerges from an independent reaction channel after a random time. This process is equivalent to an M=G=? queue [40], where M indicates a memoryless source (transcription initiation), G a general service time distribution (delay time distribution), and ? refers to the number of service channels.
In our model, the order in which initiation events enter a queue is not necessarily preserved. As Figure 1 (B) illustrates, it is possible for the initiation order to be permuted upon exit [32]. The assumption that proteins can ''skip ahead'' complicates the analysis of transient dynamics of such queues, and is essential in much of the following. While there are steps where such skipping can occur (such as protein folding), there are others for which it cannot. For instance, it is unlikely that one RNA polymerase can skip ahead of another -and similarly for ribosomes during translation off of the same transcript. Therefore, protein skipping may be more relevant in eukaryotes, where transcription and translation must occur separately, than prokaryotes, where they may occur simultaneously. However, if there is more than one copy of the gene (which is common for plasmid-based synthetic gene networks in E. coli), or more than one transcript, some skipping is likely occur. Therefore it is likely that the full results that follow are more relevant for genes of copy number greater than one.

Downstream transcriptional signaling
One purpose of transcription factors is to propagate signals to downstream target genes. Determining the dynamics and stochasticity of these signaling cascades is of both theoretical and experimental interest [41,42]. Therefore, we first examine the impact that distributed delay has on simple downstream signaling. Consider the situation depicted in Figure 1 (C), in which the product of the first gene regulates the transcription of a second gene. Using the same nomenclature as in scheme (2) we write where g 1 and g 2 are the copy numbers of the upstream and downstream genes, P 1 and P 2 are the number of functional proteins of each type, and T i is the random delay time of gene i~1,2. The transcription rate of gene 2 depends on P 1 and is given by a Hill function h(P 1 ). We consider the case in which P 1 activates g 2 (depicted in Figure 1) and the case in which P 1 represses g 2 . We now ask: If P 1 starts at zero and gene 1 is suddenly turned on, how long does it take until the signal is detected by gene 2? In other words, assume l(t)~l 0 H(t), where H(t) is the Heaviside step function. At what time t does P 1 (t) reach a level that is detectable by gene 2? In order to make the problem tractable, we assume that the Hill function is steep and switch-like, so that we can make the approximation Here h 0 is the maximum transcriptional initiation rate of g 2 and Rw0 is the threshold value of the Hill function, i.e. the number of molecules of P 1 needed for half repression (or half activation) of gene 2. The second gene therefore becomes repressed (or activated) at the time S R at which R copies of protein P 1 have been fully formed. We first examine reaction (3a). Assume that at time t~0 there are no proteins in the system. Let I 1 (t) denote the number of transcription initiation events that have occurred by time t (the arrival process of the queueing system), Q 1 (t) the number of proteins being formed at time t (the size of the queue at time t), and P 1 (t) the number of functional proteins that have been completed by time t (the exit process of the queueing system). Since the arrival process is memoryless, (I 1 (t)) t §0 is a Poisson process with constant rate l 0 for t §0. Hence, the expected value of I 1 (t) is E I 1 (t) ½ ~l 0 tH(t). The exit process, i.e. the number of fully functional proteins that have emerged from the queue, P 1 (t), is a nonhomogenous Poisson process with time-dependent rate l P1 (t)~l 0 F T1 (t), where F T1 (t) is the cumulative distribution function (CDF) of the delay time T 1 . It then follows that Inactivation (or activation) of gene 2 occurs when enough protein P 1 has accumulated to trigger a transcriptional change, according to Eq. (4a) or (4b). In other words, the random time it takes for the signal to propagate, S R , is given by S Rm inftjP 1 (t)~Rg. Trivially, S R changes by an amount identical to a change in the mean of the delay distribution. To examine the effects of randomness in delay on the signaling time, we therefore keep the mean of the delay distribution fixed, E½T 1 ~t, and vary Var½T 1 .
The probability density function of S R is given by (See Methods) Consequently, the mean and variance of the time it takes for the original signal to propagate to the downstream gene can be written as: To gain insight into the behaviors of Eqs. (6) and (7), we first examine a representative, analytically tractable example. Assume that the delay time can take on 2 discrete values, tza and t{a with equal probability. In this case, where C(x,y) is the upper incomplete gamma function. Expanding for small a, we obtain (See Methods) which is the deterministic limit. The first term is the mean delay time and the second is the average time to initiate R proteins at rate l 0 . A similar expansion for fixed R and large a gives (see panel (c) in Figure 2) It follows that for larger delay variability, the mean signaling time decreases with delay variability (See Figure 2 (A)). Indeed, Eqs. (9) and (10) form the asymptotic boundaries for the mean signaling time. The intersection of the two asymptotes at a~R=l 0 , gives an estimate of when the behavior of the system changes from the deterministic limit (for avR=l 0 ) to a regime in which increasing the variability decreases the mean signaling time (for awR=l 0 ). It follows that the deterministic approximation given by Eq. (9) is valid in an increasing range, as R grows (See Figure 2 (C)). Indeed, an asymptotic analysis of Eq. (8) shows that the corrections to Eq. (9) are approximately of size (a=R) R , and therefore rapidly decrease with R (See Methods). The bottom row of Figure 2 shows that these observations hold more generally: When T 1 is gamma distributed the mean time to produce R proteins, E½S R , is very sensitive to randomness in delay time, but only when R is small to intermediate. As expected, the densities of the times to produce R proteins, f SR , are approximately normal and independent of the delay distribution when R is large (Middle panels of Figure 2).
We therefore expect that for each fixed threshold R, E½S R is a decreasing function of the standard deviation s T1 of the delay. We have proved this to be true for symmetric delay distributions (See Methods). Intuitively, this is due to the fact that the order in which proteins enter the queue is not the same as the order in which they exit. Proteins that enter the queue before the Rth protein, but exit after the Rth protein increase S R , while the opposite is true for proteins that enter the queue after the Rth protein, and exit before it. Since only finitely many proteins enter the queue before the Rth protein, while infinitely many enter after it, the balance favors a decrease in the mean signaling time. Moreover, as delay variability increases, interchanges in exit order become more likely, and this effect becomes more pronounced. We outline the analytical argument: For each fixed time t §0, P 1 (t) is an increasing function of s T1 , hence P {1 1 (s) is decreasing function of s T1 for all s §0. Referring to Eq. (6), this implies that E½S R is a decreasing function of s T1 in the symmetric case.
In sum, mean signaling times decrease as delay variability increases (with fixed mean delay). This effect is most significant for small to moderate thresholds. We note that the decrease in mean signaling time phenomenon depends on a sufficient number of proteins entering the queue. If transcription is only active long enough for less than 2R{1 proteins to be initiated, then mean signaling time will actually increase as delay variability increases. This phenomenon is explained in the subsection of the Methods section that analyzes repressor switches.

Example: Feedforward loops
Using the above results, we now examine more complicated transcriptional signaling networks. In particular, we turn to two common feedforward loops -the type 1 coherent and the type 1 incoherent feedforward loops (FFL) [43], shown in Figure 3. Each of these networks is a transcriptional cascade resulting in the specific response of the output, gene 3. The coherent FFL generally acts as a delayed response network, while the incoherent FFL has various possible responses, such as pulsatile response [43], response time acceleration [44], and fold-change detection [45].
To examine the effect of distributed delay on these networks we assume that at t~0 gene g 1 starts transcription of protein P 1 at rate b 1 , i.e l 1 (t)~b 1 H(t). The second gene, g 2 , starts transcription after P 1 (t) reaches the threshold R 12 , so that l 2 (t)b 2 H(P 1 (t){R 12 ). For the coherent FFL, we assume that the promoter of gene g 3 acts as an AND gate so that l 3 (t)b 3 H(P 1 (t){R 13 )H(P 2 (t){R 23 ). We further assume that the promoter of g 3 in the incoherent FFL is active only in the presence of P 1 and absence of P 2 , so that we may write l 3 (t)~b 3 H(P 1 (t){R 13 )H(R 23 {P 2 (t)).
The signaling time between any two nodes i and j within the network, i.e. the random time between the initiation of transcription of gene i to the formation of a total of R ij proteins P i is denoted S ij . For each of the three pathways, the PDF of the signaling time is given by Eq. (5). In addition, because the random times S 12 and S 23 are additive (as are their variances), we can directly calculate the time at which P 2 reaches the threshold of gene 3 as S 123~S12 zS 23 . Therefore, the random time at which the coherent FFL turns on is simply given by t onm ax S 13 ,S 123 f g . Because E½S 13 and E½S 123 are decreasing functions of the delay variability, it can be expected that so is E½t on .
In contrast to the coherent FFL, the dynamics of the pulse generating incoherent FFL are less trivial. Since the repressor (P 2 ) overrides the activator (P 1 ), assuming S 123 §S 13 transcription of g 3 turns on at time S 13 and turns off at time S 123 , generating a pulse of duration Z on~S123 {S 13 . Note that E½Z on can increase or decrease as a function of the standard deviation s of the delay (see Figure 4 where s was equal for all pathways).
To see this, write E½Z on as follows: Each of the terms on the right side of Eq. (11) is the expected signaling time of a single gene (g 1 ?g 2 , g 2 ?g 3 , and g 1 ?g 3 , respectively). Consequently, E½Z on depends on s as a linear combination of 3 expected signaling time curves of the type pictured in Figure 2. The shapes of these signaling time curves determine the behavior of E½Z on as a function of s. Figure 4 shows that the behavior of the duration of the transcriptional pulse as a function of the delay variability depends on the values of each threshold within the network.

The delayed negative feedback oscillator
These observations can also be extended to networks with recurrent architectures. For instance, consider the transcriptional delayed negative feedback circuit [17], which can be described using an extension of scheme (2): where h(P) is a decreasing Hill function (i.e. P represses its own production) and m(P) is the degradation rate due to dilution and proteolysis. Mather et al. examined the oscillations produced by systems of the type described by scheme (12) when the delay T is nonrandom (degrade and fire oscillators) [17]. Starting with no proteins, P is produced at a rate governed by the Hill function h. When the level of P exceeds the midpoint of the Hill function, gene g effectively shuts down. The proteins remaining in the queue exit,  producing a spike, after which degradation diminishes P. When the protein level drops sufficiently, reaction (12a) reactivates and production of P resumes, commencing another oscillation cycle. Note that this circuit will not oscillate without delay. As a result during each oscillation the gene is turned on until its own signal reaches itself, at which time the gene is turned off [17]. Therefore, the peak height of one oscillation is determined by the length of time the gene was in the ''ON'' state. Since that time is determined by the gene's signaling time, our theory predicts that the mean peak height of the oscillations will decrease as the variability in the delay time increases. Indeed, this is exactly what our stochastic simulations show in Figure 5. This is consistent with the fact that the negative feedback circuit is dynamically similar to the g 2 {g 3 subcircuit within the incoherent FFL. Here we explicitly used a gammadistributed delay time with mean t~3min, h(P)~aC n 0 = C 0 z ð P(t)Þ n and m(P)~bP(t)zc R P(t)= R 0 zP(t) ð Þ . We can use our theory to predict the change in the peak height of the oscillator as a function of s. For a delay that is gamma-distributed, the change in signaling time as a function of s can be written as where E s ½S R is given by Eq. (6) and f (s) is the reduction in the expected signaling time. If we assume that the amount of time that protein is produced during a burst in the delayed negative feedback oscillator is also reduced by this amount, then it is possible to predict the change in the peak height accordingly. To a first approximation, if the promoter is in the ''ON'' state for a time that is f (s) less, then a total of af (s) less protein will be produced. Therefore we can write the expected peak height of the oscillator as However, due to degradation, Eq. (14) overestimates the correction to the peak height. Due to exponential degradation, only a fraction exp ({bf (s)) of the lost protein would have made it through to the peak. Also, the duration of enzymatic decay is also reduced by a time f (s). Therefore, if we assume that the enzymatic decay reaction is saturated, we need to add an amount f (s)c R to Eq. (14). This gives us a more accurate prediction of the mean peak height as Figure 5 shows that this approximation works well, even for a Hill coefficient as low as 2.

Discussion
The existence of delay in the production of protein has been known of for some time. For many systems its presence does not seriously impact performance. For example, the existence of fixed points in simple downstream regulatory networks without feedback is unaffected by delay. Delay is important if the timing of signal propagation impacts the function of the network. Delay can also change a network's dynamics. In networks with feedback, for instance, delay can result in bifurcations that are not present in the corresponding non-delayed system. The delayed negative feedback oscillator is a prime example [17]. Moreover, while the effect of delay in a single reaction may be small, it is cumulative and linearly additive in directed lines.
The intrinsic stochasticity of the reactions that create mature protein make some variation in delay time inevitable. However, we do not yet know the exact nature of this variability or the functional form of the probability density function f T . To further complicate matters, there may exist a substantial amount of extrinsic variability in the delay time -the statistics of the PDF may vary from cell to cell.
We focused on the transient dynamics of M=G=? queues in order to demonstrate the effects of distributed delay in a tractable setting. However, as mentioned earlier, M=G=? queues may not always be a good model for protein production. For genes with low copy number or few available transcripts queues with Lv? service channels (M=G=L queues) may provide a better description. For eukaryotic systems models in which transcription and translation are decoupled into separate queues may also be relevant. In addition, as protein production rates are often coupled with extrinsic factors such as growth rate and cell cycle phase, L may depend on time and on the state of the system.
The complexity of biochemical reaction networks suggests the use of networks of queues [46], and sources could be toggled on and off by other components of a reaction network. Even protein production from a single transcript may be more accurately described by a sequence of M=G=1 queues with each codon as one in a chain of service stations. In such a model ribosomes move from one codon station to the next, and are not able to skip ahead. Such models will be considered in future studies.
One further complication occurs if the burstiness of the promoter is large [47]. In the above analysis, we assumed that the initiation events of proteins were exponentially distributed in time. Since this is not necessarily the case due to the burstiness of promoters, some limits need to be put on the usefulness of the above results. Equations (9) and (10) suggest that the transition to accelerated behavior occurs when swR=l: One can think of R=l as the average time, T R , it takes to initiate R proteins, and rewrite the boundary as swT R . One can then assume that if the burstiness of the initiation events is not large, i.e. that the mean burst size is less than the signal threshold, then it does not matter what the distribution of initiation events is. In other words, as long as approximately R proteins are initiated in the time T R , and the variance of that number is not large, then Eq. (16) still holds.

Methods
Signaling time distributions associated with a single gene via queueing theory Preliminary information. We first derive the signaling time distributions for a single gene that is modeled by an M=G=?
queue. An M=G=? queue is a queueing system consisting of a memoryless arrival process (M) and infinitely many service channels (?). The service time distribution is general (G) and there exists no maximal system size. Let 1. (I(t)) t §0 denote the input (arrival) process, 2. (Q(t)) t §0 denote the queue size process, and 3. (P(t)) t §0 denote the departure (completion) process.
Thus I(t), Q(t), and P(t) are the numbers of proteins that have entered the queue, are in the queue, and have departed the queue, respectively, at time t. Note that I(t)~Q(t)zP(t) for all t §0. Suppose that (I(t)) t §0 is a nonhomogeneous Poisson process with rate function l(t). Let T denote the (random) service time and let F T denote the cumulative distribution function (CDF) of T. This is the amount of time that a protein spends in the queue after entering. If the distribution of T is absolutely continuous, let f T denote the probability density function (PDF) of T. For t §0, define Notice that E½I(t)~L(t) for all t §0.
Proposition. (transient distributions; see e.g. [40]) Let t §0. The random variables Q(t) and P(t) are Poisson with means Signaling time distributions. Let S R denote the (random) first time at which P(S R )~R, i.e. S R~m inft : P(t)~Rg. We rescale time so that the rescaled completion process is a homogeneous Poisson process with rate 1. For t §0 define P P(t) :~E½P(t)~ð t 0 F T (s)l(s)ds: Define the rescaled departure process (P P(f)) byP P(f)P ( P P {1 (f)). Let j R denote the (random) time at whichP P(j R )~R. The random time j R has a gamma distribution with PDF so S R has PDF Computing the expectation of S R , we have ~ð P P(?) We now show that if T is symmetrically distributed about its mean and l is a constant function, then for every fixed value of R, increasing the standard deviation s T of T decreases the expected signaling time.
Proposition. Assume that T is symmetrically distributed about its mean and that l is a constant function. Let R[N. The function E½S R is a decreasing function of s T .
Proof. Suppose that l:l 0 . In light of (19b), it suffices to show that for every fixed t §0, P P(t) is an increasing function of s T . We write P P(t,s T ) and F T (t,s T ) to explicitly indicate the dependence of P P and F T on s T as well as t. Fix t §0 and let c 2 wc 1 . Define t~E½T. For every 0vzƒt, we have Finally, if tw2t, then the inequality P P(t,c 2 ) § P P(t,c 1 ) follows from computation (20) and the fact that for sw2t, F T (s,s T )~1 for all relevant values of s T .
Expected value of Q(S R ). Computing the expectation of Q(S R ), we have Example -Bernoulli delay distributions. Suppose that the rate function of the input process is constant and equal to l 0 , and T is a Bernoulli random variable described by the probability measure where 0vbv1. We begin by computing E½S R . Let x 1~t {2a(1{b) and x 2~t z2ab. The CDF of T is given by l 0 (t{t), t §x 2 : The signaling time S R has PDF We compute E½S R using (19b). The inverse of P P is defined only for t §x 1 , thus Finally, we express (21) using gamma functions: where C(z)~Ð ? 0 t z{1 e {t dt and C(z,a)~Ð ? a t z{1 e {t dt. We now examine the asymptotics of E½S R in the a?0 limit. The first and second partial derivatives of E½S R with respect to a are given by Expanding E½S R for small values of a gives Using the Stirling approximation n!*n n e {n ffiffiffiffiffiffiffi ffi 2pn p we therefore obtain and therefore In particular, for b~1=2 we have In this case the correction to the deterministic limit is of order a Rz1 =R Rz3=2 .
We obtain linear large a asymptotics by noting that the first term on the right side of (22) vanishes in the a?? limit: R bl 0 zt: Figure 6 shows a comparison between these analytical results and stochastic simulations.
Example 2-Normal delay distributions. Suppose that the rate function of the input process is constant and equal to l 0 . Suppose that T is a normal random variable with mean t and standard deviation s.
The CDF of T is given by where erf is the error function. For t §0 we have Expanding P P(t) about s~0 we obtain P P(t)~0 Note that the corrections to the first terms in the expansions are exponentially small in s 2 in both regimes. We denote by P Ã the approximation for P P which omits terms exponentially small in s 2 . The signaling time PDF of S R can then be approximated by Using (19a), we have Higher moments may be obtained in a similar manner. For a repressor switch, the P 2 process and therefore the ability of gene 2 to signal downstream components depend in complex ways on the statistical properties of T 2 . We examine these complex relationships by conditioning first on S 12 and then on I 2 . Suppose that S 12~t Ã . The key observation is this: for fixed twt Ã , E½P 2 (t)jS 12~t Ã can increase or decrease with the standard deviation s T2 of T 2 . We verify this assuming T 2 is symmetrically distributed about its mean and assuming l 2 is a constant function. If the midpoint t{t Ã =2 of the time interval ½t{t Ã ,t satisfies t{t Ã =2vE½T 2 , then is an increasing function of s T 2 and therefore E½P 2 (t)jS 12~t Ã increases as s T 2 increases. By contrast, if t{t Ã =2wE½T 2 , then the integral in (24) is a decreasing function of s T 2 and therefore E½P 2 (t)jS 12~t Ã decreases as s T 2 increases. Repressive signaling can therefore qualitatively affect the response of P 2 production to changes in the variability of T 2 . We now examine the ability of gene 2 to signal downstream components by conditioning on I 2 . Let R 2Ã [N. Let S 2Ã denote the time at which P 2 first reaches level R 2Ã . The key observation is this: If we assume that gene 2 shuts off after exactly f transcription initiation events, then E½S 2Ã jf can increase or decrease as a function of s T2 . Figure 7 demonstrates this numerically for a case in which l 2 is a constant function and T 2 is symmetrically distributed. In this case, we find that E½S 2Ã jf Intuitively, this is due to the fact that the order in which proteins enter the queue is not necessarily the same as the order in which they exit. Consider again Figure 7. When the total number of transcription initiation events, f, is smaller than 19, then more proteins enter the queue before the 10th protein than after it. It is therefore more likely that a protein entering before protein 10 will exit ahead of it than that a protein entering after protein 10 will exit before it. As a result, the expected time E½S 2Ã jf increases with s T2 . When the balance favors proteins that enter the queue after protein 10, the opposite is true, and E½S 2Ã jf decreases with s T2 .
We conjecture that this trichotomy holds in general if l 2 is a constant function and T 2 is symmetrically distributed about its mean.

Stochastic simulations
Gillespie's stochastic simulation algorithm generates an exact stochastic realization for a system of N species interacting through M reactions. The state of the system is stored in the vector X , and each reaction j is characterized by a state change vector Z j and its propensity function a j : R N ?R. If the system is in state X and reaction j occurs then the system state changes to X zZ j [5].
The idea behind extending Gillespie's SSA to model distributed delay is that if a reaction is to be delayed by some amount of time then we temporarily store this reaction along with the time at which the event will occur and we only apply this reaction at the given time. We used a version of the algorithm equivalent to those described in [32,48]. Note that [48] also describes a more efficient version of the algorithm.