This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.086157v1 91/12/4350    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by van Zon, J. S. Articles by ten Wolde, P. R. Search for Related Content PubMed PubMed Citation Articles by van Zon, J. S. Articles by ten Wolde, P. R.

Diffusion of Transcription Factors Can Drastically Enhance the Noise in Gene Expression

Jeroen S. van Zon ^*, Marco J. Morelli , Sorin Tnase-Nicola  and Pieter Rein ten Wolde 

^* Division of Physics and Astronomy, Vrije Universiteit, Amsterdam, The Netherlands; and FOM Institute for Atomic and Molecular Physics, Amsterdam, The Netherlands

Correspondence: Address reprint requests to Pieter Rein ten Wolde, E-mail: tenwolde{at}amolf.nl.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
We study by Green's Function Reaction Dynamics the effect of the diffusive motion of repressor molecules on the noise in mRNA and protein levels for a gene that is under the control of a repressor. We find that spatial fluctuations due to diffusion can drastically enhance the noise in gene expression. After dissociation from the operator, a repressor can rapidly rebind to the DNA. Our results show that the rebinding trajectories are so short that, on this timescale, the RNA polymerase (RNAP) cannot effectively compete with the repressor for binding to the promoter. As a result, a dissociated repressor molecule will on average rebind many times, before it eventually diffuses away. These rebindings thus lower the effective dissociation rate, and this increases the noise in gene expression. Another consequence of the timescale separation between repressor rebinding and RNAP association is that the effect of spatial fluctuations can be described by a well-stirred, zero-dimensional, model by renormalizing the reaction rates for repressor-DNA (un) binding. Our results thus support the use of well-stirred, zero-dimensional models for describing noise in gene expression. We also show that for a fixed repressor strength, the noise due to diffusion can be minimized by increasing the number of repressors or by decreasing the rate of the open complex formation. Lastly, our results emphasize that power spectra are a highly useful tool for studying the propagation of noise through the different stages of gene expression.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
Cells process information from the outside and regulate their internal state by means of proteins and DNA that chemically and physically interact with one another. These biochemical networks are often highly stochastic because in living cells the reactants often occur in small numbers (1,2). This is particularly important in gene expression (3–10), where transcription factors are frequently present in copy numbers as low as tens of molecules per cell. Although it is generally believed that biochemical noise can be detrimental to cell function (8), it is increasingly becoming recognized that noise can also be beneficial to the organism (9). Understanding noise in gene expression is thus important for understanding cell function, and this observation has recently stimulated much theoretical and experimental work in this direction (8–10). However, the theoretical analyses usually employ the zero-dimensional chemical master equation (10–12). This approach takes into account the discrete character of the reactants and the probabilistic nature of chemical reactions. It does assume, however, that the cell is a "well-stirred" reactor, in which the particles are uniformly distributed in space at all times; the reaction rates only depend upon the global concentrations of the reactants and not upon the spatial positions of the reactant molecules. Yet, to react, reactants first have to move toward one another. They do so by diffusion, or, in the case of eukaryotes, by a combination of diffusion and active transport. Both processes are stochastic in nature and this could contribute to the noise in the network. Here, we study by computer simulation the expression of a single gene that is under the control of a repressor R in a spatially resolved model. We find that at low repressor concentration, i.e., [R] < 50 nM, the noise in gene expression is dominated by the noise arising from the diffusive motion of the repressor molecules. Our results thus show that spatial fluctuations of the reactants can be an important source of noise in biochemical networks. Our analysis also reveals that the effects of diffusion can nevertheless be described by a well-stirred model, provided that the reaction rates of repressor-DNA (un)binding are properly renormalized.

The simulations show that in gene expression significant fluctuations occur on both short and long length- and timescales. As expected from earlier work (13–15), the fluctuations on long timescales are predominantly due to protein degradation; we assume that proteins are degraded by dilution, which means that the relaxation rate of this process is on the order of 1 h. Our results, however, also elucidate an important process on much shorter length- and timescales. It is associated with the competition between repressor and RNA polymerase (RNAP) for binding to the promoter. When a repressor molecule dissociates from the DNA, it can rebind very rapidly: in our model, which neglects one-dimensional diffusion along the DNA, it can rebind on a timescale of milliseconds, or less. This timescale is much shorter than that on which the RNAP binds to the promoter, which is on the order of 0.01–0.1 s. Hence, when a repressor molecule has just dissociated, the probability that a RNAP will bind before the repressor molecule rebinds, is very small. This has two important consequences. The first is that a repressor molecule will on average rebind many times before it eventually diffuses away from the promoter and a RNAP molecule, or another repressor molecule, can bind to the promoter. This decreases the effective dissociation rate, which increases the noise in gene expression.

