1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
Modeling stochastic gene expression in growing cells Q1
David Gomez a,b,n, Rahul Marathe a,c,d, Veronika Bierbaum a, Stefan Klumpp a a
Max Planck Institute of Colloids and Interfaces, Science Park Golm, 14424 Potsdam, Germany Department of Physics, Freie Universität Berlin, Arnimallee 14, 14195 Berlin, Germany Institute for Theoretical Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany d Department of Physics, University of Pune, Ganeshkhind, 411007 Pune, India b c
H I G H L I G H T S
Volume growth is incorporated in stochastic simulations of gene expression. The bistability region diminishes via a reduced copy number of the regulatory protein. Cell volume determines the region of bistability for different noise strengths. The method can easily be applied to more complex gene circuits.
art ic l e i nf o
a b s t r a c t
Article history: Received 24 June 2013 Received in revised form 8 January 2014 Accepted 15 January 2014
Gene expression is an inherently noisy process. Fluctuations arise at many points in the expression of a gene, as all the salient reactions such as transcription, translation, and mRNA degradation are stochastic processes. The fluctuations become important when the cellular copy numbers of the relevant molecules (mRNA or proteins) are low. For regulated genes, a computational complication arises from the fact that protein synthesis rates depend on the concentrations of the transcription factors that regulate the corresponding genes. Because of the growing cell volume, such rates are effectively time-dependent. We deal with the effects of volume growth computationally using a rather simple method: the growth of the cell volume is incorporated in our simulations by stochastically adding small volume elements to the cell volume. As an application of this method we study a gene circuit with positive autoregulation that exhibits bistability. We show how the region of bistability becomes diminished by increasing the effect of noise via a reduced copy number of the regulatory protein. Cell volume determines the region of bistability for different noise strengths. The method is general and can also be applied to other cases where synthesis rates of proteins are regulated and an appropriate analytical description is difficult to achieve. & 2014 Published by Elsevier Ltd.
Keywords: Genetic circuits Stochastic gene expression Noise Cell division Volume growth Autoregulation Bistability
1. Introduction Gene expression is a complex process constituted by various interactions between proteins and DNA. Theoretical and experimental methods have been developed recently in order to understand the regulation of gene expression quantitatively (Raser and O'Shea, 2005; Longo and Hasty, 2006; Elowitz et al., 2002; So et al., 2011; Thattai Q3 and van Oudenaarden, 2001). An important issue in this quest is the appropriate mathematical description of gene expression dynamics. Mathematical models of gene regulation systems have been developed at many levels of description, including logical approaches
n
Corresponding author at: Max Planck Institute of Colloids and Interfaces,
Q2 Science Park Golm, 14424 Potsdam, Germany. Tel.: þ 49 331 567 9622. E-mail address: [email protected] (D. Gomez).
(boolean networks), deterministic dynamic systems, and stochastic approaches (McAdams and Arkin, 1997; Paulsson, 2005; Thomas, 1984). Stochastic descriptions are essential when the copy numbers of key molecules are low (Raj et al., 2006; Simpson et al., 2003; Berg, 1978; Elowitz et al., 2002), as it is often the case, or if rare events triggered by fluctuations induce large changes in the system. Low copy numbers are quite common for transcription factors in bacteria (Raj and van Oudenaarden, 2008) and can for example lead to heterogeneous timing of regulatory response to environmental signals (Boulineau et al., 2013; Megerle et al., 2008). An important example for rare stochastic events with consequences for the physiological state of the whole cell is the phenomenon of phenotype switching (Smits et al., 2006) where cells with the same genotype and under the same environmental conditions can display different stable patterns of gene expression and stochastically switch between them, for example in bacterial persistence (Balaban et al., 2004; Patra and
0022-5193/$ - see front matter & 2014 Published by Elsevier Ltd. http://dx.doi.org/10.1016/j.jtbi.2014.01.017
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
Klumpp, 2013), competence and sporulation (Dubnau and Losick, 2006; Elowitz and Leibler, 2000; Leisner et al., 2009). Stochasticity arises from several sources such as the inherent randomness of the protein synthesis reaction (Swain et al., 2002), the partitioning of molecules during cell division (Munsky et al., 2012), as well as from extrinsic fluctuations of global parameters of the cell that influence gene expression (Elowitz et al., 2002; Taniguchi et al., 2010; Paulsson, 2004) Stochastic models are often studied by computer simulations using the Gillespie algorithm (Li and Petzold, 2005; Li et al., 2008). This algorithm provides a general approach to the stochastic kinetics of biochemical reactions and other stochastic systems including ecological dynamics (Campillo and Lobry, 2012) and intracellular transport (Mueller et al., 2008). A general complication in the study of the dynamics of genetic circuits is that these circuits are coupled to the state of their host cell. Rather than simply being a bioreactor, the host cell itself is a dynamic system that grows and divides and that adapts to environmental conditions (Klumpp et al., 2009; Marathe et al., 2012a). In microbes, these processes occur on similar time scales as the dynamics of gene circuits. In stochastic simulations, the growth of the cell effectively makes all rates that depend on concentrations timedependent, as the concentration is permanently diluted by the increase of the volume of the growing cell (the number of molecules of a specific type per cell is also diluted, however not continuously by volume growth, but by cell division). In many studies these effects are not treated explicitly, but approximated by an effective degradation rate; taking them into account requires to deal with time dependent rates (Rosenfeld et al., 2002; Lu et al., 2004). In this study, we propose a simple method to take volume growth explicitly into account in stochastic simulations with the Gillespie algorithm: the basic idea is to describe volume growth as another stochastic process that adds small amounts of new volume with a certain rate. We implement both linear and exponential volume growth and also incorporate cell division and gene replication into our simulation. The method can be used either as a stochastic implementation of deterministic volume growth or as a description of truly stochastic volume growth, which contributes to the so-called ‘extrinsic noise’ (Swain et al., 2002). We apply this method to two cases, an unregulated gene and a positively autoregulated gene. The unregulated gene serves as a test case to control the additional noise generated by describing volume growth as a stochastic process. For the case with positive autoregulation, we determine the effect of noise on bistability and show that the parameter range for bistability is diminished by increasing the effect of noise via a reduced protein copy number. The method we use is general and can be applied to any situation where concentrations change continuously in time due to the increasing volume of a growing cell.
2. Model We run stochastic simulations of simple genetic circuits using the Gillespie algorithm (Gillespie, 1976; Hayot, 2008; Gillespie, 1977). To be specific, here we describe the case of one gene under the control of one transcription factor. Below, we will consider a coarse-grained description in which protein synthesis rates are dependent on the concentration of a regulatory protein. However, in general, systematic volume growth will affect all rates that depend on the concentration of one or several types of molecules. Generally unimolecular reactions such as transitions between conformational states or the spontaneous decay of a molecular complex do not depend on the volume in which the reaction takes place, so the corresponding rates are not affected by a changing
volume. Bimolecular reactions such as the formation of a complex or binding of a protein to a binding site on DNA, however, depend on the concentration of these molecules, and are thus affected by the reaction volume. The regulation of a gene involves complex processes of both types, for example protein–DNA and protein– protein binding and unbinding, initiation of transcription and translation by a DNA-bound RNA polymerase or a mRNA-bound ribosome, or protein degradation (Verma et al., 2006). These processes can either be modelled in detail or in a more coarsegrained fashion, but in either case some reactions will be dependent on the volume. Here we will use an effective description of gene regulation, where the synthesis rate of a gene product (protein or RNA) is dependent on the concentration of the corresponding regulatory protein. While our goal is a stochastic description of the dynamics of gene expression, it is simpler to formulate the corresponding deterministic model first. We describe the synthesis of the protein product of a gene by a single rate αðrÞ pg (thus we join the two steps of protein synthesis, transcription and translation into a single step and do not consider the bursty nature of stochastic gene expression (Berg, 1978; Thattai and van Oudenaarden, 2001), which however can easily be included into our approach). Here, αðrÞ is the protein synthesis rate (per gene) and g is the gene copy number. The protein synthesis rate depends on the concentration of transcription factors and other regulators, here exemplified by a single transcription factor with concentration r. Through the concentration of that transcription factor, rðtÞ ¼ RðtÞ=VðtÞ, which is obtained as the ratio of the instantaneous number R(t) of molecules of that transcription factors in the cell and the cell volume V(t), the protein synthesis rate changes as a function of time due to cell growth and synthesis of regulator molecules. Proteins are degraded with a rate β, which for simplicity we take as unregulated and thus constant. Thus, in the corresponding deterministic description, the time evolution of the protein content of the cell (P (t), in copy numbers of the molecule per cell) is given by P_ ðtÞ ¼ αðrðtÞÞgðtÞ β P:
ð1Þ
This equation is the natural starting point for a stochastic description. We note however that for deterministic descriptions, it is useful to consider the corresponding equation for the concentration pðtÞ ¼ PðtÞ=VðtÞ, which is obtained as ! αðrðtÞÞgðtÞ V_ _ ¼ βþ pðtÞ p; ð2Þ VðtÞ V Eq. (1) is valid between cell division events; at times of division the protein number is divided by two. The concentration equation, Eq. (2), can also be used beyond a division cycle, because the concentration is continuous across a division (Marathe et al., 2012a) (for exponentially growing cell cultures, V_ =V is equal to the population growth rate μ on large time scales, independent of the functional form of volume growth during the division cycle). In the following, we will consider stable proteins with negligible degradation rate, a typical case in bacteria (Nath and Koch, 1970) and thus set β ¼0. Nevertheless, the protein content is effectively diluted as the cell volume increases and as cells divide (while the concentration decreases continuously during the division cycle, the protein copy number is diluted in the discrete events of cell division). Let us now construct the stochastic model. Synthesis and degradation are described by rates and can easily be implemented in the Gillespie algorithm. In addition, we need to account for volume growth, cell division and gene replication. Volume growth leads to a time-dependence of all rates that depend on concentrations, in our case, all synthesis rates are subject to regulation by transcription factors. To deal with this time dependence, we implement the growth
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
of the cell volume as an additional stochastic reaction, which adds a small volume ΔV to the current cell volume. The concentrations of all proteins are constant between these stochastic events of volume growth. Thus all rates depend only on the copy numbers of proteins in the cell and on the instantaneous cell volume, and are constant between stochastic events. In this way, volume growth and protein synthesis are described by the simplest Gillespie algorithm without time-dependent rates. This stochastic dynamics is interrupted by cell division and by gene replication, which occur at deterministically determined times rather than at stochastic times. The implementation of volume growth and that of the two deterministic processes are described in more detail in the following two subsections. We note that the specific assumptions made here for the volume growth function and for the description of cell division are chosen to model gene expression in bacteria, but adjusting these assumptions, the method of stochastic addition of volume elements should also be applicable to eukaryotic systems. 2.1. Volume growth To deal with the computational complication of timedependent rates due to volume growth, we implement volume growth as an additional reaction into the Gillespie algorithm. In this way, our algorithm provides a systematic method for obtaining statistically correct solutions for the concentration dynamics in two general steps. First, the Gillespie algorithm iteratively generates exponentially distributed times from the sum of the protein synthesis and volume growth rates. These times set the time of the next ‘event.’ In our case, an event is either synthesis of an additional copy of the protein (or degradation of one copy, if we had not set the rate to zero) or addition of a volume element ΔV to the volume. Second, one of these events is selected stochastically with probabilities proportional to the respective rates. After the selection, the event is carried out, the time is updated to the time of the chosen event and the rates are re-calculated based on the current protein numbers and cell volume. In the following we use two different descriptions of cell volume growth, namely linear and exponential growth, two models frequently used (Cookson et al., 2009; Cooper, 2006). It is difficult to distinguish different functional forms of volume growth, thus this question has a long history of controversial discussion, see for example the recent discussion of this issue in Cooper (2006). We note, however, that new experimental techniques support exponential growth in E. coli (Mir et al., 2011; Godin, 2010). We denote the two different growth functions by the subscripts ‘lin’ and ‘exp’. The cell volume at the beginning of each cycle is called Vðt ¼ 0Þ V 0 . In a steady state of growth, the volume doubles over the course of the cell cycle time T, such that Vðt ¼ TÞ ¼ 2V 0 . Using this constraint, the cell volume within the cell cycle is given by V lin ðtÞ ¼ V 0 ð1 þ t=TÞ;
and
2.2. Gene replication and cell division In the algorithm, the stochastic reaction events ‘protein synthesis’ and ‘volume growth’ are interrupted by two deterministic processes, (i) cell division and (ii) gene replication. These are implemented as follows. Whenever the time of the next stochastic event exceeds the time of the next deterministic event, the deterministic event is executed, the time is updated and a new time for the next stochastic event is drawn. Cell division occurs at integer multiples of the doubling time T (setting time zero to a cell division time). At the division time, we divide the protein content Pðt ¼ TÞ by two and initiate a new cell cycle for one of the daughter cells. The protein content is divided either stochastically or deterministically (Fig. 1). When a cell divides deterministically, each daughter cell gets exactly half of the protein content available. For odd numbers of proteins, we take the number after division to be either ðPðt ¼ TÞ þ1Þ=2 or ðPðt ¼ TÞ 1Þ=2, each with probability of 1/2, thus still having a minimal remnant of stochasticity. When the protein content is stochastically partitioned, each protein has a probability of 1/2 to be distributed to one of the two daughter cells at the time of cell
ð4Þ
for exponential growth at time t. In the stochastic simulations, we implement volume growth by the stochastic addition of small volume elements ΔV to the cell volume. This addition occurs with rate γ and the sizes of the volume elements are given by V0 Tγ
cycle, while for the exponential volume growth, the ΔV exp depend explicitly on the instantaneous volume V(t). In both cases, the volume growth rate γ effectively defines stochastic discrete time intervals Δt with average 〈Δt〉 ¼ 1=γ in which the volume increases on an average by ΔV lin or ΔV exp . If γ is large, the discrete time intervals are short, and we arrive at a stochastic implementation of deterministic volume growth. For smaller γ, volume growth is a stochastic process that contributes to the stochasticity of gene expression.
ð3Þ
for linear growth and by ln 2t V exp ðtÞ ¼ V 0 exp T
ΔV lin ¼
3
ln 2 Tγ
ΔV exp ¼ V ðtÞ
ð5Þ
for linear and exponential growth, respectively. Note that the size of the linear volume elements ΔV lin is constant within the cell
Fig. 1. Stochastic gene expression in growing cells: at the beginning of a cell division cycle (t¼ 0), the cell has an initial volume V ðt ¼ 0Þ ¼ V 0 and protein content Pðt ¼ 0Þ, shown in red. The volume growth is described by the stochastic addition of volume elements ΔV and reaches (on an average) a final volume 2V0 at the end of the cycle. At the time of cell division T, the protein content is divided into either a stochastic or deterministic way. If the cell divides stochastically, there is a probability of 1/2 for each protein to go to one of the two daughter cells, while if the cell divides deterministically, the protein content is divided into two equal amounts. We simulate only one lineage of cells, so the cells shaded in grey are not followed in our simulations.
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
3. Results
Stochastic division and linear growth α 0 = 1 0 m in -1 (N o G e n e R e p lica tio n )
2500
α 0 = 1 0 m in -1 (G e n e R e p lica tio n ) α eff = α 0 (2 -t x /T )= 1 5 m in -1
Protein Number P
division. Note that in this work, we consider one lineage of cells, as indicated by the gray cells in Fig. 1, such that the amount of proteins refers to a single cell. Likewise, the volume, which has on an average grown from V0 to 2V0 over the cell cycle, is divided by 2 at the time of division. We note that several different descriptions of cell division are possible, in particular if volume growth is stochastic (small γ). Throughout this work we assume that division is triggered by a timer and occurs at a fixed time T after the previous division. Alternatively, division could also be triggered by cell size (or by a combination of time and size, as discussed extensively in the literature, although mostly for eukaryotic cells with distinct cell cycle phases, Smith and Martin, 1973; Tyson and Hannsgen, 1985; Nurse, 2000, for the bacterial case see, e.g. refs. Donachie and Blakely, 2003; Chien et al., 2012). In the former case, the size at division will be a fluctuating quantity whereas in the latter, the doubling time fluctuates. Both cases can be implemented in our approach and are briefly considered in the appendix. To describe gene replication, we define a replication time tx, with 0 r t x r T. For times that are smaller than tx, the gene copy number is g ¼1. When the algorithm arrives at a time t 4 t x , gene replication takes place, and we set the gene copy number to g ¼2. At times of cell division the gene copy number is set back to g ¼1. This means that the protein synthesis rate increases twofold at some time tx within the cell cycle, since from this time on two genes are transcribed simultaneously.
G ene replicates
2000
1500
1000
2000
p (molecules/μm3)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
2100
2200
2300
2400
2500
800
G e n e re p lic a te s
600
3.1. Unregulated gene expression To begin with, we run simulations for unregulated gene expression, i.e. for a constant protein synthesis rate independent of any regulator concentration. For these simulations, we consider stochastic protein partitioning during cell division and linear volume growth as given by Eq. (4), with a rather large volume growth rate γ ¼ 200 min 1. For all the simulations, the protein concentration takes 10 cell division cycles to reach its steady state. To ensure that the process is equilibrated, we extract data after a few dozen cycles. The protein copy number P(t) as well as its concentration p(t) display oscillatory behaviour due to the explicit description of the division cycle (shown in Fig. 2), as observed previously with deterministic approaches (Marathe et al., 2012a). We consider three different scenarios of protein synthesis. In the first case, the protein content increases at a constant rate α0 and gene replication is not described explicitly. Thus, the gene copy number is taken as constant, g¼ 1. This case corresponds to a gene that is replicated immediately before cell division. A constant g can also be used as a coarse-grained averaged description, as discussed below. Here, the subscript ‘zero’ indicates that protein synthesis is unregulated. Next, we consider a gene that synthesizes proteins with rate α0 and is replicated at the middle of each cell cycle, at times mT þ t x for all integers m and with t x ¼ T=2. Finally, we consider the effective synthesis rate αeff ¼ α0 ð2 t x =TÞ ¼ 1:5α0 with a constant gene copy number given by the average gene copy number for the replication time t x ¼ T=2 (Marathe et al., 2012a). In Fig. 2(A), we show the protein number P(t) for the three different descriptions. In case there is no gene replication (blue), P(t) grows uniformly with a constant slope over the whole division cycle and decreases instantaneously during division. When taking into account gene replication (red), there is a discontinuity in the slope of the protein synthesis. Before the gene replication time tx is reached, the gene copy number is given by g ¼1. For times that exceed tx, g¼ 2, i.e. two copies of the gene synthesize proteins simultaneously. Hence, the protein synthesis rate and thus the
2000
2100
2200
2300
2400
2500
Time (min) Fig. 2. Trajectories of the protein copy number (number of molecules per cell, panel A) and protein concentration (panel B) for unregulated gene expression and linear volume growth. Three different cases are shown in the blue curves, gene replication is not included in the model, in the red curves it is included explicitly (at time tx ¼30 min after each cell division) and in the black curve it is described by an effective synthesis rate αeff that corresponds to the average copy number of the gene. Notice the systematic dependence on the division cycle for the concentration in the case with explicit cell division, where for times 0 r t o t x , the volume grows faster than the protein number, but the protein number grows more rapidly after replication of the corresponding gene. No such systematic dependence is seen in the other two cases, where both the protein content and the volume growth increase on an average linearly. The parameters of the simulations are γ¼ 200 min 1, V 0 ¼ 1 μm3 and T ¼ 60 min. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
slope of P(t) increases twofold. This increase also effectively increases the synthesis rate and thus results in a larger protein amount P(t) compared to the case without gene replication. With the constant effective synthesis rate αeff (black), the protein copy number is approximately the same as in the case with explicit gene replication, but the effective protein synthesis is uniform within the cell cycle. The corresponding protein concentration p(t) for the three cases is shown in Fig. 2(B). For the cases without explicit gene replication, i.e. in the first and the third case (blue and black, respectively), the protein concentration fluctuates around a timeindependent average, i.e. the concentration is constant apart from fluctuations. Here, the corresponding rates of protein synthesis α0 and αeff are one order of magnitude smaller than the volume growth rate γ, and the protein synthesis is the main source of fluctuations (see below). In the case where the gene is replicated at time tx (red), the protein concentration varies systematically within the division cycle. Before gene replication, the volume grows faster than protein synthesis, and the protein concentration decays. After the gene is replicated, protein synthesis becomes
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 Q4111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
faster than the volume growth, and the protein concentration increases until it reaches the end of the cell cycle. The average protein concentration is equal at the beginning and the end of the cell cycle pðt ¼ 0Þ ¼ pðt ¼ TÞ, because both volume and protein copy number double (on an average) during the division cycle. It is worth noting that the constant concentration obtained with the effective synthesis rate αeff is different from the average concentration in the more detailed model with explicit gene replication. αeff is chosen such that the time dependence of the protein copy number in the more detailed model is approximately reproduced. The concentration obtained with αeff agrees with the concentration in the detailed model only at the time of cell division. Because the concentration is reduced throughout the cycle in the detailed model, the effective description via αeff leads to a concentration that is systematically larger than the average over the division cycle. Next we consider the impact of the volume growth rate γ on both the cell volume and the protein concentration. It determines the size of the volume element that is added to the volume V(t) whenever the event ‘volume growth’ occurs. For all values of γ, the volume is doubled (on an average) in every cycle, but for smaller values of γ, this is done in larger increments, and thus, fluctuations of the volume are more pronounced. This is shown in Fig. 3(A) and (B), where we plot the time dependence of the volume for three different values of γ. For the largest value, volume growth is effectively deterministic, but for the two smaller values of γ, the
γ =0.5min-1 γ =2min-1 γ =200min-1
Cell volume (μm3)
2.5
cell volume visibly increases in stochastic steps that add a relatively large volume elements to the cell volume. Since the protein concentration pðtÞ ¼ PðtÞ=VðtÞ depends effectively on the volume, fluctuations in the volume contribute to the fluctuations of the protein concentration, as shown in Fig. 3(C) and (D). To quantify the impact of the stochasticity of volume growth on the protein concentration, we systematically vary the volume growth rate in our simulations. We compare four cases, with linear and exponential volume growth, as well as stochastic and deterministic partitioning of the proteins during cell division. We first check that the variance of the protein copy number is independent of the type (linear/exponential) and of the rate of volume growth (Fig. 4A). This is expected, since in our approach, expression of an unregulated gene is independent of volume growth. Thus, the distribution of the copy number of its protein product is the same for all values of γ. For genes controlled by a transcription factor, protein synthesis is coupled to the volume growth process (via the transcription factor concentration) and thus the protein copy number distribution may exhibit a dependence on γ. Next we analyze the protein concentration, which is always dependent on the volume, so here the value of γ matters even for unregulated genes. Fig. 4(B) shows the noise strength of the protein concentration η2p ¼ δp2 =〈p〉2 (Swain et al., 2002; Marathe et al., 2012a) as a function of the volume growth rate γ. The plotted noise strength is determined as a time average, so effectively it represents
2.5
2.0
2.0
1.5
1.5
1.0
1.0 5500
5550
5600
5650
5700
5500
5550
5600
5650
5700
70 70
p (molecules/μm3)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
5
60 60 50
50
40
40
5500
5550
5600 5650 Time (min)
5700
5500
5550
5600 5650 Time (min)
5700
Fig. 3. Volume and protein concentration as functions of time for different models: (A, B) volume dynamics for the linear (A) and the exponential (B) model. For both models, three different values of the volume growth rate γ are considered, namely γ¼ 0.5, 2 and 200 min 1 (black, red and blue curves). For small γ, the volume dynamics is stochastic, while for the largest value, it is an essentially deterministic process. (C) and (D) show the corresponding time dependence of the concentration of an unregulated protein. Averages and standard deviations of the protein concentration are (C ) 〈p〉 ¼ 50:5 7 10:3 molecules/μm3 (black), 〈p〉 ¼ 51:67 8:4 molecules/μm3 (red) and 〈p〉 ¼ 51 7 7:8 molecules/μm3 (blue). (D) 〈p〉 ¼ 48:6 7 11:1 molecules/μm3 (black), 〈p〉 ¼ 49:37 8:6 molecules/μm3 (red) and 〈p〉 ¼ 50:3 7 8:1 molecules/μm3 (blue). Note that fluctuations become stronger as γ becomes smaller. The parameters of the simulations are α0 ¼ 1 min 1 , V 0 ¼ 1 μm3 and T ¼60 min. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
Protein Number Variance δP0
Stochastic division and linear volume growth Stochastic division and exponential volume growth
680 660 640 620 600
620
580 560 540
600
520 500 0
480 0
20
1
40
2
3
4
5
60
6
80
100
Protein concentration noise strength ηp2
γ (min-1)
Deterministic division and linear volume growth Deterministic division and exponential volume growth Stochastic division and linear volume growth Stochastic division and exponential volume growth
0.01
1E-3
1
10
100
γ (min-1)
Protein concentration noise strength η2p
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
1
-1
γ = 0.5 min -1 γ = 50 min
1/
0.1
0.01
3.2. Simulation time
1E-3
100
latter observation can be understood as a result of two independent sources of noise, which are stochastic synthesis of the protein and stochastic partitioning. Their contributions to the total noise are additive (Marathe et al., 2012b); in the case of deterministic partitioning, one of these sources of noise is absent, resulting in smaller noise. The difference between linear and exponential volume growth is quite small for large γ, consistent with our previous result that the deterministic contribution to the noise due to the concentration variation over the division cycle is small (Marathe et al., 2012a). For small values of γ, the volume elements added for cell growth are relatively large, and a small number of volume elements are added during a division cycle, so the volume fluctuates and contributes to the noise strength of the protein concentration. Thus, the stochastic implementation of volume growth contributes to a third source of noise, which for small γ becomes dominant. In this regime, the difference between different modes of protein partitioning is negligible. The noise generated by volume growth is due to our implementation of volume growth and should not be taken as a correct description of real stochastic processes in volume growth, thus to suppress such additional noise, the value of γ has to be chosen sufficiently large. In Fig. 4(C), we plot the protein concentration noise strength as a function of the mean protein content. For small values of 〈P〉, the noise displays the expected behaviour, η2p ¼ 1=〈P〉 (remember that we did not include bursts) (Marathe et al., 2012a). For large 〈P〉, the noise plateaus, independent of 〈P〉, but dependent on the volume growth rate γ. This behaviour is reminiscent of the transition from intrinsic to extrinsic noise as the dominant noise source seen experimentally for increasing protein abundance (Taniguchi et al., 2010; Tsuru et al., 2009), indicating that the stochasticity of volume growth in our model presents a source of extrinsic noise. However, comparison of the numerical values indicates that this contribution to the extrinsic noise is relatively small. Even for γ ¼0.5 min 1, for which the volume increases in rather distinct discrete steps that are likely unrealistic, we only obtain extrinsic noise of η2 C 0:01, which is one order of magnitude smaller than the observed extrinsic noise of ≳0:1 (Taniguchi et al., 2010; Tsuru et al., 2009).
101
102
103
104
105
Fig. 4. Effect of volume growth on protein noise. (A) The protein number variance δP 0 is independent of the volume growth rate γ. The level of δP 0 agrees with the results obtained in our previous study (Marathe et al., 2012a), where both synthesis of proteins and protein partitioning are stochastic processes. (B) The noise strength η2p of the protein concentration decays as a function of γ. (C) η2p decays as a function of the mean protein number 〈P〉. The parameters of the simulations are V 0 ¼ 1 μm3 , T ¼60 min, α0 ¼ 10 min 1 . Note that in (A), averages are taken immediately after cell division (indicated by the subscript ‘zero’), while in (B) and (C) we average over the division cycle.
averages both over different realizations of the stochastic processes of protein synthesis and partitioning (different trajectories) and over the systematic dependence during the division cycle. Here, we consider the four different combinations of volume growth and partitioning. The noise strength η2p is seen to be affected by the stochasticity of the volume growth. Only for large γ the noise strength is independent of γ. Here, η2p is different for the deterministic and stochastic partitioning. The
Simulating volume growth requires additional simulation time. To quantify this requirement, we measured the simulation time of our approach with stochastic volume growth, while varying of the volume growth rate γ and compare it with simulations with a constant volume and with discretized time. For these simulations, we used a negatively autoregulated gene rather than the unregulated gene considered so far, so that the protein synthesis rate is actually dependent on the volume. In this case, the protein synthesis rate is given by
αðpÞ ¼
α0 1 þ ðp=KÞ2
;
ð6Þ
which depends on the volume via the concentration p ¼ P/V and where we have used a Hill coefficient of 2. In the simulations with discrete time, the cell cycle is divided into intervals Δt in which the events (volume growth or synthesis of a protein) take place with probability r Δt, where r is the current rate of the event. The fixed discretization interval is chosen as Δt ¼ 1=γ for comparison with our method. Simulations without volume growth are done using the Gillespie algorithm, where cell division is included implicitly through an effective protein degradation rate β ¼ ln 2=T. The constant volume is set to the average volume of the simulations with linear volume growth, V~ ¼ 3V 0 =2. Fig. 5(A) shows that the trajectories of the protein concentration
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
are very similar in all three cases, fluctuating around 47 molecules/ μm3. Fig. 5(B) shows the dependence of the simulation time on the volume growth rate γ. The dashed line indicates the corresponding simulation time for the constant-volume case. Absolute times were determined for 104 division cycles. For small γ, all three methods perform approximately equally. If γ is increased, the simulations for both discrete time and stochastic volume growth become slower than the simulations with constant volume. The increase in simulation time with increasing γ is similar for the two methods. Overall, our stochastic growth method is slightly faster than the discrete time method (in addition, it has the advantage of being exact [through the use of the Gillespie algorithm] for all reactions that do not depend on the volume). For almost deterministic volume growth with γ ¼200 min 1 (as used in the simulations presented above), the simulation time is increased
70
Discrete time Stochastic volume growth Constant volume
p (molecules/μm3)
60
50
40
4000
4200
4400
4600
4800
5000
Time (min) 1000
Discrete time Stochastic volume growth Constant volume
100 Time (sec)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
10
1
0.1
0.1
1
10
100
1000
10000
γ (min-1) Fig. 5. Comparison of simulation times for different methods. The steady state of a negative autoregulation system was simulated with our method (red), a simple time discretization approach with time interval Δt ¼ 1=γ (black) and a constant volume model with implicit cell division (effective degradation rate β ¼ ln 2=T, blue). (A) Trajectories of the protein concentration. The average protein concentration is the same in all three cases: 〈p〉 ¼ 47:27 6:3 molecules/μm3 (black), 〈p〉 ¼ 48:47 6:1 molecules/μm3 (red) and 〈p〉 ¼ 46:6 7 4:8 molecules/μm3 (blue). (B) Simulation times for 104 division cycles as a function of the volume growth rate γ. The blue dashed line indicates the case without volume growth. The 1 parameters used are V 0 ¼ 1 μm3 , V~ ¼ 3V 0 =2, T ¼ 60 min, α0 ¼ 1 min , K¼ 100 molecules/μm3, and γ¼ 200 min 1 (in A). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
7
approximately fourfold. If the simulated system is more complex, the relative increase can be expected to be smaller, because a smaller fraction of the simulation time is used for the simulation of volume growth.
3.3. Protein synthesis with positive autoregulation As an example of a regulated gene, we apply our approach to a simple positively autoregulated gene. In this case the protein synthesis rate depends on the concentration of its product, α ¼ αðpðtÞÞ. The synthesis rate is described by an activating Hill function AðpÞ, such that the synthesis rate of the regulated protein is given by 1 p n þ αðpÞ α0 AðpÞ ¼ α0 f pK n : 1þ K
ð7Þ
Here, K is a characteristic concentration for regulation given by the dissociation constant for the binding of the protein to its corresponding operator site (Bintu et al., 2005), i.e. the binding site upstream of the promoter. f is the maximal fold-activation of the system (in the following taken to be f ¼100) and the Hill coefficient n is a measure of the cooperativity and the sensitivity of repression. The impact of autoregulation onto protein synthesis is as follows: for a concentration p that is small compared to the concentration scale K, the Hill function is approximately constant with A 1=f , and the total protein synthesis rate αðtÞ is low with α α0 =f . This rate is called the basal synthesis rate. If the concentration is of the order of K, protein synthesis is activated by its own product. For concentrations that are large compared to K, the Hill function becomes equal to one and proteins are expressed constitutively, with the maximal rate α0. The extent to which the synthesis of a gene is regulated therefore depends on the ratio between the protein expression level (controlled by α0 T=V 0 ) and the characteristic concentration K. Depending on the level of cooperativity, determined by the Hill coefficient, qualitatively different behaviours can occur (Thomas, 1984; Segel, 1993; Ferrell, 2002). For n r 1, gene expression is monostable, but for n 41, the system becomes bistable and two stable levels of the protein concentration can be found. These two levels correspond to stable steady-state solutions of the corresponding deterministic model, which are separated by an unstable steady state (Mettetal et al., 2006). In a stochastic model, transitions between these two levels occur due to fluctuations, see, e.g. Leisner et al. (2008). We set the Hill coefficient to n ¼2, so that the system becomes bistable. In the following, we consider only the realistic case of stochastic protein partitioning during cell division and assume linear volume growth. We simulate our small gene circuit with and without gene replication and vary the synthesis rate α0. Fig. 6 shows one of the hallmarks of bistability, a hysteresis loop of the protein number average 〈P 0 〉 as a function of the synthesis rate α0, for both cases (red: without explicit gene replication, blue: gene replication at tx ¼ 30 min). Here, we present hysteresis loops for two different values of the concentration scale K, i.e. KV 0 ¼ 50 (Fig. 6A) and KV 0 ¼ 1000 (Fig. 6B). Hysteresis reflects the fact that which of the steady states (high or low state) is approached depends on the initial conditions. In these simulations, we start with a single copy of the protein, i.e. P ¼1 and a synthesis rate α0 ¼ 0:1 min 1 , sufficiently small to ensure that the low state is obtained. α0 is increased in small increments. At some value of α0 the system makes a transition to the high state. At a sufficiently large value of α0, we start to decrease α0 again to observe the jump back to the low state.
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
- Without Gene Replication - With Gene Replication
300
Protein Number Average
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
8
8000
250 200
6000
150 Jump to low state
4000
100
Jump to high state
2000
50 0
0 0
1
2
3
4
5
α0 (min-1)
0
20
40
60
80
100
α0 (min-1)
Fig. 6. Hysteresis loops of the average protein number as a function of the synthesis rate α0 with (blue) and without (red) gene replication. When increasing α0, a jump from the low to the high state of 〈P 0 〉 is observed and for decreasing α0, a reverse jump at a lower value of α0. The parameters used are γ ¼200 min 1, V 0 ¼ 1 μm3 , f ¼100, T ¼ 60 min, tx ¼ 30 min, and n ¼2. The regulation threshold K is chosen to be K¼ 50 molecules/μm3 in panel A and K ¼1000 molecules/μm3 in panel B. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
In the case without gene replication, the protein concentration is low for values of α0 smaller than 2.5 for KV 0 ¼ 50. Once the protein synthesis α0 reaches values in a range of 2:5 o α0 r 2:8 min 1 , a transition from the low state to the high state occurs. For values of the protein synthesis rate α0 4 3, 〈P 0 〉 increases linearly. In the same way, when decreasing the values of α0, a transition from the high to the low protein concentration state occurs in a range 1:7 r α0 o 2 min 1 , as shown in Fig. 6(A) (red symbols). In the case with gene replication, the transition from the low to the high average protein number occurs for smaller values of α0 in the range of 1:7 min 1 o α0 r 1:8. The transition from the high to 1 the low state takes place in the range of 1:3 r α0 o 1:5 min for the lower value of K. Note that for this case, the hysteresis loop is shifted to the left, i.e. the values of α0, at which the transitions occur, are smaller. We now use our simulation approach to investigate the impact of stochasticity on bistability, a problem that was previously studied with a series expansion on the linear noise approximation (Scott et al., 2007). To this end, we perform hysteresis simulations for different values of the characteristic concentration K. From each resulting hysteresis loops, we then extract the values of the protein synthesis rate α0 at which the system crosses over from the low to the high concentration state and vice versa. An example with larger molecule numbers ðKV 0 ¼ 1000Þ is shown in Fig. 6(B). While the qualitative behaviour is the same as for the smaller value of K, for the larger K, the system is required to have a large protein number per cell in order to turn on autoregulation and exhibit a transition from the low to the high state. This protein number is obtained for a value of α0 80 min 1 (Fig. 6B). In addition, the region of bistability is wider, with an approximately 2.3-fold difference between the switching points compared to the case in Fig. 6(A) with an about 1.5-fold difference. In a deterministic description, the steady state concentrations, and thus the region of bistability, depend only on the ratio of the α0 T=V 0 , the characteristic concentration obtained by the maximal synthesis rate and concentration scale K of the autoregulation. We therefore rescale α0 with a factor T=KV 0 and plot our results as functions of the dimensionless concentration ratio α0 T=KV 0 . In the stochastic system, however, the absolute values of α0 and K matter: even if their ratio is fixed and thus regulation is set to a certain regime, smaller values of α0 and K lead to smaller protein copy numbers and thus to increased fluctuations. This point is
rather obvious if the concentration ratio is interpreted as a ratio of the protein copy numbers αT and the characteristic copy number of regulation KV, then a decrease of α and K resets the number of molecules needed for autoregulation. In Fig. 7(A) and (B), we show the bistable region as a function of the concentration ratio α0 T=ðKV 0 Þ and the characteristic regulator copy number for regulation KV0. Panel (A) shows the case without gene replication, panel (B) the case with gene replication in the middle of the division cycle. As in Fig. 6, taking gene replication into account shifts the bistable region to smaller values of α0 T=ðKV 0 Þ. This shift reflects the increase of the effective synthesis rate for a given α0 compared to the case without gene replication due to the larger average copy number of the gene. Other than that, however, the two cases show very similar results. Most importantly, we also see that the bistable region is reduced as the relevant protein copy number gets smaller. While the lower boundary of the bistable region is rather independent of KV0, the upper boundary is strongly shifted towards smaller values of α0 T=KV 0 for KV 0 ≲1000, in qualitative agreement with the approximative results of Scott et al. (2007). To analyze the effects of volume growth and of low protein copy numbers separately, we studied the case without gene replication in more detail and determined the bistability region for different values of the volume growth rate γ. Fig. 7(A) shows two cases, essentially deterministic volume growth (γ ¼ 200 min 1, black data points) and stochastic volume growth with γ ¼0.5 min 1 (red data points), together with the corresponding results for the case with a constant volume (simulated as described above for the negative autoregulation system, blue data points). The results for almost deterministic volume growth and constant volume are very similar. No bistability is found for small KV0 with KV 0 ≲10 and KV 0 ≲30. The similarity between the two cases indicates that deterministic volume growth and the resulting deterministic variation of concentrations over the cell cycle have only a small impact on the bistability. The green triangles show the boundaries of the bistable region for the fully deterministic model with deterministic volume growth, obtained by numerically integrating Eq. (1) with the auto regulatory αðpÞ and explicit cell division. The stochastic results with γ ¼200 min 1 and with constant volume are in good agreement with the deterministic limit for KV 0 ≳1000. Additionally, we also present in blue the case without explicit cell division. For the large volume growth rate and for the explicit protein degradation rate β, bistability becomes
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4. Concluding remarks γ = 200 min-1 γ = 0.5 min-1 β = Ln2/T
100000
1000
Bistable
Monostable
KV0
10000
Monostable
100
10 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 α0T/KV0 100000
1000
Bistable
Monostable
KV0
10000
Monostable
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
9
100
10
1.0
1.5
2.0
2.5
3.0
3.5
α0T/KV0 Fig. 7. Bistability range mapped as a function of the activation scale for autoregulation KV 0 and the dimensionless scale α0 T=KV 0 for genes without (A) and with (B) gene replication. Three different regions of gene expression are illustrated. Two regions are monostable and correspond to the low (left) and high (right) states of gene expression. In between, the system is bistable and the low and high states coexist. By increasing the effect of noise (via a reduction of K and all concentrations, while their relative sizes remain constant), the region of bistability becomes diminished. As K becomes larger, the bistable region approaches the deterministic limit (filled triangles). The bistability region is shown for two different values of the volume growth rate γ¼ 0.5 and γ ¼200 min 1 (black and red, respectively). Blue data shows the corresponding results for a model with constant volume and with implicit cell division described by a protein degradation rate β. Volume fluctuations effectively reduce the bistability region. In order to correct for variations, the system requires a larger characteristic copy number of regulation KV. In (B), the region for a replicated gene is shifted to smaller values of α0 T=KV 0 with respect to (A). This shift arises due to the fact that for replicated genes, the effective synthesis rate αeff 4 α0 increases the regulator concentration. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
In this paper, we have proposed a simple computational approach to incorporate the growth of the cell volume in stochastic simulations of gene expression without the need to deal with complicated time dependencies of protein synthesis rates. Such time dependent rates arise naturally for regulated genes, whose protein synthesis rates depend on the concentration of transcription factors and thus indirectly on the cell volume. The same complication arises if regulation is described by a dependence of transcription rates on the occupancy of regulatory sites in the promoter region. In that case, the binding rate of the transcription factor to the operator site is dependent on the concentration of the transcription factors. In our approach, volume growth is approximated by a stochastic process that adds small volume elements to the cell volume with a rate γ. This process is easily included in stochastic simulations with the Gillespie algorithm, where it is treated on equal footing with the protein synthesis rates. When the rate γ is chosen sufficiently large, additional noise that arises from the stochastic implementation of volume growth can be neglected and volume growth is effectively a deterministic process. For smaller values of γ, the stochasticity of growth contributed to the extrinsic noise, but the contribution is smaller than the observed extrinsic noise. We have also included explicit description of cell division and gene replication as two processes that occur at fixed times. As an application of our method, we studied a positively autoregulated gene and determined how bistability is affected by noise. We modulated the noise level by changing the copy number of the molecule while keeping all relative concentrations constant. For large copy numbers of the autoregulator, the parameter range of bistability is in good agreement with the result of the corresponding deterministic model, for smaller molecule numbers, the range of bistability shrinks. Likewise, the region of bistability becomes diminished, when volume fluctuations are a significant source of noise. The method we used here is easy to implement and rather general. It can easily be applied to more complex gene circuits with many molecular species as well as to more detailed microscopic descriptions, where for example binding of transcription factors to regulatory DNA sequences is described explicitly. It should also be useful for the study of systems where gene expression affects cell growth (e.g. for toxic gene products), and in particular for cases where bistable circuit behaviour arises from feedback between cell growth and expression (Klumpp et al., 2009). In such cases, however, one has to be careful about the implementation of the effects of growth-rate on gene expression. If the current method is naively applied, growth rate will affect gene expression only through dilution of stable proteins, while consistency with experimental data in E. coli requires that the protein synthesis rates are also dependent on the growth rate (Hintsche and Klumpp, 2013). Moreover, under conditions of such bistability, the circuit dynamics is typically coupled to the population dynamics (Patra and Klumpp, 2013) and the study of individual linages may not be sufficient.
Acknowledgments undetectable for KV 0 ≲10 and KV 0 ≲15. For the lower volume growth rate (γ ¼0.5 min 1), the bistability region is smaller, and bistability disappears already for KV 0 30, indicating that due to the additional fluctuations arising from the fluctuating volume, a larger number of molecules is required to stabilize the bistability. Thus, deterministic volume growth has a small effect on bistability, while stochasticity due to either small molecule numbers or fluctuations in volume growth reduces the region of bistability.
RM would like to acknowledge the financial support from the Department of Science and Technology, India. Appendix A. Triggers of cell division In this appendix, we briefly consider different triggers of cell division. If, as studied in the main text, cell division is triggered by
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
TimeTriggering
500
VolumeTriggering
600
400
γ =200min γ = 0.5 min-1
300
V(t=T) =V0
500
200
V(t=T) =V0
400 Count
Count
-1
300 200
100
100 0
0 1.6
1.8
2.0
2.2
2.4
40
50
60
70
80
90
600 500
V(t=T) =V(t)/2
500
V(t=T) =V(t)/2
400 Count
400 Count
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
10
300 200
300 200 100
100 0
0 1.6
1.8
2.0
2.2
2.4
40
Cell volume right before division (µm3)
50
60
70
80
90
Cell division time (min)
Fig. 8. Time- and size triggered cell division: (A and C) histograms of cell volume at time of division for time-triggered division, (B and D) histograms of the cell cycle duration for size-triggered division. In all cases, data are shown for volume growth rates of γ ¼0.5 and 200 min 1 (blue and red, respectively). In (A) and (C), the cell volume is reset to V0 at cell division, in (B) and (D), it is set to half the value before division. The parameters of the simulations are V 0 ¼ 1 μm3 and T ¼ 60 min. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
a fixed time, volume stochasticity not only has an impact on the fluctuations within the cell cycle, but also on the histories of the cells. If at the end of the cell cycle a cell is relatively small with Vðt ¼ TÞ o 2V 0 , it will give rise to smaller daughter cells, while a large cell gives rise to large daughter cells. As a measure of the variation of cell volume, we determined distributions of cell size at the time of division for two values of the volume growth rate (γ ¼0.5 and 200 min 1) and two versions of setting the volume of the daughter cells, assuming linear volume growth. The initial volume of the daughter cell is either determined as half the final volume of the mother cell, Vðt ¼ T þ 0Þ ¼ V ðt ¼ T 0Þ or it is simply reset to V0. Distributions of the volume just before division are shown for all four cases in Fig. 8(A) and (C). For the small value of γ (stochastic volume growth), the distribution of the cell volume is rather broad with a mean value of Vðt ¼ TÞ ¼ 2V 0 , but with a large variance. The cell volume at the end of the cell cycle ranges from V ðt ¼ TÞ 1:5–2:5 μm3 . Note that when the volume is reset to V0 after division, the cell has no memory of its history and the distribution is narrower than when the volume is divided by two at division. For the large value of γ (deterministic volume growth), the distribution is sharply peaked at Vðt ¼ TÞ ¼ 2V 0 in both cases. In this case, the effect of stochasticity on the history of the cells is minimal. We also simulated the case in which the division is triggered by cell size, i.e. cells divide immediately after the cell volume has reached a value 4 2V 0 . In this case, the doubling time, i.e., the duration of a cell cycle or the time between subsequent divisions,
is a fluctuating quantity. Fig. 8(B) and (D) shows histograms of the doubling time, again for two different volume growth rates and for the two versions of resting the volume. The difference between the two ways of treating the division of volume is small, because by construction the volume at the time of division is close to 2V0. For γ ¼0.5 min 1 the distribution of cycle times is quite broad, and the duration of the cell cycle can vary between 30 and 90 min. For deterministic volume growth, i.e. γ ¼200 min 1, the distribution of the duration of the cell cycle is very narrow with a clear mean value of 60 min.
References Balaban, N.Q., Merrin, J., Chait, R., Kowalik, L., Leibler, S., 2004. Bacterial persistence as a phenotypic switch. Science 305, 1622–1625. Berg, O.G., 1978. A model for the statistical fluctuations of protein numbers in a microbial population. J. Theor. Biol. 71 (4), 587–603. Bintu, L., Buchler, N.E., Garcia, H.G., Gerland, U., Hwa, T., Kondev, J., Phillips, R., 2005. Transcriptional regulation by numbers: models. Curr. Opin. Genetics Dev. 15 (2), 125–135. Boulineau, S., Tostevin, F., Kiviet, D.J., ten Wolde, P.R., Nghe, P., Tans, S.J., 2013. Single-cell dynamics reveals sustained growth during diauxic shifts. PLoS One 8 (e61686), 1–9. Campillo, F., Lobry, C., 2012. Effect of population size in a predator-prey model. Ecol. Model. 246, 1–10. Chien, A.-C., Hill, N.S., Levin, P.A., 2012. Cell size control in bacteria. Curr. Biol. 22, R340–R349. Cookson, N., Cookson, S., Tsimring, L., Hasty, J., 2009. Cell cycle-dependent variations in protein concentration. Nucl. Acids Res. 38, 2676–2681.
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
D. Gomez et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
Cooper, S., 2006. Distinguishing between linear and exponential cell growth during the division cycle: single-cell studies, cell-culture studies, and the object of cellcycle research. Theor. Biol. Med. Model. 3, 10. Donachie, W.D., Blakely, G.W., 2003. Coupling the initiation of chromosome replication to cell size in Escherichia coli. Curr. Opin. Microbiol. 6, 146–150. Dubnau, D., Losick, R., 2006. Bistability in bacteria. Mol. Microbiol. 61, 564–572. Elowitz, M.B., Leibler, S., 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403 (20), 335–338. Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S., 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. Ferrell , J.E., 2002. Self-perpetuating states in signal transduction: positive feedback, double-negative feedback and bistability. Curr. Opin. Chem. Biol. 6, 140–148. Gillespie, D.T., 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434. Gillespie, D.T., 1977. Exact stochastic simulation of coupled chemical reactions. J. Chem. Phys. 81 (25), 2340–2361. Godin, M., 2010. Using buoyant mass to measure the growth of single cells. Nat. Methods 7, 387–390. Hayot, F., 2008. Single Cell Experiments and Gillespie's Algorithm. Department of Neurology, Mount Sinai School of Medicine, New York, NY, 10029. Hintsche, M., Klumpp, S., 2013. Dilution and the theoretical description of growthrate dependent gene expression. J. Biol. Eng. 7, 22. Klumpp, S., Zhang, Z., Hwa, T., 2009. Growth rate-dependent global effects on gene expression in bacteria. Cell 139, 1366–1375. Leisner, M., Stingl, K., Frey, E., Maier, B., 2008. Stochastic switching to competence. Curr. Opin. Microbiol. 11, 553–559. Leisner, M., Kuhr, J.T., Raedler, J., Frey, E., Maier, B., 2009. Kinetics of genetic switching into the state of bacterial competence. Biophys. J. 96 (3), 1178–1188. Li, H., Cao, Y., Petzold, L.R., Gillespie, D.T., 2008. Algorithms and software for stochastic simulation of biochemical reacting systems. Biotechnol. Prog. 24, 56–61. Li, H., Petzold, L.R., 2005. Stochastic simulation of biochemical systems on the graphics processing unit. Bioinformatics 00, 1–5. Longo, D., Hasty, J., 2006. Dynamics of single-cell gene expression. Mol. Syst. Bio. 64, 1–10. Lu, T., Volfson, D., Tsimring, L., Hasty, J., 2004. Cellular growth and division in the Gillespie algorithm. Syst. Biol. 1, 121–128. Marathe, R., Bierbaum, V., Gomez, D., Klumpp, S., 2012a. Deterministic and stochastic descriptions of gene expression dynamics. J. Stat. Phys. 148, 607–626. Marathe, R., Gomez, D., Klumpp, S., 2012b. Sources of stochasticity in constitutive and autoregulated gene expression. Phys. Scr. T151, 014068. McAdams, H.H., Arkin, A., 1997. Stochastic mechanisms in gene expression. Proc. Natl. Acad. Sci. USA 94, 814–819. Megerle, J.A., Fritz, G., Gerland, U., Jung, K., Raedler, J.O., 2008. Timing and dynamics of single cell gene expression in the arabinose utilization system. Biophys. J. 95, 2103–2115. Mettetal, J.T., Muzzey, D., Pedraza, J.M., Ozbudak, E.M., van Oudernaarden, A., 2006. Predicting stochastic gene expression dynamics in single cells. Proc. Natl. Acad. Sci. USA 103 (19), 7304–7309. Mir, M., Wang, Z., Shen, Z., Bednarz, M., Bashir, R., Golding, I., Prasanth, S.G., Popescu, G., 2011. Optical measurement of cycle-dependent cell growth. Proc. Natl. Acad. Sci. USA 108, 13124–13129.
11
Mueller, M.J.I., Klumpp, S., Lipowsky, R., 2008. Motility states of molecular motors engaged in a stochastic tug-of-war. J. Stat. Phys. 133, 1059–1081. Munsky, B., Nauert, G., van Oudenaarden, A., 2012. Using gene expression noise to understand gene regulation. Science 336, 183–187. Nath, K., Koch, A.L., 1970. Protein degradation in Escherichia coli. J. Biol. Chem. 245, 2889–2900. Nurse, P., 2000. A long twentieth century of the cell cycle and beyond. Cell 100, 71–78. Patra, P., Klumpp, S., 2013. Population dynamics of bacterial persistence. PLoS One 8 (e62814). Paulsson, J., 2004. Summing up the noise in gene networks. Nature 427 (29), 415–418. Paulsson, J., 2005. Models of stochastic gene expression. Phys. Life Rev. 2, 157–175. Raj, A., van Oudenaarden, A., 2008. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell 135, 216–226. Raj, A., Peskin, C.S., Tranchina, D., Vargas, D.Y., Tyagi, S., 2006. Stochastic mRNA synthesis in mammalian cells. PloS Biol. 309, 1707–1719. Raser, J.M., O'Shea, E.K., 2005. Noise in gene expression: origins, consequences, and control. Science 309, 2010–2013. Rosenfeld, N., Elowitz, M.B., Alon, U., 2002. Negative autoregulation speeds the response times of transcription networks. J. Mol. Biol. 323, 785–793. Scott, M., Hwa, T., Ingalls, B., 2007. Deterministic characterization of stochastic gene circuits. Proc. Natl. Acad. Sci. USA 104, 7402–7407. Segel, I.H., 1993. Enzyme Kinetics. Wiley & Sons Simpson, M.L., Cox, C.D., Sayler, G.S., 2003. Frequency domain analysis of noise in autoregulated gene circuits. Proc. Natl. Acad. Sci. U.S.A. 100, 4551–4556. Smith, J.A., Martin, L., 1973. Do cells cycle? Proc. Natl. Acad. Sci. USA 70, 1263–1267. Smits, W.K., Kuipers, O.P., Veening, J., 2006. Phenotypic variation in bacteria: the role of feedback regulation. Nat. Rev. Microbiol. 4, 259–271. So, L., Ghosh, A., Zong, C., Sepulveda, L.A., Segev, R., Golding, I., 2011. General properties of transcriptional time series in Escherichia coli. Nat. Gen. 43, 554–562. Swain, P.S., Elowitz, M.B., Siggia, E.D., 2002. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc. Natl. Acad. Sci. USA 99 (20), 12795–12800. Taniguchi, Y., Choi, P.J., Li, G.W., Chen, H., Babu, M., Hearn, J., Emili, A., Xie, X.S., 2010. Quantifying Escherichia coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 329, 533–538. Thattai, M., van Oudenaarden, A., 2001. Intrinsic noise in gene regulatory networks. Proc. Natl. Acad. Sci. USA 98, 8614–8619. Thomas, R., 1984. Logical description, analysis, and synthesis of biological and other networks comprising feedback loops. Adv. Chem. Phys. 55, 247–282. Tsuru, S., Ichinose, J., Kashiwagi, A., Ying, B.-W., Kaneko, K., Yomo, T., 2009. Noisy cell growth rate leads to fluctuating protein concentration in bacteria. Phys. Biol. 6, 036015. Tyson, J.J., Hannsgen, K.B., 1985. The Distributions of cell size and generation time in a model of the cell cycle incorporating size control and random transitions. J. Theor. Biol. 113, 28–62. Verma, M., Rawool, S., Bhat, P.J., Venkatesh, K.V., 2006. Biological significance of autoregulation through steady state analysis of genetic networks. Biosystems 84, 39–48.
Please cite this article as: Gomez, D., et al., Modeling stochastic gene expression in growing cells. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.01.017i
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
