Abstract
The inherent stochastic nature of biochemical processes can drive differences in gene expression between otherwise identical cells. While celltocell variability in gene expression has received much attention, randomness in timing of events has been less studied. We investigate event timing at the singlecell level in a simple system, the lytic pathway of the bacterial virus phage λ. In individual cells, lysis occurs on average at 65 min, with an s.d. of 3.5 min. Interestingly, mutations in the lysis protein, holin, alter both the lysis time (LT) mean and variance. In our analysis, LT is formulated as the firstpassage time (FPT) for cellular holin levels to cross a critical threshold. Exact analytical formulae for the FPT moments are derived for stochastic gene expression models. These formulae reveal how holin transcription and translation efficiencies independently modulate the LT mean and variation. Analytical expressions for the LT moments are used to evaluate previously published singlecell LT data for λ phages with mutations in the holin sequence or its promoter. Our results show that stochastic holin expression is sufficient to account for the intercellular LT differences in both wildtype phages, and phage variants where holin transcription and the threshold for lysis have been experimentally altered. Finally, our analysis reveals regulatory motifs that enhance the robustness of lysis timing to cellular noise.
1. Introduction
Phenotypic variation is an intrinsic feature of biological systems. Most phenotypic variation stems from genetic or environmental variation, but some phenotypic variation can exist even for identical genotypes in the same environment. This variation often emerges from the stochastic nature of biochemical reactions within a cell. Small fluctuations in transcriptional and translational machinery can lead to large variations in important regulatory molecules, impacting cellular switches, pathways and, ultimately, cell fates [1–5].
The causes and effects of stochastic variation have been best explored in microbial systems such as bacteriophages [6–8], bacteria [9–13], yeast [14–19] and human viruses [20–23]. For example, fluctuations in gene expression make it possible for populations of isogenic λ phages in the same environment to take one of two distinct developmental fates: immediate lysis of the host cell or lysogeny as a prophage in the host genome. This developmental switch is controlled by cytoplasmic levels of bacteriophage λ protein CII [24]. When CII levels are high, the lysogenic pathway is favoured, but when CII levels are low, phages will enter the lytic pathway. While the switch is essentially stochastic, the probabilities of lysis or lysogeny can vary according to environmental conditions and the multiplicity of infection [8,25–27].
Bacteriophage λ also possesses a simpler regulatory network that determines the timing of lysis. Lysis is accomplished by a suite of four genes: S (encodes holin and antiholin), R (encodes endolysin) and Rz/Rz1 (encode spanins). Expression of these genes, and all λ structural and assembly proteins, is under the control of the λ late promoter pR′. While expression from promoter pR′ is constitutive, transcription terminates upstream of the lysis genes [28]. When the ‘decision’ is made to pursue the lytic life cycle, protein Q is produced and acts to prevent premature termination of late mRNA transcription approximately 10–15 min after lysis induction [29,30]. Following the production of the late mRNA transcript, the lysis proteins are translated by host ribosomes. Holin accumulates evenly in the Escherichia coli plasma membrane until a genetically encoded time when they spontaneously aggregate to form large holes [31]. These holes allow endolysin and the spanins to access and disrupt their targets beyond the periplasmic space [32]. Subsequently, the cell ruptures and phage progeny are released into the surrounding medium.
As hole formation and cell rupture are nearly simultaneous, lysis timing depends on the accumulation of holin in the cell membrane up to a critical threshold [33]. Dennehy and Wang observed the lysis time (LT) of single λ lysogens using a temperaturecontrolled perfusion chamber mounted on a microscope stage (see electronic supplementary material, figure S1 and movie). Their measurements show that λ lysis occurred on average at 65 min, with an s.d. of 3.5 min (figure 1). Interestingly, mutations in the holin amino acid sequence can significantly change the LT mean and variance [34]. As randomness in the onset of the late promoter expression cannot account for these variations [34,35], stochastic holin expression and accumulation may be key in determining LT differences. Ultimately, the expression and accumulation of holin in the inner membrane depend on several factors: the rate of transcription of late mRNA, the rate of holin protein translation by host ribosomes, the rate of holin insertion into the plasma membrane and the threshold holin concentration required for membrane hole formation.
The impact of these factors on LT variation is investigated using a stochastic gene expression model. LT is modelled as the firstpassage time (FPT) for cellular holin copy numbers to cross a critical threshold in the stochastic model (figure 2). Exact analytical formulae for the FPT statistical moments predict how LT mean and coefficient of variation (CV) change in response to experimental alterations of the holin transcription efficiency, translational efficiency and the threshold concentration of holin required for lysis. Some of these predictions are confirmed through analyses of LT distribution data obtained by Dennehy & Wang [34] for wildtype phage, and phages carrying mutations in either the holin sequence or the pR′ promoter. We discuss the impact of these findings on the evolution of bacteriophage lysis timing as well as the implications for the study of variation in gene expression.
2. Stochastic gene expression model formulation
We denote by x(t), the number of holin molecules inside the cell at time t, where t = 0 corresponds to the onset of holin expression. As holin does not degrade over relevant timescales [36], x(t) increases over time. The goal is to compute the FPT (time taken by a random walker to first reach a defined point) 2.1to the threshold X ≥ 1 assuming x(0) = 0. The CV^{2}, defined as 2.2is used to quantify the randomness in FPT, where the symbol denotes expected value. Let, T_{i} be independent and identically distributed random variables denoting the time interval between the i − 1th and ith mRNA transcription event from the promoter. We consider gene expression in the burst limit: transcribed mRNAs degrade instantaneously after producing a translational burst of protein molecules [37–40]. Consistent with measurements [41], we assume protein bursts are geometrically distributed with a mean translational burst size b = k_{p}/γ_{m}. Here, k_{p} denotes the protein translation rate per mRNA and γ_{m} is the mRNA degradation rate. Thus, protein count 2.3is an accumulation of independent, identical geometric bursts B_{i}, n is the number of transcription events in [0, t] and the interburst arrival time distribution is given by T_{i}.
3. Analytical calculation of firstpassage time moments
First, the minimum number of transcriptional events N required for x(t) to reach the threshold X is quantified. Let, 3.1
Then, the statistical moments for N are given by (see electronic supplementary material, appendix C) 3.2
The average number of transcriptional events required for lysis is the lysis threshold (X) divided by the mean translational burst size (b) plus one. The plus one arises from the fact that, even if the threshold is close to zero, at least one transcription event is required for x(t) ≥ X. Given N (the number of burst events) and T_{i} (the time interval between bursts), time taken to reach the threshold is 3.3Below, distribution for T_{i} is determined for two different mRNA production processes: constitutive and bursty transcription.
3.1. Constitutive transcription model
Assuming mRNA production is a Poisson process with rate k_{m}, T_{i} is an exponential random variable with mean 1/k_{m}, and (3.3) yields 3.4where is the FPT CV^{2}. For (number of transcriptional events required for lysis is large), (3.4) simplifies to 3.5Note that under similar assumptions, we obtain 3.6for deterministic arrival of bursts (i.e. ) in (3.3). Thus, for , randomness in the FPT is halved for fixed time intervals between protein bursts compared with exponentially distributed intervals. When , we obtain the Poisson limit . The results above show that for a given X, and can be independently controlled via holin translation rate (k_{p}) and transcription rate (k_{m}). Desired FPT mean and CV can be set by choosing 3.7Equation (3.5) predicts how FPT CV^{2} changes with mean as k_{m}, b or X are experimentally altered (figure 2):