The second consequence of the rapidity of the rebindings is that noise propagation during gene expression can be described by a well-stirred, zero-dimensional model. In the commonly used zero-dimensional models, chemical reactions are exponentially distributed in time. In a spatially resolved model, the distribution of repressor-DNA association times deviates markedly from Poisson statistics. Although at long timescales the distribution is exponential, at short timescales it is algebraic, due to the diffusive nature of the rebinding trajectories. However, these repressor rebindings are so fast that they do not significantly affect the dynamics of RNAP-DNA association; the latter is only affected by the repressor-DNA (un)binding dynamics at longer timescales, which obey Poisson statistics. The reason that the effect of spatial fluctuations on noise in gene expression can be described by a zero-dimensional model is thus a separation of timescales. In fact, it is conceivable that in a more realistic model of gene expression, which includes one-dimensional sliding along the DNA, the timescale of repressor rebinding is not separated from that of the RNAP dynamics. Under these conditions, the effect of spatial fluctuations might be detected in the statistics of mRNA production. However, as we discuss in the Discussion and Outlook section, the noise strength (variance) of the mRNA level can probably still be described by a zero-dimensional model, because the timescale of the spatial fluctuations, even in those more refined models, is expected to be still shorter than the typical life time of an mRNA molecule.

Because fluctuations in gene expression span orders of magnitude in length- and timescales, the simulation technique should be sufficiently detailed to resolve the events at short length- and timescales, yet also efficient enough to access the long length- and timescales. Recently, several simulation techniques have been developed for the stochastic modeling of reaction-diffusion systems (16,17). These techniques, however, do not satisfy both criteria: they either describe the system in a coarse-grained way, i.e., on the level of local concentrations rather than single particles (16,17) or are too slow to accurately model the dynamics on the long timescales (18). Our simulations have been made possible via the use of our recently developed Green's function reaction dynamics (GFRD) algorithm (19,20). GFRD is an event-driven algorithm that uses Green's functions to combine in one step the propagation of the particles in space with the reactions between them. The event-driven nature of the algorithm makes it particularly useful for problems, such as gene expression, in which the events are distributed over a wide range of length- and timescales: the algorithm takes small steps when the reactants are close to each other—such as when a repressor molecule has just dissociated from the DNA—whereas it takes large jumps in time and space when the molecules are far apart from each other, like when the repressor molecule has eventually diffused away from the promoter. The event-driven nature of GFRD makes it orders of magnitude more efficient than brute-force particle-based algorithms (20) and this has allowed us to simulate gene expression on the relevant biological timescales of hours.

Several publications (21–29) have discussed the effect of fluctuations in the binding of transcription factors to their site on the DNA (called operator) on the noise in gene expression. Most of these models are relatively simple, ignoring, for instance, production of mRNA (23–26,29). Moreover all these studies, with the exception of (24,28), ignore the role of the spatial fluctuations of the transcription factors. Our aim is to study gene expression in a biologically meaningful model. We have therefore constructed a rather detailed model, although we will also use minimal models that can be studied analytically to interpret the simulation results. The full model, which is described in the next section, contains the diffusive motion of repressor molecules, open complex formation, promoter clearance, transcription elongation, and translation (30).

In the section "Simulation results: dynamics and noise", we discuss the simulation results for both the noise in mRNA and protein level. The results reveal that for [R] < 50 nM, the noise in the spatially resolved model can be more than five times larger than the noise in the well-stirred model. We also show that a cell could minimize the effect of spatial fluctuations, either by tuning the open complex formation rate or by changing the number of repressors and their affinity for the binding site on the DNA. In the section "Simulations results: operator binding", we elucidate the origin of the enhanced noise in the spatially resolved model. In the subsequent section, we show that in the model employed here the effect of spatial fluctuations can be quantitatively described by a well-stirred model in which the reaction rates for repressor binding and unbinding are appropriately renormalized; however, as alluded to above, and as we will discuss in more detail in the last section, we expect that in a more refined model the effect of diffusion will be more complex, impeding such a simplified description. In the section "Power spectra", we discuss how the operator state fluctuations propagate through the different stages of gene expression using power spectra for the operator state, elongation complex, mRNA, and protein. The results show that these power spectra are highly useful for unraveling the dynamics of gene expression. We hope that this stimulates experimentalists to measure power spectra of not only mRNA and protein levels (31), but also of the dynamics of transcription initiation and elongation using, e.g., magnetic tweezers (32). As we argue in the last section, such experiments should make it possible to determine the importance of spatial fluctuations for the dynamics of gene expression.

	MODEL

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
Diffusive motion of repressors
We explicitly simulate the diffusive motion of the repressor molecules in space. However, since the experiments of Riggs et al. (33) and the theoretical work of Berg et al. (34), it is well known that proteins could find their target sites via a combination of one-dimensional (1D) sliding along the DNA and three-dimensional (3D) diffusion through the cytoplasm—"hopping" or "jumping" from one site on the DNA to another. This mechanism could speed up the search process and make it faster than the rate at which particles find their target by free 3D diffusion; this rate is given by k = 4D₃[R], where is the cross section, which is on the order of a protein diameter or DNA diameter, D₃ is the diffusion constant of the protein in the cytoplasm, and [R] is the concentration of the (repressor) protein. However, although it is clear that the mechanism of 3D diffusion and 1D sliding could potentially speed up the search process, whether this mechanism in living cells indeed drastically reduces the search time is still under debate (35). In this context, it is instructive to discuss the two main results of recent studies on this topic (35–40). The first is that the mean search time is given by (40)

(1)

where L is the total length of the DNA, is the average distance over which the protein slides along the DNA before it dissociates, D₁ is the diffusion constant for sliding, r is the typical mesh size in the nucleoid (the characteristic distance between two segments on the DNA (40)), and D₃ is the diffusion constant in the cytoplasm. This formula has a clear interpretation (40): ²/D₁ is the sliding time, r²/D₃ is the time spent on 3D diffusion, the sum of these terms is thus the time to perform one round of sliding and diffusion, and L/ is the total number of rounds needed to find the target. The other principal result is that the search time is minimized when the sliding distance is

