Lysogen stability is determined by the frequency of activity bursts from the fate-determining gene.

The ability of living cells to maintain an inheritable memory of their gene-expression state is key to cellular differentiation. Bacterial lysogeny serves as a simple paradigm for long-term cellular memory. In this study, we address the following question: in the absence of external perturbation, how long will a cell stay in the lysogenic state before spontaneously switching away from that state? We show by direct measurement that lysogen stability exhibits a simple exponential dependence on the frequency of activity bursts from the fate-determining gene, cI. We quantify these gene-activity bursts using single-molecule-resolution mRNA measurements in individual cells, analyzed using a stochastic mathematical model of the gene-network kinetics. The quantitative relation between stability and gene activity is independent of the fine details of gene regulation, suggesting that a quantitative prediction of cell-state stability may also be possible in more complex systems.

Introduction

The ability of living cells to maintain an inheritable memory of their gene-expression state is key to cellular differentiation (Monod and Jacob, 1961; Slack, 1991; Lawrence, 1992). A differentiated cellular state may be maintained for a long time, while at the same time allowing efficient state-switching (‘reprogramming') in response to the proper stimulus (Gurdon and Melton, 2008). However, even in the absence of external perturbation, a cell's gene-expression state may not be ‘infinitely stable' (irreversible; Lawrence, 1992). This is a consequence of the stochastic nature of all cellular reactions (Acar et al, 2005; Maheshri and O'Shea, 2007; Raj et al, 2008), which shift individual cells away from the ‘average state', and in particular may switch a cell from one state to another. A natural question then arises: how stable is a cell's gene-expression state, in the absence of an external perturbation? In other words, how long will a differentiated cell stay in the same state before spontaneously switching to an alternative one? What features of the underlying gene-regulatory network determine this stability?

The lysogenic state of an Escherichia coli cell harboring a dormant bacteriophage (prophage) lambda serves as one of the simplest examples for a stable cellular state (Ptashne, 2004, 2007; Oppenheim et al, 2005). Lysogeny is maintained by the activity of a single protein species, the lambda repressor (CI), which acts as a transcription factor to repress all lytic functions from the prophage in the E. coli cell, as well as to regulate its own production (Figure 1A–C; Ptashne, 2004). This feature of auto-regulation by the fate-determining proteins is commonly observed in systems displaying long-term cellular memory (Lawrence, 1992; Gurdon and Melton, 2008; Crews and Pearson, 2009). The lambda lysogeny system exhibits extremely high stability: spontaneous switching events occur less than once per 108 cell generations in the absence of cellular RecA activity (Little et al, 1999). At the same time, this genetic switch also exhibits fast and efficient switching in response to the appropriate stimulus, for example, damage to the bacterial genome (Oppenheim et al, 2005).

Figure 1

Maintenance of lysogeny in bacteriophage lambda. (A) The genetic circuit maintaining the lysogenic state. Cell fate is determined by a competition between two genes: the lambda repressor (cI), transcribed from the PRM promoter; and cro, produced from PR. The two gene products, CI and Cro, compete for binding to six operator sites (OR1–3, OL1–3) and mutually repress each other's transcription (Ptashne, 2004). (B) Stable states of the system. The CI-dominated lysogenic state will switch to a Cro-dominated state after a drastic decrease in the number of repressor (CI) proteins in the cell (induction). In a wild-type lysogen, Cro will activate a cascade of lytic genes leading to viral replication and cell death. In the reporter strain (NC416), the lytic pathway is blocked and thus a stable Cro-dominated state is maintained. (C) Modeling the lysogeny maintenance system. The diagram describes the features included in our stochastic simulation. A two-state transcription kinetics is assumed for both cI (PRMon and PRMoff states) and cro (PRon and PRoff states). After translation from the corresponding mRNA, CI and Cro dimers compete for binding to the OR operator sites. This process is a fast-equilibrium step. The probability of OR operator bound with CI, Cro and RNA-polymerase can be described by the grand-canonical partition function Ξ(CI2, Cro2), where s indicates the occupancy state of the operators (1–40) and i, j and k indicate the number (0, 1, 2 or 3) of the corresponding molecules bound at state s (Shea and Ackers, 1985; Darling et al, 2000). The factor μ(T) describes the reduced stability of the temperature-sensitive allele (cI857) at elevated temperatures. For more details of our stochastic simulation see Materials and methods section. (D) cI mRNA in lysogens, labeled using smFISH. Shown is an overlay of the phase-contrast and fluorescence channels. Individual cells were automatically recognized (white boundary) based on the phase-contrast image. Fluorescent foci (red) indicate the presence of cI mRNA molecules. The photon count from these foci was then used to estimate the number of mRNA molecules in each cell. The strain is wild-type lysogen MG1655(λwt). The scale bar is 2 μm. (E) cI mRNA number distribution in lysogenic cells. Images containing ∼500 cells were collected and analyzed to build the distribution of mRNA copy-number per cell. This experimental histogram was fitted to a negative binomial distribution (blue curve), parameters of which were used to calculate the transcriptional burst frequency r and burst size bTX (r=1.4±0.2, bTX=4.3±0.4, six independent experiments). The results of our stochastic simulation (red curve) are also shown for comparison. For experimental details see Materials and methods section.

