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.
Bacterial lysogeny serves as a simple paradigm for cell differentiation.
We characterize the activity of the fate‐determining genes, cI and cro, with single‐molecule resolution.
Stability of the lysogenic state is found to depend in a simple manner on the frequency of activity bursts from cI.
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).
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.
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).
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).
Source data for Figure 3B [msb201096-sup-0001-SourceData-S1.xls]
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.
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
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.
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.
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)n, 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.
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−T ln(2)/τ and 1−e−t ln(2)/τ 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.
We thank B Egan, A Raj, A van Oudenaarden, J Little, D Court, I Dodd, M Elowitz, P Ge, G Altan‐Bonnet, T Gregor and all members of the Golding lab for supplying reagents and for their advice. We thank U Alon for commenting on an earlier version of this paper. Some of the CI mutants were constructed by IG while working at the lab of Ted Cox (Princeton). Work in the Golding lab is supported by NIH grant R01GM082837, HFSP grant RGY 70/2008 and NSF Grant 082265 (PFC: Center for the Physics of Living Cells).
Author contributions: IG and CZ conceived the stability measurement project. CZ performed the majority of experiments and developed the theoretical model for lysogen stability. LhS, LAS and SOS performed additional experiments and developed analysis tools. IG, CZ, LhS, LAS and SOS wrote the paper.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Material [msb201096-sup-0001.doc]
FORTRAN source code
FORTRAN source code [msb201096-sup-0002.zip]
This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.
- Copyright © 2010 EMBO and Macmillan Publishers Limited