(2)

Under these conditions, a protein spends equal amounts of time on 3D diffusion and 1D sliding (a protein is thus half the time bound to the DNA). Equation 2 is a useful result because it shows that the average sliding distance depends upon the ratio of diffusion constants and on the typical mesh size in the nucleoid. If we now assume that D₁ and D₃ are equal (which is not obvious given that proteins bind relatively strongly to DNA—D₁ could thus very well be much smaller than D₃) and if we take the mesh size to be given by (40), where v 1 µm³ is the volume of an Escherichia coli cell and L 10³ µm, we find that 10 nm (30 bp). This corresponds to the typical diameter of a protein or DNA double helix and is thus not very large. Interestingly, recent experiments seem to confirm this: experiments from Halford et al. on restriction enzymes (EcoRV and BbcCI) with a series of DNA substrates with two target sites and varying lengths of DNA between the two sites, suggest that under the in vivo conditions, sliding is indeed limited to relatively short distances, i.e., to distances less than 50 bp (16 nm) (41,42).

Now, it should be realized that on lengthscales beyond the sliding length, the motion is essentially 3D diffusion: the sliding/hopping mechanism corresponds to 3D diffusion with a jump distance given by the sliding distance (36). Moreover, since the sliding distance is only on the order of a particle diameter, as discussed above, we have therefore decided to model the motion of the repressor molecules as 3D diffusion. But it should be remembered that on lengthscales shorter than 10–30 nm, this approach is not correct. As we discuss in the "Discussion and Outlook" section, this might have significant implications for the importance of spatial fluctuations for the noise in gene expression.

Transcription and translation
Most repressors bind to a site that (partially) overlaps with the core promoter—the binding site of the RNA polymerase (RNAP). When a repressor molecule is bound to its operator site, it prevents RNAP from binding to the promoter, thereby switching off gene expression. Only in the absence of a repressor on the operator site can RNAP bind to the promoter and initiate transcription and translation, ultimately resulting in the production of a protein. We model this by the following reaction network:

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

Equations 3 and 4 describe the competition between the binding of the repressor R and the RNAP molecules Rp to the promoter (O is the operator site). In our simulation we fix the binding site O in the center of a container with volume V = 1 µm³, comparable to the volume of a single E. coli cell. We simulate both the operator site O and the repressor molecules as spherical particles with diameter = 10 nm. The operator site O is surrounded by N_R repressor molecules that move by free 3D diffusion (see previous section) with an effective diffusion constant D = 1 µm²s^–1, as has been reported for proteins of a similar size (43). The intrinsic forward rate k_fR = 6 x 10⁹M^–1s^–1 for the repressor particles R at contact is estimated from the Maxwell-Boltzmann distribution (19). The backward rate k_bR depends on the interaction between the DNA binding site of the repressor and the operator site on the DNA and varies greatly between different operons, with stronger repressors having a lower k_bR. In our simulations, we vary k_bR between 1 and 0.01 s^–1, as discussed in more detail below. The concentration of RNAP is much higher than that of the repressor (44). Because of this, we treat the RNAP as distributed homogeneously within the cell and we do not to take diffusion of RNAP into account explicitly. Instead, RNAP associates with the promoter with a diffusion-limited rate k_fRp = 4D[Rp]. In our simulations, the concentration of free RNAP is [Rp] = 0.5 µM (44), leading to a forward rate k_fRp = 38 s^–1. Finally, the backward rate k_bRp = 0.5 is determined such that K_eq = 4D/k_bRp = 1.4 x 10⁹ M^–1 (45).

Transcription initiation is described by Eqs. 5 and 6. Before productive synthesis of RNA occurs, first the RNAP in the RNAP-promoter complex ORp unwinds approximately one turn of the promoter DNA to form the open complex ORp*. The open complex formation rate k_OC has been measured to be on the order of 0.3–3 s^–1 (32). We approximate open complex formation as an irreversible reaction. Some experiments find this step to be weakly reversible (32). However, adding a backward reaction to the model did not change the dynamics of the system in a qualitative way, as long as the backward rate is smaller than k_OC, which is in agreement with experimental results. After open complex formation, RNAP must first escape the promoter region before another RNAP or repressor can bind. Because elongation occurs at a rate of 50–100 nucleotides per second and between 30 and 60 nucleotides must be cleared by RNAP before the promoter is accessible, a waiting time of t_clear = 1s is required before another binding can occur. Since promoter clearance consists of many individual elongation events that obey Poisson statistics individually, we model the step as one with a fixed time delay t_clear, not as a Poisson process with rate 1/t_clear.

