Revisiting Bistability in the Lysis/Lysogeny Circuit of Bacteriophage Lambda

The lysis/lysogeny switch of bacteriophage lambda serves as a paradigm for binary cell fate decision, long-term maintenance of cellular state and stimulus-triggered switching between states. In the literature, the system is often referred to as “bistable.” However, it remains unclear whether this term provides an accurate description or is instead a misnomer. Here we address this question directly. We first quantify transcriptional regulation governing lysogenic maintenance using a single-cell fluorescence reporter. We then use the single-cell data to derive a stochastic theoretical model for the underlying regulatory network. We use the model to predict the steady states of the system and then validate these predictions experimentally. Specifically, a regime of bistability, and the resulting hysteretic behavior, are observed. Beyond the steady states, the theoretical model successfully predicts the kinetics of switching from lysogeny to lysis. Our results show how the physics-inspired concept of bistability can be reliably used to describe cellular phenotype, and how an experimentally-calibrated theoretical model can have accurate predictive power for cell-state switching.

Use of the temperature-sensitive allele: To tune CI protein level in the cell, we use the temperaturesensitive allele cI857 (Hershey, 1971). In this mutant, CI bears a GCA to ACA mutation at bp 199 (Lieb, 1981).
The precise way by which temperature dependence arises in this mutant is not known. Phenomenologically, two alternative scenarios have been invoked to describe the effect of temperature on CI in the cell: (1) active degradation of CI protein, at a temperature-dependent rate ( ) (Isaacs et al, 2003;Villaverde et al, 1993); (2) loss of functionality, such that only a fraction ( ) of CI proteins is active, i.e. able to bind DNA and regulate transcription (Zong et al, 2010;Gaitanaris et al, 1994). Importantly, both scenarios result in the same steady-state promoter activity levels if the following mapping is made: ( ) = ( ) ( ( ) ⁄ + ( )), where is the growth rate. In this work, we assume the second scenario, i.e. that the effect of temperature is to change the fraction of active CI molecules in the cell. This choice is motivated by three considerations. First, CI857 has been found to be capable of renaturing in the cell (Gaitanaris et al, 1994), an observation inconsistent with irreversible degradation. Second, a model assuming a temperature-dependent fraction of active CI successfully predicts the stability of cI857 lysogens at different temperatures (Zong et al, 2010). Third, the increased degradation scenario was incapable of quantitatively reproducing the delayed switching kinetics observed in our experiments ( Figure 3B-D).
Predicting the steady-states of the lambda switch: The steady-states are de�ined by the requirement that CI production is balanced by CI loss through dilution where ( ) is P RM activity (the amount of CI produced per generation time) and is the mean concentration of CI molecules per cell. The rate of loss is 1, since units are chosen to be generation times (cell interdivision time/ln (2)).
In our reporter strain, equation (1) must be modi�ied due to the fact that only a fraction of CI proteins, ( ), is active.
Where we have used the transformation → so that now represents the concentration of active CI molecules. This expression may be rewritten: An additional feature of our strain is the presence of the �luorescent reporter. Its steady-state is de�ined by the expression: where is the cell fluorescence and 0 is a normalization constant. We can therefore equate = ( ) at steady-state using the appropriate normalization.
We can also rewrite equation (3) as: This last relation allows us to construct the regulatory curve of the PRM promoter with respect to CI from experimental fluorescence data, i.e. plot ( ) vs. as we do in Fig. 2A. Conversely, given a theoretical fit for the dose response function ( ) we can obtain a theoretical prediction of the steady-states of the system ( ) vs. ( ), since, at steady-state, This theoretical curve is shown as the solid line in Fig. 2B, where we have plotted the predicted steady-states using the theoretical fit for the dose-response obtained in Fig. 2A.
Error propagation: To propagate errors resulting from the fitting of experimental data, we first calculated the prediction bounds and confidence intervals for the fit of the PRM versus active CI data ( Fig. 2A) to a Hill function, and for the fit of PRM activity variance versus mean to a power law (Fig. S6B,). This was done using Matlab's curve fitting toolbox. The prediction bounds represent the inference from the data that the true functional value of the fit lies within the prediction bounds with 95% certainty. To then estimate the prediction bounds for the PRM steady-states, the individual functions representing the upper and lower bounds of the PRM(CI) curve were transformed analogously to the nominal fit.
To estimate the parameter bounds for the stochastic model (see Supplemental Information Below), 100 random parameter sets were chosen from within the confidence intervals of both experimental fits above, subject to the further requirement that the mean standard error (MSE) of the fit be less than twice that of the fit given by Matlab's curve fitting algorithm. This requirement was imposed in order to closely match the prediction bounds of the individual fits. The prediction bounds for the delay time and mean switching time in our stochastic model represent the minimum and maximum values of those observables obtained through this random sampling.

