- Email: [email protected]
Evolving kinetics of gene expression in stochastic environments
Computational Biology and Chemistry 37 (2012) 11–16
Contents lists available at SciVerse ScienceDirect
Computational Biology and Chemistry journal homepage: www.elsevier.com/locate/compbiolchem
Evolving kinetics of gene expression in stochastic environments Antti Häkkinen, Andre S. Ribeiro ∗ Laboratory of Biosystem Dynamics, Computational Systems Biology Research Group, Department of Signal Processing, Tampere University of Technology, Finland
a r t i c l e
i n f o
Article history: Received 4 November 2011 Received in revised form 18 January 2012 Accepted 14 February 2012 Keywords: Gene expression Stochastic dynamics Genetic operators Fluctuating environments Transcriptional bursting
a b s t r a c t Recent studies have shown that the in vivo dynamics of RNA numbers in bacteria is regulated, to a great extent, by the kinetics of rate limiting steps in transcription. Strong evidence suggests that the kinetics of these steps is sequence dependent. We investigate the selective advantages of rate limiting steps of differing kinetics. For that, we model the kinetics of expression of a gene responsible for promoting cell division at the expense of resources in the environment in individual cells of a population. We model mutations that affect the kinetics of the rate limiting steps and selective pressure in various environmental conditions. Depletion of resources leads to cell death. We find that small changes in the evolutionary constraints can favor widely different noise levels in RNA and protein numbers. Increasing the cost in nutrients for division favors noisier expression. The results provide a better understanding of why different genes differ in the kinetics of production of RNA and proteins. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Different genes have evolved to express at different rates and with different dynamical patterns (McClure, 1985). Genome wide observations of protein numbers in individual Escherichia coli cells have shown that the cell to cell diversity varies widely from one protein to the next and can range from sub- to super-Poissonian (Taniguchi et al., 2010). This diversity is likely present at the RNA level as well, since protein numbers tend to follow RNA numbers in bacteria (Kaern et al., 2005). Evidence suggests that most RNA molecules exist in low copy numbers in organisms such as E. coli, since most genes are expressed one to a few times during the cell lifetime (Bernstein et al., 2002). Due to this, and because protein numbers follow RNA numbers, small fluctuations in RNA numbers are likely to have downstream effects on the phenotype (Choi et al., 2010). Organisms may therefore have evolved mechanisms to stringently control the stochasticity in gene expression. In general, essential genes appear to have expression levels that are less noisy than that of genes associated with, for example, response to perturbations (Fraser et al., 2004; Taniguchi et al., 2010). The expression rate of a gene is determined by several factors and varies for different environmental conditions. On one hand,
∗ Corresponding author at: Department of Signal Processing, Tampere University of Technology, P.O. Box 553, 33101 Tampere, Finland. Tel.: +358 3 3115 4975; fax: +358 3 3115 4989. E-mail addresses: antti.hakkinen@tut.fi (A. Häkkinen), andre.ribeiro@tut.fi (A.S. Ribeiro). 1476-9271/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiolchem.2012.02.002
genes need to express at a rate fast enough to produce the minimum number of proteins required to execute a task. On the other hand, the rate of production ought to be limited to avoid waste of energy by the cell. In addition, expression rates are likely to be limited by other constraints such as toxicity beyond certain concentrations (Eckert and Beck, 1989). There is evidence suggesting that one core regulatory process of fluctuations in protein and RNA temporal levels is the multi-step process of transcription initiation (Ribeiro et al., 2010; Kandhavelu et al., 2011). This process has, in general, one to a few rate limiting steps, that is, steps whose duration determine the rate of transcription. Therefore, these steps have a strong effect on the strength of fluctuations in RNA and, consequently, protein numbers in bacteria. Other events also affecting fluctuations can occur at later stages, such as sequence dependent pauses in transcriptional elongation (Rajala et al., 2010; Ribero et al., 2010). Current techniques do not allow in vivo measurements of the durations of the underlying steps of transcription initiation. However, it is possible to measure the intervals between the completion of consecutive transcription events (Kandhavelu et al., 2011). This, along with recent measurement techniques of protein numbers (Yu et al., 2006), has allowed the development of accurate stochastic models of gene expression (Zhu et al., 2007; Ribero, 2010). The first stochastic models of genetic circuits assumed gene expression as an instantaneous process (McAdams and Arkin, 1998). The dynamics was driven by the Stochastic Simulation Algorithm (SSA, see Gillespie, 1976). Several recent studies (Bratsun et al., 2005; Roussel and Zhu, 2006; Ribeiro et al., 2006) have shown that the time for gene expression to be completed once initiated is not negligible in the dynamics of genetic circuits. This duration,
12
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
while sequence dependent, also has a random nature, varying even between transcription events of the same gene (Kandhavelu et al., 2011) and thus, is also a source of stochasticity. The delayed Stochastic Simulation Algorithm (DSSA, see Roussel and Zhu, 2006) was developed as an extension of the original SSA to allow the use of delayed reaction products besides the elementary reactions with exponential statistics. This is useful to simulate genetic circuits, since some of the underlying processes consist of sequential events, exhibiting Gaussian-like rather than exponential-like statistics. Since the kinetics of these underlying events are unknown (e.g. in the case of transcription initiation), the explicit modeling is not be possible. If one is to model evolutionary processes of genetic circuits, the evolvability of these rate limiting steps in transcription and translation need to be considered. Finally, genes and genetic circuits are not isolated systems. Rather, they have evolved to function within particular ranges of environmental conditions and, to some extent, to cope with stochastic fluctuations in the environment. The stochastic nature of the environments also needs to be accounted for when simulating and studying how genetic circuits evolve. Several works have focused on the evolutionary dynamics of genes or small genetic circuits, when subject to stochastic environments (see e.g. Cohen, 1966; Levins, 1968; Schaffer, 1974). Recently, this topic has been revisited in the context of how initially isogenic populations may adapt to fluctuating environments (Lachmann and Jablonka, 1996; Kussell et al., 2005; Kussell and Leibler, 2005). One of the interesting findings from these studies is that the stochastic switching between two phenotypes may be used by cells as a survival strategy when subject to fluctuating environments (Acar et al., 2008). Another study investigated how genetic switches can better adjust their switching dynamics provided that they possess a sensing mechanism of the environment state, such as the detection of the amount of toxins or nutrients in the environment at a given moment (Ribeiro, 2008). In these models, gene expression delays are either not present, or modeled as constant from one event to the next. Additionally, the cell populations are constrained to fixed sizes. Here we study the dynamics and evolutionary pathways of a delayed stochastic model of gene expression in cell populations of variable size subject to stochastic environments. For that, we model the evolution of the dynamics of a gene responsible for expressing a protein required for cell division in an environment with limited, stochastic supply of nutrients required for protein expression. We allow mutations and selection to act on the kinetics of the rate limiting steps of gene expression, and test whether the environmental pressure drives the evolution of different kinetics. To test this hypothesis, we developed an extension for SGN Sim (Ribeiro and Lloyd-Price, 2007), that allows modeling evolutionary patterns of genetic circuits of populations variable in size and the fluctuating environments they are subject to. The features of this extension are described in the next section.
interpreter. A simulation involves simulating a set of entities, who are independent and whose dynamics are governed by the DSSA, and scheduling events that manipulate these entities to the event queue. The entities can be created, destroyed, and mutated dynamically, and the times at which they occur are drawn from arbitrary distributions, either using the built-in random number generation facilities or ones interfaced with Lua. The simulation of the entities uses DSSA, which is an extension of the SSA (Gillespie, 1976). This extension allows non-Markovian dynamics, by attaching an arbitrarily distributed time delay before the release of each of the reaction products. The DSSA simulation engine is based on that of SGN Sim (Ribeiro and Lloyd-Price, 2007), which implements the algorithm as specified by Roussel and Zhu (2006). The simulation algorithm is a form of the logarithmic direct method (Li and Petzold, 2006), modified for DSSA. The entities can be created based on a previously created entity, or based on a system description. The system descriptions are written in SGN Sim format (Ribeiro and Lloyd-Price, 2007). Moreover, the entities can be modified immediately after their creation, such that schemes like reproduction by division can be implemented. Mutations can be applied by queueing mutating operators in the event queue, with time intervals drawn from arbitrary distributions, affecting the molecular counts or the parameters of the model. The scheme of mutations, such as to which parameters the mutations affect, and to what extent they are mutated are left to be specified by the user. It is also possible model more complex schemes of that introduced genetic variation, for example, by exchanging the parameters values between separate entities. At any given point in time, a subpopulation of the simulated entities can be selected, based on arbitrary combination of features in their state, such as concentration of some molecules. Combined with the operators of reproduction or destruction, this allows modeling the process of natural selection, based on an arbitrary fitness function. Moreover, the events that can be scheduled are not limited to the creation, destruction, and modification of the entities. It is possible to schedule events that are arbitrary Lua functions; however, this has the disadvantage that the parallel simulations must be synchronized. In cases where it is applicable, such as in the case of cell division, this synchronization is not necessary and can be avoided by combining the non-locking built-in operations. The built-in Lua interpreter is used to parse the description of the evolutionary process of the simulator. This description consist of queuing the events to the event queues, and to set up parameters such as the length of the simulation. The events can be queued as one-shot events, possibly being queued again at the occurrence of some event, or as recurring events. The time intervals between these events can be drawn from random number generators of various distributions, or specified arbitrarily. It can also be used to execute arbitrary pieces of Lua code.
2.2. Model 2. Methods 2.1. Simulation implementation To study the evolutionary phenomena in genetic networks, we implement a simulator that combines simulation of stochastic molecular kinetics with evolutionary operators. The simulator is written in C++and can be compiled on Microsoft Windows and typical Unix-like systems, such as Linux. Multithreading is exploited to accelerate simulations on hardware that supports parallel execution. The simulator is available upon request. The simulator consists of a DSSA simulation engine, a parallel event queueing system (priority queue), and a built-in Lua
We study the evolutionary patterns of the kinetics of expression of a gene responsible for the production of a protein necessary for cell division. The expression is triggered by the presence of nutrients in the cell. This is a common process in real cells. For example, evidence suggests that, when starving, the gene expression patterns of the cells change, the nature of change depending on which nutrients are lacking (Peterson et al., 2005). Also, starving cells commonly inhibit DNA replication (Joseleau-Petit et al., 1999). For example, under amino acid starvation, DNA replication is inhibited due to the inhibition of DnaA expression (Wang and Levin, 2009), a protein required in multiple copies to initiate the DNA replication (Mott and Berger, 2007).
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
13
Fig. 1. A pictorial overview of the model. Nutrients are created in the environment and distributed to the cells. The cells require nutrients to express proteins. Accumulation of these proteins is required for cell division. Also depicted are the kinetic rate constants of the several processes in the model, described below.
To mimic these mechanisms and their evolution, in the model proposed here the expression of the protein is triggered by the presence of nutrients in the cell. However, cells cannot deplete all the nutrients in the division process, as this leads to death. The size of the cell population is dynamic and all individuals compete for dividing during their limited lifetime, , in an environment where nutrients are generated by a stochastic process with a constant rate and distributed uniformly to the cells. The simulations of the evolutionary process aim at determining the optimal set of parameters values for both the mean and the strength of the fluctuations in protein levels, as a function of the environmental conditions. Also, we study how the optima change with the environmental settings. A pictorial description of the model is shown in Fig. 1. The model consists of a set of chemical interactions that are implemented in each individual cell. In the reactions, i denotes the cell index, to make a distinction between the variables that are not shared by the individuals. Nutrients ni are transported to a cell by the following reaction equation (RE): kn N(t)−1
∅ −→ ni ,
(1)
where kn is the total rate at which novel nutrients units become available, and N(t) is the cell population size at moment t. Consequently, reaction (1) is such that, the larger is the population, the less nutrients ni will be available for each cell. The expression of proteins uses nutrient units as substrate (RE (2)). Proteins and nutrients are subject to degradation via reactions (3) and (4), respectively: ∞
Pi + ni −→Pi (i (t)) + pi dp
pi −→∅ dn
ni −→∅
(2) (3)
,
(4)
where Pi denotes the promoter region of the gene, and pi denotes the protein. The rate of RE (2) is set to infinity, such that the production of proteins is only limited by the absence of nutrients and by the time delay i .
The time delay i (t) determines the time between the production of two consecutive proteins. This delay accounts for the underlying processes during transcription initiation (McClure, 1985) and recent evidence suggests that these processes are the most relevant in determining the rate of RNA production in bacteria (Kandhavelu et al., 2011). In prokaryotes, evidence suggests that the rate of protein production generally tends to follow that of the RNA production (Kaern et al., 2005). We set the delay i (t) ∼ Gamma(˛i (t), ˇi (t)) to be gamma distributed with shape parameter ˛i (t) and rate parameter ˇi (t). The gamma distribution allows modeling a wide range of dynamics by varying its parameters ˛i (t) and ˇi (t). Depending on the values of the parameters, this distribution allows a sub-Poissonian (˛i (t) > 1) dynamics of protein production, e.g. resultant of a sequential rate-limiting mechanism of gene expression (Kandhavelu et al., 2011), or a super-Poissonian (˛i (t) < 1) dynamics, which results from models such as the on–off model (Peccoud and Ycart, 1995). Mutations to the parameters ˛i (t) and ˇi (t) allow the kinetics of protein production to change. The mutations are performed with exponentially distributed time intervals with a rate of km , in a linear scale with step size of m, with uniform probability to either increase or decrease. Finally, when a threshold T in the number of pi is reached, the cell is able to divide. In this process, two daughter cells are formed, who split equally pi and ni of the mother cell, and inherit the values of ˛i (t) and ˇi (t), which regulate the kinetics of gene expression. 3. Results In optimal growth conditions the division time of E. coli is 20 min. Thus, we set the maximum cell lifetime = 100 min. Given the parameters described below, we observed that for the initial parameter values, the average division time was approximately 40 min. We choose initial shape of ˛i (0) = 1, which corresponds to Poissonian production dynamics of proteins. Initially, the expression and degradation rates are set such that the mean protein level is 5. For that we set the production rate ˇi (0) = 5 min−1 . This value,
14
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
while high, is within realistic intervals, since e.g. the lac promoter can initialize an RNA production every three seconds (Kennell and Riezman, 1977), while the lar promoter produces one RNA per 15 min under full induction (Golding et al., 2005; Kandhavelu et al., 2011). Next, we set the environmental conditions, that is, the amount of available nutrients. We want to compare how the gene expression dynamics evolves in two different conditions, namely, when many nutrients are required for dividing and when very few nutrients are required. However, we also need to have populations of comparable size in the two simulations, since the number of cells determines e.g. how many mutations overall can occur in a generation. To have populations of similar size, we scale the amount of nutrients in the environment by the amount of nutrients needed for division. Thus, we set kn = 20T min−1 and dn = 0.1 min−1 . Finally, regarding the gene expression mechanism, we need to set a mutation rate that determines the speed at which the cell population can adapt to the environment. As mentioned, ˛i and ˇi are subject to changes due to mutations. We found that the speed by which these parameters change does not, within a wide range of values, affect the final outcome of the selection process. Due to that, for simplicity, we set the mutation rate to km = 1 min−1 and step size m = 0.02, arbitrarily, as these values permit reaching near optimal values in a reasonable number of generations. In the results presented, the initial population size is set to N(0) = 10. We investigated several values of initial population sizes, namely N(0) ∈ {1, 2, 5, 10, 20, 50, 100, 200, 500}, and found no significant differences in their long term behavior. This is because, in the beginning of each simulation, the population size is first reduced to a small value due to starvation, after which the population size starts growing toward a limiting value as the cells become better adapted to the environment. In the first generation, the promoter is initially set to be unoccupied ([Pi ] ← 0), and the protein concentration is set to reflect the expected protein number given infinite source of food ([pi ] ← ˛i (0)−1 ˇi (0) dp−1 ). Moreover, the initial number of nutrients [ni ] in the cell is set to zero. We study the evolved noise levels in gene expression of the model for two values of threshold, namely, T ← 1 and T ← 10 corresponding to low and high values of threshold, respectively. Each model was simulated up to 500 000 min to verify that the solution had arrived near the steady state. Each condition was simulated 100 times to guarantee that the result was repeatable. Note that the parameters in the model will never converge since the mutation rates are fixed. In Fig. 2, we present examples of the time evolution of the distribution of the shape parameter ˛i (t) in populations subject to each of the threshold values. This parameter determines the coefficient of variation of the gamma distribution (which equals ˛−1/2 ), and thus, is a measure of the variability of intervals between consecutive gene expression events. For any level of variability, an arbitrary rate of production can be attained by evolving ˇi (t) appropriately. As shown in the figure, when subject to a low threshold, the cell populations opt to evolve low noise in gene expression (i.e. lower noise than that of the Poisson process), whereas in the case of high threshold they favor more noisy expression dynamics. In either case, the initial condition (˛ = 1, that is, Poissonian dynamics) is less fit, implying that there is evolutionary pressure to favor nonPoissonian dynamics under these conditions. The evolved parameter values allow for higher probability of threshold-passing than a Poisson distribution with equivalent mean level would. This results in higher probability of reproduction with comparable level of nutrients consumption. In the case of low threshold, the nutrients can be accumulated to prevent starvation, while competence can be reached with more precisely timed productions of the proteins. In the case of high threshold, a
Fig. 2. Time evolution of the distribution of the scale parameter ˛(t), in an environment with low threshold for competence (upper panel), and in an environment with high threshold (lower panel). Different levels of gray represent quantiles of the population, with median denoted by the white dashed line, and maximum and minimum denoted by the black dashed lines.
different strategy is preferred. Namely, it is favorable to have noisier dynamics of gene expression, allowing a sufficiently large fraction of the population to survive by producing proteins in bursts (one of which will cross the threshold), while coping with the starvation by having lower mean levels. These resulting protein distributions are shown in Fig. 3, along with the expected distribution if gene expression was Poissonian, that is, for ˛ = 1, and of equal rate. While visually these distributions do not appear to be drastically different to the Poisson distribution, the evolution of the shape parameter ˛ in both cases (see Fig. 3) shows that the non-Poissonian distributions provide significant advantage in comparison with the Poissonian, since the cells with the former are completely displaced the cells with the latter. Also, we performed a Chi-squared test to obtain the probability that such distributions or more extreme values (p-value) would be obtained from the Poisson distributions with equivalent means. In both cases the p-value was found to be smaller than the precision of
Fig. 3. Example distributions of protein concentration in cells, at moment 50 000 min, in an environment with low threshold for competence (upper panel), and in an environment with high threshold for competence (lower panel). Solid lines with round markers represent Poisson distributions with equivalent production rate and the vertical dashed lines represent the threshold.
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
Fig. 4. Distribution of the scale parameter ˛(50 000) as a function of the threshold for competence (T). Different levels of gray represent quantiles of the population, with the median denoted by the white dashed line, and the maximum and minimum denoted by the black dashed lines.
the calculations (as performed in MATLAB software) implying that the distributions are likely not to be Poissonian. For both values of threshold, more probability mass of the distribution is placed in the region that passes the threshold than in the Poissonian case. This provides each cell with greater probability of crossing the threshold for competence, and consequently allowing reproduction. Additionally, in both cases, the bulk of the distribution lies below the threshold. Since the rate of gene expression is proportional to the consumption of nutrients, it is favorable to decrease the rate of the protein production to conserve the available nutrients and prevent starvation, as long as the competence can still be reached with a probability high enough. In the case of low threshold, the probability mass has been shifted to the bar corresponding to having one protein per cell. The rest of the tail is reduced, since producing extra proteins is not advantageous. Consequently, the overall variance of this distribution is reduced in comparison with the Poisson distribution. Due to the constraints on the shape of the distribution of the intervals between proteins productions in this model, the shape of the tail can have a restricted form, preventing its complete elimination. Similarly, the effects to the shape of the distribution are also visible in the case of high threshold. However, since the threshold is shifted, the distribution evolves to a different shape. Again, the bulk of the mass is placed below the threshold, making the expected consumption of nutrients small. This is accompanied with a fat right tail, that allow strong fluctuations in protein levels by which the cell can transiently reach competence. Finally, we tested how changing the threshold T affects to the evolved value of the shape parameter ˛. For this, we simulated the model for several values of threshold and recorded the distribution of ˛s in the end of the simulation. The results are presented in Fig. 4. From these results, with the current set of parameters values, for threshold values smaller than 3 the sub-Poissonian gene expression dynamics is favored. For threshold values greater than 3 the superPoissonian dynamics is favored. These results are in agreement with the previous simulations for the two threshold values (T = 1 and T = 10). 4. Conclusions Genome-wide studies of protein numbers in E. coli have shown that these differ widely in mean and in cell to cell diversity
15
(Taniguchi et al., 2010). Recent studies have associated these differences to, among other reasons, differences in the kinetics of the rate limiting steps during gene expression (Kandhavelu et al., 2011), which have been found to exist, for example, at the level of transcription initiation (McClure, 1985). Here we explored how different kinetics of gene expression may evolve, pressured by the needs to adapt to stochasticity in environmental conditions associated to the process that the protein is involved in. For that, we investigated the evolutionary dynamics of a gene responsible for cell division in a stochastic environment where nutrients are scarce and fluctuate in numbers. We found that small changes in the environmental constraints can drive a population of cells to evolve widely diverse strategies. Moreover, we showed that small changes in the kinetics of gene expression, such as tuning the strength of the fluctuations in the protein levels, can result in significant selective advantage. It is of interest to note that in both environmental conditions tested, selection leads to a far-from-Poissonian gene expression kinetics. This is in agreement with recent genome wide measurements of cell to cell diversity in protein numbers (Taniguchi et al., 2010). We hypothesize that the differing stochasticity of the processes that the various proteins are involved in is likely to be a major driver of different evolutionary pathways taken by the genes, in what concerns the kinetics of their expression. For example, processes fundamental to the survival of cells that require very few proteins may have evolved stringent control of stochasticity so as to guarantee that all cells have a one to few proteins at any given moment. Other processes, requiring large number of proteins, or occurring very seldom in time, may have led to the evolution of more noisier processes of protein production. In the future, it would be of interest to extend this study to the level of small genetic circuits associated with specific functions, such as tracking time or involved in decision-making. Acknowledgement This work was supported by the Academy of Finland. References Acar, M., Mettetal, J.T., van Oudenaarden, A., 2008. Stochastic switching as a survival strategy in fluctuating environments. Nat. Gen. 40 (4), 471–475. Bernstein, J.A., Khodursky, A.B., Lin, P.-H., Lin-Chao, S., Cohen, S.N., 2002. Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. Proc. Natl. Acad. Sci. U.S.A. 99 (15), 9697–9702. Bratsun, D., Volfson, D., Tsimring, L.S., Hasty, J., 2005. Delay-induced stochastic oscillations in gene regulation. Proc. Natl. Acad. Sci. U.S.A. 102 (41), 14593–14598. Choi, P.J., Xie, X.S., Shaknovich, E.I., 2010. Stochastic switching in gene networks can occur by a single-molecule event or many molecular steps. J. Mol. Biol. 396 (1), 230–244. Cohen, D., 1966. Optimizing reproduction in a randomly varying environment. J. Theor. Biol. 12 (1), 119–129. Eckert, B., Beck, C.F., 1989. Overproduction of transposon Tn10-encoded tetracycline resistance protein results in cell death and loss of membrane potential. J. Bacteriol. 171 (6), 3557–3559. Fraser, H.B., Hirsh, A.E., Giaever, G., Kumm, J., Eisen, M.B., 2004. Noise minimization in eukaryotic gene expression. PLoS Biol. 2 (6), e137. Gillespie, D.T., 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22 (4), 403–434. Golding, I., Paulsson, J., Zawilski, S.M., Cox, E.C., 2005. Real-time kinetics of gene activity in individual bacteria. Cell 123 (6), 1025–1036. Joseleau-Petit, D., Vinella, D., D’ari, R., 1999. Metabolic alarms and cell division in Escherichia coli. J. Bacteriol. 181 (1), 9–14. Kaern, M., Elston, T.C., Blake, W.J., Collins, J.J., 2005. Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet. 6 (6), 451–464. Kandhavelu, M., Mannerstrom, H., Gupta, A., Hakkinen, A., Lloyd-Price, J., Yli-Harja, O., Ribeiro, A.S., 2011. In vivo kinetics of transcription initiation of the lar promoter in Escherichia coli: evidence for a sequential mechanism with two rate-limiting steps. BMC Syst. Biol. 5 (1), 149. Kennell, D., Riezman, H., 1977. Transcription and translation initiation frequencies of the Escherichia coli lac operon. J. Mol. Biol. 114 (1), 1–21. Kussell, E., Kishony, R., Balaban, N.Q., Leibler, S., 2005. Bacterial persistence: a model of survival in changing environments. Genetics 169 (4), 1807–1814.
16
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
Kussell, E., Leibler, S., 2005. Phenotypic diversity, population growth, and information in fluctuating environments. Science 309 (5743), 2075–2078. Lachmann, M., Jablonka, E., 1996. The inheritance of phenotypes: an adaptation to fluctuating environments. J. Theor. Biol. 181 (1), 1–9. Levins, R., 1968. Evolution in Changing Environments: Some Theoretical Explorations. Princeton University Press, Princeton, NJ. Li, H., Petzold, L., 2006. Logarithmic Direct Method for Discrete Stochastic Simulation of Chemically Reacting Systems. Tech. Rep., Department of Computer Science, University of California, Santa Barbara. McAdams, H.H., Arkin, A., 1998. Simulations of prokaryotic genetic circuits. Annu. Rev. Biophys. Biomol. Struct. 27 (1), 199–224. McClure, W.R., 1985. Mechanism and control of transcription initiation in prokaryotes. Annu. Rev. Biochem. 54 (1), 171–204. Mott, M.L., Berger, J.M., 2007. DNA replication initiation: mechanisms and regulation in bacteria. Nat. Rev. Microbiol. 5 (5), 343–354. Peccoud, J., Ycart, B., 1995. Markovian modelling of gene product synthesis. Theor. Popul. Biol. 48 (2), 222–234. Peterson, C.N., Mandel, M.J., Silhavy, T.J., 2005. Escherichia coli starvation diets: essential nutrients weigh in distinctly. J. Bacteriol. 187 (22). Rajala, T., Hakkinen, A., Healy, S., Yli-Harja, O., Ribeiro, A.S., 2010. Effects of transcriptional pausing on gene expression dynamics. PLoS Comput. Biol. 6 (3), e1000704. Ribeiro, A.S., 2008. Dynamics and evolution of stochastic bistable gene networks with sensing in fluctuating environments. Phys. Rev. E 78 (6), 061902. Ribeiro, A.S., Hakkinen, A., Mannerstrom, H., Lloyd-Price, J., Yli-Harja, O., 2010. Effects of the promoter open complex formation on gene expression dynamics. Phys. Rev. E 81 (1), 011912.
Ribeiro, A.S., Lloyd-Price, J., 2007. SGN Sim, a stochastic genetic networks simulator. Bioinformatics 23 (6), 777–779. Ribeiro, A.S., Zhu, R., Kauffman, S.A., 2006. A general modeling strategy for gene regulatory networks with stochastic dynamics. J. Comp. Biol. 13 (9), 1630–1639. Ribero, A.S., 2010. Stochastic and delayed stochastic models of gene expression and regulation. Math. Biosci. 223 (1), 1–11. Ribero, A.S., Hakkinen, A., Healy, S., Yli-Harja, O., 2010. Dynamical effects of transcriptional pause-prone sites. Comp. Biol. Chem. 34 (3), 143–148. Roussel, M.R., Zhu, R., 2006. Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression title. Phys. Biol. 3 (4), 274–284. Schaffer, W.M., 1974. Optimal reproductive efforts in fluctuating environments. Am. Nat. 108 (964), 783–790. Taniguchi, Y., Choi, P.J., Li, G.-W., Chen, H., Babu, M., Hearn, J., Emili, A., Xie, X.S., 2010. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 329 (5991), 533–538. Wang, J.D., Levin, P.A., 2009. Metabolism, cell growth and the bacterial cell cycle. Nat. Rev. Microbiol. 7 (11), 822–827. Yu, J., Xiao, J., Run, X., Lao, K., Xie, X.S., 2006. Probing gene expression in live cells, one protein molecule at a time. Science 311 (5767), 1600–1603. Zhu, R., Ribeiro, A.S., Salahub, D., Kauffman, S.A., 2007. Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models. J. Theor. Biol. 246 (4), 725–745.
Contents lists available at SciVerse ScienceDirect
Computational Biology and Chemistry journal homepage: www.elsevier.com/locate/compbiolchem
Evolving kinetics of gene expression in stochastic environments Antti Häkkinen, Andre S. Ribeiro ∗ Laboratory of Biosystem Dynamics, Computational Systems Biology Research Group, Department of Signal Processing, Tampere University of Technology, Finland
a r t i c l e
i n f o
Article history: Received 4 November 2011 Received in revised form 18 January 2012 Accepted 14 February 2012 Keywords: Gene expression Stochastic dynamics Genetic operators Fluctuating environments Transcriptional bursting
a b s t r a c t Recent studies have shown that the in vivo dynamics of RNA numbers in bacteria is regulated, to a great extent, by the kinetics of rate limiting steps in transcription. Strong evidence suggests that the kinetics of these steps is sequence dependent. We investigate the selective advantages of rate limiting steps of differing kinetics. For that, we model the kinetics of expression of a gene responsible for promoting cell division at the expense of resources in the environment in individual cells of a population. We model mutations that affect the kinetics of the rate limiting steps and selective pressure in various environmental conditions. Depletion of resources leads to cell death. We find that small changes in the evolutionary constraints can favor widely different noise levels in RNA and protein numbers. Increasing the cost in nutrients for division favors noisier expression. The results provide a better understanding of why different genes differ in the kinetics of production of RNA and proteins. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Different genes have evolved to express at different rates and with different dynamical patterns (McClure, 1985). Genome wide observations of protein numbers in individual Escherichia coli cells have shown that the cell to cell diversity varies widely from one protein to the next and can range from sub- to super-Poissonian (Taniguchi et al., 2010). This diversity is likely present at the RNA level as well, since protein numbers tend to follow RNA numbers in bacteria (Kaern et al., 2005). Evidence suggests that most RNA molecules exist in low copy numbers in organisms such as E. coli, since most genes are expressed one to a few times during the cell lifetime (Bernstein et al., 2002). Due to this, and because protein numbers follow RNA numbers, small fluctuations in RNA numbers are likely to have downstream effects on the phenotype (Choi et al., 2010). Organisms may therefore have evolved mechanisms to stringently control the stochasticity in gene expression. In general, essential genes appear to have expression levels that are less noisy than that of genes associated with, for example, response to perturbations (Fraser et al., 2004; Taniguchi et al., 2010). The expression rate of a gene is determined by several factors and varies for different environmental conditions. On one hand,
∗ Corresponding author at: Department of Signal Processing, Tampere University of Technology, P.O. Box 553, 33101 Tampere, Finland. Tel.: +358 3 3115 4975; fax: +358 3 3115 4989. E-mail addresses: antti.hakkinen@tut.fi (A. Häkkinen), andre.ribeiro@tut.fi (A.S. Ribeiro). 1476-9271/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiolchem.2012.02.002
genes need to express at a rate fast enough to produce the minimum number of proteins required to execute a task. On the other hand, the rate of production ought to be limited to avoid waste of energy by the cell. In addition, expression rates are likely to be limited by other constraints such as toxicity beyond certain concentrations (Eckert and Beck, 1989). There is evidence suggesting that one core regulatory process of fluctuations in protein and RNA temporal levels is the multi-step process of transcription initiation (Ribeiro et al., 2010; Kandhavelu et al., 2011). This process has, in general, one to a few rate limiting steps, that is, steps whose duration determine the rate of transcription. Therefore, these steps have a strong effect on the strength of fluctuations in RNA and, consequently, protein numbers in bacteria. Other events also affecting fluctuations can occur at later stages, such as sequence dependent pauses in transcriptional elongation (Rajala et al., 2010; Ribero et al., 2010). Current techniques do not allow in vivo measurements of the durations of the underlying steps of transcription initiation. However, it is possible to measure the intervals between the completion of consecutive transcription events (Kandhavelu et al., 2011). This, along with recent measurement techniques of protein numbers (Yu et al., 2006), has allowed the development of accurate stochastic models of gene expression (Zhu et al., 2007; Ribero, 2010). The first stochastic models of genetic circuits assumed gene expression as an instantaneous process (McAdams and Arkin, 1998). The dynamics was driven by the Stochastic Simulation Algorithm (SSA, see Gillespie, 1976). Several recent studies (Bratsun et al., 2005; Roussel and Zhu, 2006; Ribeiro et al., 2006) have shown that the time for gene expression to be completed once initiated is not negligible in the dynamics of genetic circuits. This duration,
12
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
while sequence dependent, also has a random nature, varying even between transcription events of the same gene (Kandhavelu et al., 2011) and thus, is also a source of stochasticity. The delayed Stochastic Simulation Algorithm (DSSA, see Roussel and Zhu, 2006) was developed as an extension of the original SSA to allow the use of delayed reaction products besides the elementary reactions with exponential statistics. This is useful to simulate genetic circuits, since some of the underlying processes consist of sequential events, exhibiting Gaussian-like rather than exponential-like statistics. Since the kinetics of these underlying events are unknown (e.g. in the case of transcription initiation), the explicit modeling is not be possible. If one is to model evolutionary processes of genetic circuits, the evolvability of these rate limiting steps in transcription and translation need to be considered. Finally, genes and genetic circuits are not isolated systems. Rather, they have evolved to function within particular ranges of environmental conditions and, to some extent, to cope with stochastic fluctuations in the environment. The stochastic nature of the environments also needs to be accounted for when simulating and studying how genetic circuits evolve. Several works have focused on the evolutionary dynamics of genes or small genetic circuits, when subject to stochastic environments (see e.g. Cohen, 1966; Levins, 1968; Schaffer, 1974). Recently, this topic has been revisited in the context of how initially isogenic populations may adapt to fluctuating environments (Lachmann and Jablonka, 1996; Kussell et al., 2005; Kussell and Leibler, 2005). One of the interesting findings from these studies is that the stochastic switching between two phenotypes may be used by cells as a survival strategy when subject to fluctuating environments (Acar et al., 2008). Another study investigated how genetic switches can better adjust their switching dynamics provided that they possess a sensing mechanism of the environment state, such as the detection of the amount of toxins or nutrients in the environment at a given moment (Ribeiro, 2008). In these models, gene expression delays are either not present, or modeled as constant from one event to the next. Additionally, the cell populations are constrained to fixed sizes. Here we study the dynamics and evolutionary pathways of a delayed stochastic model of gene expression in cell populations of variable size subject to stochastic environments. For that, we model the evolution of the dynamics of a gene responsible for expressing a protein required for cell division in an environment with limited, stochastic supply of nutrients required for protein expression. We allow mutations and selection to act on the kinetics of the rate limiting steps of gene expression, and test whether the environmental pressure drives the evolution of different kinetics. To test this hypothesis, we developed an extension for SGN Sim (Ribeiro and Lloyd-Price, 2007), that allows modeling evolutionary patterns of genetic circuits of populations variable in size and the fluctuating environments they are subject to. The features of this extension are described in the next section.
interpreter. A simulation involves simulating a set of entities, who are independent and whose dynamics are governed by the DSSA, and scheduling events that manipulate these entities to the event queue. The entities can be created, destroyed, and mutated dynamically, and the times at which they occur are drawn from arbitrary distributions, either using the built-in random number generation facilities or ones interfaced with Lua. The simulation of the entities uses DSSA, which is an extension of the SSA (Gillespie, 1976). This extension allows non-Markovian dynamics, by attaching an arbitrarily distributed time delay before the release of each of the reaction products. The DSSA simulation engine is based on that of SGN Sim (Ribeiro and Lloyd-Price, 2007), which implements the algorithm as specified by Roussel and Zhu (2006). The simulation algorithm is a form of the logarithmic direct method (Li and Petzold, 2006), modified for DSSA. The entities can be created based on a previously created entity, or based on a system description. The system descriptions are written in SGN Sim format (Ribeiro and Lloyd-Price, 2007). Moreover, the entities can be modified immediately after their creation, such that schemes like reproduction by division can be implemented. Mutations can be applied by queueing mutating operators in the event queue, with time intervals drawn from arbitrary distributions, affecting the molecular counts or the parameters of the model. The scheme of mutations, such as to which parameters the mutations affect, and to what extent they are mutated are left to be specified by the user. It is also possible model more complex schemes of that introduced genetic variation, for example, by exchanging the parameters values between separate entities. At any given point in time, a subpopulation of the simulated entities can be selected, based on arbitrary combination of features in their state, such as concentration of some molecules. Combined with the operators of reproduction or destruction, this allows modeling the process of natural selection, based on an arbitrary fitness function. Moreover, the events that can be scheduled are not limited to the creation, destruction, and modification of the entities. It is possible to schedule events that are arbitrary Lua functions; however, this has the disadvantage that the parallel simulations must be synchronized. In cases where it is applicable, such as in the case of cell division, this synchronization is not necessary and can be avoided by combining the non-locking built-in operations. The built-in Lua interpreter is used to parse the description of the evolutionary process of the simulator. This description consist of queuing the events to the event queues, and to set up parameters such as the length of the simulation. The events can be queued as one-shot events, possibly being queued again at the occurrence of some event, or as recurring events. The time intervals between these events can be drawn from random number generators of various distributions, or specified arbitrarily. It can also be used to execute arbitrary pieces of Lua code.
2.2. Model 2. Methods 2.1. Simulation implementation To study the evolutionary phenomena in genetic networks, we implement a simulator that combines simulation of stochastic molecular kinetics with evolutionary operators. The simulator is written in C++and can be compiled on Microsoft Windows and typical Unix-like systems, such as Linux. Multithreading is exploited to accelerate simulations on hardware that supports parallel execution. The simulator is available upon request. The simulator consists of a DSSA simulation engine, a parallel event queueing system (priority queue), and a built-in Lua
We study the evolutionary patterns of the kinetics of expression of a gene responsible for the production of a protein necessary for cell division. The expression is triggered by the presence of nutrients in the cell. This is a common process in real cells. For example, evidence suggests that, when starving, the gene expression patterns of the cells change, the nature of change depending on which nutrients are lacking (Peterson et al., 2005). Also, starving cells commonly inhibit DNA replication (Joseleau-Petit et al., 1999). For example, under amino acid starvation, DNA replication is inhibited due to the inhibition of DnaA expression (Wang and Levin, 2009), a protein required in multiple copies to initiate the DNA replication (Mott and Berger, 2007).
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
13
Fig. 1. A pictorial overview of the model. Nutrients are created in the environment and distributed to the cells. The cells require nutrients to express proteins. Accumulation of these proteins is required for cell division. Also depicted are the kinetic rate constants of the several processes in the model, described below.
To mimic these mechanisms and their evolution, in the model proposed here the expression of the protein is triggered by the presence of nutrients in the cell. However, cells cannot deplete all the nutrients in the division process, as this leads to death. The size of the cell population is dynamic and all individuals compete for dividing during their limited lifetime, , in an environment where nutrients are generated by a stochastic process with a constant rate and distributed uniformly to the cells. The simulations of the evolutionary process aim at determining the optimal set of parameters values for both the mean and the strength of the fluctuations in protein levels, as a function of the environmental conditions. Also, we study how the optima change with the environmental settings. A pictorial description of the model is shown in Fig. 1. The model consists of a set of chemical interactions that are implemented in each individual cell. In the reactions, i denotes the cell index, to make a distinction between the variables that are not shared by the individuals. Nutrients ni are transported to a cell by the following reaction equation (RE): kn N(t)−1
∅ −→ ni ,
(1)
where kn is the total rate at which novel nutrients units become available, and N(t) is the cell population size at moment t. Consequently, reaction (1) is such that, the larger is the population, the less nutrients ni will be available for each cell. The expression of proteins uses nutrient units as substrate (RE (2)). Proteins and nutrients are subject to degradation via reactions (3) and (4), respectively: ∞
Pi + ni −→Pi (i (t)) + pi dp
pi −→∅ dn
ni −→∅
(2) (3)
,
(4)
where Pi denotes the promoter region of the gene, and pi denotes the protein. The rate of RE (2) is set to infinity, such that the production of proteins is only limited by the absence of nutrients and by the time delay i .
The time delay i (t) determines the time between the production of two consecutive proteins. This delay accounts for the underlying processes during transcription initiation (McClure, 1985) and recent evidence suggests that these processes are the most relevant in determining the rate of RNA production in bacteria (Kandhavelu et al., 2011). In prokaryotes, evidence suggests that the rate of protein production generally tends to follow that of the RNA production (Kaern et al., 2005). We set the delay i (t) ∼ Gamma(˛i (t), ˇi (t)) to be gamma distributed with shape parameter ˛i (t) and rate parameter ˇi (t). The gamma distribution allows modeling a wide range of dynamics by varying its parameters ˛i (t) and ˇi (t). Depending on the values of the parameters, this distribution allows a sub-Poissonian (˛i (t) > 1) dynamics of protein production, e.g. resultant of a sequential rate-limiting mechanism of gene expression (Kandhavelu et al., 2011), or a super-Poissonian (˛i (t) < 1) dynamics, which results from models such as the on–off model (Peccoud and Ycart, 1995). Mutations to the parameters ˛i (t) and ˇi (t) allow the kinetics of protein production to change. The mutations are performed with exponentially distributed time intervals with a rate of km , in a linear scale with step size of m, with uniform probability to either increase or decrease. Finally, when a threshold T in the number of pi is reached, the cell is able to divide. In this process, two daughter cells are formed, who split equally pi and ni of the mother cell, and inherit the values of ˛i (t) and ˇi (t), which regulate the kinetics of gene expression. 3. Results In optimal growth conditions the division time of E. coli is 20 min. Thus, we set the maximum cell lifetime = 100 min. Given the parameters described below, we observed that for the initial parameter values, the average division time was approximately 40 min. We choose initial shape of ˛i (0) = 1, which corresponds to Poissonian production dynamics of proteins. Initially, the expression and degradation rates are set such that the mean protein level is 5. For that we set the production rate ˇi (0) = 5 min−1 . This value,
14
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
while high, is within realistic intervals, since e.g. the lac promoter can initialize an RNA production every three seconds (Kennell and Riezman, 1977), while the lar promoter produces one RNA per 15 min under full induction (Golding et al., 2005; Kandhavelu et al., 2011). Next, we set the environmental conditions, that is, the amount of available nutrients. We want to compare how the gene expression dynamics evolves in two different conditions, namely, when many nutrients are required for dividing and when very few nutrients are required. However, we also need to have populations of comparable size in the two simulations, since the number of cells determines e.g. how many mutations overall can occur in a generation. To have populations of similar size, we scale the amount of nutrients in the environment by the amount of nutrients needed for division. Thus, we set kn = 20T min−1 and dn = 0.1 min−1 . Finally, regarding the gene expression mechanism, we need to set a mutation rate that determines the speed at which the cell population can adapt to the environment. As mentioned, ˛i and ˇi are subject to changes due to mutations. We found that the speed by which these parameters change does not, within a wide range of values, affect the final outcome of the selection process. Due to that, for simplicity, we set the mutation rate to km = 1 min−1 and step size m = 0.02, arbitrarily, as these values permit reaching near optimal values in a reasonable number of generations. In the results presented, the initial population size is set to N(0) = 10. We investigated several values of initial population sizes, namely N(0) ∈ {1, 2, 5, 10, 20, 50, 100, 200, 500}, and found no significant differences in their long term behavior. This is because, in the beginning of each simulation, the population size is first reduced to a small value due to starvation, after which the population size starts growing toward a limiting value as the cells become better adapted to the environment. In the first generation, the promoter is initially set to be unoccupied ([Pi ] ← 0), and the protein concentration is set to reflect the expected protein number given infinite source of food ([pi ] ← ˛i (0)−1 ˇi (0) dp−1 ). Moreover, the initial number of nutrients [ni ] in the cell is set to zero. We study the evolved noise levels in gene expression of the model for two values of threshold, namely, T ← 1 and T ← 10 corresponding to low and high values of threshold, respectively. Each model was simulated up to 500 000 min to verify that the solution had arrived near the steady state. Each condition was simulated 100 times to guarantee that the result was repeatable. Note that the parameters in the model will never converge since the mutation rates are fixed. In Fig. 2, we present examples of the time evolution of the distribution of the shape parameter ˛i (t) in populations subject to each of the threshold values. This parameter determines the coefficient of variation of the gamma distribution (which equals ˛−1/2 ), and thus, is a measure of the variability of intervals between consecutive gene expression events. For any level of variability, an arbitrary rate of production can be attained by evolving ˇi (t) appropriately. As shown in the figure, when subject to a low threshold, the cell populations opt to evolve low noise in gene expression (i.e. lower noise than that of the Poisson process), whereas in the case of high threshold they favor more noisy expression dynamics. In either case, the initial condition (˛ = 1, that is, Poissonian dynamics) is less fit, implying that there is evolutionary pressure to favor nonPoissonian dynamics under these conditions. The evolved parameter values allow for higher probability of threshold-passing than a Poisson distribution with equivalent mean level would. This results in higher probability of reproduction with comparable level of nutrients consumption. In the case of low threshold, the nutrients can be accumulated to prevent starvation, while competence can be reached with more precisely timed productions of the proteins. In the case of high threshold, a
Fig. 2. Time evolution of the distribution of the scale parameter ˛(t), in an environment with low threshold for competence (upper panel), and in an environment with high threshold (lower panel). Different levels of gray represent quantiles of the population, with median denoted by the white dashed line, and maximum and minimum denoted by the black dashed lines.
different strategy is preferred. Namely, it is favorable to have noisier dynamics of gene expression, allowing a sufficiently large fraction of the population to survive by producing proteins in bursts (one of which will cross the threshold), while coping with the starvation by having lower mean levels. These resulting protein distributions are shown in Fig. 3, along with the expected distribution if gene expression was Poissonian, that is, for ˛ = 1, and of equal rate. While visually these distributions do not appear to be drastically different to the Poisson distribution, the evolution of the shape parameter ˛ in both cases (see Fig. 3) shows that the non-Poissonian distributions provide significant advantage in comparison with the Poissonian, since the cells with the former are completely displaced the cells with the latter. Also, we performed a Chi-squared test to obtain the probability that such distributions or more extreme values (p-value) would be obtained from the Poisson distributions with equivalent means. In both cases the p-value was found to be smaller than the precision of
Fig. 3. Example distributions of protein concentration in cells, at moment 50 000 min, in an environment with low threshold for competence (upper panel), and in an environment with high threshold for competence (lower panel). Solid lines with round markers represent Poisson distributions with equivalent production rate and the vertical dashed lines represent the threshold.
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
Fig. 4. Distribution of the scale parameter ˛(50 000) as a function of the threshold for competence (T). Different levels of gray represent quantiles of the population, with the median denoted by the white dashed line, and the maximum and minimum denoted by the black dashed lines.
the calculations (as performed in MATLAB software) implying that the distributions are likely not to be Poissonian. For both values of threshold, more probability mass of the distribution is placed in the region that passes the threshold than in the Poissonian case. This provides each cell with greater probability of crossing the threshold for competence, and consequently allowing reproduction. Additionally, in both cases, the bulk of the distribution lies below the threshold. Since the rate of gene expression is proportional to the consumption of nutrients, it is favorable to decrease the rate of the protein production to conserve the available nutrients and prevent starvation, as long as the competence can still be reached with a probability high enough. In the case of low threshold, the probability mass has been shifted to the bar corresponding to having one protein per cell. The rest of the tail is reduced, since producing extra proteins is not advantageous. Consequently, the overall variance of this distribution is reduced in comparison with the Poisson distribution. Due to the constraints on the shape of the distribution of the intervals between proteins productions in this model, the shape of the tail can have a restricted form, preventing its complete elimination. Similarly, the effects to the shape of the distribution are also visible in the case of high threshold. However, since the threshold is shifted, the distribution evolves to a different shape. Again, the bulk of the mass is placed below the threshold, making the expected consumption of nutrients small. This is accompanied with a fat right tail, that allow strong fluctuations in protein levels by which the cell can transiently reach competence. Finally, we tested how changing the threshold T affects to the evolved value of the shape parameter ˛. For this, we simulated the model for several values of threshold and recorded the distribution of ˛s in the end of the simulation. The results are presented in Fig. 4. From these results, with the current set of parameters values, for threshold values smaller than 3 the sub-Poissonian gene expression dynamics is favored. For threshold values greater than 3 the superPoissonian dynamics is favored. These results are in agreement with the previous simulations for the two threshold values (T = 1 and T = 10). 4. Conclusions Genome-wide studies of protein numbers in E. coli have shown that these differ widely in mean and in cell to cell diversity
15
(Taniguchi et al., 2010). Recent studies have associated these differences to, among other reasons, differences in the kinetics of the rate limiting steps during gene expression (Kandhavelu et al., 2011), which have been found to exist, for example, at the level of transcription initiation (McClure, 1985). Here we explored how different kinetics of gene expression may evolve, pressured by the needs to adapt to stochasticity in environmental conditions associated to the process that the protein is involved in. For that, we investigated the evolutionary dynamics of a gene responsible for cell division in a stochastic environment where nutrients are scarce and fluctuate in numbers. We found that small changes in the environmental constraints can drive a population of cells to evolve widely diverse strategies. Moreover, we showed that small changes in the kinetics of gene expression, such as tuning the strength of the fluctuations in the protein levels, can result in significant selective advantage. It is of interest to note that in both environmental conditions tested, selection leads to a far-from-Poissonian gene expression kinetics. This is in agreement with recent genome wide measurements of cell to cell diversity in protein numbers (Taniguchi et al., 2010). We hypothesize that the differing stochasticity of the processes that the various proteins are involved in is likely to be a major driver of different evolutionary pathways taken by the genes, in what concerns the kinetics of their expression. For example, processes fundamental to the survival of cells that require very few proteins may have evolved stringent control of stochasticity so as to guarantee that all cells have a one to few proteins at any given moment. Other processes, requiring large number of proteins, or occurring very seldom in time, may have led to the evolution of more noisier processes of protein production. In the future, it would be of interest to extend this study to the level of small genetic circuits associated with specific functions, such as tracking time or involved in decision-making. Acknowledgement This work was supported by the Academy of Finland. References Acar, M., Mettetal, J.T., van Oudenaarden, A., 2008. Stochastic switching as a survival strategy in fluctuating environments. Nat. Gen. 40 (4), 471–475. Bernstein, J.A., Khodursky, A.B., Lin, P.-H., Lin-Chao, S., Cohen, S.N., 2002. Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. Proc. Natl. Acad. Sci. U.S.A. 99 (15), 9697–9702. Bratsun, D., Volfson, D., Tsimring, L.S., Hasty, J., 2005. Delay-induced stochastic oscillations in gene regulation. Proc. Natl. Acad. Sci. U.S.A. 102 (41), 14593–14598. Choi, P.J., Xie, X.S., Shaknovich, E.I., 2010. Stochastic switching in gene networks can occur by a single-molecule event or many molecular steps. J. Mol. Biol. 396 (1), 230–244. Cohen, D., 1966. Optimizing reproduction in a randomly varying environment. J. Theor. Biol. 12 (1), 119–129. Eckert, B., Beck, C.F., 1989. Overproduction of transposon Tn10-encoded tetracycline resistance protein results in cell death and loss of membrane potential. J. Bacteriol. 171 (6), 3557–3559. Fraser, H.B., Hirsh, A.E., Giaever, G., Kumm, J., Eisen, M.B., 2004. Noise minimization in eukaryotic gene expression. PLoS Biol. 2 (6), e137. Gillespie, D.T., 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22 (4), 403–434. Golding, I., Paulsson, J., Zawilski, S.M., Cox, E.C., 2005. Real-time kinetics of gene activity in individual bacteria. Cell 123 (6), 1025–1036. Joseleau-Petit, D., Vinella, D., D’ari, R., 1999. Metabolic alarms and cell division in Escherichia coli. J. Bacteriol. 181 (1), 9–14. Kaern, M., Elston, T.C., Blake, W.J., Collins, J.J., 2005. Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet. 6 (6), 451–464. Kandhavelu, M., Mannerstrom, H., Gupta, A., Hakkinen, A., Lloyd-Price, J., Yli-Harja, O., Ribeiro, A.S., 2011. In vivo kinetics of transcription initiation of the lar promoter in Escherichia coli: evidence for a sequential mechanism with two rate-limiting steps. BMC Syst. Biol. 5 (1), 149. Kennell, D., Riezman, H., 1977. Transcription and translation initiation frequencies of the Escherichia coli lac operon. J. Mol. Biol. 114 (1), 1–21. Kussell, E., Kishony, R., Balaban, N.Q., Leibler, S., 2005. Bacterial persistence: a model of survival in changing environments. Genetics 169 (4), 1807–1814.
16
A. Häkkinen, A.S. Ribeiro / Computational Biology and Chemistry 37 (2012) 11–16
Kussell, E., Leibler, S., 2005. Phenotypic diversity, population growth, and information in fluctuating environments. Science 309 (5743), 2075–2078. Lachmann, M., Jablonka, E., 1996. The inheritance of phenotypes: an adaptation to fluctuating environments. J. Theor. Biol. 181 (1), 1–9. Levins, R., 1968. Evolution in Changing Environments: Some Theoretical Explorations. Princeton University Press, Princeton, NJ. Li, H., Petzold, L., 2006. Logarithmic Direct Method for Discrete Stochastic Simulation of Chemically Reacting Systems. Tech. Rep., Department of Computer Science, University of California, Santa Barbara. McAdams, H.H., Arkin, A., 1998. Simulations of prokaryotic genetic circuits. Annu. Rev. Biophys. Biomol. Struct. 27 (1), 199–224. McClure, W.R., 1985. Mechanism and control of transcription initiation in prokaryotes. Annu. Rev. Biochem. 54 (1), 171–204. Mott, M.L., Berger, J.M., 2007. DNA replication initiation: mechanisms and regulation in bacteria. Nat. Rev. Microbiol. 5 (5), 343–354. Peccoud, J., Ycart, B., 1995. Markovian modelling of gene product synthesis. Theor. Popul. Biol. 48 (2), 222–234. Peterson, C.N., Mandel, M.J., Silhavy, T.J., 2005. Escherichia coli starvation diets: essential nutrients weigh in distinctly. J. Bacteriol. 187 (22). Rajala, T., Hakkinen, A., Healy, S., Yli-Harja, O., Ribeiro, A.S., 2010. Effects of transcriptional pausing on gene expression dynamics. PLoS Comput. Biol. 6 (3), e1000704. Ribeiro, A.S., 2008. Dynamics and evolution of stochastic bistable gene networks with sensing in fluctuating environments. Phys. Rev. E 78 (6), 061902. Ribeiro, A.S., Hakkinen, A., Mannerstrom, H., Lloyd-Price, J., Yli-Harja, O., 2010. Effects of the promoter open complex formation on gene expression dynamics. Phys. Rev. E 81 (1), 011912.
Ribeiro, A.S., Lloyd-Price, J., 2007. SGN Sim, a stochastic genetic networks simulator. Bioinformatics 23 (6), 777–779. Ribeiro, A.S., Zhu, R., Kauffman, S.A., 2006. A general modeling strategy for gene regulatory networks with stochastic dynamics. J. Comp. Biol. 13 (9), 1630–1639. Ribero, A.S., 2010. Stochastic and delayed stochastic models of gene expression and regulation. Math. Biosci. 223 (1), 1–11. Ribero, A.S., Hakkinen, A., Healy, S., Yli-Harja, O., 2010. Dynamical effects of transcriptional pause-prone sites. Comp. Biol. Chem. 34 (3), 143–148. Roussel, M.R., Zhu, R., 2006. Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression title. Phys. Biol. 3 (4), 274–284. Schaffer, W.M., 1974. Optimal reproductive efforts in fluctuating environments. Am. Nat. 108 (964), 783–790. Taniguchi, Y., Choi, P.J., Li, G.-W., Chen, H., Babu, M., Hearn, J., Emili, A., Xie, X.S., 2010. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 329 (5991), 533–538. Wang, J.D., Levin, P.A., 2009. Metabolism, cell growth and the bacterial cell cycle. Nat. Rev. Microbiol. 7 (11), 822–827. Yu, J., Xiao, J., Run, X., Lao, K., Xie, X.S., 2006. Probing gene expression in live cells, one protein molecule at a time. Science 311 (5767), 1600–1603. Zhu, R., Ribeiro, A.S., Salahub, D., Kauffman, S.A., 2007. Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models. J. Theor. Biol. 246 (4), 725–745.