Equations 7–11 describe the dynamics of mRNA and protein numbers. After clearing the promoter region, RNAP starts elongation of the transcript T. As for clearance, the elongation step is modeled as a process with a fixed time delay t_elon = 30 s, corresponding to an elongation rate of 50–100 nucleotides per second and a 1500 bp gene. When an mRNA M is formed, it can degrade with a rate k_dm. Here, the mRNA degradation rate is determined by fixing the average mRNA concentration in the unrepressed state, as described below. Furthermore, an mRNA molecule can form an mRNA-ribosome complex M_ribo and start translation. We assume that b = 5 proteins are produced on average from a single mRNA molecule (7), so that the start of translation occurs at a rate k_ribo = bk_dm. After a fixed time delay, t_trans = 30 s, a protein P is produced. The mRNA is available for ribosome binding immediately after the start of translation. Due to the delay in protein production, M can start to be degraded, while the mRNA-ribosome complex M_ribo is still present; M thus represents the mRNA leader region rather than the entire mRNA molecule. Finally, the protein P degrades at a rate k_dp, which is determined by the requirement that the average protein concentration in the unrepressed state has a desired value, as we describe now.

We vary the free parameters in the reaction network described in Eqs. 3–11, N_R, k_bR, k_dm, k_dp, in the following way: first, we choose the concentration of mRNA and protein in the absence of repressor molecules. In this case, tuning of the concentrations is most straightforward by adjustment of the mRNA and protein decay rates k_dm and k_dp. For the above reaction network one can show that the average mRNA number N_M and protein number N_P is given by

(12)

(13)

where K₁ = k_fRp/(k_bRp + k_OC), K₂ = k_fR/k_bR, K₃ = k_OCt_clear, K₄ = k_OC/k_dm, and K₅ = k_ribo/k_dp are equilibrium constants, V is the volume of the cell, and N_R is the total number of repressors. The unrepressed state corresponds to N_R = 0. In our simulations, we fix the mRNA and protein numbers in the unrepressed state at N_M = 50 and N_P = 2 x 10⁵. The mRNA and protein decay rates then follow straightforwardly from Eqs. 12 and 13: the mRNA degradation rate is k_dm = 0.019 s^–1 (46) and the protein degradation rate is k_dp = 2.4 x 10^–4s^–1; the latter corresponds to protein degradation by dilution with a cell cycle time of 1 h.

Next, we determine by what factor these concentrations should decrease in the repressed state. This can be done by changing the number of repressors, N_R, and the repressor backward rate, k_bR. We define the repression level f as the transcription initiation rate in the absence of repressors, divided by the initiation rate in the repressed state (47). For a repression level f, the concentration of mRNA and proteins in the repressed state is a fraction 1/f of the concentration in the unrepressed state and it follows that

(14)

Thus, a fixed repression level f does not specify a unique combination of N_R and k_bR: increasing the number of repressors twofold, while also increasing the repressor backward rate by the same factor, gives the same repression level. This means that the cell can control mRNA and protein levels in the repressed state either by having a large number of repressors that stay on the DNA for a short time or by having a small number of repressors, possibly even one, that stay on the DNA for a long time. Even though it is conceivable that the latter is preferable for economic reasons, there is no difference between the two extremes in terms of the average gene expression. In our simulations, we vary N_R and k_bR, but use a fixed repression level f = 100. Consequently, in the repressed state, on average N_M = 0.5 and N_P = 200.

Lastly, we would like to emphasize that, while all reaction rates were, as much as possible, taken from experiments, it should be realized that the measured rates might not be very precise. However, we believe that this does not affect the main conclusions of our work.

	SIMULATION TECHNIQUE

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
We simulate the above reaction network using GFRD (19,20). GFRD is an event-driven algorithm, which combines in one step the propagation of the particles in space with the reactions between them. The main idea is to determine at each iteration of the simulation, a maximum time step, such that only single particles or pairs of particles have to be considered. For these particles, the Smoluchowski equation (48) can be solved exactly using Green's functions. For each single particle, the Green's function is just the Gaussian distribution function p₁(r, t|r₀, t₀), which yields the probability that, given that the particle is at point r₀ at time t₀, it is at position r at a later time t. For each pair of particles, two Green's functions are obtained: one for their center-of-mass, and one for their interparticle vector r; the latter, p₂(r, t|r₀, t₀), yields the probability that the interparticle vector r₀ at time t₀ becomes r at a later time t. Importantly, the interparticle Green's function does not only take into account the diffusion of the particles, but also the reactions between them. This makes it possible to derive for each pair of particles the propensity function q(t|r₀), which yields the probability that the pair will react for the first time at time t, given that the particles were separated by a distance r₀ initially. The propensity functions, then, can be used to set up an event-driven algorithm, quite analogous to kinetic Monte Carlo algorithms for zero-dimensional master equations, such as the Gillespie algorithm (49). The event-driven nature allows GFRD to make large jumps in time and space when the particles are far apart from each other, making it up to five orders of magnitude more efficient than brute-force Brownian dynamics. For details of the algorithm, in particular on how the Green's functions and the propensity functions are derived, we refer to Refs. (19,20).