— is inversely related to for perturbations in the lysis threshold. More specifically, increasing X will increase but decrease .

— For alterations in translational efficiency k_{p} (and hence, the mean translational burst size b), increase in comes with a decrease in .

— Changes in the transcriptional efficiency k_{m} alter without affecting .
Next, a more sophisticated mRNA production process is considered, where the promoter switches between different transcriptional states.
3.2. Stochastic promoter switching model
Promoter transitions between active and inactive states have been well documented in various genes [42–46]. In the stochastic formulation, the promoter resides in active and inactive states for exponentially distributed time intervals with means 1/k_{off} and 1/k_{on}, respectively, where k_{on} and k_{off} represent transition rates (figure 3, top right). From the active state, mRNA production is a Poisson process with rate k_{m}. As before, we assume transcribed mRNAs degrade instantaneously after producing a geometrically distributed burst of protein molecules of mean size b. Let, T_{i} be independent and identically distributed random variables denoting the time interval between the i − 1th and ith transcriptional event from the active state. Then, the mean and variance for the interburst time interval is given by 3.8for all i ≥ 1 (see electronic supplementary material, appendix D). Using (3.3), (3.8) and assuming , we obtain 3.9aand 3.9b
As expected, (3.9) reduces to (3.5) for constitutive transcription, i.e. and promoter is always in the active state.
In the absence of promoter switching, decreasing transcriptional efficiency is predicted to increase without affecting (figure 2). However, inclusion of stochastic promoter switching complicates the relationship between and . Qualitative versus trends for perturbations in the transcriptional efficiency (measured by either k_{m}, k_{on} or k_{off}) are shown in figure 3. Interestingly, modulating the transcriptional efficiency can have different effects on depending on the underlying parameter that is being altered. For example, decreasing k_{on} or k_{m} will both increase , but the former is predicted to increase , while the latter decreases (figure 3). As in the previous section, FPT CV^{2} changes inversely with the mean for alterations in X or b.
4. Measured lysis time variation consistent with model predictions
If LT variation is indeed a result of noisy holin expression, then randomness in event timing should vary with mean as predicted by (3.5) or (3.9), i.e. should decrease with increasing as the lysis threshold is perturbed. Before we test this prediction, quantification of and from singlecell LT measurements is discussed.
4.1. Connecting lysis time to firstpassage time
LT can be related to the FPT as LT = FPT + T_{pR′}, where T_{pR′} is the time for the onset of transcription from the late promoter. As T_{pR′} is shown to be uncorrelated with FPT [35], 4.1where σ_{I} denotes the s.d. in , and is taken to be 15 min [29,30]. As T_{pR′} variation only accounts for a small fraction of LT variability [34,35], we assume . Nonzero values of do not change the versus trend described below. For wildtype λ (see electronic supplementary material, appendix B), 4.2which implies from (4.1), 4.3As we ignore T_{pR′} variation in these calculations, the reported values for are an upper bound. Assuming constitutive transcription from the wildtype pR’ and a lysis threshold X of 1500 holin molecules [36], we obtain the following estimates for the translational burst size and the transcription rate from (3.5) and (4.3) 4.4
Noisy holin expression with these values is sufficient to generate intercellular LT differences similar to data (figure 4). These parameter values are consistent with previous estimates of holin transcription/translation efficiencies measured using bulk assays [36].
4.2. Mutations in the holin sequence
Mutations in the holin amino acid sequence can both increase and decrease LT mean and s.d. compared with wildtype [34]. These mutations may alter holin selfaffinity or membrane insertion rate, thus affect the total cellular holin threshold required for lysis. For example, mutant holin with a stronger membrane binding affinity than the wildtype will have a larger fraction of the expressed holin present in the plasma membrane. Consequently, the critical holin concentration required for hole formation will be reached with less total holin expression. In our model, this corresponds to having a lower threshold X, and hence, a shorter LT compared with the wildtype.
Singlecell LT measurements for six holin amino acid sequence mutations (see electronic supplementary material, appendix B) reveal that shorter LTs are associated with higher a CV^{2} (figure 5). Thus, as predicted by a stochastic model of holin expression, results show an inverse relationship between and for changing X.
It is important to point out that the inverse correlation between and in figure 5 is inconsistent with a model where LT differences arise due to celltocell variations in holin production rate. Consider an alternative model, where holin molecules are produced at a rate k, in which k has a certain distribution across the cell population. Celltocell differences in k would reflect intercellular variation in the abundances of host ribosomes and/or RNA polymerases, which arise due to differences in cell volume, cell cycle or stochastic expression of these enzymes. As holin levels build up at different rates inside individual cells, some cells would lyse earlier than others generating a LT distribution similar to data. However, in this case, 4.5aand 4.5band is independent of X. Hence, is predicted to be uncorrelated with for alterations in X. The inverse trend between and obtained from experimental measurements (figure 5) suggests that LT differences are due to stochasticity in holin expression rather than difference in holin production rates between cells.
4.3. Mutations in the late promoter
Previous work measured LT variation for four λ variants carrying mutations in the late promoter pR′ TATA box region [34]. Only two of these mutants showed a significant decrease in promoter strength, which increased both the LT mean and CV^{2} (see electronic supplementary material, appendix B). For example, 4.6for promoter mutant SYP208 [34]. Note a twofold increase in , and a sixfold increase in , compared with their respective values in wildtype λ (see (4.3)). The observed increase in cannot be explained by a reduced transcription rate k_{m} in the constitutive model, as in this case CV is invariant of k_{m} (figure 2). However, as discussed below, these results can be explained by a model that includes random promoter transitions between active and inactive states (figure 3).
Results in §3.2 show that in the stochastic promoter switching model (figure 3, top right), reducing transcription from the active state (decreasing k_{m}) increases but decreases . By contrast, enhancing the stability of the inactive state (decreasing k_{on}) increases both and . Finally, for decreasing activestate stability (increasing k_{off}), exhibits an inverted Ushaped profile with increasing . Thus, an increase in both and in response to a reduced holin transcription is consistent with promoter mutations decreasing k_{on} or increasing k_{off} in the promoter switching model (figure 3). A recent study has shown that promoter mutations in the TATA box region primarily modulate the transcription burst frequency k_{on} [48]. Assuming constitutive transcription from the wildtype pR′ , mutations that significantly decrease k_{on} will induce additional promoter switching noise in holin expression. Consequently, weaker nonconstitutive late promoter expression makes cell lysis longer, and more unpredictable compared with wildtype λ.
5. Discussion
The timing of cellular events, such as mitosis, apoptosis, developmental differentiation and reversion from latency, frequently depends on the concentrations of essential regulator molecules. Stochastic gene expression causes variation in regulator molecule concentrations even for identical cells induced at the same time. Consequently, molecular threshold concentrations will be reached at different times in different cells, giving rise to significant variation in the timing of critical cellular events across a homogeneous cell population. The holin lysis system of bacteriophage λ provides a simple system to study how event timing is regulated at the singlecell level. The main contributions of this work are as follows:
(i) LT was determined by the FPT for cellular holin levels to reach a critical threshold. Exact analytical formulae for the FPT mean and CV^{2} were derived for different stochastic gene expression models. These formulae are given by (3.4) and (3.9) for the constitutive transcription model and the stochastic promoter switching model, respectively.
(ii) FPT derivations reveal that for a fixed lysis threshold, LT CV^{2} is minimized by reducing the translation burst size b. As discussed later, mechanisms for reducing b are embedded in λ's lytic pathway in order to ensure precision in lysis timing.
(iii) Using (3.4) and (3.9), predictions on how biologically meaningful parameters (transcription rate, translation rate, lysis threshold, etc.) impact LT statistics were generated. These predictions are summarized in figure 2 (for the constitutive transcription model) and figure 3 (for the stochastic promoter switching model).
(iv) Analysis of data from [34] confirms model predictions. In particular, when the threshold for lysis is altered through holin sequence mutations, increases in mean LT are associated with decreases in LT CV^{2}. By contrast, when holin transcription is altered through holin promoter mutations, increases in mean LT are associated with increases in CV^{2}.
Below, we discuss the implications of some of these findings in further detail.
5.1. Independent control of lysis time mean and variation
Formulae for the FPT moments reveal how transcriptional/translational efficiencies independently modulate the mean and variation in lysis timing. For constitutive transcription, FPT CV^{2} increases linearly with the translational burst size b, but is independent of the transcription efficiency (figure 2). Precision in the FPT can be achieved by setting b sufficiently small. When , the limit of noise suppression for a given lysis threshold X is attained: . This limit corresponds to a Poisson arrival process for holin production. Once b is chosen for a desired , holin transcription efficiency can be adjusted to have any mean FPT.
Note that b = k_{p}/γ_{m} is proportional to the translation efficiency k_{p}, where γ_{m} is the mRNA degradation rate. In λ, all lysis and structural proteins are translated from the same polycistronic mRNA transcript. Hence, reducing b by making k_{p} small is a more viable option than increasing γ_{m}, as the latter would affect the translational burst sizes of structural proteins. S gene, which encodes holin, is reported to have a low translational efficiency [36]. As discussed later, holin translation bursts are further reduced through production of a second protein, antiholin, from the S mRNA. Given the direct relationship between and b, mechanisms in place for minimizing b are essential for buffering randomness in lysis timing.
Our results on the FPT are complementary to previous work quantifying protein copy number variations in stochastic gene expression models. These studies have also shown that a reduced translation efficiency is critical for minimizing stochastic fluctuations in protein levels [13,38]. However, there are qualitative differences in the formulae for the FPT and protein copy number statistical moments. For example, the FPT CV^{2} is proportional to the translation efficiency k_{p}. By contrast, for , the steadystate protein copy number CV^{2} is independent of k_{p} [13,38,49].
5.2. An incoherent feedforward circuit minimizes lysis time stochasticity
The λ S gene encodes two proteins with opposing function: holin and antiholin [50–52]. These proteins are produced in 2 : 1 ratio, i.e. for every two holins, there is one antiholin molecule. These two proteins differ only in that antiholin has an extra two residues, methionine and lysine, at the Nterminus. The positively charged lysine prevents antiholin from accessing holin's characteristic three transmembrane domain configuration in the E. coli inner membrane [53]. Experimental evidence suggests that antiholin binds holin and sequesters it in the cytoplasm, preventing it from participating in hole formation [54,55]. This creates an incoherent feedforward circuit in λ's lytic pathway [56].
Preliminary analysis provides some clues on the functioning of this feedforward circuit. Let, holin and antiholin translation from the mRNA be Poisson processes with rates k_{a} and k_{ah}, respectively. Then, the probability of producing i holin and j antiholin molecules from the mRNA in time interval [0, t] will be 5.1
We denote by x and y, the total number of holin and antiholin molecules synthesized by an individual mRNA before it decays. Random variables x and y would be correlated, because an mRNA that lives longer than average by random chance would create more of both holin and antiholin. Assuming each mRNA lives for an exponentially distributed time interval with mean 1/γ_{m}, the joint probability distribution of x and y is 5.2Consider holin–antiholin interaction in a limiting case: antiholin rapidly binds to holin making holin inactive in an irreversible reaction. Then the active holin translational burst size (number of active holin made per mRNA) would be x − y (if x > y) or 0 (if x ≤ y). Equation (5.2) yields the following mean translational burst size b for active holin 5.3aand 5.3bwhich decreases with increasing antiholin production rate k_{ah}. Recall from (3.5), that for a fixed X and , is minimized by reducing b (see equation (3.5)). Thus, antiholin may enhance precision in lysis timing by lowering b, and making the active holin production process less bursty. Removing antiholin (k_{ah} = 0) would increase b, which is predicted to decrease but increase (equation (3.5)). Indeed, λ variant IN62 (no antiholin expression) lyses 10 min faster but with a twofold higher compared with wildtype [34]. Our results show that as in other systems [57–59], antiholinmediated incoherent feedforward loop is a key regulatory motif in λ's lytic pathway for buffering LT from stochastic gene expression.
5.3. Evolution of lysis timing
In lifehistory theory, the timing of lysis is viewed as a tradeoff between present and future reproduction. Assuming a linear rate of progeny assembly [60], postponing lysis increases progeny number, but delays infection of new hosts. In other words, increasing fecundity also increases generation time. Several studies have theoretically and experimentally examined the optimal timing of lysis [60–66]. These studies reach similar conclusions: when host density and quality are high, early lysing genotypes will be selectively favoured as they will have more grandchildren than competing genotypes. Conversely, poor host quality and low host densities should favour longer LTs.
We expect phage fitness will be maximized if the LT CV is minimized around the optimum LT for any particular environment. Indeed, we find that wildtype λ has a lower LT CV and a higher fitness than all laboratorygenerated mutant strains [34,60]. Furthermore, λ's lytic pathway exhibits the ingredients required for precision in timing: low S gene translational efficiency [36] and antiholin synthesis to further reduce the active holin translational burst size. At first, production of antiholin seems energy inefficient as it inactivates holin, translated from the same mRNA. Our analysis and data suggest that antiholin may reduce LT variations and is reminiscent of the classical tradeoff between energy efficiency and noise buffering in cellular systems [67,68].
Our results predict how the mean LT should evolve under selection, while minimizing LT CV. In principle, mean LT can be modulated through mutations in the following regions: the pR′ promoter sequence, the holin leader sequence and the holin coding sequence. Because of the polycistronic nature of the pR′ promoter sequence, we expect it can be ruled out. Presumably, mutations in pR′ will have undesirable pleotropic effects on the expression of λ's structural and regulatory proteins. We predict that selection for shorter LTs (i.e. when host quality and density is high) will favour mutations that increase translational efficiency, and not those that lower the threshold concentration of holin required for lysis (X). From (3.5), we see that, when translational burst size (b) is small, CV^{2} ≈ 1/X, and more sensitive to changes in X rather than b. So decreasing mean LT by lowering X will increase CV^{2}. By contrast, selection for longer LTs will favour mutations that increase X because this will decrease LT CV^{2}, while mutations decreasing translational efficiency will not substantially change CV^{2}. These predictions can be tested by genetic sequencing of λ phage following experimental evolution under high or low host densities.
5.4. Burst size implications of lysis timing variation
Given that progeny are assembled at a linear rate, we expect that differences in LT will be reflected by differences in progeny burst size (total progeny produced in a single infection). Consequently, we should observe a direct relationship between LT stochasticity and burst size stochasticity. Seminal work by Delbrück showed that phage burst size can vary significantly from celltocell, even after controlling for cell size [69]. Interestingly, measured viral burst size CV is an order of magnitude larger than the LT CV (see electronic supplementary material, appendix E). Recent singlecell measurements also show a broad viral burst size distribution for vesicular stomatitis virus (VSV), an animal virus that can occasionally infect humans. Data show a 300fold celltocell difference in the number of VSV progeny produced per cell, and these differences cannot be attributed to viral genetic variation [70,71]. Phage λ provides an ideal system to study how viral burst size variations arise from LT variations and noisy expression of structural genes. Understanding how nongenetic factors affect viral fitness at the singlecell level has important implications for health and disease.
6. Summary
In summary, celltocell LT differences arise due to the stochastic transcription and translation of the lysis timekeeper, holin. We derive analytical formulae for the mean and CV^{2} of the time it takes for holin to reach a critical concentration in the inner membrane. These formulae are used to generate predictions regarding the impact of biologically relevant parameters, such as transcription and translation rates, on LT mean and CV^{2}. We evaluate these predictions using previously published data on singlecell variation in LT [34]. Our analysis confirms that stochastic holin expression can account for variation in LT. Furthermore, we find that the holin lysis system contains several features designed to minimize noise.
Funding statement
J.J.D. was supported by grants from the Professional Staff Congress of the City University of New York and the National Science Foundation (Division of Environmental Biology Awards 0804039 and 1148879, and Division of Molecular and Cellular Biosciences Award 0918199). A.S. was supported by the National Science Foundation grant no. DMS1312926, University of Delaware Research Foundation (UDRF) and Oak Ridge Associated Universities (ORAU).
Acknowledgements
We thank Andrew Rutenberg, Gillian Ryan, IngNang Wang, PakWing Fok, Pavol Bokes, Ryan Zurakowski, Khem Ghusinga and Ryland Young for discussion on some of the ideas contained herein.
 Received February 6, 2014.
 Accepted March 14, 2014.
 © 2014 The Author(s) Published by the Royal Society. All rights reserved.