Title: The Effect of Transcription and Translation Initiation Frequencies on the Stochastic Fluctuations in Prokaryotic Gene Expression
Abstract: The kinetics of prokaryotic gene expression has been modelled by the Monte Carlo computer simulation algorithm of Gillespie, which allowed the study of random fluctuations in the number of protein molecules during gene expression. The model, when applied to the simulation of LacZ gene expression, is in good agreement with experimental data. The influence of the frequencies of transcription and translation initiation on random fluctuations in gene expression has been studied in a number of simulations in which promoter and ribosome binding site effectiveness has been changed in the range of values reported for various prokaryotic genes. We show that the genes expressed from strong promoters produce the protein evenly, with a rate that does not vary significantly among cells. The genes with very weak promoters express the protein in "bursts" occurring at random time intervals. Therefore, if the low level of gene expression results from the low frequency of transcription initiation, huge fluctuations arise. In contrast, the protein can be produced with a low and uniform rate if the gene has a strong promoter and a slow rate of ribosome binding (a weak ribosome binding site). The implications of these findings for the expression of regulatory proteins are discussed. The kinetics of prokaryotic gene expression has been modelled by the Monte Carlo computer simulation algorithm of Gillespie, which allowed the study of random fluctuations in the number of protein molecules during gene expression. The model, when applied to the simulation of LacZ gene expression, is in good agreement with experimental data. The influence of the frequencies of transcription and translation initiation on random fluctuations in gene expression has been studied in a number of simulations in which promoter and ribosome binding site effectiveness has been changed in the range of values reported for various prokaryotic genes. We show that the genes expressed from strong promoters produce the protein evenly, with a rate that does not vary significantly among cells. The genes with very weak promoters express the protein in "bursts" occurring at random time intervals. Therefore, if the low level of gene expression results from the low frequency of transcription initiation, huge fluctuations arise. In contrast, the protein can be produced with a low and uniform rate if the gene has a strong promoter and a slow rate of ribosome binding (a weak ribosome binding site). The implications of these findings for the expression of regulatory proteins are discussed. RNA polymerase ribosome binding site transcription initiation frequency Stochastic fluctuations in gene expression have drawn the attention of several research groups. It has been postulated (1McAdams H.H. Arkin A. Proc. Natl. Acad. Sci. U. S. A. 1997; 94: 814-819Crossref PubMed Scopus (1464) Google Scholar, 2Arkin A. Ross J. McAdams H.H. Genetics. 1998; 149: 1633-1648PubMed Google Scholar) and observed (3Newlands S. Levitt L.K. Robinson C.S. Karpf C.A.B. Hodgson V.R.M. Wade R.P. Hardeman E.C. Genes Dev. 1998; 12: 2748-2758Crossref PubMed Scopus (150) Google Scholar) that proteins are produced in short "bursts" occurring at random time intervals rather than in a continuos manner. This expression pattern arises from the fact that elementary processes such as polymerase binding and open complex formation involve small number of molecules and therefore show a wide distribution of reaction times. The fluctuations in a number of protein molecules resulting from stochastic gene expression are of special importance if the concentration of the protein constitutes the message that regulates expression of another gene. In this case, one may expect that the activation pattern of the regulated gene is also random, to some extent. The stochastic nature of gene expression and other biochemical processes involving small numbers of molecules may explain the variations observed in isogenic populations of bacterial cells. Examples are individual chemotactic responses (4Levin M.D. Morton-Firth C.J. Abouhamad W.N. Bourret R.B. Bray D. Biophys. J. 1998; 74: 175-181Abstract Full Text Full Text PDF PubMed Scopus (82) Google Scholar), the distribution of generation times observed in Escherichia coli cells (5Tyson J.J. Hannsgen K.B. J. Theor. Biol. 1985; 113: 29-62Crossref PubMed Scopus (37) Google Scholar), and the individual cell responses observed during induction of lactose (6Maloney P.C. Rotman B. J. Mol. Biol. 1973; 73: 77-91Crossref PubMed Scopus (75) Google Scholar) and arabinose (7Siegele D.A. Hu J.C. Proc. Natl. Acad. Sci. U. S. A. 1997; 94: 8168-8172Crossref PubMed Scopus (348) Google Scholar) operons at subsaturating inducer concentrations. Stochastic effects have also been observed in engineered genetic networks in E. coli (8Elowitz M.B. Leibler S. Nature. 2000; 403: 335-338Crossref PubMed Scopus (3375) Google Scholar). Taking into account the stochastic phenomena described above, it is clear that the cell must contain various mechanisms that ensure deterministic regulation of cellular processes. The check points in eukaryotic cell cycle regulation (9Alberts B. Bray D. Lewis J. Raft M. Roberts K. Watson J.D. Molecular Biology of the Cell. 3rd Ed. Garland Publishing, London and New York1994Google Scholar) are examples of such a mechanism. Computer simulations implementing various Monte Carlo algorithms proved to be useful in studying stochastic processes in biochemical reaction networks. This is especially true for several prokaryotic systems for which, due to their relative simplicity, a large number of quantitative measurements are available. Examples are the kinetic analysis of developmental pathway bifurcation in phage λ (2Arkin A. Ross J. McAdams H.H. Genetics. 1998; 149: 1633-1648PubMed Google Scholar), computer simulation of the chemotactic signaling pathway (4Levin M.D. Morton-Firth C.J. Abouhamad W.N. Bourret R.B. Bray D. Biophys. J. 1998; 74: 175-181Abstract Full Text Full Text PDF PubMed Scopus (82) Google Scholar), and mechanistic modeling of the all or none effect in lactose operon regulation (10Carrier T.A. Keasling J.D. J. Theor. Biol. 1999; 201: 25-36Crossref PubMed Scopus (43) Google Scholar). The goal of this article is to study the kinetic mechanisms that are responsible for the stochastic phenomena of gene expression. Because both RNA and protein elongation rates are approximately equal for all genes, different levels of gene expression are the result of varying frequencies of transcription and translation initiations. Therefore, we examined the influence of transcription and translation initiation frequencies on the pattern of gene expression. The model of a prokaryotic gene was built using kinetic parameters of LacZ determined by Kennell and Riezman (11Kennell D. Riezman H. J. Mol. Biol. 1977; 114: 1-21Crossref PubMed Scopus (238) Google Scholar). The model is consistent with the experimental results presented in their work. We have systematically varied the frequencies of transcription and translation initiation reactions, keeping other parameters of the model constant. Our results indicate that the genes expressed from strong promoters (high frequency of transcription initiation) produce proteins with constant velocity and low variation in the number of molecules. A decrease in the frequency of transcription initiation causes a noisy pattern of gene expression, i.e. the protein is produced in bursts at random time intervals. In contrast, a decrease in the frequency of translation initiation lowers the speed of protein synthesis but does not lead to noisy expression patterns. We present the kinetic model of single prokaryotic gene expression. The transcription and translation processes are represented by the collection of kinetic equations that are numerically integrated with the use of the Monte Carlo approach of Gillespie (12Gillespie D.T. J. Phys. Chem. 1977; 81: 2340-2361Crossref Scopus (7224) Google Scholar, 13Gillespie D.T. Phys. A. 1992; 188: 404-425Crossref Scopus (774) Google Scholar). We assume that the cells in which the model gene is expressed are in the exponential growth phase. All of the experimental data cited in our work concern the cells in this growth phase. Moreover, under these conditions, the concentrations of free RNA polymerase and ribosomes available for the single gene are most probably kept constant by the delicate balance between many coupled reactions (14Bremer H. Dennis P.P. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 1553-1569Google Scholar, 15McClure W.R Annu. Rev. Biochem... 1985; 54: 171-204Google Scholar, 16Neidhardt F.C. J. Bacteriol. 1999; 181: 7405-7408Crossref PubMed Google Scholar) such as synthesis of RNA polymerase and ribosomes, binding of RNA polymerase and ribosomes to other expressing genes, nonspecific binding of RNA polymerase to DNA, and its idling at pause sites. A detailed model of all the processes that buffer concentrations of free RNA polymerase and ribosomes would require taking into account the biosynthesis of all the components of gene expression machinery and their binding within all the transcription units in the cell. Most of the parameters of such a complex model would be unknown and impossible to estimate. Instead, we included in the model randomly changing pools of the free RNA polymerase and ribosomes. Because these pools result from many small contributions of other processes, their distributions should be Gaussian. The mean values of these distributions have been set close to the estimations presented by other authors. The fluctuations in the RNA polymerase and ribosome concentrations described by the standard deviations of the distributions are unknown. Therefore, we checked the sensitivity of our results to these parameters. Although many research groups are developing detailed kinetic models of various elements of the gene expression machinery (e.g. transcription initiation (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar)), the in vivo values of the rate constants of many elementary reactions remain unknown. In contrast, many experimental works describe the kinetics of prokaryotic gene expression (mainly of the LacZ gene in E. coli) in terms of simple effective parameters such as the frequencies of transcription and translation initiation, the rate of RNA polymerase and ribosome movement, mRNA half-life, etc. (11Kennell D. Riezman H. J. Mol. Biol. 1977; 114: 1-21Crossref PubMed Scopus (238) Google Scholar,18Yarchuk O. Jacques N. Guillerez J. Dreyfus M. J. Mol. Biol. 1992; 226: 581-596Crossref PubMed Scopus (147) Google Scholar). In our model, we adjust unknown values of elementary reactions to reproduce measured values of transcription and translation inititiation frequencies, mRNA level, and the protein synthesis rate in vivo. We assume that effective parameters obtained in this way implicitly include reactions involving nucleotides, translation initiation factors, bivalent ions, and other compounds necessary for proper function of the gene expression machinery. The summary of the kinetic model and its parameters for the LacZ gene is given in Table I. Details of the model and simulation algorithm are described below.Table ISummary of the kinetic model and its parameters for the LacZ geneReaction/processRate constantsCommentsThe pool of available RNA polymerase The pool of available ribosomesGaussian distribution with μ = 35 and δ = 3.5 molecules Gaussian distribution with μ = 350 and δ = 35 moleculesOrder of magnitude of the mean is set according to experimental estimations (14Bremer H. Dennis P.P. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 1553-1569Google Scholar, 17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). Sensitivity of the results to the value of δ has been tested.RNA polymerase binding and dissociationAssociation rate: 108m−1 s−1 Dissociation rate: 10 s−1Association rate is set according to experimental estimations (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). Dissociation rate is set to reproduce transcription initiation frequency measured by Kennell & Riezman (11Kennell D. Riezman H. J. Mol. Biol. 1977; 114: 1-21Crossref PubMed Scopus (238) Google Scholar).Closed complex isomerization1 s−1According to experimental estimations (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar)Promoter clearance and RBS synthesis1 s−1Estimated using the rate of polymerase movement and the lengths of RBS and promoter sequencesRibosome binding and dissociationAssociation rate: 108m−1s−1 Dissociation rate: 2.25 s−1Association rate is set to the order of magnitude of diffusion limited aggregation. Dissociation rate is set to reproduce proper translation initiation frequency.RBS clearance0.5 s−1According to experimental estimations (19Draper D.E. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 902-908Google Scholar)mRNA degradation0.3 s−1Set close to transcription initiation frequency. Values at least as high as transcription initiation frequency are necessary to assure that the mRNA level is lowered by the decrease in protection of RBS by bound ribosomes (as observed by Yarchuk et al. (18Yarchuk O. Jacques N. Guillerez J. Dreyfus M. J. Mol. Biol. 1992; 226: 581-596Crossref PubMed Scopus (147) Google Scholar)).Elongation of protein chain0.015 s−1According to commonly accepted value of ribosome movement rate (15 amino acids/s) and the length of β-galactosidase.Protein degradation6.42 × 10−5According to β-galactosidase half-life of 3 h measured by Bergquist and Truman (22Berquist P.L. Truman P. Mol. Gen. Genet. 1978; 164: 105-108Crossref PubMed Scopus (8) Google Scholar). Open table in a new tab Transcription initiation is represented in our model as a two-step process: P+RNAP→P_RNAPEquation 1 P_RNAP→P+RNAPEquation 2 P_RNAP→TrRNAPEquation 3 where P represents the promoter region of the gene,RNAP represents RNA polymerase, and TrRNAPrepresents the transcribing RNA polymerase. Reactions in Equations 1and 2 describe reversible RNA polymerase binding; the reaction in Equation 3 represents isomerization of closed binary complex to open binary complex. The second order rate constant of RNAP1 binding has been set to 108m−1s−1, which is the order of magnitude of RNAP DNA binding determined in in vitro assays (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). The values of closed complex to open complex isomerization rates reported in the literature are of the order of 1 s−1. The rate of reaction (Equation 3) was set to this value. If the rate constants of reactions in Equations 1 and 3 are set to the values listed above, the rate of RNAP dissociation must be of the order of 10 s−1 to reproduce a transcription initiation frequency close to the value of 0.3 s−1reported by Kennell and Riezman (11Kennell D. Riezman H. J. Mol. Biol. 1977; 114: 1-21Crossref PubMed Scopus (238) Google Scholar). Moreover, the binding constant resulting from the rates of reactions in Equations 1 and 2is 107m−1, which is in reasonable agreement with values estimated according to in vitro measurements (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). To clear the promoter region, active RNA polymerase must move 30 to 60 nucleotides (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). Taking into account that the rate of polymerase movement is about 40 nucleotides/s, this step takes about 1 s (17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). The length of the mRNA chain that is synthesized during this time corresponds roughly to the length of the leader region containing the ribosome binding site (RBS). Therefore, the synthesis of the RBS and promoter clearance occur at approximately the same rate of 1 s−1. It is irrelevant for the conclusions of our model whether the RBS appears slightly earlier or later than promoter clearance. Therefore, we modeled these two processes by the single first order reaction with a rate constant of 1 s−1: TrRNAP→RBS+P+ElRNAPEquation 4 where ElRNAP denotes polymerase elongating a given mRNA molecule. The elongation of the mRNA is not modeled. We assume that the RNA polymerase completes synthesis of the mRNA molecule with a rate sufficient to allow for a ribosome movement rate of 15 amino acids/s. In prokaryote transcription, translation and mRNA degradation are tightly coupled. Following the experimental results of Yarchuk et al. (18Yarchuk O. Jacques N. Guillerez J. Dreyfus M. J. Mol. Biol. 1992; 226: 581-596Crossref PubMed Scopus (147) Google Scholar), we have assumed the following interdependence between translation and mRNA decay: (i) RNase E and ribosomes compete for the RBS; (ii) if RNase E binds to the RBS faster than the ribosome, it degrades mRNA in the 5′ to 3′ direction and does not interfere with the movement of the ribosomes that have been already bound; and (iii) every ribosome that successfully binds to the RBS completes translation of the protein. The above assumptions have been modeled by the following reactions: Ribosome+RBS→RibRBSEquation 5 RibRBS→RBS+RibosomeEquation 6 RibRBS→ElRib+RBSEquation 7 RBS→decayEquation 8 where Ribosome represents the pool of free ribosomes,RibRBS represents the ribosome binding site protected by the bound ribosome, and ElRib denotes the ribosome elongating protein chain. The reaction in Equation 8 models decay of unprotected RBS. Parameters of the reactions in Equations Equation 5, Equation 6, Equation 7, Equation 8 are more difficult to estimate than those concerning transcription initiation. The second order rate constant for ribosome binding was set to 108m−1s−1, which is of the order of diffusion limited macromolecule binding. The rate of the RBS clearance step in translation initiation (the reaction in Equation 7) was set to 0.5 s−1, which is close to experimental estimates (19Draper D.E. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 902-908Google Scholar). Parameters of the reactions in Equations 6 and 8 have been adjusted to reproduce the experimentally known rate of protein synthesis and ribosome spacing of the LacZ gene. Another important constraint is provided by experiments of Yarchuk et al. (18Yarchuk O. Jacques N. Guillerez J. Dreyfus M. J. Mol. Biol. 1992; 226: 581-596Crossref PubMed Scopus (147) Google Scholar). The authors introduced mutations in the RBS region of the LacZ gene that resulted in a 200-fold variation in the expression level. A decrease in the protein synthesis rate has been followed by a nearly equivalent change in the stationary mRNA level. This means that for the less effective ribosome binding sites, the rate of mRNA decay exceeded or at least was no lower than the transcription initiation frequency. Otherwise, the significant mRNA level would be observed even in the case of RBS practically unprotected by ribosomes. We have set the rate of the reaction in Equation 8 to 0.3 s−1, which is equal to the transcription initiation frequency. Therefore, the rate of ribosome dissociation has to be set to 2.25 s−1 to reproduce the rate of protein synthesis measured by Kenell and Riezman. The elongation of the protein chain has been modeled as a single first order reaction: RibRBS→proteinEquation 9 We have used the commonly accepted rate of ribosome movement of 15 amino acids/s. Therefore, the rate of the reaction in Equation 9 has been set to 15/Ls−1, where L denotes the length of protein chain. For the 1024-amino acid-long product of theLacZ gene, the value of this parameter is 0.015. This treatment of the elongation reaction is simplified with respect to other Monte Carlo models of prokaryotic gene expression (1McAdams H.H. Arkin A. Proc. Natl. Acad. Sci. U. S. A. 1997; 94: 814-819Crossref PubMed Scopus (1464) Google Scholar, 20Carrier A.T. Keasling J.D. J. Theor. Biol. 1997; 189: 195-209Crossref PubMed Scopus (26) Google Scholar). In these publications, the addition of a single amino acid by a single ribosome has been simulated as a separate reaction, and ribosome overlap has been prevented by excluded volume constraints. It will be shown in the following sections that despite neglecting fine details of ribosome movement and using the average rate of protein chain synthesis instead, the model is able to reproduce the experimental data concerning expression of the LacZ gene. Simplification of the elongation reaction made the simulation computationally fast, thus we could collect good statistics from long timescale trajectories for a large number of parameter combinations. Protein degradation has been simulated by the following reaction: Protein→decayEquation 10 Although protein degradation rates are variable, most of the native proteins have half-lives of the order of hours (21Goldberg A.L. Dice J.F. Annu. Rev. Biochem. 1974; 43: 835-869Crossref PubMed Scopus (519) Google Scholar). In our simulations, we have applied the value of 3 h that has been measured for wild type β-galactosidase in E. coli cells (22Berquist P.L. Truman P. Mol. Gen. Genet. 1978; 164: 105-108Crossref PubMed Scopus (8) Google Scholar). The numbers of available RNAP and ribosome molecules were generated from the Gaussian distribution. In the case of RNAP, the mean of the distribution was set to 35 molecules. This value is close to estimations made by other authors according to experimental data (14Bremer H. Dennis P.P. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 1553-1569Google Scholar, 15McClure W.R Annu. Rev. Biochem... 1985; 54: 171-204Google Scholar, 17Record T.M. Reznikoff W.S. Craig M.L. McQuade K.L. Schlax P.J. Escherichia coli and Salmonella. 2nd Ed. ASM Press, Washington, D. C.1996: 792-821Google Scholar). Because the number of ribosome molecules present in the cell is approximately an order of magnitude larger than the number of RNAP molecules, the mean of the ribosome pool distribution has been set to 350. We are not able to estimate the S.D. of the above-mentioned pools. We have repeated all the experiments with the standard deviations equal to 10% and 50% of the means to study the influence of the fluctuations in RNAP and ribosome pools on the gene expression pattern. In the following sections, results will be presented for the standard deviations of the pools equal to 10% of the means, unless it is stated explicitly that standard deviations are set to 50% of the means. To compute the time evolution of the system described above, we have used the Monte Carlo algorithm of Gillespie (12Gillespie D.T. J. Phys. Chem. 1977; 81: 2340-2361Crossref Scopus (7224) Google Scholar). This method has a sound physical basis and has already been proven useful in the simulation of prokaryotic gene expression (1McAdams H.H. Arkin A. Proc. Natl. Acad. Sci. U. S. A. 1997; 94: 814-819Crossref PubMed Scopus (1464) Google Scholar,2Arkin A. Ross J. McAdams H.H. Genetics. 1998; 149: 1633-1648PubMed Google Scholar). For the convenience of the reader, a brief description of the Gillespie algorithm is given below. The rigorous derivation of this computational protocol can be found elsewhere (12Gillespie D.T. J. Phys. Chem. 1977; 81: 2340-2361Crossref Scopus (7224) Google Scholar, 13Gillespie D.T. Phys. A. 1992; 188: 404-425Crossref Scopus (774) Google Scholar). Gillespie's algorithm is formulated to numerically study systems composed of a large number of chemical species whose interactions are described by rate equations. At each step of simulation, two random numbers are generated. The first one is used to determine which reaction in the system will happen; the second one determines the waiting time for this event. The probability of a given reaction to be chosen is proportional to the product of its stochastic rate constant and numbers of substrate molecules. The distribution of the waiting times is given by the following equation: P(τ,μ)=aμexp(−a0τ)Equation 11 where aμ=hμcμ;a0=∑aμEquation 12 P(τ, μ) is the probability that the waiting time for the reaction is τ, provided that reaction μ was chosen to happen; h μ is the product of the numbers of substrate molecules in reaction μ; and c μ is the stochastic rate constant of this reaction. After the reaction and its waiting time are determined, one elementary reaction is executed (e.g. if the reaction is A + B→AB, one A molecule and one B molecule will be removed from the system, and one AB molecule will be added). The time of the simulation is then increased by τ, and the next simulation step is executed. The stochastic rate constant of the reaction is defined as the probability that the elementary reaction will happen in the infinitesimally short time interval. It can be calculated from the rate constant using simple relations. For the first order reactions, both constants are equal. For the second order reaction, the stochastic rate constant is equal to the rate constant divided by the volume of the system. It has been rigorously proven that for an infinitely large number of molecules, the Markov Chain generated by the Monte Carlo scheme given above converges to the same trajectory of the system as the deterministic integration of differential equations (12Gillespie D.T. J. Phys. Chem. 1977; 81: 2340-2361Crossref Scopus (7224) Google Scholar, 13Gillespie D.T. Phys. A. 1992; 188: 404-425Crossref Scopus (774) Google Scholar). The advantage of the Gillespie algorithm over the differential equation approach is that it also remains exact for arbitrary low numbers of molecules in the system. This is extremely important in the simulation of biochemical systems, where some processes involve single molecules (e.g. one molecule of promoter in transcription initiation). Randomly fluctuating pools of available ribosomes and RNA polymerase have been modeled in the following way. Before executing every step of the Gillespie algorithm, the number of RNAP molecules and ribosomes has been generated from Gaussian distributions with the means and standard deviations defined in the previous section. Then, the values ofa μ (see Equation 12) for every reaction were calculated, and the step of the Gillespie algorithm has been executed. Generation of the number of substrate molecules from the appropriate probability distribution has already been applied by Arkinet al. (2Arkin A. Ross J. McAdams H.H. Genetics. 1998; 149: 1633-1648PubMed Google Scholar) in the Gillespie algorithm simulation to model rapid equilibrium phenomena in regulatory protein binding. We have applied the Gillespie algorithm to the model described in the previous sections. At the beginning of the simulation, only a single promoter (P) was present in the system. The numbers of all other molecules, except random pools of RNAP and ribosomes, were set to 0. Simulation was pursued until it reached the generation time ofE. coli, and then cell division was simulated,i.e. simulation was continued with the system containing one promoter molecule and all other numbers of molecules divided by 2. The generation time was set to 35 min. We have neglected the fact that during DNA replication, two copies of the bacterial gene can be expressed. Linear volume change of the growing cell has been assumed. Before every step of the Gillespie algorithm, stochastic rate constants of second order reactions have been recomputed by dividing the appropriate rate constants by the current volume of the cell. Simulation of several generations is important in the simulation of constitutive gene expression. As shown below, after starting our simulation with 0 protein molecules, the level of gene expression equilibrates after a few generations. We have found that 10 generations are sufficient to equilibrate the system. For each simulation, we computed 100 independent trajectories of the system covering 10 generations each. We have used our own implementation of the Gillespie algorithm. The software has been tested against examples given in the original work (12Gillespie D.T. J. Phys. Chem. 1977; 81: 2340-2361Crossref Scopus (7224) Google Scholar). To check the validity of the model, we performed the simulation with the parameters describing theLacZ gene. Fig. 1 Apresents the gene expression pattern spanning 10 bacterial generations. Standard deviation of RNAP and ribosome pools has been set to 3.5 and 35, respectively (10% of the mean values). To compare these data with the results of Kennell and Riezman (11Kennell D. Riezman H. J. Mol. Biol. 1977; 114: 1-21Crossref PubMed Scopus (238) Google Scholar), we have analyzed the changes in the number of protein molecules in the first generation because the experimental data covered the time scale of about 30 min after induction. The transcription initiation frequency obtained in this simulation was equal to 0.258 s−1, close to that of 0.3 s−1 reported by Kennell and Riezman for the LacZ gene. As seen in Fig.1 B, after about 500 s, the number of protein molecules begins