As discussed above, only the operator site O and the repressor particles R are simulated in space. All other reactions are assumed to occur homogeneously within the cell and are simulated according to the well-stirred model (49) or with fixed time delays for reaction steps involving elongation. A few modifications with respect to the algorithm described in van Zon and ten Wolde (19,20) are implemented to improve simulation speed. First, we neglect excluded volume interactions between repressor particles mutually because the concentration of repressor is very low. This means that the only potential reaction pairs we consider are operator-repressor pairs. Secondly, we use periodic boundary conditions instead of a reflecting boundary, which leads to a larger average time step. Because the operator site O is both small compared to the volume of the cell and is far removed from the cell boundary, this has no effect on the dynamics of the system. Finally, because the repressor backward rate k_bR is rather small, the operator site can be occupied by a repressor for a time long compared to the average simulation time step. If the repressor is bound to the operator site longer than a time L²/6D, where L is the length of the sides of our container, the other repressor molecules diffuse on average from one side of the box to the other. Consequently, when the repressor eventually dissociates from the operator site, the other repressor molecules have lost all memory of their positions at the time of repressor binding. Here, when a repressor will dissociate after a time longer than L²/6D, we do not propagate the other repressors with GFRD, but we only update the master equation and fixed delay reactions. We update the positions of the free repressors at the moment that the operator site becomes accessible again, by assigning each free repressor molecule a random position in the container; the dissociated repressor is put at contact with the operator site. We see no noticeable difference between this scheme and results obtained by the full GFRD algorithm described previously (19,20).

To obtain accurate statistics, especially upon notoriously difficult quantities such as noise and power spectra, very long simulations were performed. A total number of 24 simulations were performed, one for each combination of parameter values (N_R, k_OC). A single simulation took on average 24 h of CPU time on a Pentium IV 3.0 GHz processor.

	SIMULATION RESULTS: DYNAMICS AND NOISE

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
To study the effect of spatial fluctuations on the repression of genes, we simulate the reaction network described in Eqs. 3–11 both by GFRD—thus explicitly taking into account the diffusive motion of the repressor particles—and according to the well-stirred model, where the repressor particles are assumed to be homogeneously distributed in space and the dynamics depends only on the concentration of repressor. In Fig. 1 we show the behavior of mRNA and protein numbers for a system with open complex formation rate k_OC = 30 s^–1 and with varying numbers of repressors N_R. We keep the repression factor fixed at f = 100 so that with increasing N_R the repressor backward rate k_bR is also increased, i.e., repressor particles are bound to the DNA for a shorter time.

View larger version (32K): [in this window] [in a new window]   FIGURE 1  Dynamics of mRNA and protein numbers in the repressed state for different number of repressors NR. The number of mRNA and protein molecules is shown for simulations with GFRD (black line) and according to the zero-dimensional master equation, the well-stirred model (WS, gray line). In the GFRD simulation, diffusion of repressor particles is explicitly included. (a, b) NR = 5. (c, d) NR = 20. (e, f) NR = 80. In general, there is a dramatic difference in dynamics due to the spatial fluctuations of the repressor molecules. This difference becomes more pronounced as the number of repressors decreases. However, we find that in all cases NM = 0.5 and NP = 200, on average. The rate of open complex formation kOC = 30 s–1 and, when NR varies, the repression strength is kept constant at f = 100 by changing kbR (see Eq. 14).

When we increase the number of repressors N_R and change k_bR in such a way that the repression level f remains constant, we find that both for GFRD and the well-stirred model the fluctuations in mRNA and protein number decrease. In the absence of spatial fluctuations this effect is minor, but for GFRD this decrease is sharp: for a large number of repressors, the burst in mRNA become both weaker and more frequent. This in turn leads to smaller peaks and shorter periods of exponential decay in protein numbers. In fact, as N_R is increased both approaches converge to the same behavior. At around N_R 100, the dynamics of the protein number is similar for the well-stirred model and the spatially resolved model. The same happens for mRNA number when N_R 500.

In Fig. 2, we quantify the noise in mRNA and protein number, defined as standard deviation divided by the mean, while we change the number of repressors, N_R. As we keep the amount of repression fixed at f = 100, we simultaneously vary the backward rate k_bR according to Eq. 14. When all parameters are the same, the noise for the GFRD simulation, including the diffusive motion of the repressors, is always larger than the noise for the well-stirred model, where the diffusive motion is ignored. In both cases, the noise decreases when the number of repressors is increased and the repressor backward rate becomes larger. This is consistent with the mRNA and protein tracks shown in Fig. 1. We also investigated the effect of changing the open complex formation rate k_OC. In nature, this rate can be tuned by changing the basepair composition of the promoter region on the DNA. When we change k_OC, we change the mRNA decay rate k_dm so that the average mRNA and protein concentrations remain unchanged (see section "Transcription and translation"). We find that when k_OC is lowered, the fluctuations in mRNA and protein levels are sharply reduced. When k_OC is much larger than the RNAP backward rate k_bRp = 0.5 s^–1, almost every RNAP binding to the promoter DNA will result in transcription of an mRNA. For k_OC smaller than k_bRp, RNAP binding will lead to transcription only infrequently. As a consequence, the operator filters out part of the fluctuations in RNAP binding due to the diffusive motion of the repressor particles, leading to the decrease in noise observed in Fig. 2. This shows that the open complex formation rate plays a considerable role in controlling noise in gene expression.

View larger version (19K): [in this window] [in a new window]   FIGURE 2  Noise in (a) mRNA number and (b) protein number as a function of the number of repressors NR and for constant repression factor f = 100. Data obtained by GFRD simulation is shown for kOC = 0.3(), 3(), and 30(*) s–1. Noise levels for the well-stirred model (WS) are shown as gray lines and those for the well-stirred model with reaction rates renormalized according to Eqs. 16 and 17 (the renormalized well-stirred (RWS) model) are shown as black lines, both for kOC = 0.3 (solid lines), 3 (dashed lines), and 30 (dotted lines) s–1. Only when the reaction rates are properly renormalized does the noise in the well-stirred model agree well with the noise in the GFRD simulations, which include the effect of diffusion. (Insets) Noise levels as a function of kOC for NR = 5. Symbols indicate results for GFRD and lines are results for the chemical master equation with renormalized reaction rates (RWS model).

	SIMULATIONS RESULTS: OPERATOR BINDING

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

