- Email: [email protected]
Reconstructing the single-cell-level behavior of a toggle switch from population-level measurements
FEBS Letters 582 (2008) 1067–1072
Reconstructing the single-cell-level behavior of a toggle switch from population-level measurements Hirokazu Tozakia,b,*, Tetsuya J. Kobayashic,d,*, Hiroyuki Okanob, Ryo Yamamotoe, Kazuyuki Aiharaa,e, Hidenori Kimurab,* b
a ERATO Aihara Complexity Modelling Project, JST, 45-18 Oyama, Shibuya-ku, Tokyo 151-0065, Japan Bio-Mimetic Control Research Center, The Institute of Physical and Chemical Research, RIKEN, 2271-130 Anagahora, Shimoshidami, Moriyama-ku, Nagoya, Aichi Prefecture 463-0003, Japan c Research Fellow of the Japan Society for the Promotion of Science, Japan d Center for Developmental Biology, The Institute of Physical and Chemical Research, RIKEN, 2-2-3 Minatojima Minamimachi, Chuouku, Kobe-shi, Hyougo Prefecture 650-0047, Japan e Graduate school of Information Science and Technology, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Received 15 October 2007; revised 5 February 2008; accepted 22 February 2008 Available online 4 March 2008 Edited by Paul Bertone
Abstract Single-cell-level behaviors of cells are typically inferred from ensemble measurements. However, such inferences implicitly assume a biological version of ergodicity: the percentage of cells in a state is identical to the probability to find a cell in that state. While the ergodicity does not always hold, it has been rarely tested. Here, we reveal that the ergodicity does not necessarily hold even for simple toggle switches and that apparent stabilities of the switches are due to a balance between single-cell-level biased stabilities and growth rates differences. Therefore, verification of the ergodicity and reconstructing single-cell-level behaviors are crucial for understanding intracellular systems. 2008 Federation of European Biochemical Societies. Published by Elsevier B.V. All rights reserved. Keywords: Intracellular switch; Stochastic flipping; Nonlinear master equation; Growth difference; Flow cytometry
1. Introduction Genetically clonal cells show substantial non-genetic individuality [1,2], which indicates the importance of the characterization of intracellular systems from a single-cell-level stochastic viewpoint [3–6]. To infer such single-cell behaviors, ensemble measurements by flow cytometry or image cytometry have been widely used [3,4,6] because of the technical difficulties in single-cell time-lapse measurements [5]. Despite the popularity of the ensemble measurements, inferences from the ensemble measurements implicitly assume a biological version of ergodicity: namely, the percentage of the cell population in a particular state is identical to the probability to find a single cell in that state. Only when this biological ergodicity holds, stochastic behavior of a cell can be inferred directly from a population-level measurement of cells at the stationary state. However, there are many biological sit* Corresponding authors, H. Tozaki for experiment and T.J. Kobayashi for theory. E-mail addresses: [email protected] (H. Tozaki), tetsuya@ mail.crmind.net (T.J. Kobayashi), [email protected] (H. Kimura).
uations under which this assumption no longer holds. For example, when the growth rate of each cell depends on its state, then the stochastic change of the state leads to the heterogeneity of growth rates in a population of cells. This heterogeneity in growth rate biases the ratios of cells in different states in a population. In addition, the assumption does not hold when intercellular interactions exist between cells [7,8] or underlying single-cell dynamics is complicated [9,10] (see Supplementary Information S1). Thus, the validation of the biological ergodicity is crucial when we infer single-cell behavior from the population-level observation. In addition, when the assumption does not apply to a population-level experiment, the real behavior of a single cell should be reconstructed from ensemble data with mathematical models. Despite this substantial influence of the growth-rate difference to the validity of the biological ergodicity, the growth-rate difference seems to be underestimated compared with inter-cellular interactions that attract much attention, and thus the problem of biological ergodicity has been rarely tested especially when no intercellular interaction is expected.
2. Materials and methods 2.1. Strains, plasmids, and media The E. coli strain JM2300 (CGSC strain 5002) and the toggle plasmids pTAK131, 132 and pIKE107 were used. When the growth rates of low- and high-state cells were measured, JM2300 cells carrying the pBluescript plasmid were used as control cells. The LB (Miller) medium with 100 lg ml1 ampicillin was used in all experiments. 2.2. Induction and long-term culture of toggle switches Cells were firstly cultured at 37 C for about 12 h, and then induced toggle switches to either low or high state as described previously [11]. We call the state in which GFP is expressed dominantly the high state and we call the other state the low state. After induction to high or low states, cells were cultured for 50 h at 32 C. The cells were kept in the logarithmic growth phase by dilution every 5 h. At 5 h and 50 h the cells were measured by using a Becton-Dickinson FACSCalibur. 2.3. Measuring transition between low and high states After induction to the high or low states, the high- or low-state cells were diluted to appropriate densities and cultured for several hours to obtain a logarithmic growth phase. The high or low state cells were
0014-5793/$34.00 2008 Federation of European Biochemical Societies. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.febslet.2008.02.057
1068
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
then mixed at various ratios and cultured for 5 h at 32 C except pIKE107 which was cultured at 37 C. Changes of ratios of the highand low-state cells in 5-h cultures were measured by flow cytometry. 2.4. Measuring growth rates JM2300 cells carrying pBluescript were used as control cells. After culturing for about 12 h and dilution to an appropriate density, they were cultured for several hours in order to obtain a logarithmic growth phase. The cells carrying the toggle switches, on the other hand, were, after induction to the high or low state, diluted to appropriate densities and cultured for several hours. Then the control and the high- or lowstate cells in a logarithmic growth phase were mixed and cultured for 5 h at 32 C except pIKE107 which was cultured at 37 C. Changes of ratios of the control and high- or low-state cells in 5-h cultures were measured by flow cytometry.
3. Results and discussion 3.1. Testing the biological ergodicity by a mathematical model and mixing experiment As a bench mark system for testing the biological ergodicity, we use the toggle switch [11] (Fig. 1A) because of its simplicity. In addition, it may show stochastic transition between two states [12], and has no specific intercellular interaction in contrast to natural switches. We constructed the following mathematical model describing stochastic transitions between two states (high and low) under a condition that the biological ergodicity holds (see Supplementary Information S2 for the derivation): dP l ¼ k lh P l þ k hl P h ; dt dP h ¼ k hl P h þ k lh P l : dt
ð1Þ
Here Pl and Ph are, respectively, the percentages of cells in the low and high states, and kl-h and kh-l, respectively, show the frequencies of flipping from low to high and high to low. Solving these equations analytically yields k lh k lh P h ðtÞ ¼ eðk lh þkhl Þt þ P h ð0Þ þ : k lh þ k hl k lh þ k hl ð2Þ Eq. (2) means that Ph(t) depends linearly on the initial percentage of the high-state cells, Ph(0). To test this linear relation experimentally, the toggle plasmid pTAK131 was used. We mixed low- and high-state cells in various proportions, cultured the mixed cells for 5 h, and measured the changes in the percentage of the high-state cells (Fig. 1B). The percentage of the high-state cells decreased for most of initial conditions, and the percentage of the high-state cells after 5-h incubation did not linearly depend on the initial percentage. The similar non-linear relations were also observed in other toggle switches (see Supplementary Information S7). This result is not only quantitatively but also qualitatively inconsistent with the theoretical prediction of Eq. (2), suggesting that the biological ergodicity does not hold. 3.2. Modification of model predicts the existence of growth difference This inconsistency may be due to difference in growth rates or to intercellular interactions [7,8]. The former is more probable for our system because the toggle switches are not designed to interact with each other intercellularly. To see
whether the difference of growth rates can explain the inconsistency or not, we constructed a modified mathematical model in which both the stochastic transition of the switch and the cell growth are incorporated as follows (see Supplementary Information S3 for the derivation): dP l ¼ k lh P l þ k hl P h þ ðgl < g >ÞP l ; dt dP h ¼ k hl P h þ k lh P l þ ðgh < g >ÞP h ; dt
ð3Þ
where Ægæ represents the average growth rate of the whole population defined as Ægæ = glPl + ghPh. This equation can be solved analytically as P l ðtÞ exp½ðK þ GÞtðP l ð0Þ P h ð0ÞÞT ð4Þ ; ¼ P h ðtÞ ð1 1Þ exp½ðK þ GÞtðP l ð0Þ P h ð0ÞÞT k lh k hl where K ¼ and G ¼ diagð gl gh Þ. The k lh k hl numerical simulation of this equation showed that the model can qualitatively explain the non-linear relationship between the initial percentage of cells and that after 5-h incubation (Fig. 1C). 3.3. Experimental verification of the growth-rate difference To experimentally prove and quantify the growth-rate difference predicted by the modified model, we compared the growth rate of low-state cells with that of high-state cells by using flow cytometry. We mixed high- or low-state cells with control cells, cultured them for 5 h, and observed the change of cell proportion (Fig. 2). The ratio of high- and low-state cells was almost constant for 5 h compared with the change in the ratio of control and toggle-carrying cells, indicating that this dynamics is dominated not by inter-state but by inter-cellular growth difference. Since no switching occurs between the control cells and the cells with toggle switches, the change in the percentages of high- or low-state cells and control cells can be approximately described with the following mathematical model: dpi ¼ ðgi hgiÞpi ; dt dpc ¼ ðgc hgiÞpc ; dt
ð5Þ
where c is the index of control cells, and i is either l or h depending on whether the toggle switch is in the low state or the high state. These equations can be solved analytically as h i c ð0Þ =s. ðgc gi Þ ¼ log 11=p 11=pc ðsÞ By using this equation, the growth-rate difference of the control cells and the cells with toggle switches can be estimated from the experimentally observed percentages of the control cells at t = 0 and t = s. The estimated values were gc gh = 0.069 h1 and gc gl = 0.066 h1. As the difference of these values, we have gl gh = 0.135 h1, which indicates that the low-state cells grow faster than the high-state cells and that the difference is less than 10% of the absolute growth rate of the low- and high-state cells. The growth differences were also observed in the other toggle switches, and verified by optical density (O.D.) (see Supplementary Informations S7 and S5 for more detail). Because the high state has lower growth rate than the other for all plasmid tested, we think that
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
1069
Fig. 1. Toggle switch and its transient behaviors at the population level. (A) The scheme of the toggle switch pTAK131. The pTAK131 plasmid includes the Lac repressor (lacI) and the Ptrc-2 promoter as one promoter-repressor pair and includes the temperature-sensitive lambda repressor (CIts) and the PLs1con promoter as the other pair [11]. (B) Population-level dynamics of the toggle switch. Ph(t) at t = 5 h as a function of Ph(0). The individual points represent each experimental data, and blue curves correspond to the numerical solution of Eq. (4). The blue curve was obtained with klh = 103 h1, khl = 108 h1, and gl gh = 0.135 h1. The distributions of the low- and high-state cells at t = 0 and t = 5 h are shown for individual experimental data. The dashed diagonal line representing Ph(t) = Ph(0) is shown for readability. The reproducibility and generality of this result was tested in Supplementary Information S6 and S7. (C) Simulated population-level dynamics of the toggle switch with growth-rate difference. Ph(t) at t = 5 h as a function of Ph(0) was calculated with Eq. (3) for different values of growth-rate difference, Dg = gl gh, and a fixed value of the flipping frequencies. The red and blue curves respectively represent the numerical simulations for negative and positive values of the growth-rate difference, which ranges from 0.5 h1 to 0.5 h1.
the weak toxicity of GFP is the most plausible origin of the growth difference. While the existence of growth rate difference was verified, it is still unknown whether this very small growthrate difference can induce the qualitative discrepancy between the theoretical prediction (Eq. (2)) and the experimental data (Fig. 1B).
3.4. Quantitative agreement between the experiment and the mathematical model To further verify whether the experimentally observed growth-rate difference can consistently and quantitatively explain the mixing experiment or not, we re-calculated Eq. (3) by assigning the experimentally estimated growth-rate. Then
1070
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
Fig. 2. Growth rates of low- and high-state cells (A and B) The percentages of the low-state cells and control cells (cells carrying pBluescript) when mixed (A) and after being cultured for 5 h (B). (C and D) The percentages of the high-state cells and control cells (cells carrying pBluescript) when mixed (C) and after being cultured for 5 h (D).
we found that the relationship between Ph(t) and Ph(0) can be reproduced not only qualitatively but also quantitatively with Eq. (3) for certain values of kl-h and kh-l as shown in Fig. 1B with a blue curve. We also found that the agreement is good enough when kl-h and kh-l are less than 102 h1 (see Fig. S1 in the Supplementary Information). The similar results were also obtained for the other toggle switches tested (see Supplementary Information S7). This quantitative consistency between the mathematical model and the experimental data indicates that the inconsistency between the simplest model of Eqs. (1) and (2) and the experimental data can be attributed to growth-rate difference rather than unknown intercellular interaction between switches. 3.5. Estimation of switching parameters Even though growth-rate difference exists, the biological ergodicity may hold. Identification of the other parameters, flipping frequencies of the switch and subsequent reconstruction of the single-cell behavior are indispensable for validating the biological ergodicity. However, the quantification of the very small flipping frequencies by direct time-lapse measurement is technically difficult because the estimated upper bound of the flipping frequencies, 102 h1, indicates that we need to
monitor quickly growing single-cells for at least 4 days. We employed Eq. (3) and the experimental data on the long-term population dynamics of the cells to further narrow down the estimated ranges of kl-h and kh-l. When we cultured cells carrying toggle switches induced to the low state, most of the cells stayed in the low state for a long time (Fig. S3). However, a small number of the high-state cells appeared soon as shown in Fig. S3. Since the growth rate of the low-state cells is higher than that of the high-state cells, this change can be attributed not to the growth-rate difference but to the flipping of the switch. Furthermore, after a small number of high-state cells appeared, the percentage of high-state cells was kept almost constant (Fig. S3), indicating that the population is in the stationary state that corresponds to limtfi1(Pl(t),Ph(t)) in Eq. (5). From the final percentage of the high-state cells observed in Fig. S3 we estimated that klh @ 103 h1 (see Supplementary Informations S8 and Fig. S4). On the other hand, we confirmed that toggle-carrying cells being induced completely to the high state can sustain the high state for more than 50 h (Fig. S5). It is difficult to distinguish whether the breakdown of the high state is attributed to a stochastic flipping of the switch or to contamination of the low-state cells that escaped from been inducted to high state.
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
1071
Fig. 3. Reconstructed behavior of the toggle switch (A) The numerically calculated Ph(t) for different initial conditions. The red and blue curves respectively represent the dynamics of Ph(t) in a cell population (Eq. (1)) and a single cell (Eq. (3)) where klh = 103 h1, khl = 108 h1, and gl gh = 0.135 h1. (B) Schematic diagram of the behavior of the toggle switch. Cells in the low state grow fast and switch to the high state frequently, whereas cells in the high state grow slowly and switch to the low state infrequently. (C) Schematic diagram of the behavior of the toggle switch inferred under the assumption that the biological ergodicity holds. Cells in the low state switch to the high state infrequently, whereas cells in the high state switch to the low state frequently.
Because the growth rate of the low-state cells is higher than that of the high-state cells, contamination of low-state cells shorten the duration of the high state but never lengthen it. Thus the result that the high state can be sustained for more than 50 h defines the upper bound of the switching frequency from high to low as khl < 108 h1 (see Supplementary Informations S8 and Fig. S6). 3.6. Reconstruction of the single-cell-level dynamics of the switch Using Eq. (4) and the estimated values of parameters gl, gh, kl-h, and kh-l, we theoretically reproduced in silico the behavior of the toggle switch in a cell population (Fig. 3A, red line) and the behavior of the toggle switch in a single cell (Fig. 3A, blue line). The results elucidated that the behavior of a single cell and a cell population can differ significantly. At the singlecell-level the low state is relatively less stable than the high state (Fig. 3A, blue line), while at the population level the low state is apparently more stable than the high state (Fig. 3A, red line). In the case of the pTAK131 toggle plasmid, the apparent relative stability at the population level is established by the balance between the growth-rate difference and biased relative stability of the two states at the single-cell-level (Fig. 3B and Supplementary Information S9).
4. Discussion 4.1. Reconstruction revealed the underlying single-cell-level behavior of the switch Our analysis demonstrated that the real single-cell-level behavior of an intracellular switch can be completely different
from the behavior inferred under the assumption that the biological ergodicity holds (Fig. 3C). Specifically, the growth-rate difference whose value is less than the 10% of absolute growth rate can induce the substantial discrepancy between the apparent relative stability at the population-level and the real relative stability at the single-cell-level. Therefore, a criterion to assess the influence of growth-rate difference is important to quickly check the possible discrepancy. A theoretical analysis proves that the maximum discrepancy is determined by the ratio of the growth-rate difference and the sum of the flipping frequencies (see Supplementary Information S10 for more detail). When the growth-rate difference is less than the twice of the sum of the flipping frequencies, the maximum discrepancy is limited to be one fourth of their ratio. Thus, the break of the biological ergodicity by growthrate difference is not negligible when the expected time scale of the flipping of a switch is much less than the absolute growth rate of cells. 4.2. Biological implications of biological ergodicity Since bistability at the population-level has been observed not only in toggle switches but also in many intracellular systems [3,11,13–17] the verification of the biological ergodicity and the reconstruction of single-cell-level behavior are also crucial to fully understand the molecular mechanisms and biological functions of intracellular bistability in immune cells and ES-cells. Since the mixing experiment can be conducted to these cells without serious technological difficulties, the combined use of the mixing experiments and the mathematical model can contribute to reveal the mechanism and strategy of differentiation underlying the population-level bistability.
1072 Acknowledgements: We are grateful to Dr. J.J. Collins (Boston University) for the toggle switch plasmids. We thank Drs. H. Kobayashi, R. Wang, K. Tsumoto, G. Kurosawa, and Mr. Y. Matsuura for their experimental instructions, valuable advice and critical reading. This work was supported by Research Fellowships of JSPS for Young Scientists (17-4169) and Grant-in-Aid for Scientific Research on Priority Areas – System study on higher-order brain functions – from MEXT of Japan (17022012).
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.febslet.2008. 02.057. References [1] Kaern, M., Elston, T.C., Blake, W.J. and Collins, J.J. (2005) Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet. 6, 451–464. [2] Raser, J.M. and OÕShea, E.K. (2005) Noise in gene expression: origins, consequences, and control. Science 309, 2010–2013. [3] Acar, M., Becskei, A. and van Oudenaarden, A. (2005) Enhancement of cellular memory by reducing stochastic transitions. Nature 435, 228–232. [4] Blake, W.J., Kaern, M., Cantor, C.R. and Collins, J.J. (2003) Noise in eukaryotic gene expression. Nature 422, 633–637. [5] Elowitz, M.B., Levine, A.J., Siggia, E.D. and Swain, P.S. (2002) Stochastic gene expression in a single cell. Science 297, 1183–1186. [6] Ozbudak, E.M., Thattai, M., Lim, H.N., Shraiman, B.I. and Van Oudenaarden, A. (2004) Multistability in the lactose utilization network of Escherichia coli. Nature 427, 737–740.
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072 [7] Cox, C.D. et al. (2003) Analysis of noise in quorum sensing. Omics 7, 317–334. [8] Goryachev, A.B., Toh, D.J., Wee, K.B., Lee, T., Zhang, H.B. and Zhang, L.H. (2005) Transition to quorum sensing in an Agrobacterium population: a stochastic model. PLoS Comput. Biol. 1, e37. [9] Ma, L., Wagner, J., Rice, J.J., Hu, W., Levine, A.J. and Stolovitzky, G.A. (2005) A plausible model for the digital response of p53 to DNA damage. Proc. Natl. Acad. Sci. USA 102, 14266–14271. [10] Ribeiro, A.S. (2007) Dynamics of a two-dimensional model of cell tissues with coupled stochastic gene networks. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 76, 051915. [11] Gardner, T.S., Cantor, C.R. and Collins, J.J. (2000) Construction of a genetic toggle switch in Escherichia coli. Nature 403, 339–342. [12] Morishita, Y. and Aihara, K. (2004) Noise-reduction through interaction in gene expression and biochemical reaction processes. J. Theor. Biol. 228, 315–325. [13] Becskei, A., Seraphin, B. and Serrano, L. (2001) Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion. Embo J. 20, 2528–2535. [14] Isaacs, F.J., Hasty, J., Cantor, C.R. and Collins, J.J. (2003) Prediction and measurement of an autoregulatory genetic module. Proc. Natl. Acad. Sci. USA 100, 7714–7719. [15] Kramer, B.P., Viretta, A.U., Daoud-El-Baba, M., Aubel, D., Weber, W. and Fussenegger, M. (2004) An engineered epigenetic transgene switch in mammalian cells. Nat. Biotechnol. 22, 867– 870. [16] Weinberger, L.S., Burnett, J.C., Toettcher, J.E., Arkin, A.P. and Schaffer, D.V. (2005) Stochastic gene expression in a lentiviral positive-feedback loop: HIV-1 Tat fluctuations drive phenotypic diversity. Cell 122, 169–182. [17] Maamar, H. and Dubnau, D. (2005) Bistability in the Bacillus subtilis K-state (competence) system requires a positive feedback loop. Mol. Microbiol. 56, 615–624.
Reconstructing the single-cell-level behavior of a toggle switch from population-level measurements Hirokazu Tozakia,b,*, Tetsuya J. Kobayashic,d,*, Hiroyuki Okanob, Ryo Yamamotoe, Kazuyuki Aiharaa,e, Hidenori Kimurab,* b
a ERATO Aihara Complexity Modelling Project, JST, 45-18 Oyama, Shibuya-ku, Tokyo 151-0065, Japan Bio-Mimetic Control Research Center, The Institute of Physical and Chemical Research, RIKEN, 2271-130 Anagahora, Shimoshidami, Moriyama-ku, Nagoya, Aichi Prefecture 463-0003, Japan c Research Fellow of the Japan Society for the Promotion of Science, Japan d Center for Developmental Biology, The Institute of Physical and Chemical Research, RIKEN, 2-2-3 Minatojima Minamimachi, Chuouku, Kobe-shi, Hyougo Prefecture 650-0047, Japan e Graduate school of Information Science and Technology, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Received 15 October 2007; revised 5 February 2008; accepted 22 February 2008 Available online 4 March 2008 Edited by Paul Bertone
Abstract Single-cell-level behaviors of cells are typically inferred from ensemble measurements. However, such inferences implicitly assume a biological version of ergodicity: the percentage of cells in a state is identical to the probability to find a cell in that state. While the ergodicity does not always hold, it has been rarely tested. Here, we reveal that the ergodicity does not necessarily hold even for simple toggle switches and that apparent stabilities of the switches are due to a balance between single-cell-level biased stabilities and growth rates differences. Therefore, verification of the ergodicity and reconstructing single-cell-level behaviors are crucial for understanding intracellular systems. 2008 Federation of European Biochemical Societies. Published by Elsevier B.V. All rights reserved. Keywords: Intracellular switch; Stochastic flipping; Nonlinear master equation; Growth difference; Flow cytometry
1. Introduction Genetically clonal cells show substantial non-genetic individuality [1,2], which indicates the importance of the characterization of intracellular systems from a single-cell-level stochastic viewpoint [3–6]. To infer such single-cell behaviors, ensemble measurements by flow cytometry or image cytometry have been widely used [3,4,6] because of the technical difficulties in single-cell time-lapse measurements [5]. Despite the popularity of the ensemble measurements, inferences from the ensemble measurements implicitly assume a biological version of ergodicity: namely, the percentage of the cell population in a particular state is identical to the probability to find a single cell in that state. Only when this biological ergodicity holds, stochastic behavior of a cell can be inferred directly from a population-level measurement of cells at the stationary state. However, there are many biological sit* Corresponding authors, H. Tozaki for experiment and T.J. Kobayashi for theory. E-mail addresses: [email protected] (H. Tozaki), tetsuya@ mail.crmind.net (T.J. Kobayashi), [email protected] (H. Kimura).
uations under which this assumption no longer holds. For example, when the growth rate of each cell depends on its state, then the stochastic change of the state leads to the heterogeneity of growth rates in a population of cells. This heterogeneity in growth rate biases the ratios of cells in different states in a population. In addition, the assumption does not hold when intercellular interactions exist between cells [7,8] or underlying single-cell dynamics is complicated [9,10] (see Supplementary Information S1). Thus, the validation of the biological ergodicity is crucial when we infer single-cell behavior from the population-level observation. In addition, when the assumption does not apply to a population-level experiment, the real behavior of a single cell should be reconstructed from ensemble data with mathematical models. Despite this substantial influence of the growth-rate difference to the validity of the biological ergodicity, the growth-rate difference seems to be underestimated compared with inter-cellular interactions that attract much attention, and thus the problem of biological ergodicity has been rarely tested especially when no intercellular interaction is expected.
2. Materials and methods 2.1. Strains, plasmids, and media The E. coli strain JM2300 (CGSC strain 5002) and the toggle plasmids pTAK131, 132 and pIKE107 were used. When the growth rates of low- and high-state cells were measured, JM2300 cells carrying the pBluescript plasmid were used as control cells. The LB (Miller) medium with 100 lg ml1 ampicillin was used in all experiments. 2.2. Induction and long-term culture of toggle switches Cells were firstly cultured at 37 C for about 12 h, and then induced toggle switches to either low or high state as described previously [11]. We call the state in which GFP is expressed dominantly the high state and we call the other state the low state. After induction to high or low states, cells were cultured for 50 h at 32 C. The cells were kept in the logarithmic growth phase by dilution every 5 h. At 5 h and 50 h the cells were measured by using a Becton-Dickinson FACSCalibur. 2.3. Measuring transition between low and high states After induction to the high or low states, the high- or low-state cells were diluted to appropriate densities and cultured for several hours to obtain a logarithmic growth phase. The high or low state cells were
0014-5793/$34.00 2008 Federation of European Biochemical Societies. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.febslet.2008.02.057
1068
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
then mixed at various ratios and cultured for 5 h at 32 C except pIKE107 which was cultured at 37 C. Changes of ratios of the highand low-state cells in 5-h cultures were measured by flow cytometry. 2.4. Measuring growth rates JM2300 cells carrying pBluescript were used as control cells. After culturing for about 12 h and dilution to an appropriate density, they were cultured for several hours in order to obtain a logarithmic growth phase. The cells carrying the toggle switches, on the other hand, were, after induction to the high or low state, diluted to appropriate densities and cultured for several hours. Then the control and the high- or lowstate cells in a logarithmic growth phase were mixed and cultured for 5 h at 32 C except pIKE107 which was cultured at 37 C. Changes of ratios of the control and high- or low-state cells in 5-h cultures were measured by flow cytometry.
3. Results and discussion 3.1. Testing the biological ergodicity by a mathematical model and mixing experiment As a bench mark system for testing the biological ergodicity, we use the toggle switch [11] (Fig. 1A) because of its simplicity. In addition, it may show stochastic transition between two states [12], and has no specific intercellular interaction in contrast to natural switches. We constructed the following mathematical model describing stochastic transitions between two states (high and low) under a condition that the biological ergodicity holds (see Supplementary Information S2 for the derivation): dP l ¼ k lh P l þ k hl P h ; dt dP h ¼ k hl P h þ k lh P l : dt
ð1Þ
Here Pl and Ph are, respectively, the percentages of cells in the low and high states, and kl-h and kh-l, respectively, show the frequencies of flipping from low to high and high to low. Solving these equations analytically yields k lh k lh P h ðtÞ ¼ eðk lh þkhl Þt þ P h ð0Þ þ : k lh þ k hl k lh þ k hl ð2Þ Eq. (2) means that Ph(t) depends linearly on the initial percentage of the high-state cells, Ph(0). To test this linear relation experimentally, the toggle plasmid pTAK131 was used. We mixed low- and high-state cells in various proportions, cultured the mixed cells for 5 h, and measured the changes in the percentage of the high-state cells (Fig. 1B). The percentage of the high-state cells decreased for most of initial conditions, and the percentage of the high-state cells after 5-h incubation did not linearly depend on the initial percentage. The similar non-linear relations were also observed in other toggle switches (see Supplementary Information S7). This result is not only quantitatively but also qualitatively inconsistent with the theoretical prediction of Eq. (2), suggesting that the biological ergodicity does not hold. 3.2. Modification of model predicts the existence of growth difference This inconsistency may be due to difference in growth rates or to intercellular interactions [7,8]. The former is more probable for our system because the toggle switches are not designed to interact with each other intercellularly. To see
whether the difference of growth rates can explain the inconsistency or not, we constructed a modified mathematical model in which both the stochastic transition of the switch and the cell growth are incorporated as follows (see Supplementary Information S3 for the derivation): dP l ¼ k lh P l þ k hl P h þ ðgl < g >ÞP l ; dt dP h ¼ k hl P h þ k lh P l þ ðgh < g >ÞP h ; dt
ð3Þ
where Ægæ represents the average growth rate of the whole population defined as Ægæ = glPl + ghPh. This equation can be solved analytically as P l ðtÞ exp½ðK þ GÞtðP l ð0Þ P h ð0ÞÞT ð4Þ ; ¼ P h ðtÞ ð1 1Þ exp½ðK þ GÞtðP l ð0Þ P h ð0ÞÞT k lh k hl where K ¼ and G ¼ diagð gl gh Þ. The k lh k hl numerical simulation of this equation showed that the model can qualitatively explain the non-linear relationship between the initial percentage of cells and that after 5-h incubation (Fig. 1C). 3.3. Experimental verification of the growth-rate difference To experimentally prove and quantify the growth-rate difference predicted by the modified model, we compared the growth rate of low-state cells with that of high-state cells by using flow cytometry. We mixed high- or low-state cells with control cells, cultured them for 5 h, and observed the change of cell proportion (Fig. 2). The ratio of high- and low-state cells was almost constant for 5 h compared with the change in the ratio of control and toggle-carrying cells, indicating that this dynamics is dominated not by inter-state but by inter-cellular growth difference. Since no switching occurs between the control cells and the cells with toggle switches, the change in the percentages of high- or low-state cells and control cells can be approximately described with the following mathematical model: dpi ¼ ðgi hgiÞpi ; dt dpc ¼ ðgc hgiÞpc ; dt
ð5Þ
where c is the index of control cells, and i is either l or h depending on whether the toggle switch is in the low state or the high state. These equations can be solved analytically as h i c ð0Þ =s. ðgc gi Þ ¼ log 11=p 11=pc ðsÞ By using this equation, the growth-rate difference of the control cells and the cells with toggle switches can be estimated from the experimentally observed percentages of the control cells at t = 0 and t = s. The estimated values were gc gh = 0.069 h1 and gc gl = 0.066 h1. As the difference of these values, we have gl gh = 0.135 h1, which indicates that the low-state cells grow faster than the high-state cells and that the difference is less than 10% of the absolute growth rate of the low- and high-state cells. The growth differences were also observed in the other toggle switches, and verified by optical density (O.D.) (see Supplementary Informations S7 and S5 for more detail). Because the high state has lower growth rate than the other for all plasmid tested, we think that
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
1069
Fig. 1. Toggle switch and its transient behaviors at the population level. (A) The scheme of the toggle switch pTAK131. The pTAK131 plasmid includes the Lac repressor (lacI) and the Ptrc-2 promoter as one promoter-repressor pair and includes the temperature-sensitive lambda repressor (CIts) and the PLs1con promoter as the other pair [11]. (B) Population-level dynamics of the toggle switch. Ph(t) at t = 5 h as a function of Ph(0). The individual points represent each experimental data, and blue curves correspond to the numerical solution of Eq. (4). The blue curve was obtained with klh = 103 h1, khl = 108 h1, and gl gh = 0.135 h1. The distributions of the low- and high-state cells at t = 0 and t = 5 h are shown for individual experimental data. The dashed diagonal line representing Ph(t) = Ph(0) is shown for readability. The reproducibility and generality of this result was tested in Supplementary Information S6 and S7. (C) Simulated population-level dynamics of the toggle switch with growth-rate difference. Ph(t) at t = 5 h as a function of Ph(0) was calculated with Eq. (3) for different values of growth-rate difference, Dg = gl gh, and a fixed value of the flipping frequencies. The red and blue curves respectively represent the numerical simulations for negative and positive values of the growth-rate difference, which ranges from 0.5 h1 to 0.5 h1.
the weak toxicity of GFP is the most plausible origin of the growth difference. While the existence of growth rate difference was verified, it is still unknown whether this very small growthrate difference can induce the qualitative discrepancy between the theoretical prediction (Eq. (2)) and the experimental data (Fig. 1B).
3.4. Quantitative agreement between the experiment and the mathematical model To further verify whether the experimentally observed growth-rate difference can consistently and quantitatively explain the mixing experiment or not, we re-calculated Eq. (3) by assigning the experimentally estimated growth-rate. Then
1070
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
Fig. 2. Growth rates of low- and high-state cells (A and B) The percentages of the low-state cells and control cells (cells carrying pBluescript) when mixed (A) and after being cultured for 5 h (B). (C and D) The percentages of the high-state cells and control cells (cells carrying pBluescript) when mixed (C) and after being cultured for 5 h (D).
we found that the relationship between Ph(t) and Ph(0) can be reproduced not only qualitatively but also quantitatively with Eq. (3) for certain values of kl-h and kh-l as shown in Fig. 1B with a blue curve. We also found that the agreement is good enough when kl-h and kh-l are less than 102 h1 (see Fig. S1 in the Supplementary Information). The similar results were also obtained for the other toggle switches tested (see Supplementary Information S7). This quantitative consistency between the mathematical model and the experimental data indicates that the inconsistency between the simplest model of Eqs. (1) and (2) and the experimental data can be attributed to growth-rate difference rather than unknown intercellular interaction between switches. 3.5. Estimation of switching parameters Even though growth-rate difference exists, the biological ergodicity may hold. Identification of the other parameters, flipping frequencies of the switch and subsequent reconstruction of the single-cell behavior are indispensable for validating the biological ergodicity. However, the quantification of the very small flipping frequencies by direct time-lapse measurement is technically difficult because the estimated upper bound of the flipping frequencies, 102 h1, indicates that we need to
monitor quickly growing single-cells for at least 4 days. We employed Eq. (3) and the experimental data on the long-term population dynamics of the cells to further narrow down the estimated ranges of kl-h and kh-l. When we cultured cells carrying toggle switches induced to the low state, most of the cells stayed in the low state for a long time (Fig. S3). However, a small number of the high-state cells appeared soon as shown in Fig. S3. Since the growth rate of the low-state cells is higher than that of the high-state cells, this change can be attributed not to the growth-rate difference but to the flipping of the switch. Furthermore, after a small number of high-state cells appeared, the percentage of high-state cells was kept almost constant (Fig. S3), indicating that the population is in the stationary state that corresponds to limtfi1(Pl(t),Ph(t)) in Eq. (5). From the final percentage of the high-state cells observed in Fig. S3 we estimated that klh @ 103 h1 (see Supplementary Informations S8 and Fig. S4). On the other hand, we confirmed that toggle-carrying cells being induced completely to the high state can sustain the high state for more than 50 h (Fig. S5). It is difficult to distinguish whether the breakdown of the high state is attributed to a stochastic flipping of the switch or to contamination of the low-state cells that escaped from been inducted to high state.
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072
1071
Fig. 3. Reconstructed behavior of the toggle switch (A) The numerically calculated Ph(t) for different initial conditions. The red and blue curves respectively represent the dynamics of Ph(t) in a cell population (Eq. (1)) and a single cell (Eq. (3)) where klh = 103 h1, khl = 108 h1, and gl gh = 0.135 h1. (B) Schematic diagram of the behavior of the toggle switch. Cells in the low state grow fast and switch to the high state frequently, whereas cells in the high state grow slowly and switch to the low state infrequently. (C) Schematic diagram of the behavior of the toggle switch inferred under the assumption that the biological ergodicity holds. Cells in the low state switch to the high state infrequently, whereas cells in the high state switch to the low state frequently.
Because the growth rate of the low-state cells is higher than that of the high-state cells, contamination of low-state cells shorten the duration of the high state but never lengthen it. Thus the result that the high state can be sustained for more than 50 h defines the upper bound of the switching frequency from high to low as khl < 108 h1 (see Supplementary Informations S8 and Fig. S6). 3.6. Reconstruction of the single-cell-level dynamics of the switch Using Eq. (4) and the estimated values of parameters gl, gh, kl-h, and kh-l, we theoretically reproduced in silico the behavior of the toggle switch in a cell population (Fig. 3A, red line) and the behavior of the toggle switch in a single cell (Fig. 3A, blue line). The results elucidated that the behavior of a single cell and a cell population can differ significantly. At the singlecell-level the low state is relatively less stable than the high state (Fig. 3A, blue line), while at the population level the low state is apparently more stable than the high state (Fig. 3A, red line). In the case of the pTAK131 toggle plasmid, the apparent relative stability at the population level is established by the balance between the growth-rate difference and biased relative stability of the two states at the single-cell-level (Fig. 3B and Supplementary Information S9).
4. Discussion 4.1. Reconstruction revealed the underlying single-cell-level behavior of the switch Our analysis demonstrated that the real single-cell-level behavior of an intracellular switch can be completely different
from the behavior inferred under the assumption that the biological ergodicity holds (Fig. 3C). Specifically, the growth-rate difference whose value is less than the 10% of absolute growth rate can induce the substantial discrepancy between the apparent relative stability at the population-level and the real relative stability at the single-cell-level. Therefore, a criterion to assess the influence of growth-rate difference is important to quickly check the possible discrepancy. A theoretical analysis proves that the maximum discrepancy is determined by the ratio of the growth-rate difference and the sum of the flipping frequencies (see Supplementary Information S10 for more detail). When the growth-rate difference is less than the twice of the sum of the flipping frequencies, the maximum discrepancy is limited to be one fourth of their ratio. Thus, the break of the biological ergodicity by growthrate difference is not negligible when the expected time scale of the flipping of a switch is much less than the absolute growth rate of cells. 4.2. Biological implications of biological ergodicity Since bistability at the population-level has been observed not only in toggle switches but also in many intracellular systems [3,11,13–17] the verification of the biological ergodicity and the reconstruction of single-cell-level behavior are also crucial to fully understand the molecular mechanisms and biological functions of intracellular bistability in immune cells and ES-cells. Since the mixing experiment can be conducted to these cells without serious technological difficulties, the combined use of the mixing experiments and the mathematical model can contribute to reveal the mechanism and strategy of differentiation underlying the population-level bistability.
1072 Acknowledgements: We are grateful to Dr. J.J. Collins (Boston University) for the toggle switch plasmids. We thank Drs. H. Kobayashi, R. Wang, K. Tsumoto, G. Kurosawa, and Mr. Y. Matsuura for their experimental instructions, valuable advice and critical reading. This work was supported by Research Fellowships of JSPS for Young Scientists (17-4169) and Grant-in-Aid for Scientific Research on Priority Areas – System study on higher-order brain functions – from MEXT of Japan (17022012).
Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.febslet.2008. 02.057. References [1] Kaern, M., Elston, T.C., Blake, W.J. and Collins, J.J. (2005) Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet. 6, 451–464. [2] Raser, J.M. and OÕShea, E.K. (2005) Noise in gene expression: origins, consequences, and control. Science 309, 2010–2013. [3] Acar, M., Becskei, A. and van Oudenaarden, A. (2005) Enhancement of cellular memory by reducing stochastic transitions. Nature 435, 228–232. [4] Blake, W.J., Kaern, M., Cantor, C.R. and Collins, J.J. (2003) Noise in eukaryotic gene expression. Nature 422, 633–637. [5] Elowitz, M.B., Levine, A.J., Siggia, E.D. and Swain, P.S. (2002) Stochastic gene expression in a single cell. Science 297, 1183–1186. [6] Ozbudak, E.M., Thattai, M., Lim, H.N., Shraiman, B.I. and Van Oudenaarden, A. (2004) Multistability in the lactose utilization network of Escherichia coli. Nature 427, 737–740.
H. Tozaki et al. / FEBS Letters 582 (2008) 1067–1072 [7] Cox, C.D. et al. (2003) Analysis of noise in quorum sensing. Omics 7, 317–334. [8] Goryachev, A.B., Toh, D.J., Wee, K.B., Lee, T., Zhang, H.B. and Zhang, L.H. (2005) Transition to quorum sensing in an Agrobacterium population: a stochastic model. PLoS Comput. Biol. 1, e37. [9] Ma, L., Wagner, J., Rice, J.J., Hu, W., Levine, A.J. and Stolovitzky, G.A. (2005) A plausible model for the digital response of p53 to DNA damage. Proc. Natl. Acad. Sci. USA 102, 14266–14271. [10] Ribeiro, A.S. (2007) Dynamics of a two-dimensional model of cell tissues with coupled stochastic gene networks. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 76, 051915. [11] Gardner, T.S., Cantor, C.R. and Collins, J.J. (2000) Construction of a genetic toggle switch in Escherichia coli. Nature 403, 339–342. [12] Morishita, Y. and Aihara, K. (2004) Noise-reduction through interaction in gene expression and biochemical reaction processes. J. Theor. Biol. 228, 315–325. [13] Becskei, A., Seraphin, B. and Serrano, L. (2001) Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion. Embo J. 20, 2528–2535. [14] Isaacs, F.J., Hasty, J., Cantor, C.R. and Collins, J.J. (2003) Prediction and measurement of an autoregulatory genetic module. Proc. Natl. Acad. Sci. USA 100, 7714–7719. [15] Kramer, B.P., Viretta, A.U., Daoud-El-Baba, M., Aubel, D., Weber, W. and Fussenegger, M. (2004) An engineered epigenetic transgene switch in mammalian cells. Nat. Biotechnol. 22, 867– 870. [16] Weinberger, L.S., Burnett, J.C., Toettcher, J.E., Arkin, A.P. and Schaffer, D.V. (2005) Stochastic gene expression in a lentiviral positive-feedback loop: HIV-1 Tat fluctuations drive phenotypic diversity. Cell 122, 169–182. [17] Maamar, H. and Dubnau, D. (2005) Bistability in the Bacillus subtilis K-state (competence) system requires a positive feedback loop. Mol. Microbiol. 56, 615–624.