Estimation of burst parameters from single-cell data:
The mean protein burst size ( ) and burst frequency ( ) were extracted from experimental data by calculating the dependence of the cell-to-cell variance 2 in PRM activity on the mean activity . This was done as follows: the single-cell fluorescence data from experiments at all temperatures ( Fig. 2C) were combined, and then binned based on the level of active CI . The variance and mean of red cell fluorescence were then approximated using a power law relation, (Fig. S6B), with = 0.68 ± 0.14 and = 5.3 ± 3.5 −1 in molecular units (see below).
This allowed us to estimate the burst size and frequency for a given level of x, using the relation = 2 2 and = 2 , (Shahrezaei & Swain, 2008;Friedman et al, 2006;Mehta et al, 2008) (Fig. S6C). The parameter was rescaled from fluorescence units to molecule number, using the estimated value of ~300 total CI molecules/cell (Levine et al, 1979;Reichardt & Kaiser, 1971).
A stochastic model for the lambda switch: We assume that the switching between lysogeny and lysis is governed by the stochastic kinetics of CI copy numbers. We therefore model these kinetics according to three assumptions (shown schematically in Fig. S6A): 1. CI is produced from PRM in geometrically distributed bursts of mean size b and frequency per generation a, both of which are functions of the instantaneous level of active CI (see Estimation of burst parameters from single-cell data).

CI is diluted by cell division.
3. The onset of lysis occurs once CI level in the cell drops below a certain threshold.
These assumptions can be captured using a master equation, which specifies how the population structure (in terms of CI copy number) evolves in time (Kampen, 1981;Mehta et al, 2008). Specifically, the master equation is written as: where is the probability of having n CI molecules/cell. is a matrix coefficient describing the transition rate from state m to state n, and is completely defined by the following 4 rules: (1) The value of for = + is ( ) × ( ; � is a geometric distribution with success probability (2) The value of for = − 1 is . CI is diluted in proportion to growth rate (which is 1 in the units chosen) and copy-number.
(3) The value of for < − 1 is 0. CI copy-number cannot instantly decrease by more than 1.
(4) In order to conserve probability ≡ − ∑ ≠ i.e. the flux of probability into all other states from state is the flux out of .
By using the experimentally determined burst parameters (see Estimation of burst parameters from single-cell data), we were able to fully parameterize the master equation above. The solution of the master equation at arbitrary times was found by numerical matrix exponentiation: ( , ) = ( , 0), utilizing the expm function of Matlab, with a cutoff applied for very high values of (>2.5x the amount of CI produced per generation if the gene is constantly on). The initial distribution ( , 0) was obtained by evolving a negative binomial distribution for 100 generations, under the initial condition parameters (100% of CI being active.) Finally, in order to determine the fraction of lysogenic cells left at any given time, a threshold value of active CI was chosen based on a value that corresponds to the experimentally chosen partition between states (see Estimating pleotropic effects of changes in temperature: In analyzing our results, we make the simplifying assumption that the effect of changing temperature is only to vary the fraction of active CI in the cell. Other physiological effects resulting from the change in temperature are neglected. To numerically test the validity of this approximation, we refer to the results of (Herendeen et al, 1979), who examined the change in abundance of 100 E. coli proteins in the temperature range 13.5-46℃. They found that the majority of proteins measured displayed <20% variation in abundance from 30-42℃, while >95% displayed less than 60% variation. Assuming this is a representative sample, and that this variation is linear in temperature change, we can determine the limits of what is likely to occur in the case of CI (Fig. S3). P RM activity was allowed to vary by 30% from the nominal fit (solid line in Fig. 2A and Fig. S3A) when temperature is changed from 34℃ to 40℃ (effectively the worst-case scenario, since the temperature range is ½ that of the 30-42℃ range of the study). This variation is shown in the two dashed lines in Fig. S3A. This variability was then propagated to the predicted steady states of the system shown in Fig. 2B where the dashed lines are again the estimated limits of temperature variation. As seen in Fig. S3B, the estimated limits of temperature variation are smaller than the experimental uncertainty in the fit (shaded region in Fig. 2A Fig. 2A)              B. Figure S7