View larger version (39K): [in this window] [in a new window]   FIGURE 3  Dynamics of repressor binding for a repression factor of f = 2 and NR = 5. (a) The OR-bias for GFRD (black line) and the well-stirred model (gray line). The OR-bias is defined as the fraction of time a repressor is bound to the operator site in the last 50 s. When the diffusive motion of repressor molecules is included (black line), the OR-bias switches between periods where repressors are continuously bound to or absent from the DNA for long times. (b, c) Time trace of the occupancy of the operator site by repressor molecules. When OR = 1 a repressor is bound to the operator site and OR = 0 indicates either a free operator site or one with RNAP bound. For the GFRD simulations, an initial binding is followed by several rapid rebindings, whereas for the well-stirred model binding and rebinding is much more unstructured. Note that here, for reasons of clarity, f = 2 instead of f = 100 as used in the text and Figs. 1 and 2.

The occurrence of rapid rebindings is intimately related to the nature of diffusion. When diffusion and the positions of the reactants are ignored all dynamics is based only on the average concentration of the reactants. As a consequence, when in this approach a repressor dissociates from the operator site, the probability of rebinding depends only on the concentration of repressor in the cell. On the level of actual positions of the reactants, this amounts to placing the repressor at a random position in the container. The situation is very different for the GFRD approach, where the positions of the reactants are taken into account. After a dissociation from the operator site, the repressor particle is placed at contact with the operator site. Because of the close proximity of the repressor to its binding site, it has a high probability of rapidly rebinding to, and only a small probability of diffusing away from, the binding site. At the same time, when the repressor eventually diffuses away from the operator site, the probability that the same, or more likely, another repressor diffuses to and binds the operator site is much smaller than the probability of binding in the well-stirred model, as will be shown quantitatively in the section "Two-step kinetic scheme". This results in the behavior observed in Fig. 3 B.

It can now be understood that the bursts in mRNA production correspond to the prolonged absence of repressor from the operator site compared to the well-stirred model. Especially for low repressor concentrations, these periods of absence can be long enough that the concentration of mRNA reaches values comparable to those in the unrepressed state for brief periods of time. When a repressor binds to the operator site, due to the rapid rebindings it will remain bound effectively for a time much longer than the mRNA lifetime, leading to long periods where mRNA is absent in the cell. This shows that under these conditions spatial fluctuations and not stochastic chemical kinetics are the dominant contribution to the noise in mRNA and protein numbers in the repressed state.

	TWO-STEP KINETIC SCHEME

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
In the current model, the average repressor concentration profile is uniform. It is therefore natural to investigate to what extent the effect of diffusion on the repressor dynamics can be described by a zero-dimensional, well-stirred, model, via the following two-step kinetic scheme (50,51):

(15)

The first step in Eq. 15 describes the diffusion of repressor to the operator site resulting in the encounter complex OR, with the rates k₊ and k_– depending on the diffusion coefficient D and the size of the particles. The next step describes the subsequent binding of repressor to the DNA. In this case the rates are related to the microscopic rates defined in Eq. 3. When the encounter complex is assumed to be in steady state, the two-step kinetic scheme can be mapped onto the reaction described in Eq. 3, but with effective rate constants k'_fR = k₊k_a/(k_– + k_a) and k'_bR = k_–k_d/(k_– + k_a) (50). The two-step kinetic scheme should yield the same average concentrations as the scheme in Eq. 3, so that the equilibrium constant K = k_a/k_d = k'_fR/k'_bR = k_fR/k_bR, where k_fR and k_bR are the reaction rates defined in Eq. 3.

It is possible to express the effective rate constants k'_fR and k'_bR in terms of the microscopic rate constants k_fR and k_bR. For the setup used here, where a single operator O is surrounded by a homogeneous distribution of repressor R, the rate k₊ follows from the solution of the steady-state diffusion equation with a reactive boundary condition with rate k = k_a at contact (48,51) and is given by the diffusion-limited reaction rate k_D = 4D. The rates k_– and k_a depend on the exact definition of the encounter complex OR. It is natural to identify the rate k_d with the intrinsic dissociation rate k_bR, thus k_d = k_bR. From these expressions for k₊ and k_d and the requirement that the equilibrium constant should remain unchanged, one finds that k_a/k_– = k_fR/k_D. Using this result one obtains k'_fR = k_Dk_fR/(k_D + k_fR) and k'_bR = k_Dk_bR/(k_D + k_fR).

These renormalized rate constants have a clear interpretation. For the effective forward rate it follows, for instance, that: 1/k'_fR = 1/k_D + 1/k_fR; that is, on average, the time required for repressor binding is given by the time needed to diffuse toward the operator plus the time for a reaction to occur when the repressor is in contact with the operator site (51). The effective backward rate has a similar interpretation. The probability that after dissociation the repressor diffuses away from the operator site and never returns is given by , where S_irr(t, r₀) is the irreversible survival probability for two reacting particles (52). Using that , the expression for k'_bR can be written as : that is, the effective dissociation rate is the microscopic dissociation rate multiplied by the probability that after dissociation the repressor escapes from the operator site (51).

