Stochastic Delay Accelerates Signaling in Gene Networks
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 Gillespies 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–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–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–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 network 15 ., Delay can also affect the stochastic nature of gene expression , and the relation between the two can be subtle and complex 28–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 , , is not constant but instead a random variable ., If denotes the probability density function ( PDF ) of , this situation can be described deterministically by an integro-delay differential equation 32 of the form ( 1 ) where is a positive definite state vector of protein concentrations , and 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 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–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 ., 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 ( indicating a nearly deterministic completion time ) to near ( 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 reactions may impact the delay time distribution in complicated ways ., Therefore , we do not try to derive the actual shape of , 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 ( 2 ) where is the gene , and transcription is initiated at rate , which can depend explicitly on both time and protein number , ., After initiation , it takes a random time , , 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 ., In the setting of scheme ( 2 ) , first assume that does not depend on , and protein formation is initiated according to a memoryless process with rate ., A fully formed protein enters the population a random time 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 queue 40 , where indicates a memoryless source ( transcription initiation ) , 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 ., 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 ( 3a ) ( 3b ) where and are the copy numbers of the upstream and downstream genes , and are the number of functional proteins of each type , and is the random delay time of gene ., The transcription rate of gene 2 depends on and is given by a Hill function ., We consider the case in which activates ( depicted in Figure 1 ) and the case in which represses ., We now ask: If 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 , where is the Heaviside step function ., At what time does 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 ( 4a ) ( 4b ) Here is the maximum transcriptional initiation rate of and is the threshold value of the Hill function , i . e . the number of molecules of needed for half repression ( or half activation ) of gene, 2 . The second gene therefore becomes repressed ( or activated ) at the time at which copies of protein have been fully formed ., We first examine reaction ( 3a ) ., Assume that at time there are no proteins in the system ., Let denote the number of transcription initiation events that have occurred by time ( the arrival process of the queueing system ) , the number of proteins being formed at time ( the size of the queue at time ) , and the number of functional proteins that have been completed by time ( the exit process of the queueing system ) ., Since the arrival process is memoryless , is a Poisson process with constant rate for ., Hence , the expected value of is ., The exit process , i . e . the number of fully functional proteins that have emerged from the queue , , is a nonhomogenous Poisson process with time-dependent rate , where is the cumulative distribution function ( CDF ) of the delay time ., It then follows that Inactivation ( or activation ) of gene 2 occurs when enough protein 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 , , is given by ., Trivially , 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 , , and vary ., The probability density function of is given by ( See Methods ) ( 5 ) Consequently , the mean and variance of the time it takes for the original signal to propagate to the downstream gene can be written as: ( 6 ) ( 7 ) 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 discrete values , and with equal probability ., In this case , ( 8 ) where is the upper incomplete gamma function ., Expanding for small , we obtain ( See Methods ) ( 9 ) which is the deterministic limit ., The first term is the mean delay time and the second is the average time to initiate proteins at rate ., A similar expansion for fixed and large gives ( see panel ( c ) in Figure 2 ) ( 10 ) 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 , gives an estimate of when the behavior of the system changes from the deterministic limit ( for ) to a regime in which increasing the variability decreases the mean signaling time ( for ) ., It follows that the deterministic approximation given by Eq ., ( 9 ) is valid in an increasing range , as grows ( See Figure 2 ( C ) ) ., Indeed , an asymptotic analysis of Eq ., ( 8 ) shows that the corrections to Eq ., ( 9 ) are approximately of size , and therefore rapidly decrease with ( See Methods ) ., The bottom row of Figure 2 shows that these observations hold more generally: When is gamma distributed the mean time to produce proteins , , is very sensitive to randomness in delay time , but only when is small to intermediate ., As expected , the densities of the times to produce proteins , , are approximately normal and independent of the delay distribution when is large ( Middle panels of Figure 2 ) ., We therefore expect that for each fixed threshold , is a decreasing function of the standard deviation 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 protein , but exit after the protein increase , while the opposite is true for proteins that enter the queue after the protein , and exit before it ., Since only finitely many proteins enter the queue before the 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 , is an increasing function of , hence is decreasing function of for all ., Referring to Eq ., ( 6 ) , this implies that is a decreasing function of 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 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 ., 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 gene starts transcription of protein at rate , i . e ., The second gene , , starts transcription after reaches the threshold , so that ., For the coherent FFL , we assume that the promoter of gene acts as an AND gate so that ., We further assume that the promoter of in the incoherent FFL is active only in the presence of and absence of , so that we may write ., The signaling time between any two nodes and within the network , i . e . the random time between the initiation of transcription of gene to the formation of a total of proteins is denoted ., For each of the three pathways , the PDF of the signaling time is given by Eq ., ( 5 ) ., In addition , because the random times and are additive ( as are their variances ) , we can directly calculate the time at which reaches the threshold of gene as ., Therefore , the random time at which the coherent FFL turns on is simply given by ., Because and are decreasing functions of the delay variability , it can be expected that so is ., In contrast to the coherent FFL , the dynamics of the pulse generating incoherent FFL are less trivial ., Since the repressor ( ) overrides the activator ( ) , assuming transcription of turns on at time and turns off at time , generating a pulse of duration ., Note that can increase or decrease as a function of the standard deviation of the delay ( see Figure 4 where was equal for all pathways ) ., To see this , write as follows: ( 11 ) Each of the terms on the right side of Eq ., ( 11 ) is the expected signaling time of a single gene ( , , and , respectively ) ., Consequently , depends on as a linear combination of expected signaling time curves of the type pictured in Figure, 2 . The shapes of these signaling time curves determine the behavior of as a function of ., 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 ., 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 ) : ( 12a ) ( 12b ) where is a decreasing Hill function ( i . e . represses its own production ) and 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 is nonrandom ( degrade and fire oscillators ) 17 ., Starting with no proteins , is produced at a rate governed by the Hill function ., When the level of exceeds the midpoint of the Hill function , gene effectively shuts down ., The proteins remaining in the queue exit , producing a spike , after which degradation diminishes ., When the protein level drops sufficiently , reaction ( 12a ) reactivates and production of 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 genes 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 sub-circuit within the incoherent FFL ., Here we explicitly used a gamma-distributed delay time with mean , and ., We can use our theory to predict the change in the peak height of the oscillator as a function of ., For a delay that is gamma-distributed , the change in signaling time as a function of can be written as ( 13 ) where is given by Eq ., ( 6 ) and 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 less , then a total of less protein will be produced ., Therefore we can write the expected peak height of the oscillator as ( 14 ) However , due to degradation , Eq ., ( 14 ) overestimates the correction to the peak height ., Due to exponential degradation , only a fraction of the lost protein would have made it through to the peak ., Also , the duration of enzymatic decay is also reduced by a time ., Therefore , if we assume that the enzymatic decay reaction is saturated , we need to add an amount to Eq ., ( 14 ) ., This gives us a more accurate prediction of the mean peak height as ( 15 ) Figure 5 shows that this approximation works well , even for a Hill coefficient as low as ., 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 networks 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 ., 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 queues in order to demonstrate the effects of distributed delay in a tractable setting ., However , as mentioned earlier , queues may not always be a good model for protein production ., For genes with low copy number or few available transcripts queues with service channels ( 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 , 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 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 ( 16 ) One can think of as the average time , , it takes to initiate proteins , and rewrite the boundary as ., 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 proteins are initiated in the time , and the variance of that number is not large , then Eq ., ( 16 ) still holds ., Gillespies stochastic simulation algorithm generates an exact stochastic realization for a system of species interacting through reactions ., The state of the system is stored in the vector , and each reaction is characterized by a state change vector and its propensity function ., If the system is in state and reaction occurs then the system state changes to 5 ., The idea behind extending Gillespies 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 .
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 .
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 .
systems biology, stochastic processes, mathematics, theoretical biology, regulatory networks, synthetic biology, biology, computational biology, signaling networks, molecular biology, genetics and genomics, probability theory