Home
Search
Collections
Journals
About
Contact us
My IOPscience
Transcriptional dynamics with time-dependent reaction rates
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Biol. 12 016015 (http://iopscience.iop.org/1478-3975/12/1/016015) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 198.91.37.2 This content was downloaded on 21/02/2015 at 19:22
Please note that terms and conditions apply.
Phys. Biol. 12 (2015) 016015
doi:10.1088/1478-3975/12/1/016015
PAPER
RECEIVED
Transcriptional dynamics with time-dependent reaction rates
10 September 2014 REVISED
Shubhendu Nandi and Anandamohan Ghosh
16 January 2015
Indian Institute of Science Education and Research Kolkata Mohanpur 741246, India
ACCEPTED FOR PUBLICATION
E-mail: [email protected] and [email protected]
16 January 2015 PUBLISHED
Keywords: transcriptional regulation, non-Poissonian dynamics, information transfer, time-dependent model
11 February 2015
Abstract Transcription is the first step in the process of gene regulation that controls cell response to varying environmental conditions. Transcription is a stochastic process, involving synthesis and degradation of mRNAs, that can be modeled as a birth–death process. We consider a generic stochastic model, where the fluctuating environment is encoded in the time-dependent reaction rates. We obtain an exact analytical expression for the mRNA probability distribution and are able to analyze the response for arbitrary time-dependent protocols. Our analytical results and stochastic simulations confirm that the transcriptional machinery primarily act as a low-pass filter. We also show that depending on the system parameters, the mRNA levels in a cell population can show synchronous/asynchronous fluctuations and can deviate from Poisson statistics.
1. Introduction Transcription is the process of synthesis of RNAs from segments of DNA and is the first step in gene expression. Coarse-grained models of gene transcription often replace the complex biochemical reactions occurring during transcription with a single reaction step. In the simplest model of gene transcription, the process of mRNA synthesis and its subsequent degradation is modeled as a birth–death process [1] and mRNAs are seen to follow Poisson statistics [2]. This scenario is successful in constitutive gene expression seen mostly in housekeeping genes [3–5], where the gene is always in an active state and transcription occurs at random with an average rate equal to the ratio of the transcription and degradation rates. A refinement to this basic model is the ON–OFF model [6, 7] in which the gene toggles between an active (ON) and inactive (OFF) state with random transition probabilities. The distribution of mRNAs in this case tends to a negative binomial distribution in the limit when the transcription rate is large compared to the rate of degradation, and mRNA production is seen to be more burst-like. The ON–OFF model has been successful in explaining bursty transcriptional dynamics in facultative genes which turn on only when needed, but due to the intrinsic ON–OFF random switching in a population of cells, this © 2015 IOP Publishing Ltd
model is incapable of producing a synchronous response as is often observed in a multitude of cell activities such as the stress response of cells subject to a heat-shock [4]. The main focus of the present study is to consider a general model that incorporates the explicit time dependence of mRNA transcription and degradation rates. We consider a birth–death process with timedependent rates of transcription and degradation and obtain an exact analytical expression for P (m, t ), the probability of finding m mRNAs at time t. We find that P (m, t ) follows a non-homogeneous Poisson distribution (NHPD) which offers a considerably richer dynamics than its homogeneous counterpart. We can now consider arbitrary environmental fluctuations, captured in our model by the time-dependent variation of transcription rates. In particular we consider periodic fluctuations and demonstrate how the cell can perform simple regulatory tasks by acting as a lowpass filter. We validate our exact results with numerical simulations and also compute mutual information (MI) to illustrate the efficiency of transcriptional coding. We show that the synchronous behavior in a cell population is a natural outcome in our time-dependent model and the system being driven out of equilibrium can lead to interesting statistical inference such as differences in single cell dynamics and population averages. In section 2, we present the time-dependent model and in section 3, we discuss our results.
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
k on
2. Model
I⟶A k off
In a basic two-step model of gene transcription, a single gene which is always in an active state, is assumed to be transcribed at a constant rate kt, and the synthesized mRNAs are assumed to be degraded at a constant rate kd [2]. The ‘birth’ and ‘death’ of an mRNA is modeled by a pair of first order kinetic reactions: kt
A⟶I kt
A⟶A+M kd
M ⟶ ∅.
Here kon and koff are the rates of gene activation and inactivation respectively. In the steady state, one can derive exact analytic expressions for the variance and mean number of mRNAs and the Fano factor can be written as
D⟶D+M kd
M ⟶ ∅,
(1)
where mRNAs are denoted by M, an active gene by D and ∅ denotes the degradation products of mRNAs. Genes are often present in just one or two copies per cell and for our choice of reaction rates, mRNAs are typically on the order of tens of molecules. Consequently, the process of gene transcription is fundamentally stochastic and the rate constants kt and kd in equation (1) take on the meaning of probabilities per unit time. Since the mRNA expression levels are stochastic variables, a relevant quantity of interest is the probability mass function P (m, t ) which gives the probability of finding m mRNAs at a given time t. The birth–death process, in equation (1), can be cast as a master equation and can be solved for P (m, t ) which follows a Poisson distribution. The model, equation (1), henceforth referred to as the homogeneous Poisson (HP) model, predicts that mRNA production events occur at random with a constant average rate. In the steady state, the rate is equal to the average number of mRNAs, given by μ = kt k d and the process being Poissonian, the variance σ 2 = μ. A measure of the fluctuations is typically given by the Fano factor, F = σ 2 μ, and for the HP model, F = 1. However, recent experiments have shown that transcription events occur in irregular ‘bursts’ characterized by short and intense periods of mRNA production interspersed between longer periods of little or no activity [7–11]. The degree of burstiness can be quantified by the Fano factor F with referral to the Poisson process (F = 1). A Fano factor of F > 1 indicates super-Poissonian or bursty behavior whereas a Fano factor of F < 1 indicates sub-Poissonian or antibursty behavior. Since, the HP model cannot account for bursty behavior, random ON–OFF switching of the gene is introduced so that instead of being active at all times, the gene is assumed to switch stochastically between active and inactive states. In this model [6], henceforth referred to as the ON–OFF model, transcription is only possible when the gene is in the active state (A) and there is no transcription when the gene is in the inactive state (I). The following set of reactions describe the process: 2
(2)
F=1+
kt koff
( kon + koff )( kon + koff + kd )
.
(3)
It can be seen from equation (3) that F ⩾ 1. In the limit k off = 0 , the ON–OFF model reduces to the HP model. It must be pointed out however that the ON– OFF model does not show sub-Poissonian or antibursty behavior since F can never be less than 1. 2.1. Time-dependent model In both models discussed so far, the rates kt and kd are treated as constants. However, cells are constantly exposed to fluctuating environmental signals. Often, the signals activate complex biochemical signal transduction pathways that result in a chemical modification of specific transcription factors [12]. In other systems, the signaling molecules may simply enter the cell and directly bind the transcription factor. Usually the signal causes a physical change in the shape of the transcription factor protein causing it to assume an active molecular state. The transcription factor in its active state binds the promoter region of a gene and increases the rate of transcription [13]. Thus environmental signals that fluctuate in time, causes fluctuations in the rate of transcription, and consequently the transcription rate kt should in general be timedependent. The ON–OFF model indirectly addresses this issue by switching the gene on and off in concert with the fluctuating signals. But it does not allow transcription to proceed when the gene is in the inactive state although many genes are known to have non-zero basal expression levels and show bursty activity as well [3]. Also, the ON–OFF model only switches between two ON–OFF states, but recent experimental data on eve stripe 2 expression in living Drosophila embryos support the idea of a multi-state transcription model [14]. Furthermore, none of the existing models of gene transcription are able to account for arbitrary time variations in rate kinetics kt , k d . We model the synthesis and degradation of mRNAs by the pair of reactions k t (t )
D⟶D+M k d (t )
M ⟶ ∅,
(4)
where the transcription and degradation rates kt and kd are now both functions of time. The symbols D, M and
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 1. (a) A schematic representation of the process of transcription of an active gene with time dependent propensities. For m mRNAs at time t, the propensities for transcription and degradation of an mRNA are given by kt(t) and mk d (t ) respectively. (b) Results of the stochastic simulation algorithm (SSA) fitted with the probability distribution P (m, t ) obtained from our model. The histogram was plotted using binned data from 1000 realizations of the SSA over two complete periods of the input signal. Inset: functional form of the input signal kt(t) which is a square-wave oscillating between k1 and k2 with k d = 0.05 (arb. time)−1.
∅ have the same meaning as in the HP model. A schematic of the process is shown in figure 1(a). The probability P (m, t ) of finding m mRNAs at time t obeys the master equation: ∂P ( m , t ) = kt (t ) [P (m − 1, t ) − P (m , t ) ] ∂t + k d (t ) [(m + 1) P (m + 1, t ) − mP (m , t ) ] .
(6)
⎛ dμ ⎞ ∂ϕ ∂ϕ + uϕ ⎜ + k d μ − kt ⎟ + k d u = 0. (7) ⎝ dt ⎠ ∂t ∂u
If we choose the function μ such that it satisfies dμ + k d μ − kt = 0 dt
(8)
∂ϕ ∂ϕ + kd u = 0. ∂t ∂u
(9)
The function μ can be obtained by solving equation (8) giving t
t
k d dτ
∫
= (1 + u)m0 .
kt e ∫0
τ
k d ds
dτ .
(10)
0
We can solve equation (9) by employing the method of characteristics obtaining 3
(12)
Since F = ϕ (u, t )e μ(t ) u and μ (0) = 0, it follows from equation (12) that F (u , 0) = ϕ (u , 0) = ψ (u) = (1 + u)m0 .
(13)
We note that the choice of ϕ (u , t ) = ψ (U ) = (1 + U )m0
(14)
satisfies the boundary condition given by equation (12) and is therefore a solution to equation (7). Assuming that we start with zero initial number of mRNAs such that m0 = 0, the solution to equation (6) with timedependent kt and kd can be written as F (z , t ) = e μ(z − 1).
By definition, 1 ∂m P (m, t ) = m ! ∂z m F (z , t )
z=0
(15)
satisfies P (m, t ) which generates an
NHPD: P (m , t ) =
then equation (7) reduces to
∑ m
δm, m0 z m = z
(5)
A trial-solution of the form F = ϕ (u, t )e μ(t ) u upon substitution in equation (6) yields
μ (t ) = e ∫0
t
m m0
⎛ ∂F ⎞ ∂F = u ⎜ kd − kt F ⎟ . ⎝ ∂u ⎠ ∂t
(11)
where ψ is any arbitrary function of U = ue− ∫0 kd dτ . If there are initially m0 mRNAs then the generating function F satisfies F (z , 0) = ∑P (m , 0) z m =
2.2. Exact solution By defining the generating function F (z , t ) as ∑m z mP (m, t ) and u = z − 1, we can convert equation (5) into a first order partial differential equation:
−
ϕ (u , t ) = ψ (U ),
1 −μ m e μ m!
(16)
with a time-dependent function μ (t ) given by equation (10). In the following sections we consider variations in time-dependent rates and study an ensemble of genetically identical cells and assume that each cell has exactly one copy of the gene. We also compare our exact results with numerical simulations utilizing a stochastic simulation algorithm (SSA) (see appendix A) which is a modification of Gillespie’s original algorithm [15, 16]. We choose, for example, kt(t) to be a square-wave alternating between two
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 2. Transcriptional time-series generated for two different time-periods (a) T = 10k d−1 and (b) T = k d−1, corresponding to low and high frequency oscillations in kt respectively and k d = 0.05 (arb. time)−1. For each time-period, two separate realizations of the SSA are shown. The simulation results are fitted with the theoretical curves (blue) which correspond to the average number of mRNAs obtained by integrating equation (17). The corresponding kt are also shown on top. (c) Plot of the mutual information (MI) of two time-series data (the transcription rate and the corresponding mRNA levels) for different values of kd, as the time-period T of the input signal kt is varied. Each curve is averaged over 100 realizations.
constant transcription rates k1 and k2 and verify that equation (16) gives the correct distribution. The result is plotted in figure 1(b).
3. Results and discussion 3.1. Low pass filtering Low pass filtering is a common survival strategy employed by gene regulatory networks and has been experimentally observed in synthetic transcriptional cascades [17] and in the metabolic system of yeasts subject to a periodically varying external glucose source [18]. The ability of gene transcription networks to respond reliably to a slowly varying signal, while ignoring fast fluctuations, confers a distinct fitness advantage. Maintaining a regulatory network that can interpret an external signal and act accordingly has its price. Regulation is only feasible when the benefit exceeds this cost of having a regulatory system in place. It is well known that short pulses in the input signal, whose duration is less than a critical value, can have a deleterious effect on growth [13]. The reason is that ultrashort pulses lead to a reduction in fitness because the gene products do not have sufficient time to reach levels in which the accumulated benefit exceeds the cost of production. Therefore, it is beneficial for an organism to filter out rapid fluctuations in the input signal. We show using a minimal stochastic model of gene transcription how this may be achieved in a living cell. To probe the dynamic response of our model, we first choose kt(t) to be a square-wave alternating between two constant transcription rates k1 and k2 with a time-period T and assume the degradation probability to be a constant kd. Thus equation (10) reduces to t
μ (t ) = e
−k d t
∫
kt (τ)e k d τdτ ,
(17)
0
4
where μ (t ) is the average number of mRNAs at time t. To see the temporal evolution of μ (t ), we then integrate equation (17) for two different time periods corresponding to low and high frequency oscillations in kt(t). The theoretical curves for each time-period are compared to two separate realizations of the SSA and the results are plotted in figure 2. We find that for high frequency oscillations in kt, the transcription network is insensitive to the variations in the transcription rate and behaves as a low-pass filter. Noting the formal equivalence of equation (8) with equations that describe driven filter circuits, we understand this result by comparing the ‘response time’ of the transcription network given by k d−1 with the period T of the oscillatory transcription rate. Just as in electronic filter circuits, if we squeeze in too many periods so that T is small, the network filters out the high frequency fluctuations in kt(t). The number of transcripts (output) encode the temporal variation of the reaction rates (input). The efficiency of the filtering property can be quantified by the rate of information transfer between the input and the output signal. A quantitative measure of the information transfer can be obtained by estimating the MI (see appendix B) between the transcription rate kt and the corresponding mRNA levels obtained from the SSA. When the system filters out the high-frequency fluctuations in the signal, the corresponding MI is expected to decrease. As the frequency decreases and the system becomes sensitive to the oscillations in kt, the MI increases reflecting this change. As seen in figure 2(c), when the period T is small, the information transfer is close to zero signifying that the system filters out the fluctuations in the transcription rate, and as T becomes large, the MI gradually increases and becomes close to 1. We also study the effect of different degradation rates kd on the information transfer and in figure 2(c),
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
we show the corresponding MI curves for three different values of kd corresponding to kd = 0.025, kd = 0.05, and kd = 0.1 (measured in arbitrary units of time). We see that as kd increases, the MI saturates faster. According to our previous analysis, this is to be expected since k d−1 is the ‘time-constant’ in our transcription filter. As kd increases, the system becomes more responsive to high frequency fluctuations in kt and the MI curves get steeper. Thus with the increase in kd, the sensitivity of the transcriptional network to fluctuations in kt increases. The effect of time-dependent degradation rates also has a substantial bearing on transcriptional dynamics. Our preliminary results with kd chosen to be a square-wave oscillating between k d1 = 0.005 and k d2 = 0.05 in phase with kt (also measured in the same units) shows some interesting features when we vary the time-period Td of kd keeping the period Tt of kt fixed at Tt = 200. The results are shown in figure 3. We see that when the ratio of the time-periods ρ = Td Tt equals 1, the burstiness in transcription disappears and mRNA expression levels resemble constitutive transcription. This is because the oscillations in both rates increase and decrease in tandem which help maintain an approximately constant transcript level. However, when ρ deviates from 1, there is a reappearance of bursts and the burst-frequency is roughly proportional to ρ. The dynamics can be more intricate if kt and kd have a non-zero relative phase. Further detailed studies need to be undertaken to elucidate the complex interplay between kt and kd. 3.2. Synchronous/asynchronous response in a cell population In order to test the system response for noisy signals, we next generate kt(t) so that it consists of a sequence of pulses in which the pulse durations d are drawn from an exponential distribution P (d) = λ exp (−λd). We explore two different signal modalities: one in which every copy of the system has a different uncorrelated transcription rate but the pulse widths for each signal are drawn from the same exponential distribution, and the other in which all copies of the system transcribe with identical transcription rates. The results of the simulation for the two different modalities are plotted in figure 4. In scenario (a), where the transcription rates are specific to each copy of the system with the only commonality being the average pulse duration, we see that the individual time-series are bursty but there is considerable heterogeneity in the burst timings (figures 4 (a1)–(a3)). This is due to the lack of correlation between the individual transcription rates and both the time-averaged and ensemble-averaged values of the Fano factor deviate from unity. This is somewhat similar to the ON–OFF model where switching of the gene on and off is seen as a stochastic process and the switching events between different copies of the 5
system are uncorrelated and we expect the ensemble and time averaged Fano factors to be equal. We see that in scenario (b), where all copies have the same kt, individual realizations of the SSA show synchronized bursting. The time-averaged Fano factor for each realization is super-Poissonian showing that the time-series is bursty but the ensemble-averaged Fano factor at each time instant is close to 1 (figures 4 (b1)–(b2)). This shows that although each copy of the system shows burstiness in time, at any given instant, the number of mRNAs produced is fairly homogeneous across the population. This is expected as the underlying distribution is an NHPD, as predicted by our model, and at each instant the average number of mRNAs has a simple Poisson distribution. The timeaveraged and ensemble-averaged values of the Fano factor do not match because the system is repeatedly perturbed and does not have time to relax.
4. Conclusions Sensory transcription networks are designed to respond to dynamically changing environments. We offer a new conceptual framework in which the environmental fluctuations are encoded in the transcription rate which varies in time. It is well known that the fluctuations in gene expression networks can have both extrinsic and intrinsic components [19, 20]. The time-dependent modulation of the transcription rate may arise from propagation of changes in upstream signaling as in fluctuations in activator concentrations in regulatory networks [21, 22] or it may be a result of intrinsic switching of the gene between ON–OFF states in the absence of any genetic regulation or external signaling [23]. It becomes difficult to deconvolve experimentally the contribution of each component [24] and in our model any change in the environment, extrinsic and/or intrinsic, is encoded in the transcription rate kt. We obtain an exact analytic solution of the chemical master equation resulting from a birth–death process of mRNA production and degradation with time-dependent kt , k d . Since the system is inherently stochastic, the reaction rates are taken as probabilities per unit time. We find that the probability mass density P (m, t ) follows an NHPD and we validate our result by stochastic simulations using a modified Gillespie algorithm [25]. For periodic variations in kt, we find that the system behaves as a low pass filter and the information transfer depends on the period of the signal kt and the response time k d−1 such that when T is small relative to k d−1, fluctuations in kt are filtered out and the MI becomes close to zero. We also test the system response for noisy signals consisting of a train of pulses in which the pulse widths are sampled from an exponential distribution. The two different signal modalities that we explore indicate that the system can have quite different responses
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 3. (a) Temporal variation in kt with period Tt = 200 . (b1), (c1) and (d1) are the temporal variations in kd with periods Td = 200 , Td = 400 and Td = 100 , respectively, and (b2), (c2) and (d2) are the corresponding transcript levels showing two separate realizations of the SSA for each kd (orange and green). The theoretical curves (shown in blue) are obtained by integrating equation (10) and time is in arbitrary units.
based on the temporal correlations between the transcription rates of different copies of the system. Our simulations show that a population of identical cells with the same transcription rate shows synchronized bursting and at any instant, the average number of mRNAs produced is Poisson distributed. This corresponds to genes that faithfully transmit fluctuations in the environment, such as the concentration of a given hormone [5]. Such synchrony in mRNA levels in cell populations are often observed during cellular stress response [4, 19]. However, if the transcription rates across the population lack temporal correlation, the synchrony in the burst timings is lost and the ensemble-averaged Fano factor deviates from the Poisson 6
case. This can be viewed as an extension of the ON– OFF model for finite times and with non-zero basal expression levels which result in phenotypic heterogeneity. With the advent of new experimental techniques, it is now possible to study transcription dynamics at the population level as well as in single cells [3]. Our work suggests that transcription is a nonequilibrium process and we need to exercise caution in quantifying statistical fluctuations, since depending on the burst synchrony/asynchrony, the time-averaged and the ensemble-averaged Fano factor may or may not be equal. Recent experimental studies have provided insight into the origins of burstiness where chromatin
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 4. (a1) Transcriptional time-series generated by the SSA where transcription rate kt has exponentially distributed pulse-widths with parameter λ = 0.5k d and kt stochastically fluctuates between k1 k d = 2 and k 2 k d = 16 . (a2) Another realization of the SSA with a kt which is temporally uncorrelated to the kt used in (a1), but has the same parameters λ, k1, k2 and kd. (a3) Ensemble-averaged Fano factor computed at different time-instants. (b1) Transcriptional time-series generated for an ensemble-wide transcription rate kt with a parameter λ = 0.5k d stochastically fluctuating between k1 k d = 2 and k 2 k d = 16 . Two different realizations of the SSA are shown. (b2) Ensemble-averaged Fano factor computed at different time-instants. In all of (a1–a3) and (b1–b2), mRNA time-series generated by the SSA are fitted with the theoretical value of the average mRNA number (blue) obtained by integrating equation (17). The ensemble Fano factor fluctuates in time around a mean value indicated by the red dotted lines. (c) Table summarizing the values of the time-averaged and ensemble-averaged Fano factors for the two different scenarios a and b.
remodeling [26] and DNA supercoiling [27] have been implicated as possible molecular mechanisms underlying transcriptional bursts. Here we consider a coarse-grained stochastic model and show how the transcriptional machinery may display bursty dynamics depending on the time-scale of fluctuation of the time-dependent transcription rate. Our study indicates that there exist a possible connection between bursty mRNA dynamics and information transfer due to transcription. Simple population dynamics models which treat transcription as a birth–death process can account for burstiness by introducing stochastic switching of the gene [6], but these models cannot account for sub-Poissonian or anti-bursty behavior which is also a recurring theme in diverse areas of biology [28–32]. The phenomenon of transcriptional antibursts has recently been observed experimentally in individual cells of E.coli., driven by the lar promoter, under weak and medium induction [32]. Simulations of our time-dependent model also indicate low Fano factors in certain time-scales which may be relevant in experimental situations but it remains to be ascertained if it is a true statistical phenomenon. 7
Further detailed studies need to be undertaken to investigate these issues in realistic gene regulatory networks.
Appendix A. Stochastic simulation algorithm The standard Gillespie algorithm [15, 16] is used to simulate chemical systems where the reaction propensities are time-independent. Since kt(t) in our model varies with time, we use a modified version of Gillespie’s first-reaction method to account for the time dependent transcription rates. If a time-dependent reaction rate k(t) is a piece-wise linear function with a finite number of discontinuities, following [25], we treat the discontinuities as detailed below. Let k(t) be discontinuous at a point tn so that k(t) satisfies ⎧ k< (t ) when t < tn k (t ) = ⎨ ⎩ k> (t ) when t > tn.
(A.1)
To simulate a discontinuous change in the rate, we change k(t) according to equation (A.1) when the sum
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
of the current simulation time and the next putative reaction time pass tn. k(t) is changed from k < (t ) to k > (t ) and the time is reset to t = tn . Consequently, our algorithm proceeds as follows. (1)Initialize the number of mRNAs. Set time t = 0. (2)Calculate the propensity for each chemical reaction. (3)For each reaction i generate a putative next reaction time τi . (4)Let τμ , corresponding to reaction μ, be the minimum τi . (5)Let tn be the time for the next discontinuous change in kt(t). If t + τμ < tn, then change the number of mRNAs appropriately for the occurrence of reaction μ. Change t to t + τμ. If t + τμ > tn, then change kt(t) accordingly. Set t = tn . (6)Go to step 2. To simulate arbitrary functions, we need to approximate the given function by a series of step functions or a piece-wise linear function and proceed as outlined above.
Appendix B. Mutual information The MI of two discrete random variables X and Y is a measure of the reduction in uncertainty about one random variable given knowledge of the other and can be defined as [33]: MI (X ; Y ) =
⎛ p (x , y ) ⎞ ⎟ , (B.1) ⎝ p (x ) p (y ) ⎠
∑ ∑ p (x, y) log ⎜ y∈Y x∈X
where p (x , y) is the joint probability distribution function of X and Y and p(x) and p(y) are the marginal probability distribution functions of X and Y respectively. If the base of the logarithm is 2, (the usual choice and the one we have used here), MI is measured in bits. Alternatively, MI of X and Y can be expressed as MI (X ; Y ) = H (X ) − H (X Y ),
(B.2)
where H(X) is the marginal entropy of X and H ((X ∣ Y )) is the conditional entropy of X conditioned on Y. Using the above definition, we generate the MI of two time-series, the signal kt(t) and the response m(t), where kt(t) and m(t) are the transcription rate and mRNA level respectively at time t.
References [1] van Kampen N G 2007 Stochastic Processes in Physics and Chemistry (Amsterdam: North-Holland, Elsevier) [2] Thattai M and van Oudenaarden A 2001 Intrinsic noise in gene regulatory networks Proc. Natl Acad. Sci. USA 98 8614–9
8
[3] Sanchez A and Golding I 2013 Genetic determinants and cellular constraints in noisy gene expression Science 342 1188–93 [4] Weake V M and Workman J L 2010 Inducible gene expression: diverse regulatory mechanisms Nat. Rev. Genetics 11 426–37 [5] Lionnet T and Singer R H 2012 Transcription goes digital EMBO Rep. 13 313–21 [6] Peccoud J and Ycart B 1995 Markovian modeling of geneproduct synthesis Theor. Population Biol. 48 222–34 [7] Suter D M, Molina N, Gatfield D, Schneider K, Schibler U and Naef F 2011 Mammalian genes are transcribed with widely different bursting kinetics Science 332 472–4 [8] Golding I, Paulsson J, Zawilski S M and Cox E C 2005 Realtime kinetics of gene activity in individual bacteria Cell 123 1025–36 [9] Chubb J R, Trcek T, Shenoy S M and Singer R H 2006 Transcriptional pulsing of a developmental gene Curr. Biol. 16 1018–25 [10] So L-H, Ghosh A, Zong C, Sepúlveda L A, Segev R and Golding I 2011 General properties of transcriptional time series in Escherichia coli Nat. Genetics 43 554–60 [11] Raj A, Peskin C S, Tranchina D, Vargas D Y and Tyagi S 2006 Stochastic mRNA synthesis in mammalian cells PLoS Biol. 4 1707–19 [12] Gomperts B D, Kramer I M and Tatham P E R 2002 Signal Transduction vol 18 (New York: Academic) [13] Alon U 2007 An Introduction to Systems Biology: Design Principles of Biological Circuits vol 10 (Boca Raton, FL: CRC) [14] Bothma J P, Garcia H G, Esposito E, Schlissel G, Gregor T and Levine M 2014 Dynamic regulation of eve stripe 2 expression reveals transcriptional bursts in living Drosophila embryos Proc. Natl Acad. Sci. USA 7 1–6 [15] Gillespie D T 1976 A general method for numerically simulating the stochastic time evolution of coupled chemical reactions J. Comput. Phys. 22 403–34 [16] Gillespie D T 1977 Exact stochastic simulation of coupled chemical reactions J. Phys. Chem. 81 2340–61 [17] Hooshangi S, Thiberge S and Weiss R 2005 Ultrasensitivity and noise propagation in a synthetic transcriptional cascade Proc. Natl Acad. Sci. USA 102 3581–6 [18] Bennett M R, Pang W L, Ostroff N A, Baumgartner B L, Nayak S, Tsimring L S and Hasty J 2008 Metabolic gene regulation in a dynamically changing environment Nature 454 1119–22 [19] Raj A and van Oudenaarden A 2008 Nature, Nurture, or chance: stochastic gene expression and its consequences Cell 135 216–26 [20] Tkačik G, Gregor T and Bialek W 2008 The role of input noise in transcriptional regulation PLoS One 3 e2774 [21] Elowitz M B and Leibler S 2000 A synthetic oscillatory network of transcriptional regulators Nature 403 335–8 [22] Locke J C W, Young J W, Fontes M, Jimenez M J H and Elowitz M B 2011 Stochastic pulse regulation in bacterial stress response Science 334 366–9 [23] Harper C V et al 2011 Dynamic analysis of stochastic transcription cycles PLoS Biol. 9 e1000607 [24] Hilfinger A and Paulsson J 2011 Separating intrinsic from extrinsic fluctuations in dynamic biological systems Proc. Natl Acad. Sci. USA 108 12167–72 [25] Shahrezaei V, Ollivier J F and Swain P S 2008 Colored extrinsic fluctuations and stochastic gene expression Mol. Syst. Biol. 4 196 [26] Voss T C and Hager G L 2014 Dynamic regulation of transcriptional states by chromatin and transcription factors Nat. Rev. Genetics 15 69–81 [27] Chong S, Chen C, Ge H and Xie X S 2014 Mechanism of transcriptional bursting in bacteria Cell 158 314–26 [28] Gur M, Beylin A and Snodderly D M 1997 Response variability of neurons in primary visual cortex (V1) of alert monkeys J. Neurosci. 17 2914–20 [29] Kara P, Reinagel P and Reid R C 2000 Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons Neuron 27 635–46
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
[30] Maimon G and Assad J A 2009 Beyond poisson: increased spiketime regularity across primate parietal cortex Neuron 62 426–40 [31] Scaglione A, Moxon K A, Aguilar J and Foffani G 2011 Trialto-trial variability in the responses of neurons carries information about stimulus location in the rat whisker thalamus Proc. Natl Acad. Sci. USA 108 14956–61
9
[32] Kandhavelu M, Häkkinen A, Yli-Harja O and Ribeiro A S 2012 Single-molecule dynamics of transcription of the lar promoter Phys. Biol. 9 026004 [33] Shannon C E 1948 A mathematical theory of communication Bell Syst. Tech. J. 27 379–423
Search
Collections
Journals
About
Contact us
My IOPscience
Transcriptional dynamics with time-dependent reaction rates
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Biol. 12 016015 (http://iopscience.iop.org/1478-3975/12/1/016015) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 198.91.37.2 This content was downloaded on 21/02/2015 at 19:22
Please note that terms and conditions apply.
Phys. Biol. 12 (2015) 016015
doi:10.1088/1478-3975/12/1/016015
PAPER
RECEIVED
Transcriptional dynamics with time-dependent reaction rates
10 September 2014 REVISED
Shubhendu Nandi and Anandamohan Ghosh
16 January 2015
Indian Institute of Science Education and Research Kolkata Mohanpur 741246, India
ACCEPTED FOR PUBLICATION
E-mail: [email protected] and [email protected]
16 January 2015 PUBLISHED
Keywords: transcriptional regulation, non-Poissonian dynamics, information transfer, time-dependent model
11 February 2015
Abstract Transcription is the first step in the process of gene regulation that controls cell response to varying environmental conditions. Transcription is a stochastic process, involving synthesis and degradation of mRNAs, that can be modeled as a birth–death process. We consider a generic stochastic model, where the fluctuating environment is encoded in the time-dependent reaction rates. We obtain an exact analytical expression for the mRNA probability distribution and are able to analyze the response for arbitrary time-dependent protocols. Our analytical results and stochastic simulations confirm that the transcriptional machinery primarily act as a low-pass filter. We also show that depending on the system parameters, the mRNA levels in a cell population can show synchronous/asynchronous fluctuations and can deviate from Poisson statistics.
1. Introduction Transcription is the process of synthesis of RNAs from segments of DNA and is the first step in gene expression. Coarse-grained models of gene transcription often replace the complex biochemical reactions occurring during transcription with a single reaction step. In the simplest model of gene transcription, the process of mRNA synthesis and its subsequent degradation is modeled as a birth–death process [1] and mRNAs are seen to follow Poisson statistics [2]. This scenario is successful in constitutive gene expression seen mostly in housekeeping genes [3–5], where the gene is always in an active state and transcription occurs at random with an average rate equal to the ratio of the transcription and degradation rates. A refinement to this basic model is the ON–OFF model [6, 7] in which the gene toggles between an active (ON) and inactive (OFF) state with random transition probabilities. The distribution of mRNAs in this case tends to a negative binomial distribution in the limit when the transcription rate is large compared to the rate of degradation, and mRNA production is seen to be more burst-like. The ON–OFF model has been successful in explaining bursty transcriptional dynamics in facultative genes which turn on only when needed, but due to the intrinsic ON–OFF random switching in a population of cells, this © 2015 IOP Publishing Ltd
model is incapable of producing a synchronous response as is often observed in a multitude of cell activities such as the stress response of cells subject to a heat-shock [4]. The main focus of the present study is to consider a general model that incorporates the explicit time dependence of mRNA transcription and degradation rates. We consider a birth–death process with timedependent rates of transcription and degradation and obtain an exact analytical expression for P (m, t ), the probability of finding m mRNAs at time t. We find that P (m, t ) follows a non-homogeneous Poisson distribution (NHPD) which offers a considerably richer dynamics than its homogeneous counterpart. We can now consider arbitrary environmental fluctuations, captured in our model by the time-dependent variation of transcription rates. In particular we consider periodic fluctuations and demonstrate how the cell can perform simple regulatory tasks by acting as a lowpass filter. We validate our exact results with numerical simulations and also compute mutual information (MI) to illustrate the efficiency of transcriptional coding. We show that the synchronous behavior in a cell population is a natural outcome in our time-dependent model and the system being driven out of equilibrium can lead to interesting statistical inference such as differences in single cell dynamics and population averages. In section 2, we present the time-dependent model and in section 3, we discuss our results.
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
k on
2. Model
I⟶A k off
In a basic two-step model of gene transcription, a single gene which is always in an active state, is assumed to be transcribed at a constant rate kt, and the synthesized mRNAs are assumed to be degraded at a constant rate kd [2]. The ‘birth’ and ‘death’ of an mRNA is modeled by a pair of first order kinetic reactions: kt
A⟶I kt
A⟶A+M kd
M ⟶ ∅.
Here kon and koff are the rates of gene activation and inactivation respectively. In the steady state, one can derive exact analytic expressions for the variance and mean number of mRNAs and the Fano factor can be written as
D⟶D+M kd
M ⟶ ∅,
(1)
where mRNAs are denoted by M, an active gene by D and ∅ denotes the degradation products of mRNAs. Genes are often present in just one or two copies per cell and for our choice of reaction rates, mRNAs are typically on the order of tens of molecules. Consequently, the process of gene transcription is fundamentally stochastic and the rate constants kt and kd in equation (1) take on the meaning of probabilities per unit time. Since the mRNA expression levels are stochastic variables, a relevant quantity of interest is the probability mass function P (m, t ) which gives the probability of finding m mRNAs at a given time t. The birth–death process, in equation (1), can be cast as a master equation and can be solved for P (m, t ) which follows a Poisson distribution. The model, equation (1), henceforth referred to as the homogeneous Poisson (HP) model, predicts that mRNA production events occur at random with a constant average rate. In the steady state, the rate is equal to the average number of mRNAs, given by μ = kt k d and the process being Poissonian, the variance σ 2 = μ. A measure of the fluctuations is typically given by the Fano factor, F = σ 2 μ, and for the HP model, F = 1. However, recent experiments have shown that transcription events occur in irregular ‘bursts’ characterized by short and intense periods of mRNA production interspersed between longer periods of little or no activity [7–11]. The degree of burstiness can be quantified by the Fano factor F with referral to the Poisson process (F = 1). A Fano factor of F > 1 indicates super-Poissonian or bursty behavior whereas a Fano factor of F < 1 indicates sub-Poissonian or antibursty behavior. Since, the HP model cannot account for bursty behavior, random ON–OFF switching of the gene is introduced so that instead of being active at all times, the gene is assumed to switch stochastically between active and inactive states. In this model [6], henceforth referred to as the ON–OFF model, transcription is only possible when the gene is in the active state (A) and there is no transcription when the gene is in the inactive state (I). The following set of reactions describe the process: 2
(2)
F=1+
kt koff
( kon + koff )( kon + koff + kd )
.
(3)
It can be seen from equation (3) that F ⩾ 1. In the limit k off = 0 , the ON–OFF model reduces to the HP model. It must be pointed out however that the ON– OFF model does not show sub-Poissonian or antibursty behavior since F can never be less than 1. 2.1. Time-dependent model In both models discussed so far, the rates kt and kd are treated as constants. However, cells are constantly exposed to fluctuating environmental signals. Often, the signals activate complex biochemical signal transduction pathways that result in a chemical modification of specific transcription factors [12]. In other systems, the signaling molecules may simply enter the cell and directly bind the transcription factor. Usually the signal causes a physical change in the shape of the transcription factor protein causing it to assume an active molecular state. The transcription factor in its active state binds the promoter region of a gene and increases the rate of transcription [13]. Thus environmental signals that fluctuate in time, causes fluctuations in the rate of transcription, and consequently the transcription rate kt should in general be timedependent. The ON–OFF model indirectly addresses this issue by switching the gene on and off in concert with the fluctuating signals. But it does not allow transcription to proceed when the gene is in the inactive state although many genes are known to have non-zero basal expression levels and show bursty activity as well [3]. Also, the ON–OFF model only switches between two ON–OFF states, but recent experimental data on eve stripe 2 expression in living Drosophila embryos support the idea of a multi-state transcription model [14]. Furthermore, none of the existing models of gene transcription are able to account for arbitrary time variations in rate kinetics kt , k d . We model the synthesis and degradation of mRNAs by the pair of reactions k t (t )
D⟶D+M k d (t )
M ⟶ ∅,
(4)
where the transcription and degradation rates kt and kd are now both functions of time. The symbols D, M and
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 1. (a) A schematic representation of the process of transcription of an active gene with time dependent propensities. For m mRNAs at time t, the propensities for transcription and degradation of an mRNA are given by kt(t) and mk d (t ) respectively. (b) Results of the stochastic simulation algorithm (SSA) fitted with the probability distribution P (m, t ) obtained from our model. The histogram was plotted using binned data from 1000 realizations of the SSA over two complete periods of the input signal. Inset: functional form of the input signal kt(t) which is a square-wave oscillating between k1 and k2 with k d = 0.05 (arb. time)−1.
∅ have the same meaning as in the HP model. A schematic of the process is shown in figure 1(a). The probability P (m, t ) of finding m mRNAs at time t obeys the master equation: ∂P ( m , t ) = kt (t ) [P (m − 1, t ) − P (m , t ) ] ∂t + k d (t ) [(m + 1) P (m + 1, t ) − mP (m , t ) ] .
(6)
⎛ dμ ⎞ ∂ϕ ∂ϕ + uϕ ⎜ + k d μ − kt ⎟ + k d u = 0. (7) ⎝ dt ⎠ ∂t ∂u
If we choose the function μ such that it satisfies dμ + k d μ − kt = 0 dt
(8)
∂ϕ ∂ϕ + kd u = 0. ∂t ∂u
(9)
The function μ can be obtained by solving equation (8) giving t
t
k d dτ
∫
= (1 + u)m0 .
kt e ∫0
τ
k d ds
dτ .
(10)
0
We can solve equation (9) by employing the method of characteristics obtaining 3
(12)
Since F = ϕ (u, t )e μ(t ) u and μ (0) = 0, it follows from equation (12) that F (u , 0) = ϕ (u , 0) = ψ (u) = (1 + u)m0 .
(13)
We note that the choice of ϕ (u , t ) = ψ (U ) = (1 + U )m0
(14)
satisfies the boundary condition given by equation (12) and is therefore a solution to equation (7). Assuming that we start with zero initial number of mRNAs such that m0 = 0, the solution to equation (6) with timedependent kt and kd can be written as F (z , t ) = e μ(z − 1).
By definition, 1 ∂m P (m, t ) = m ! ∂z m F (z , t )
z=0
(15)
satisfies P (m, t ) which generates an
NHPD: P (m , t ) =
then equation (7) reduces to
∑ m
δm, m0 z m = z
(5)
A trial-solution of the form F = ϕ (u, t )e μ(t ) u upon substitution in equation (6) yields
μ (t ) = e ∫0
t
m m0
⎛ ∂F ⎞ ∂F = u ⎜ kd − kt F ⎟ . ⎝ ∂u ⎠ ∂t
(11)
where ψ is any arbitrary function of U = ue− ∫0 kd dτ . If there are initially m0 mRNAs then the generating function F satisfies F (z , 0) = ∑P (m , 0) z m =
2.2. Exact solution By defining the generating function F (z , t ) as ∑m z mP (m, t ) and u = z − 1, we can convert equation (5) into a first order partial differential equation:
−
ϕ (u , t ) = ψ (U ),
1 −μ m e μ m!
(16)
with a time-dependent function μ (t ) given by equation (10). In the following sections we consider variations in time-dependent rates and study an ensemble of genetically identical cells and assume that each cell has exactly one copy of the gene. We also compare our exact results with numerical simulations utilizing a stochastic simulation algorithm (SSA) (see appendix A) which is a modification of Gillespie’s original algorithm [15, 16]. We choose, for example, kt(t) to be a square-wave alternating between two
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 2. Transcriptional time-series generated for two different time-periods (a) T = 10k d−1 and (b) T = k d−1, corresponding to low and high frequency oscillations in kt respectively and k d = 0.05 (arb. time)−1. For each time-period, two separate realizations of the SSA are shown. The simulation results are fitted with the theoretical curves (blue) which correspond to the average number of mRNAs obtained by integrating equation (17). The corresponding kt are also shown on top. (c) Plot of the mutual information (MI) of two time-series data (the transcription rate and the corresponding mRNA levels) for different values of kd, as the time-period T of the input signal kt is varied. Each curve is averaged over 100 realizations.
constant transcription rates k1 and k2 and verify that equation (16) gives the correct distribution. The result is plotted in figure 1(b).
3. Results and discussion 3.1. Low pass filtering Low pass filtering is a common survival strategy employed by gene regulatory networks and has been experimentally observed in synthetic transcriptional cascades [17] and in the metabolic system of yeasts subject to a periodically varying external glucose source [18]. The ability of gene transcription networks to respond reliably to a slowly varying signal, while ignoring fast fluctuations, confers a distinct fitness advantage. Maintaining a regulatory network that can interpret an external signal and act accordingly has its price. Regulation is only feasible when the benefit exceeds this cost of having a regulatory system in place. It is well known that short pulses in the input signal, whose duration is less than a critical value, can have a deleterious effect on growth [13]. The reason is that ultrashort pulses lead to a reduction in fitness because the gene products do not have sufficient time to reach levels in which the accumulated benefit exceeds the cost of production. Therefore, it is beneficial for an organism to filter out rapid fluctuations in the input signal. We show using a minimal stochastic model of gene transcription how this may be achieved in a living cell. To probe the dynamic response of our model, we first choose kt(t) to be a square-wave alternating between two constant transcription rates k1 and k2 with a time-period T and assume the degradation probability to be a constant kd. Thus equation (10) reduces to t
μ (t ) = e
−k d t
∫
kt (τ)e k d τdτ ,
(17)
0
4
where μ (t ) is the average number of mRNAs at time t. To see the temporal evolution of μ (t ), we then integrate equation (17) for two different time periods corresponding to low and high frequency oscillations in kt(t). The theoretical curves for each time-period are compared to two separate realizations of the SSA and the results are plotted in figure 2. We find that for high frequency oscillations in kt, the transcription network is insensitive to the variations in the transcription rate and behaves as a low-pass filter. Noting the formal equivalence of equation (8) with equations that describe driven filter circuits, we understand this result by comparing the ‘response time’ of the transcription network given by k d−1 with the period T of the oscillatory transcription rate. Just as in electronic filter circuits, if we squeeze in too many periods so that T is small, the network filters out the high frequency fluctuations in kt(t). The number of transcripts (output) encode the temporal variation of the reaction rates (input). The efficiency of the filtering property can be quantified by the rate of information transfer between the input and the output signal. A quantitative measure of the information transfer can be obtained by estimating the MI (see appendix B) between the transcription rate kt and the corresponding mRNA levels obtained from the SSA. When the system filters out the high-frequency fluctuations in the signal, the corresponding MI is expected to decrease. As the frequency decreases and the system becomes sensitive to the oscillations in kt, the MI increases reflecting this change. As seen in figure 2(c), when the period T is small, the information transfer is close to zero signifying that the system filters out the fluctuations in the transcription rate, and as T becomes large, the MI gradually increases and becomes close to 1. We also study the effect of different degradation rates kd on the information transfer and in figure 2(c),
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
we show the corresponding MI curves for three different values of kd corresponding to kd = 0.025, kd = 0.05, and kd = 0.1 (measured in arbitrary units of time). We see that as kd increases, the MI saturates faster. According to our previous analysis, this is to be expected since k d−1 is the ‘time-constant’ in our transcription filter. As kd increases, the system becomes more responsive to high frequency fluctuations in kt and the MI curves get steeper. Thus with the increase in kd, the sensitivity of the transcriptional network to fluctuations in kt increases. The effect of time-dependent degradation rates also has a substantial bearing on transcriptional dynamics. Our preliminary results with kd chosen to be a square-wave oscillating between k d1 = 0.005 and k d2 = 0.05 in phase with kt (also measured in the same units) shows some interesting features when we vary the time-period Td of kd keeping the period Tt of kt fixed at Tt = 200. The results are shown in figure 3. We see that when the ratio of the time-periods ρ = Td Tt equals 1, the burstiness in transcription disappears and mRNA expression levels resemble constitutive transcription. This is because the oscillations in both rates increase and decrease in tandem which help maintain an approximately constant transcript level. However, when ρ deviates from 1, there is a reappearance of bursts and the burst-frequency is roughly proportional to ρ. The dynamics can be more intricate if kt and kd have a non-zero relative phase. Further detailed studies need to be undertaken to elucidate the complex interplay between kt and kd. 3.2. Synchronous/asynchronous response in a cell population In order to test the system response for noisy signals, we next generate kt(t) so that it consists of a sequence of pulses in which the pulse durations d are drawn from an exponential distribution P (d) = λ exp (−λd). We explore two different signal modalities: one in which every copy of the system has a different uncorrelated transcription rate but the pulse widths for each signal are drawn from the same exponential distribution, and the other in which all copies of the system transcribe with identical transcription rates. The results of the simulation for the two different modalities are plotted in figure 4. In scenario (a), where the transcription rates are specific to each copy of the system with the only commonality being the average pulse duration, we see that the individual time-series are bursty but there is considerable heterogeneity in the burst timings (figures 4 (a1)–(a3)). This is due to the lack of correlation between the individual transcription rates and both the time-averaged and ensemble-averaged values of the Fano factor deviate from unity. This is somewhat similar to the ON–OFF model where switching of the gene on and off is seen as a stochastic process and the switching events between different copies of the 5
system are uncorrelated and we expect the ensemble and time averaged Fano factors to be equal. We see that in scenario (b), where all copies have the same kt, individual realizations of the SSA show synchronized bursting. The time-averaged Fano factor for each realization is super-Poissonian showing that the time-series is bursty but the ensemble-averaged Fano factor at each time instant is close to 1 (figures 4 (b1)–(b2)). This shows that although each copy of the system shows burstiness in time, at any given instant, the number of mRNAs produced is fairly homogeneous across the population. This is expected as the underlying distribution is an NHPD, as predicted by our model, and at each instant the average number of mRNAs has a simple Poisson distribution. The timeaveraged and ensemble-averaged values of the Fano factor do not match because the system is repeatedly perturbed and does not have time to relax.
4. Conclusions Sensory transcription networks are designed to respond to dynamically changing environments. We offer a new conceptual framework in which the environmental fluctuations are encoded in the transcription rate which varies in time. It is well known that the fluctuations in gene expression networks can have both extrinsic and intrinsic components [19, 20]. The time-dependent modulation of the transcription rate may arise from propagation of changes in upstream signaling as in fluctuations in activator concentrations in regulatory networks [21, 22] or it may be a result of intrinsic switching of the gene between ON–OFF states in the absence of any genetic regulation or external signaling [23]. It becomes difficult to deconvolve experimentally the contribution of each component [24] and in our model any change in the environment, extrinsic and/or intrinsic, is encoded in the transcription rate kt. We obtain an exact analytic solution of the chemical master equation resulting from a birth–death process of mRNA production and degradation with time-dependent kt , k d . Since the system is inherently stochastic, the reaction rates are taken as probabilities per unit time. We find that the probability mass density P (m, t ) follows an NHPD and we validate our result by stochastic simulations using a modified Gillespie algorithm [25]. For periodic variations in kt, we find that the system behaves as a low pass filter and the information transfer depends on the period of the signal kt and the response time k d−1 such that when T is small relative to k d−1, fluctuations in kt are filtered out and the MI becomes close to zero. We also test the system response for noisy signals consisting of a train of pulses in which the pulse widths are sampled from an exponential distribution. The two different signal modalities that we explore indicate that the system can have quite different responses
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 3. (a) Temporal variation in kt with period Tt = 200 . (b1), (c1) and (d1) are the temporal variations in kd with periods Td = 200 , Td = 400 and Td = 100 , respectively, and (b2), (c2) and (d2) are the corresponding transcript levels showing two separate realizations of the SSA for each kd (orange and green). The theoretical curves (shown in blue) are obtained by integrating equation (10) and time is in arbitrary units.
based on the temporal correlations between the transcription rates of different copies of the system. Our simulations show that a population of identical cells with the same transcription rate shows synchronized bursting and at any instant, the average number of mRNAs produced is Poisson distributed. This corresponds to genes that faithfully transmit fluctuations in the environment, such as the concentration of a given hormone [5]. Such synchrony in mRNA levels in cell populations are often observed during cellular stress response [4, 19]. However, if the transcription rates across the population lack temporal correlation, the synchrony in the burst timings is lost and the ensemble-averaged Fano factor deviates from the Poisson 6
case. This can be viewed as an extension of the ON– OFF model for finite times and with non-zero basal expression levels which result in phenotypic heterogeneity. With the advent of new experimental techniques, it is now possible to study transcription dynamics at the population level as well as in single cells [3]. Our work suggests that transcription is a nonequilibrium process and we need to exercise caution in quantifying statistical fluctuations, since depending on the burst synchrony/asynchrony, the time-averaged and the ensemble-averaged Fano factor may or may not be equal. Recent experimental studies have provided insight into the origins of burstiness where chromatin
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
Figure 4. (a1) Transcriptional time-series generated by the SSA where transcription rate kt has exponentially distributed pulse-widths with parameter λ = 0.5k d and kt stochastically fluctuates between k1 k d = 2 and k 2 k d = 16 . (a2) Another realization of the SSA with a kt which is temporally uncorrelated to the kt used in (a1), but has the same parameters λ, k1, k2 and kd. (a3) Ensemble-averaged Fano factor computed at different time-instants. (b1) Transcriptional time-series generated for an ensemble-wide transcription rate kt with a parameter λ = 0.5k d stochastically fluctuating between k1 k d = 2 and k 2 k d = 16 . Two different realizations of the SSA are shown. (b2) Ensemble-averaged Fano factor computed at different time-instants. In all of (a1–a3) and (b1–b2), mRNA time-series generated by the SSA are fitted with the theoretical value of the average mRNA number (blue) obtained by integrating equation (17). The ensemble Fano factor fluctuates in time around a mean value indicated by the red dotted lines. (c) Table summarizing the values of the time-averaged and ensemble-averaged Fano factors for the two different scenarios a and b.
remodeling [26] and DNA supercoiling [27] have been implicated as possible molecular mechanisms underlying transcriptional bursts. Here we consider a coarse-grained stochastic model and show how the transcriptional machinery may display bursty dynamics depending on the time-scale of fluctuation of the time-dependent transcription rate. Our study indicates that there exist a possible connection between bursty mRNA dynamics and information transfer due to transcription. Simple population dynamics models which treat transcription as a birth–death process can account for burstiness by introducing stochastic switching of the gene [6], but these models cannot account for sub-Poissonian or anti-bursty behavior which is also a recurring theme in diverse areas of biology [28–32]. The phenomenon of transcriptional antibursts has recently been observed experimentally in individual cells of E.coli., driven by the lar promoter, under weak and medium induction [32]. Simulations of our time-dependent model also indicate low Fano factors in certain time-scales which may be relevant in experimental situations but it remains to be ascertained if it is a true statistical phenomenon. 7
Further detailed studies need to be undertaken to investigate these issues in realistic gene regulatory networks.
Appendix A. Stochastic simulation algorithm The standard Gillespie algorithm [15, 16] is used to simulate chemical systems where the reaction propensities are time-independent. Since kt(t) in our model varies with time, we use a modified version of Gillespie’s first-reaction method to account for the time dependent transcription rates. If a time-dependent reaction rate k(t) is a piece-wise linear function with a finite number of discontinuities, following [25], we treat the discontinuities as detailed below. Let k(t) be discontinuous at a point tn so that k(t) satisfies ⎧ k< (t ) when t < tn k (t ) = ⎨ ⎩ k> (t ) when t > tn.
(A.1)
To simulate a discontinuous change in the rate, we change k(t) according to equation (A.1) when the sum
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
of the current simulation time and the next putative reaction time pass tn. k(t) is changed from k < (t ) to k > (t ) and the time is reset to t = tn . Consequently, our algorithm proceeds as follows. (1)Initialize the number of mRNAs. Set time t = 0. (2)Calculate the propensity for each chemical reaction. (3)For each reaction i generate a putative next reaction time τi . (4)Let τμ , corresponding to reaction μ, be the minimum τi . (5)Let tn be the time for the next discontinuous change in kt(t). If t + τμ < tn, then change the number of mRNAs appropriately for the occurrence of reaction μ. Change t to t + τμ. If t + τμ > tn, then change kt(t) accordingly. Set t = tn . (6)Go to step 2. To simulate arbitrary functions, we need to approximate the given function by a series of step functions or a piece-wise linear function and proceed as outlined above.
Appendix B. Mutual information The MI of two discrete random variables X and Y is a measure of the reduction in uncertainty about one random variable given knowledge of the other and can be defined as [33]: MI (X ; Y ) =
⎛ p (x , y ) ⎞ ⎟ , (B.1) ⎝ p (x ) p (y ) ⎠
∑ ∑ p (x, y) log ⎜ y∈Y x∈X
where p (x , y) is the joint probability distribution function of X and Y and p(x) and p(y) are the marginal probability distribution functions of X and Y respectively. If the base of the logarithm is 2, (the usual choice and the one we have used here), MI is measured in bits. Alternatively, MI of X and Y can be expressed as MI (X ; Y ) = H (X ) − H (X Y ),
(B.2)
where H(X) is the marginal entropy of X and H ((X ∣ Y )) is the conditional entropy of X conditioned on Y. Using the above definition, we generate the MI of two time-series, the signal kt(t) and the response m(t), where kt(t) and m(t) are the transcription rate and mRNA level respectively at time t.
References [1] van Kampen N G 2007 Stochastic Processes in Physics and Chemistry (Amsterdam: North-Holland, Elsevier) [2] Thattai M and van Oudenaarden A 2001 Intrinsic noise in gene regulatory networks Proc. Natl Acad. Sci. USA 98 8614–9
8
[3] Sanchez A and Golding I 2013 Genetic determinants and cellular constraints in noisy gene expression Science 342 1188–93 [4] Weake V M and Workman J L 2010 Inducible gene expression: diverse regulatory mechanisms Nat. Rev. Genetics 11 426–37 [5] Lionnet T and Singer R H 2012 Transcription goes digital EMBO Rep. 13 313–21 [6] Peccoud J and Ycart B 1995 Markovian modeling of geneproduct synthesis Theor. Population Biol. 48 222–34 [7] Suter D M, Molina N, Gatfield D, Schneider K, Schibler U and Naef F 2011 Mammalian genes are transcribed with widely different bursting kinetics Science 332 472–4 [8] Golding I, Paulsson J, Zawilski S M and Cox E C 2005 Realtime kinetics of gene activity in individual bacteria Cell 123 1025–36 [9] Chubb J R, Trcek T, Shenoy S M and Singer R H 2006 Transcriptional pulsing of a developmental gene Curr. Biol. 16 1018–25 [10] So L-H, Ghosh A, Zong C, Sepúlveda L A, Segev R and Golding I 2011 General properties of transcriptional time series in Escherichia coli Nat. Genetics 43 554–60 [11] Raj A, Peskin C S, Tranchina D, Vargas D Y and Tyagi S 2006 Stochastic mRNA synthesis in mammalian cells PLoS Biol. 4 1707–19 [12] Gomperts B D, Kramer I M and Tatham P E R 2002 Signal Transduction vol 18 (New York: Academic) [13] Alon U 2007 An Introduction to Systems Biology: Design Principles of Biological Circuits vol 10 (Boca Raton, FL: CRC) [14] Bothma J P, Garcia H G, Esposito E, Schlissel G, Gregor T and Levine M 2014 Dynamic regulation of eve stripe 2 expression reveals transcriptional bursts in living Drosophila embryos Proc. Natl Acad. Sci. USA 7 1–6 [15] Gillespie D T 1976 A general method for numerically simulating the stochastic time evolution of coupled chemical reactions J. Comput. Phys. 22 403–34 [16] Gillespie D T 1977 Exact stochastic simulation of coupled chemical reactions J. Phys. Chem. 81 2340–61 [17] Hooshangi S, Thiberge S and Weiss R 2005 Ultrasensitivity and noise propagation in a synthetic transcriptional cascade Proc. Natl Acad. Sci. USA 102 3581–6 [18] Bennett M R, Pang W L, Ostroff N A, Baumgartner B L, Nayak S, Tsimring L S and Hasty J 2008 Metabolic gene regulation in a dynamically changing environment Nature 454 1119–22 [19] Raj A and van Oudenaarden A 2008 Nature, Nurture, or chance: stochastic gene expression and its consequences Cell 135 216–26 [20] Tkačik G, Gregor T and Bialek W 2008 The role of input noise in transcriptional regulation PLoS One 3 e2774 [21] Elowitz M B and Leibler S 2000 A synthetic oscillatory network of transcriptional regulators Nature 403 335–8 [22] Locke J C W, Young J W, Fontes M, Jimenez M J H and Elowitz M B 2011 Stochastic pulse regulation in bacterial stress response Science 334 366–9 [23] Harper C V et al 2011 Dynamic analysis of stochastic transcription cycles PLoS Biol. 9 e1000607 [24] Hilfinger A and Paulsson J 2011 Separating intrinsic from extrinsic fluctuations in dynamic biological systems Proc. Natl Acad. Sci. USA 108 12167–72 [25] Shahrezaei V, Ollivier J F and Swain P S 2008 Colored extrinsic fluctuations and stochastic gene expression Mol. Syst. Biol. 4 196 [26] Voss T C and Hager G L 2014 Dynamic regulation of transcriptional states by chromatin and transcription factors Nat. Rev. Genetics 15 69–81 [27] Chong S, Chen C, Ge H and Xie X S 2014 Mechanism of transcriptional bursting in bacteria Cell 158 314–26 [28] Gur M, Beylin A and Snodderly D M 1997 Response variability of neurons in primary visual cortex (V1) of alert monkeys J. Neurosci. 17 2914–20 [29] Kara P, Reinagel P and Reid R C 2000 Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons Neuron 27 635–46
Phys. Biol. 12 (2015) 016015
S Nandi and A Ghosh
[30] Maimon G and Assad J A 2009 Beyond poisson: increased spiketime regularity across primate parietal cortex Neuron 62 426–40 [31] Scaglione A, Moxon K A, Aguilar J and Foffani G 2011 Trialto-trial variability in the responses of neurons carries information about stimulus location in the rat whisker thalamus Proc. Natl Acad. Sci. USA 108 14956–61
9
[32] Kandhavelu M, Häkkinen A, Yli-Harja O and Ribeiro A S 2012 Single-molecule dynamics of transcription of the lar promoter Phys. Biol. 9 026004 [33] Shannon C E 1948 A mathematical theory of communication Bell Syst. Tech. J. 27 379–423