For diffusion-limited reactions, such as the reaction considered here, we have that . Now, the renormalized rate constants reduce to:

(16)

(17)

In Fig. 2, we compare the noise profiles for the GFRD algorithm with those obtained by a simulation of the well-stirred model, where instead of the microscopic rates k_fR and k_fB we use the renormalized rates from Eqs. 16 and 17. Surprisingly, we find complete agreement. One of the main reasons why this is unexpected, is that for the master equation the time between events is Poisson-distributed, whereas after a dissociation the time to the next rebinding is distributed according to a power-law distribution when diffusion is taken into account (52). The reason that this power-law behavior of rebinding times is not of influence on the noise profile, is that the timescale of rapid rebinding is much shorter than any of the other relevant timescales in the network. Specifically, rebinding times are so short that the probability that a RNAP will bind before a rebinding is negligible. As a consequence, the transcription network is not at all influenced by the brief period the operator site is accessible before rebinding: for the transcription machinery the series of consecutive rebindings, albeit distributed algebraically in time individually, is perceived as a single event. And on much longer timescales, when a repressor diffuses in from the bulk toward the operator site, the distribution of arrival times is expected to be Poissonian, because on these timescales the repressors are distributed homogeneously in the bulk. This is succinctly summarized in Fig. 4, which shows the distribution of association times. It is seen that at short timescales, the association events are algebraically distributed in time—these arise from the rapid rebindings—whereas at long timescales, they are distributed exponentially in time. For comparison, we also show the distribution of the repressor-DNA association times in the well-stirred model, with appropriately renormalized rate constants for repressor (un)binding (Eqs. 16 and 17). As expected, the number of association events is much smaller at short timescales, but follows the same distribution as that of the spatially resolved model at long timescales. As described quantitatively in the section "Noise propagation", the rate constant for the exponential relaxation is given not only by the diffusion-limited rate of repressor-DNA association, but also by the RNAP promoter occupancy.

View larger version (14K): [in this window] [in a new window]   FIGURE 4  Distribution of repressor-DNA association times for the spatially resolved model (GFRD, solid black line) and for the well-stirred, zero-dimensional model in which the rate constants for repressor-DNA (un)binding are given by the intrinsic (un)binding rates divided by the number of repressor rebindings in the absence of RNAP (see Eqs. 18 and 19)—the renormalized well-stirred model (RWS, dot-dashed line). It is seen that for the spatially resolved model, at short times the distribution follows a power law, while at long times it is exponential. Moreover, at long times, the distribution of the spatially resolved model (GFRD) follows that of the renormalized well-stirred model (RWS). To a good approximation, the relaxation rate at long times in both the spatially resolved model and in the renormalized well-stirred model is given by k'fR[RT][O]', where k'fR is the effective repressor association rate (see Eq. 18), [RT] is the total repressor concentration, and [O]' is the probability that the promoter is not occupied by RNAP, given that it is not occupied by repressor (see Eq. 28). The dip in the distribution at t 0.1–1.0 s–1 is due to the competition with the RNAP for binding to the promoter. The model parameters are f = 100, kOC = 30 s–1, NR = 5.

(18)

(19)

In words, after an initial binding the repressor spends N_RB times longer on the DNA than expected on the basis of the microscopic backward rate, as it rebinds on average N_RB times. Because the average occupancy should not change, the forward rate should be renormalized in the same way. In conclusion, in this model the effects of diffusion can be properly described by a well-stirred model when the reaction rates are renormalized by the average number of rebindings.

	POWER SPECTRA

TOP ABSTRACT INTRODUCTION MODEL SIMULATION TECHNIQUE SIMULATION RESULTS: DYNAMICS AND... SIMULATIONS RESULTS: OPERATOR... TWO-STEP KINETIC SCHEME POWER SPECTRA DISCUSSION AND OUTLOOK APPENDIX: COMPUTING POWER... ACKNOWLEDGEMENTS REFERENCES

 
In this section, we study how the noise due to the stochastic dynamics of the repressor molecules propagates through the different steps of gene expression for both the spatially resolved model and the well-stirred model. This analysis will also provide further insight into why the well-stirred model with renormalized rate constants for the (un)binding of the repressor molecules works so well.

In biochemical networks, the noise in the output signal depends upon the noise in the biochemical reactions that constitute the network, the so-called intrinsic noise, and on the noise in the input signal, called extrinsic noise (6,14,53–56). In our case, the output signal is the protein concentration, whereas the input signal is provided by the repressor concentration. The intrinsic noise arises from the biochemical reactions that constitute the transcription and translation steps. Moreover, we consider the noise in the protein concentration that is due to the (un)binding of the RNAP to (from) the DNA to be part of the intrinsic noise. The extrinsic noise is provided by the fluctuations in the binding of the repressor to the operator, i.e., in the state OR. Because the total repressor concentration, [R_T] = [R] + [OR], is constant, the extrinsic noise is also given by the fluctuations in the concentration of unbound repressor.

The noise properties of biochemical networks are most clearly elucidated via the power spectra of the time traces of the copy numbers of the components. Recently, we have shown that if the fluctuations in the input signal are uncorrelated with the noise in the biochemical reactions that constitute the processing network, the power spectrum of the output signal is given (56)