The lambda system has been well characterized in terms of the regulatory circuitry that creates the stable lysogenic state. Specifically, the regulation of the two key promoters, PRM (producing CI) and PR (which initiates the lytic cascade at low repressor levels) has been mapped as a function of CI and Cro (the ‘anti-repressor') concentrations (Dodd et al, 2001; Ptashne, 2004; Figure 1A–C). A thermodynamic model using grand-canonical ensemble has been used to describe the occupancy of the operator sites controlling promoter activities and the corresponding protein levels (Shea and Ackers, 1985; Darling et al, 2000; Dodd et al, 2004; Anderson and Yang, 2008).

To predict the stability of the lysogenic state, characterization of the steady state has to be accompanied by quantification of the stochastic dynamics of gene activity. Recent studies have demonstrated that the production of both mRNA (Golding et al, 2005) and proteins (Cai et al, 2006; Yu et al, 2006) exhibit intermittent, non-Poissonian kinetics. Such ‘bursty' gene activity has been previously suggested to affect the switching of cellular states (Kaufmann et al, 2007; Schultz et al, 2007; Choi et al, 2008; Mehta et al, 2008; Gordon et al, 2009). Below, we characterize in detail the stochastic kinetics of gene activity in our system, in particular the frequency of activity bursts from the promoter PRM, which maintains the lysogenic state. Knowing this frequency allows us, in turn, to make a direct prediction of the stability of the lysogenic state.

Results

Single-molecule-resolution characterization of gene activity in a lysogen

Gene activity in individual cells was characterized using single-molecule fluorescence in situ hybridization (smFISH; Raj et al, 2008). We first quantified the statistics of cI mRNA numbers in a stable lysogen (MG1655(λwt) at 37°C, Figure 1D and E). The observed mRNA statistics displayed a variance-to-mean ratio larger than 1 (σ2/〈n〉=5.3±0.4, six independent experiments, ∼500 cells per experiment), indicating non-Poissonian kinetics for mRNA production.

mRNA number statistics were analyzed in the framework of a two-state model for transcription (Golding et al, 2005; Raj et al, 2006; Shahrezaei and Swain, 2008; Zenklusen et al, 2008; Figure 1C; Supplementary Figure 3). The gene is assumed to switch stochastically between ‘on' and ‘off' states, and mRNA is produced only in the ‘on' state. The resulting time-series of mRNA production is intermittent or ‘bursty' (Golding et al, 2005; Chubb et al, 2006; Raj et al, 2006; Zenklusen et al, 2008). The measured mRNA copy-number distribution allowed us to estimate the average transcriptional burst size bTX (number of mRNA molecules produced at each bursting event) and the average number of burst events r per mRNA lifetime τmRNA. The lifetime of cI mRNA (and similarly cro mRNA, see Materials and methods section) was measured using quantitative RT–PCR after inhibition of transcription with rifampicin (Bernstein et al, 2002). Together, these measurements allowed us to estimate kon=r/τmRNA, the rate of switching the gene ‘on' in the two-state model (see Materials and methods section for detailed derivations). Thus, based on the combined smFISH and mRNA lifetime experiments, we were able to estimate the average burst size and the burst frequency (i.e. frequency of activity events) of cI transcription. In the case of MG1655(λwt) at 37°C, we found a frequency of 1.4±0.2 events per min with an average burst size of 4.3±0.4 transcripts per event (six independent experiments).

Next, we extended the survey of system behavior by quantifying gene activity in the reporter strain NC416 (Svenningsen et al, 2005). The reporter strain carries a temperature-sensitive allele, cI857 (Hershey, 1971; Hecht et al, 1983). In this allele, a single mutation in the cI gene leads to decreased structural stability of the repressor protein at higher temperatures, and thus to a temperature-sensitive phenotype of the lysogenic state. The reporter strain contains the complete PRM/PR circuitry, but not the lytic genes. Therefore, cells do not die after switching occurs; instead, switched cells enter a Cro-dominated state (Svenningsen et al, 2005; Figure 1A and B).

We measured the copy-number distribution of cI and cro mRNA at different temperatures between 30 and 40°C (Figure 2A and B; two independent experiments; ∼500 cells per experiment). One can observe the expected transition from cI dominance (lysogeny) at low temperatures to cro dominance at higher temperatures. Both mRNA species exhibited the typical negative binomial statistics, indicating a bursty mode of transcription from both PRM and PR promoters throughout the temperature range. Each of the promoters maintained an approximately constant burst size when active (4.1±0.5 for cI and 1.7±0.5 for cro).

Figure 2

Tuning the state of the lysogeny system. (A) Two-color smFISH was used for counting cI (red) and cro (green) mRNA in the reporter strain NC416. Sample images of cells grown at different temperatures are shown in the top row. Images containing ∼500 cells at each temperature were processed to yield the mRNA number distributions in the population (bottom two rows). The observed distributions (red/green circles) were fitted with negative binomial distributions (red/green lines) and compared to results of the stochastic simulation (blue lines). For details of smFISH and the stochastic simulation see Materials and methods section. (B) mRNA numbers as a function of temperature. Diamonds represent experimental data; dashed lines represent simulation data. cI and cro mRNA in NC416 strain (cro+) are plotted against temperature. The measured cI mRNA level for wild-type lysogen at 37°C is also plotted. (C) Comparison of cI mRNA and protein levels. Protein levels at each temperature were quantified using immunofluorescence (see Materials and methods section). Each protein data point represents the mean of >100 cells. mRNA levels are from smFISH data. Both mRNA and protein levels were normalized by the value at 30°C. Error bars correspond to the s.e.m. values obtained from two independent experiments. The population-averaged mRNA and protein levels were found to be proportional to each other (correlation coefficient is 0.99).

We note that, at the population level, mRNA numbers were directly reflected in the protein level. We measured the CI protein level at different temperatures using quantitative immunofluorescence (see Materials and methods section). As shown in Figure 2C, the population-averaged levels of cI mRNA and CI protein were proportional to each other in the range: 30−40°C (R2=0.99).

Measurement of lysogenic stability

We quantified stability using the ‘switching rate' (S), the probability of switching from lysogeny to lysis in one cell generation (S is actually the switching rate per ∼1.4 cell generations; see Materials and methods section). S was measured experimentally in a fully functional lysogen (Figure 3). As host, we used a RecA-deficient strain, JL5902 (Little et al, 1999), because in wild-type RecA background the stability is masked by frequent spontaneous activation of the cell's SOS response (Little et al, 1999; see Figure 3A). We estimated S based on the number of free phages in an exponentially growing culture of lysogens (Little et al, 1999). Specifically, S≈φ/BM, where φ is the number of free phages in the culture, B is the number of bacterial cells and M is the average number of phages released per cell lysis (∼200 at 30°C and 40°C; see Materials and methods section). It is noteworthy that a constant switching rate S implies a constant ratio of free phages to bacteria during cell growth. Our data suggests that this is indeed the case (see Supplementary Figure 9). We measured S values for the temperature-sensitive prophage (cI857) in the temperature range 28–36°C (Figure 3A). The observed S values covered approximately eight orders of magnitude. We also conducted measurements for wild-type prophage, in both RecA+ and RecA− backgrounds, and observed very little change in S over the temperature range (Figure 3A), suggesting that the changes to repressor activity in the cI857 allele dominate over all other temperature-dependent effects (Ryals et al, 1982; Farewell and Neidhardt, 1998).

Figure 3

Stability of the lysogenic state. (A) Measured values of spontaneous switching-rate per cell generation (S) for temperature-sensitive (cI857) prophages in both RecA+ and RecA– hosts (red and blue triangles, respectively), and wild-type prophage (red and blue squares, respectively). For experimental details see Materials and methods section. (B) The relation between lysogen stability and PRM activity. Plotted is the measured switching rate (S) as a function of the number of activity events from PRM in one protein lifetime (R), for the wild-type lysogen (red circle), cI857 at different temperatures (blue triangles), and mutants in both the cI gene (white triangles) and the PRM promoter (white squares). The simulation data includes a 20% error estimate on the effective CI activity (shaded blue area; see Materials and methods section). The points fall close to the theoretical prediction given by S=exp(−R) (solid black line), as compared to the prediction of two alternative hypotheses: non-bursty (Poissonian) gene activity (dashed black line); or protein production from an efficient ribosome-binding site yielding bCI=bLacZ (dot-dashed black line). Source data is available for this figure at www.nature.com/msb.

Stability is determined by the frequency of activity bursts from PRM

When examining the relation between gene activity and lysogen stability (Figure 3B), we observed a simple exponential dependence of the switching rate S on the frequency of activity bursts from the PRM promoter. Specifically, S was well-described by the expression S=exp(−R)=exp(−konτ/ln 2), where kon is the rate of transcription bursts and τ is the cell doubling time. Both parameters were measured in experiment. For the temperature-sensitive allele, R is further multiplied by a factor μ(T), which describes the decreased fraction of active CI proteins at increased temperatures. The value for μ(T) was calculated using a comparison of the measured mRNA levels to the predictions of the stochastic simulation (see Materials and methods section; Supplementary Figure 4). Our estimation of μ(T) also agrees well with previous experimental data (Villaverde et al, 1993; Isaacs et al, 2003).

The exponential dependence found above can be intuitively understood using the following simple model: we assume that CI molecules are produced from the PRM promoter following discrete bursts of cI mRNA, and that the occurrence of the transcription-burst events obeys Poissonian statistics (Golding et al, 2005; Friedman et al, 2006). The average frequency of events is kon. Thus, the probability of the cell NOT producing any cI mRNA (and therefore repressor proteins) for a duration t is P0(t)=exp(−kont). Next, we note that the mean lifetime of CI proteins, due to cell growth and division, is τ/ln 2 (where τ is the cell-doubling time). This timescale is much longer than the lifetime of mRNA, as measured directly (see Supplementary Figure 1). By plugging the protein lifetime into the expression for P0(t), we obtain the probability of not producing new CI for the whole lifetime of the protein: exp(−konτ/ln 2)=exp(−R). This is just the behavior observed in experiment for the spontaneous switching rate. Thus, the loss of lysogeny (switching) occurs if no new cI mRNA (and the downstream proteins) is produced for the mean lifetime of CI protein, which is ∼1.4 cell generations.

The expression R=konτ/ln 2 can also be written as R=〈CI〉/bCI, where 〈CI〉 is the average number of repressor protein per cell and bCI is the average number of proteins produced after one transcription burst. bCI=bTXbTL, where bTX is the average number of cI mRNAs produced per transcription burst and bTL is the average number of CI proteins made from one mRNA molecule during its lifetime. For CI, we estimated bCI to be ∼20 (see Materials and methods section).

As an additional test for the agreement between the predicted switching rate and experiments, we plot in Figure 3B the predicted switching rate under two alternative hypotheses (while maintaining 〈CI〉 unchanged): the first one is that proteins are produced in a Poissonian manner, one-at-a-time (no burstiness), in which case one expects S=exp(−〈CI〉). As shown in Figure 3B, the predicted switching rate in this case is orders of magnitude lower than the experimental data. The strong discrepancy demonstrates the critical role of bursty mRNA and protein production in limiting the stability of the lysogenic state (Mehta et al, 2008). A second case we test is what would happen if the cI gene had an efficient ribosome-binding site, for example, as that of the lacZ gene (Shean and Gottesman, 1992). In such a case, ∼six-fold more proteins would be produced from each mRNA (Kennell and Riezman, 1977). As seen in Figure 3B, the predicted switching rates in this case are considerably higher than those found in experiment. Thus, the relatively inefficient translation of the cI transcript increases the stability of the lysogenic state.

Finally, we asked whether the quantitative relation between the number of activity bursts and the spontaneous switching rate can be demonstrated in other alleles of PRM-cI beyond the cI857 case. The lambda lysogeny system has been studied for many years, and many mutations have been created, targeting multiple features of the lysis/lysogeny switch (see e.g. Little et al, 1999; Atsumi and Little, 2004, 2006a, 2006b; Michalowski et al, 2004; Michalowski and Little, 2005). However, using PRM-cI mutants to test the gene-activity/stability relation presents the following challenge: unlike cI857, most mutants do not provide a continuous spectrum of lysogen stability but instead only a single stability phenotype. In most cases, the stability phenotype obtained is one of two: close to the wild-type value (and to the sensitivity limit of measuring lysogen stability, due to the appearance of unstable mutants (Little et al, 1999; Aurell and Sneppen, 2002; see Figure 3B) or the lysogen is ‘too unstable', such that lysogenization of the host cell fails in the first place. As seen in Table I below, out of 18 alleles examined, 12 exhibited a switching rate very close to wild-type and four were ‘too unstable'. Only two gave an intermediate switching rate (see Figure 3B). Despite these limitations, it can be seen in Figure 3B that the measured stability of the mutant lysogens was consistent with the theoretical prediction. Taken together with our previous data set, these results further support our key observation that the stability of the lysogenic state depends in a simple manner on the number of activity bursts from the fate-determining gene, and that this relation is general, that is, it holds even when the promoter and gene-coding region are modified, as demonstrated by the mutants.

Table 1

PRM-cI mutants used for analyzing lysogen stability

Name	Genotype	PRM/CI phenotype	Reference	Source	Burst size, b	S [generations−1]a	R [generations−1]b
aS, switching rate (per 1.4 cell generation), measured in JL5902 (RecA−). Mean and s.e. for each strain were calculated from two to eight independent experiments.
bR, frequency of transcription bursts (per 1.4 cell generation), calculated using the formula R=rτ/τmRNA ln 2. τ is the cell doubling time at 37°C, τmRNA is the cI mRNA lifetime (∼4 min for NC416 strain, see Materials and methods section) and r is the frequency of transcription bursts (per mRNA lifetime), obtained from to the single-cell mRNA data. R, mean and s.e. for each strain were calculated from two to four independent smFISH experiments, except those of NP7, NP8 and NP11 (one smFISH experiment only).
cNL, non-lysogenizable.
dFor the temperature-sensitive allele cI857, switching rate data is from the λIG2903 lysogen, whereas burst size is from the NC416 reporter strain (Svenningsen et al, 2005).
λIG831	Wild type	Wild type		Lab stocks	4.3±0.4	9.0 × 10−9±6.4 × 10−9	16±1
λIG2504	cI T88CTAT → TGT	Covalent dimerization	Sauer et al, 1990	Lab stocks	4.3±0.5	5.8 × 10−9±5.9 × 10−9	19±3
λIG1006	2 × TAG (amber) at AA102GAGTAC → TAGTAG	N-terminus only	Reidhaar-Olson and Sauer, 1988; Sauer et al, 1979	Lab stocks	NLc	NL	NL
λIG04061	cI E34KGAA → AAA	Tight operator binding	Sauer et al, 1990	Lab stocks	0.3±0.3	2.6 × 10−3±1.4 × 10−3	3.3±3.3
λIG15052	cI L18FCTT → TTT	Misfolding	Lim and Sauer, 1989	Lab stocks	NL	NL	NL
λIG04062	cI Q44YCAG → TAT	Weak operator binding	Sauer et al, 1990	Lab stocks	NL	NL	NL
λIG15051	cI Q44LCAG → CTG	Weak operator binding	Sauer et al, 1990	Lab stocks	NL	NL	NL
λIG28061	cI V36IGTC → ATC	Misfolding	Lim and Sauer, 1989	Lab stocks	3.6±0.5	1.4 × 10−8±1.1 × 10−8	12±4
λIG28062	cI L18V CTT → GTT	Misfolding	Lim and Sauer, 1989	Lab stocks	4.0±0.6	7.5 × 10−8±3.4 × 10−8	16±2
λJL8525	cI D38N	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	4.0±0.5	3.7 × 10−4±7.5 × 10−5	9±1
λNP2	cI D38NPRM-35 TAGA → GCTGPRM-10 GATT → TATT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	4.5±0.6	2.0 × 10−8±6.4 × 10−9	18±4
λNP3	cI D38NPRM-35 TAGA → CATTPRM-10 GATT → GAAT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	0.9±0.4	1.0 × 10−8±2.2 × 10−9	14±1
λNP4	cI D38NPRM-35 TAGA → CCTTPRM-10 GATT → CCAT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	1.4±0.5	1.7 × 10−8±2.3 × 10−9	16±1
λNP5	cI D38NPRM-35 TAGA → CTAAPRM-10 GATT → GAAT	No positive autoregulation, reduced PRM activity	Michalowski and Little, 2005	J. Little	2.1±0.5	3.0 × 10−8±1.4 × 10−8	15±4
λNP6	cI D38NPRM-35 TAGA → CCCAPRM-10 GATT → TGAT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	1.8±0.3	1.5 × 10−8±2.6 × 10−9	15±2
λNP7	cI D38NPRM-35 TAGA → TACCPRM-10 GATT → TACT	No positive autoregulation, reduced PRM activity	Michalowski and Little, 2005	J. Little	0.5±0.1	2.3 × 10−8±2.8 × 10−9	25±7
λNP8	cI D38NPRM-35 TAGA → GTGTPRM-10 GATT → GTAT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	2.8±0.6	2.8 × 10−8±1.0 × 10−8	14±3
λNP10	cI D38NPRM-35 TAGA → CCAAPRM-10 GATT → GAAT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	1.1±0.4	1.5 × 10−8±4.7 × 10−9	16±3
λNP11	cI D38NPRM-35 TAGA → CTCAPRM-10 GATT → TGCT	No positive autoregulation; reduced PRM activity	Michalowski and Little, 2005	J. Little	1.4±0.7	3.1 × 10−8±2.0 × 10−9	21±3
λIG2903 d	cI A66TGCA → ACA	Temperature-sensitive	Hecht et al, 1983	Lab stocks
		30°C			4.7±0.6	5.8 × 10−9±1.5 × 10−9	17±3
		32°C			4.8±0.7	1.1 × 10−7±7.3 × 10−8	16±2
		34°C			4.1±0.7	2.1 × 10−6±1.6 × 10−6	10±2
		36°C			6.5±0.3	1.7 × 10−2±8.0 × 10−3	4±1

Discussion

We have shown that the stability of a bacterial lysogen is given by the simple expression exp([number of gene-activity events in τ]), where τ is the relevant time-scale for maintaining the lysogenic state. We note that exponential switching probabilities have been previously predicted in theoretical models of gene regulatory circuits, assuming weak-noise conditions (Bialek, 2001; Roma et al, 2005; Mehta et al, 2008). In particular, a few past studies applied thermal barrier-crossing (‘Kramers problem') approaches to the problem of cellular state-switching (Bialek, 2001; Aurell and Sneppen, 2002; Mehta et al, 2008). However, it is important to note that when a system is controlled by pure dynamical rules (rather than thermal fluctuations), converting the master equation into a stochastic differential equation (e.g. Langevin) becomes challenging, especially when large fluctuations affect the outcome—as in the case of the lysogen. Barrier-crossing analysis is also restricted by whether an effective potential and an appropriate reaction coordinate can be identified or not, and by the question of how to define the effective ‘temperature' (Lu et al, 2006).

In our study, we avoided these theoretical challenges by taking the following approach: the rate of switching was first calculated by Gillespie simulation of the master equation describing the gene network. An intuitive understanding of the results from both experiments and simulation was obtained using the argument of survival probability; switching occurs if there are not enough CI molecules to maintain lysogeny. This event only happens at the rare events that no CI is made for a specific period of time. The ‘survival' probability naturally explains the exponential behavior seen in experiment and simulations, with the key parameters being the promoter burst frequency and the protein lifetime.

At first glance, it may seem surprising that a simple expression captures the behavior of a real-life, naturally evolved system in which the stability is believed to be an important phenotype (Little et al, 1999). Specifically, the lambda lysogeny circuit has long served as a paradigm for the intricacy and precision of gene regulation (Ptashne, 2006; Court et al, 2007), in which the proper state of the system depends on the finely-tuned balance between the affinities of CI and Cro to their six DNA targets (OR1–3 and OL1–3). In contrast, our findings suggest that the stability of a genetic switch can be estimated simply based on the rate of gene activity, thus the intricacy is absent in the expression describing the stability of the switch.

We note that in line with the observation that lysogen stability is insensitive to many system parameters, there is a body of work from the last decade, mainly from the Little lab (Little et al, 1999; Atsumi and Little, 2004, 2006a, 2006b; Michalowski et al, 2004; Michalowski and Little, 2005), pointing to the robust performance of the lambda lysogeny switch even when the underlying gene circuit is modified. For example, it has been shown that a stable lysogenic state could be maintained when the relative affinities of the operator sites to CI and Cro were reversed (Little et al, 1999), when the positive autoregulation by CI was deleted (Michalowski and Little, 2005), when PRM was made stronger or weaker (Michalowski et al, 2004), and even when Cro and CI were replaced by the Lac and Tet repressors, respectively (Atsumi and Little, 2006b). Although only semi-quantitative in nature, these studies suggest that the genetic circuitry found in lambda is not unique, and many alternative systems can maintain a stable lysogenic state. According to our findings here, the critical element is whether the new design can produce the proper rate of gene activity.

Finally, it is tempting to contemplate the possible relevance of our results concerning the stability of cellular states to higher systems, in which the ability of cells to maintain an inheritable memory of their gene-expression state is key to cellular differentiation (Monod and Jacob, 1961; Gurdon and Melton, 2008). Admittedly, the maintenance of bacterial lysogeny does not exhibit the complexity of cell differentiation in higher eukaryotic systems. However, even though a range of additional mechanisms have a role in cellular memory in the higher systems (Burrill and Silver, 2010), the fundamental feature of autoregulation by the fate-determining protein seems to be a central element (Slack, 1991; Lawrence, 1992; Gurdon and Melton, 2008; Crews and Pearson, 2009). We thus look forward to investigating the stability of epigenetic states in higher systems.

Materials and methods

Strains and media

Phage and bacterial strains used in this work are listed in Supplementary Table 1. Cell growth and manipulation followed standard protocols (Sambrook and Russell, 2001).

Two main phage strains, λIG831 (‘λwt') and λIG2903 (‘λts'), were used in this study. Both were kanamycin-resistant due to a cassette inserted at the bor gene (a gift from Barry Egan, Adelaide University). This antibiotic marker was used to select lysogens after infection. λwt carried the wild-type cI allele, whereas λts carried the temperature-sensitive allele, cI857. λts was constructed following the method of Dodd et al (2001). Briefly, mutations were introduced using the Quickchange kit (Stratagene) into a plasmid carrying the required fragment of the lambda genome. Cells carrying the mutated plasmid were infected with λ434 and the mutant sequences were crossed onto the lambda chromosome. The immλ recombinants were selected among the progeny as plaques on λ434 lysogens. Recombination led to the replacement of imm434 region with lambda sequences.

We infected RecA+ strain MG1655 and RecA− strain JL5902 (a gift from John Little, University of Arizona) with λwt and λts to form four types of lysogen: (RecA+, cIwt), (RecA+, cI857), (RecA−, cIwt) and (RecA−, cI857).

Two reporter strains were used, NC416 and NC417 (a gift from Don Court, NCI Fredrick) (Svenningsen et al, 2005). They were both kanamycin-resistant. The immunity region of the reporter strains is shown in Figure 1A.

Single-molecule fluorescence in situ hybridization (smFISH)

Our protocol is based on the ‘Star-FISH' method developed by Raj et al (2008). Modifications were made to adapt the protocol to E. coli.

Sample preparation

An overnight culture was diluted 103-fold into 10 ml LB, and incubated on a shaker. When the optical density (OD600) reached ∼0.4, cell culture was cooled in an ice-water bath and the cells were collected by centrifugation at 4°C. The cells were washed by resuspending in 1 ml ice-cold 1 × PBS. The cell suspension was centrifuged again and the supernatant was removed. The pellet was gently resuspended in 1 ml 1 × PBS with 3.7% formaldehyde. The solution was gently mixed at room temperature for 15 min. Cells were pelleted by centrifugation and washed twice with 1 ml 1 × PBS. The pellet was resuspended in 1 ml 70% ethanol and mixed gently on a nutator at room temperature for at least 1 h to permeabilize cell membrane.

Hybridization procedure

The cells in 70% ethanol were pelleted by centrifugation. The pellet was resuspended in 1 ml wash solution (40% formamide, 2 × SSC), then set aside for a few minutes. Meanwhile, 2 μl cI probes and 2 μl cro probes were added to 50 μl hybridization solution (10% dextran sulfate, 2 mM vanadyl-ribonucleoside complex, 0.02% RNase-free BSA, 50 μg E. coli tRNA, 2 × SSC, 40% formamide). Probes were mixed well with the viscous hybridization solution by vortexing. Next, the cells in wash solution were pelleted by centrifugation. After the supernatant was removed, the pellet was resuspended in the hybridization solution, and then incubated at 30°C overnight. After hybridization, 5 μl of the hybridization mix was washed with 100 μl wash solution. Cells were washed twice with 100 μl wash solution. The pellet was then resuspended in 20 μl 2 × SSC solution for imaging under the microscope.

DNA oligonucleotide probes were designed using the program (www.singlemoleculefish.com) from Arjun Raj (MIT), and were ordered with 3′-end amine modification from Biosearch Technologies. The probes for hybridizing cI mRNA were labeled with TAMRA (Invitrogen), and the probes for hybridizing cro mRNA were labeled with Cy5 (GE Health Life Science). The labeling process followed the procedures described previously (Raj et al, 2008). The sequences of the probes are given in Supplementary information.

Microscopy and image analysis

After hybridization and washing (above), 1–2 μl sample was placed between a piece of 1.5% agarose gel (∼1 mm thick) and a glass coverslip. Still images of cells were acquired using an inverted epifluorescence microscope (Eclipse TE2000-E, Nikon) and a cooled CCD camera (Cascade 512, Photometrics). A × 100 oil immersion phase objective was used in conjunction with a × 2.5 lens in front of the camera. The microscope and camera were controlled by the Metamorph software (Molecular Devices). Phase contrast images were taken with 100-ms exposure time. TexasRed filter set (Nikon) was used for fluorescence imaging of cI mRNA, with 400 ms exposure time and nine-slice z-stacks with 120-nm spacing. Cy5 filter set (Nikon) was used for fluorescence imaging of cro mRNA, with 8-s exposure time. For each sample, multiple images were taken and ∼500 cells were analyzed from those images.

The images were automatically analyzed using home-made MATLAB programs. Cell recognition was done using the Schnitzcell program (a gift from Michael Elowitz (Caltech)). Spot recognition was performed with a home-made program developed in our lab. To localize the foci, we scanned for the local maxima in both x- and y-directions. The spot at the local maximum in both directions was taken as the center of one fluorescent spot. As multiple mRNA can be in proximity of one fluorescence spot, the fluorescence intensity of the spot was measured. The histogram of the spot fluorescence intensity was plotted and we identified the position of the primary peak as the intensity of a single mRNA. The total mRNA in a cell was estimated by summing the estimated number of mRNA from all the spots in the cell. We then fit the total mRNA per cell data to a negative binomial distribution with parameters r and p0 to obtain estimates of the transcription kinetics parameters: r (burst frequency) and bTX (burst size that is obtained from p0, see below). We report the mean and s.e.m. for these parameters, from two to six independent experiments in each case.

Extraction of transcription parameters from smFISH data

Transcription is described using a phenomenological two-state model for gene activity (Supplementary Figure 3; Peccoud and Ycart, 1995; Golding et al, 2005; Golding and Cox, 2006; Raj et al, 2006; Shahrezaei and Swain, 2008). In this model, the gene stochastically switches between ‘on' and ‘off' states, and mRNA is produced only when the gene is ‘on'. The mRNA distribution at the steady state can be analytically solved for this model (Raj et al, 2006; Shahrezaei and Swain, 2008). When koff≫krd, that is, when the dwell time at the ‘on' state is much shorter than the mRNA lifetime, the mRNA statistics follows the negative binomial (NB) distribution:

where n is the number of mRNA molecules per cell. The parameter p0 equals koff/(koff+kTX), and describes the probability of the gene turning from the ‘on' state to the ‘off' state. The parameter r denotes the average number of gene activity events (i.e. the number of times that the gene has been ‘on') during one mRNA lifetime. The mRNA molecules in a cell in general come from multiple gene-activity events, and the time window in which we can detect them is the lifetime of mRNA. In the special case that r is equal to 1, we get P(n, r=1, p0)=p0(1−p0), which is simply the geometric distribution and it describes the probability of transcribing n mRNAs when the gene has been ‘on' once.

The measured mRNA distribution allowed us to estimate the average transcriptional burst size bTX (number of mRNA molecules produced at each bursting event): bTX=kTX/koff=(1−p0)/p0, as well as the average number of burst events r per mRNA lifetime τmRNA.

The expression for kon is derived as follows. The average number of mRNA per cell follows

where 〈Pon〉=kon/(kon+koff)≈kon/koff and krd is the decay rate of mRNA. At the steady state,

therefore the switching-on rate of the gene is

where we used the relation 〈mRNA〉=bTXr. As r is known from smFISH data, we can calculate the rate kon once we know mRNA lifetime τmRNA. Standard RNA lifetime assay (Bernstein et al, 2002) was carried out as described below. The decay curves of cI mRNA and cro mRNA are shown in Supplementary Figure 1. Thus, both kon and the burst statistics: bTX for transcription are obtained. The transcription rate kTX is set at bTXkoff.

Immunofluorescence

Reagents, incubation times and temperatures were based on established immunofluorescence protocols in bacteria (Maddock and Shapiro, 1993; Harry et al, 1995; Raj et al, 2008). Unless otherwise stated, all steps were performed at room temperature. Cells were grown overnight at 30°C, diluted 1:1000 in 15 ml of LB medium, and grown at the experimentally desired temperature (30–40°C) until cultures reached OD600∼0.4. Cells were collected by centrifugation at 4°C and resuspended in 1 ml of fixation solution (1 × PBS, 3.7% formaldehyde), transferred to an Eppendorf tube, and mixed gently for 30 min. Cells were then washed twice with 1 × PBS, resuspended in 300 μl water, ethanol added to reach a final concentration of 70% and mixed gently for 1 h. The suspension was centrifuged and permeabilized by adding 100 μl lysozyme solution (10 μg/ml lysozyme, 50 mM glucose, 25 mM Tris–HCl, 10 mM EDTA (pH 8.0)) and left for 10 min. 1 ml of 1 × PBS was then added, the cells collected, and washed twice in 1 × PBS. Before antibody labeling, two blocking steps were performed to decrease non-specific binding. Cells were incubated in 100 μl blocking solution (1 × PBS, 2% BSA, 0.05% Tween 20) for 30 min, and washed twice with 1 × PBS. Then, cells were resuspended in 100 μl of blocking solution containing 1:100 rabbit anti-CI serum (a gift from Ian Dodd, University of Adelaide; Dodd et al, 2001) and incubated for 1 h. Cells were collected by centrifugation, washed with 1 × PBS, blocking solution added and incubated for 30 min. Cells were collected by centrifugation and resuspended in 100 μl of blocking solution containing 1:100 mouse anti-rabbit secondary antibody labeled with Cy3 dye (Jackson Immunochemicals) and incubated for 1 h in the dark. The suspension was then diluted in 1 ml 1 × PBS, kept for 5 min and then washed twice in 1 × PBS. The immunolabeled cells were resuspended in 1 × PBS (5–10 ml), placed on a #1 slide, covered by a thin 1.5% agarose gel pad in 1 × PBS and imaged under the microscope.

Quantification of protein levels

To quantify the relative levels of CI in the NC416 reporter strain, a fiducial marker (FM) strain was constructed by transforming NC416 with a pTrc99A plasmid expressing GFP. In each experiment, three strains: MG1655 (wt), FM and NC416 were grown as described above. FM was grown overday with 0.1 mM IPTG to induce the production of GFP. FM and MG1655 were grown overday at 30°C, while independent cultures of NC416 were grown at 30, 32, 34, 36, 38 and 40°C. Before fixation, ∼7.5 ml of MG1655 or NC416 cultures were mixed with the same volume of FM and the immunofluorescence protocol performed. After imaging, cells were recognized using automated procedures in MATLAB (see Microscopy and Image Analysis). At each temperature, the mean CI expression level was calculated as the ratio of the mean cell intensity of the Cy3 signal of NC416 (or MG1655) to the one from FM. The Cy3 level obtained from MG1655 was used as baseline and subtracted from the NC416 levels.

RNA lifetime measurement using quantitative RT–PCR (qRT–PCR)

RNA decay experiments (Bernstein et al, 2002) were carried out as follows. A volume of 20 ml culture was grown to OD600∼0.4, and divided into two flasks at t=0. In one of the flasks, rifampicin was added to a final concentration of 500 μg/ml to inhibit transcription. The culture in the other flask, without rifampicin, served as the reference. We continued to incubate the cultures with shaking. A volume of 1.5 ml of culture was taken from each flask every 3 min and immediately mixed with RNAprotect Bacterial Reagent (Qiagen) to stabilize cellular RNA. This was continued until 12 min after adding rifampicin. RNA was then extracted from the RNA-stabilized samples using RNeasy Mini Kit (Qiagen). After RNA was extracted, the total RNA concentration and the purity were estimated from the UV absorption spectrum (Cary 100, Varian). All RNA samples were diluted to a final concentration of 50 μg/ml for the qRT–PCR experiment. qRT–PCR was performed on a MiniOpticon System (Bio-Rad) using HotStart-IT SYBR Green One-step qRT–PCR Master Mix kit (USB). The primers for cI mRNA were: 5′-AGGACGCACGTCGCCTTAAAG-3′ and 5′-GCCATTAAATAAAGCACCAACGCCTG-3′. The primers for cro mRNA were: 5′-GCTTTGGGCAAACCAAGACAGCTAAA-3′ and 5′-GCTTTACCTCTTCCGCATAAACGCTTC-3′.

The relative amount of target mRNA was measured for samples at different times and the decay rates were extracted from the exponential fit of experimental data (Supplementary Figure 1). The decay rates of cI mRNA in the lysogen MG1655(λwt) and the reporter system NC416 were measured at 30°C; the decay rate of cro mRNA in the reporter system NC416 was measured at 40°C. The lifetime values of cI mRNA and cro mRNA are 4.0±0.8 and 2.8±0.4 min, respectively.

Measuring spontaneous lysis rate

The measurement of lysogen stability is based on the protocol of Little et al (1999). Lysogen culture was grown overnight in LBGM (LB with 1 mM MgSO4 and 0.2% glucose). The overnight culture was centrifuged and resuspended in 1 ml LBGM. The washed culture was diluted 103-fold in 15 ml LBGM and grown at the desired temperature (between 28 and 36°C). During the exponential growth phase (OD600∼0.1–0.5), 1 ml culture was taken. OD600 value was immediately measured and 15 μl chloroform was quickly added to the sample and mixed by inverting 10 times. The samples were then stored at 4°C for phage titration later.

Phages were titered using the following procedures (Little et al, 1999). Indicator strain (LE392) was grown overnight in LBMM (LB with 1 mM MgSO4 and 0.2% maltose). The overnight culture was diluted 100-fold in LBMM, and then grown to OD600∼0.5. The cells were pelleted by centrifugation at 4°C, 1000 g for 10 min. The pellet was resuspended to one-tenth original volume in 1 × PBS with 10 mM MgSO4 and stored at 4°C. A volume of 10 or 100 μl of the chloroform-treated solution or its dilution series was added to 100 μl host culture and mixed well. The tube was incubated at 37°C for 15 min, to allow phages to attach to cells. Meanwhile, NZYM (1% NZ amine, 0.5% yeast extract, 0.5% NaCl, 0.2% MgSO4, 0.75% agar (pH 7.0)) top agar was melted and kept at 48°C. After infection, the bacteria/phage solution was mixed with 3 ml NZYM top agar and plated onto dry and pre-warmed NZYM plates. Plaques became visible after ∼12-h incubation at 37°C. The number of bacterial cells, B, was estimated based on OD600 value. The average number of phages released per lysis, M, was measured separately, as described in the next section. The spontaneous induction rate was calculated based on the formula derived below.

Derivation of the spontaneous induction rate

For the case of exponentially growing culture, the number of cells B at time t is:

where τ is the generation time and B0 is a constant reflecting the initial conditions. The rate of increase in phage concentration is the product of the three factors: the spontaneous induction rate I, the number of phages released at each lytic event M, and the bacterial concentration B.

Note that B is calculated at time t−T, where T (the ‘latent period' Little et al, 1999) is the time interval between induction and lysis. We can then write the expression for the free phage concentration in the lysogen culture at time t:

The ratio of free phage to bacteria in culture can be written as:

This expression converges within a few cell generations to a constant value. In addition, both the generation time and the latent period are of the order ∼1 h, and thus e− and 1−e− are close to ∼1. We then find that the ratio of free phage to bacteria in a growing culture of lysogens is approximately given by the product of the spontaneous switching rate per generation and the number of phages released at lysis:

Note that this simple model predicts that the free phage-to-bacteria ratio will be approximately constant during exponential growth. As seen in Supplementary Figure 9, we found this to be the case. In our study, we used this measured ratio to calculate the spontaneous induction rate, using: S=Iτ/ln 2≈φ/BM (per 1.4 cell generation).

Measuring the number of phages released at cell lysis

We used the following protocol (Little et al, 1999) to measure the average number of phages released per cell lysis. Bacterial culture (lysogen with temperature-sensitive CI) was grown in LBGM overnight at 30°C. The overnight culture was diluted 103-fold in 20 ml LB and incubated at 30°C. At OD600∼0.4, 10 μl culture was taken for titering bacterial concentration. The rest of the culture was immediately transferred to 40°C and incubated for 15 min to induce lysis. The culture was then divided and moved to the studied temperatures (30°C and 40°C in this case). A volume of 1 ml culture was taken out every 15 min, gently mixed with 15 μl chloroform to kill the bacteria, and stored at 4°C. This procedure was continued until the culture became clear, indicating the complete lysis of bacterial cells. The phage concentration of each sample was determined using the titering assay described above. As shown in Supplementary Figure 2, the phage concentration reached plateau after all bacterial cells had lysed. The number of phages released per lysis was calculated as the ratio between the phage concentrations at the plateau to the bacterial cell concentration determined from colony counting.

Stochastic simulation of the lysis/lysogeny switch

Transcription parameters estimated from smFISH data, combined with known parameters of protein translation (Kennell and Riezman, 1977; Shean and Gottesman, 1992), allowed us to construct a fully calibrated stochastic model for the maintenance of lysogeny. We used the Gillespie algorithm (Gillespie, 1977) to simulate the stochastic kinetics of the system shown in Figure 1C. In the algorithm, the reaction probabilities are obtained from the kinetic rate parameters (summarized in Supplementary Table 2). Each gene under consideration (cI or cro, chemical rates for each species are denoted by superscripts) can be in the inactive (‘off') or active (‘on') state. When the gene is in the inactive state, it switches to the active state at a rate kon=PRNAPfon. In this study, PRNAP is the probability that an RNA polymerase is bound to the promoter, calculated from a grand canonical ensemble that takes into account all occupancy configurations of CI and Cro dimers at the operator sites (Ackers et al, 1982; Shea and Ackers, 1985; Anderson and Yang, 2008). fon is defined as the rate of transcription initiation when CI or Cro are bound at the proper operator to enhance transcription. We assume fon is constant. At the level of ensemble average, we expect that 〈PRNAP〉fon as calculated from simulation to equal kon as measured from smFISH experiment. We estimated fon using the following procedure: we first use the experimentally derived kon value as the ansatz for fon, used in the first round of simulation. We then calculate 〈PRNAP〉 and adjust fon so that 〈PRNAP〉fon becomes closer to the experimental kon value. We repeat this procedure till 〈PRNAP〉fon agrees well with the experimental kon value. By the above procedure, we are able to tune our simulation model so that its mRNA profile matches with experimentally measured mRNA profile. For the temperature sensitive strain NC416, we performed the above procedure for cI using the experimental mRNA profile at 30°C and for cro using the mRNA profile at 40°C.

In the simulation run, we update the value of PRNAP after each time step, and calculate the new kon for each time step. When the gene is in the active state, it can switch back to the inactive state at rate koff, or be transcribed to give a new mRNA at rate kTX. To estimate koff and kTX, we note that kTX=bTXkoff. bTX is the transcriptional burst size, estimated from smFISH experiments for each gene. We choose koff to be an arbitrary fast rate, assuming the dwell time in the active state is very short. Each existing mRNA is translated into protein monomers at rate kTL and is degraded at rate krd. The translation rate of cI mRNA is taken to be six-fold lower than that of lacZ mRNA (Kennell and Riezman, 1977; Shean and Gottesman, 1992), or 0.02 s−1. The translation rate of cro mRNA into Cro monomers is assumed to be the same as that of lacZ mRNA. The degradation rates of cI and cro mRNA were estimated from qRT–PCR as described above. CI dimerization was assumed to be a fast process, such that CI monomers and dimers are in equilibrium at all times (Johnson et al, 1980). Cro dimers association and dissociation rates were kascro and kdiscro, respectively (Jia et al, 2005). The resulting number of CI dimers discounted by μ(T) (see below) and the number of Cro dimers were fed into the calculation of PRNAP to update the promoter activity. Cell growth rate k0 was measured directly. All mRNA and protein species were binomially partitioned upon cell division.

As a first indication of its validity, the average CI protein level predicted by the model (∼300 molecules per cell in a lysogen) is in good agreement with the known value (Ackers et al, 1982). Note that no parameters were adjusted to obtain this estimate. Beyond previous models, our model accurately captures the fluctuations of gene activity that set a limit on the stability of the lysogenic state because the underlying burst processes are described in the model and calibrated by experimental data.

Modeling the cI857 allele

To capture the effect of temperature in our model, we introduced the parameter μ(T), which denotes the fraction of active repressor proteins at a given temperature. For wild-type CI, μ(T)=1. For cI857, μ(T) decreases as temperature is increased, indicating only partial functionality of repressor (Supplementary Figure 4). The exact molecular mechanism underlying the reduced activity of repressor at higher temperatures is unknown: CI may be degraded at a higher rate (Isaacs et al, 2003); alternatively, only a fraction of CI molecules may remain active (Oppenheim and Noff, 1975). Importantly, however, these details are not critical for our analysis. We identified a simple sigmoidal shape for μ(T) based on the average number of cI and cro mRNA over the complete range of temperatures examined (Supplementary Figure 4). μ(T) relation thus allows a mapping between the simulation and experiments. As a test for the validity of this description, the same curve for μ(T) can be used to predict cI activity in a cro– reporter system (NC417; Svenningsen et al, 2005). This was done by eliminating Cro in the theoretical model. The good agreement between theory and experiment (Supplementary Figure 5) again supports the validity of our model.

Estimating the switching rate

In the simulation, the rate of switching events was determined by enumerating the events in which CI numbers reach zero for the first time. ‘Brute-force' sampling was performed when the rate was fast enough (>0.0001 per generation). However, when the rate becomes smaller and the switching events become extremely rare, brute-force sampling is inefficient for generating reliable statistics. Instead, we used the Forward Flux Sampling (FFS) method (Allen et al, 2005). In this method, the frequency of rare events is sampled based on the initial flux and the probabilities of the flux reaching the successive interfaces along the reaction coordinate. We compared the results from FFS with those from brute-force sampling under the same conditions, and they agree with each other very well.

Estimating the number of CI proteins generated from each transcription burst, bCI

We can write bCI=bTXbTL, in which bTX is the average number of cI mRNAs produced per transcription burst and bTL is the average number of CI proteins made from one mRNA molecule during its lifetime. bTX was estimated from the single-cell mRNA data and found to be ≈4 (see Figure 2A). To estimate bTL, we note that the translation rate of cI mRNA is about six-fold lower than that of lacZ mRNA (Shean and Gottesman, 1992), or 0.02 s−1 (Kennell and Riezman, 1977). With cI mRNA lifetime measured to be ≈4 min (see Supplementary Figure 1), we obtain bTL≈5 and thus bCI≈20.

Gene activity in a Cro− reporter strain

We used a Cro− reporter strain, NC417 (Svenningsen et al, 2005) as a control for the validity of our stochastic model, as the difference between Cro+ and Cro− has been measured before (Schubert et al., 2007). We performed the same smFISH experiments, as well as qRT–PCR, for this strain. The strain carries the cro27 allele, bearing a missense mutation. It also carries the same temperature-sensitive cI allele as in the Cro+ reporter strain (NC416). The experiments were carried out between 30 and 40°C in 2°C intervals. Our stochastic model was extended to describe the Cro– strain, by simply setting the initial number of Cro proteins to zero and the transcription rate from promoter PR to zero. The mRNA distributions and the averaged level from experiments and simulations are summarized in Supplementary Figure 5. The good agreement between theory and experiment in the Cro− strain further strengthens our confidence in the validity of the stochastic model, and in particular the mapping between physical temperature and CI activity μ. It is noteworthy that we used μ(T) from the results of the Cro+ strain (Supplementary Figure 4).

Comparing results from smFISH, qRT–PCR and LacZ activity assay

We compared cI mRNA levels determined using smFISH with the results from qRT–PCR. The comparison between the two methods is shown in Supplementary Figure 6. The good agreement between the smFISH and qRT–PCR data (correlation coefficient ∼0.95) indicates that the quantitative measurement of mRNA using smFISH is reliable.

We also compared the mean level of cI mRNA obtained from smFISH experiments with the promoter activity reported in the literature, measured using β-galactosidase (LacZ) assay (Schubert et al, 2007). In those experiments, CI was expressed from an inducible promoter and the activity of PRM was assayed using LacZ as a reporter. The comparison is shown in Supplementary Figure 7 and our smFISH data is consistent with the LacZ assay results.

Temperature does not significantly affect the number of phages released at cell lysis

Figure 3A shows that temperature has a very small effect on the spontaneous lysis rate. We also measured the number of phages released per lysis at 30 and 40°C for both RecA− strain (JL5902) and RecA+ strain (MG1655). As shown in Supplementary Figure 2, there is no significant difference in the number of phages released per lysis between 30 and 40°C for both strains.

‘Continuous phase transition'-like behavior

As shown in the right panel of Supplementary Figure 8, there is no bimodality in the mRNA distribution at the transition temperature 36.5°C. This behavior suggests a continuous phase transition-like behavior. As the temperature increases, the CI-dominant state (left panel of Supplementary Figure 8) shifts gradually toward the direction of Cro-dominant state. At the transition temperature, the different cells ‘spread' across the region between the CI-dominant state and the Cro-dominant state.

FORTRAN source code

FORTRAN source code