(20)

Here, S_P() is the power spectrum of the output signal, the protein concentration. The spectrum S_in() denotes the intrinsic noise of the processing network; it is defined as the noise in the output signal in the absence of noise in the input signal. Here, the intrinsic noise is due to the biochemical reactions of transcription and translation. The spectrum S_ex() is the power spectrum of the input signal, which, in this case, is given by the noise in the concentration of unbound repressor: S_ex() = S_R(); because the total repressor concentration is constant, this power spectrum is also directly related to that of the repressor-bound state of the operator, S_OR(). The function g() is a transfer function, which indicates how fluctuations in the input signal are transmitted toward the output signal. If the extrinsic noise is uncorrelated with the intrinsic noise, then g() is an intrinsic quantity that only depends upon properties of the processing network, and not upon properties of the incoming signal (56). However, for the network studied here, the noise in the input signal is not uncorrelated with the intrinsic noise (56). As we have shown recently, this means that Eq. 20 is not strictly valid (56); the extrinsic contribution to the power spectrum of the output signal can no longer be factorized into a function that only depends upon intrinsic properties of the network, g(), and one that only depends upon the input signal, S_ex(). This relation is nevertheless highly instructive. Indeed, Eq. 20 could be interpreted as a heuristic definition of the transfer function g().

The diffusive motion of the repressor molecules impede an analytical evaluation of the power spectrum for the extrinsic noise. Moreover, whereas power spectra can be calculated analytically for linear reaction networks (57), the delays in transcription resulting from promoter clearance and elongation, preclude the derivation of an analytical expression for the power spectrum of the intrinsic noise. We have, therefore, obtained the power spectra S_P(), S_ex(), and S_in(), directly from the time traces of the copy numbers. The power spectrum of a component X is given by , where is the Fourier transform of the concentration X(t) of component X. Conventional FFT algorithms are not convenient for computing the power spectra, because our signals vary over a wide range of timescales. We therefore adopted a novel and efficient procedure, which is described in the Appendix. This procedure should prove useful for computing the power spectra of time traces of copy numbers of species in biochemical networks, as obtained by kinetic Monte Carlo simulations.

As indicated above, the intrinsic noise, S_in(), is defined as the noise in the output signal in the absence of fluctuations in the input signal. To determine the intrinsic contribution to the noise in the protein concentration, we discarded the (un)binding reaction of the repressor to the DNA (Eq. 3), while rescaling the rate k_bRp for the dissociation reaction of the RNAP from the DNA (Eq. 4) in such a way that the average concentration of the protein P remains unchanged. This eliminates the extrinsic noise arising from the repressor dynamics, thereby allowing us to obtain the intrinsic noise of the reactions in Eqs. 4–11. The rescaled backward rate is given by

(21)

For the interpretation of the power spectra of the mRNA and protein concentration, as discussed below, it is instructive to recall the power spectrum of a linear birth-and-death process,

(22)

with rate constants k and µ. For the interpretation of the spectra of repressor binding to the DNA, it is useful to recall the spectrum of a two-state model,

(23)

with rate constants k₁ and k₂. For both models, the power spectrum is a Lorentzian function of the form (58):

(24)

For the birth-and-death process, the variance in the concentration of A, ², is k/µ, whereas for the two-state system, the variance ² in the occupancy n is n(1 – n); the decay rate in the two-state model is µ = k₁ + k₂. The corner frequency µ (in both models) yields the timescale on which fluctuations relax back to steady state. We also note that the noise strength ² is given by the integral of the power spectrum S(): . The noise strength is thus dominated by those frequencies at which the power spectrum is largest.

In the next subsection, we discuss the effect of spatial fluctuations on the noise in gene expression and explain why a well-stirred model with renormalized rate constants for repressor (un)binding can capture its effect. In the subsequent section, we discuss how the noise is propagated through the different stages of gene expression.

Spatial fluctuations
In Fig. 5, we show the power spectra for the input and output signals, for both the spatially resolved model and the well-stirred model with renormalized rate constants for repressor (un)binding (see previous section). We recall that the output signal is the protein concentration, whereas the input signal is the concentration of unbound repressor (the extrinsic noise). Fig. 5 also shows the power spectrum of the intrinsic noise. This is the noise in the protein concentration (the output signal), when the noise in the input signal resulting from the repressor dynamics has been eliminated by the procedure outlined above. The power spectra have been obtained in a parameter regime where the diffusing repressors have a large effect on the noise: k_OC = 30 s^–1, N_R = 5 (see Fig. 2).

View larger version (21K): [in this window] [in a new window]   FIGURE 5  Power spectra of the repressor and protein concentrations for f = 100, kOC = 30 s–1, NR = 5. Data are shown both for the renormalized well-stirred model (RWS) with reaction rates renormalized according to Eqs. 16–17, and for GFRD, taking into account the spatial fluctuations of the repressor molecules explicitly. Also shown is the power spectrum of the intrinsic noise, which is the power spectrum of the protein concentration in the absence of fluctuations in the repressor concentration (extrinsic noise). For large , the repressor spectrum (extrinsic noise) differs between the well-stirred and the spatially resolved model. However, this difference does not appear in the power spectra of the protein concentration. The inset shows the frequency-dependent gain g() (see Eq. 20).