Stochastic Kinetic Analysis of Developmental Pathway Bifurcation in Phage -Infected Escherichia coli Cells

Corresponding author: Harley H. McAdams, Department of Developmental Biology, Stanford University School of Medicine, Stanford, CA 94305., mcadams{at}cmgm.stanford.edu (E-mail).

Communicating editor: R. S. HAWLEY

	ABSTRACT

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

Fluctuations in rates of gene expression can produce highly erratic time patterns of protein production in individual cells and wide diversity in instantaneous protein concentrations across cell populations. When two independently produced regulatory proteins acting at low cellular concentrations competitively control a switch point in a pathway, stochastic variations in their concentrations can produce probabilistic pathway selection, so that an initially homogeneous cell population partitions into distinct phenotypic subpopulations. Many pathogenic organisms, for example, use this mechanism to randomly switch surface features to evade host responses. This coupling between molecular-level fluctuations and macroscopic phenotype selection is analyzed using the phage lysis-lysogeny decision circuit as a model system. The fraction of infected cells selecting the lysogenic pathway at different phage:cell ratios, predicted using a molecular-level stochastic kinetic model of the genetic regulatory circuit, is consistent with experimental observations. The kinetic model of the decision circuit uses the stochastic formulation of chemical kinetics, stochastic mechanisms of gene expression, and a statistical-thermodynamic model of promoter regulation. Conventional deterministic kinetics cannot be used to predict statistics of regulatory systems that produce probabilistic outcomes. Rather, a stochastic kinetic analysis must be used to predict statistics of regulatory outcomes for such stochastically regulated systems.

IN MCADAMS and ARKIN 1997 , we analyzed properties of a representative single bacterial genetically coupled link, that is, a configuration where one promoter controls a gene whose protein product regulates another promoter. In that analysis, an integrated molecular-level model of the mechanisms controlling gene transcription and translation was developed, and the expected time pattern of protein production from the controlled gene was investigated using the stochastic formulation of chemical kinetics (GILLESPIE 1977 , GILLESPIE 1992B ). The results suggested that the stochastic fluctuations of the reaction rates of gene expression reactions can produce a highly erratic time pattern of protein production in each individual cell and a wide diversity of protein concentrations across a cell population at any instant of time (MCADAMS and ARKIN 1997 ). ("Stochastic" is used here in the technical sense of "arising from a random process.")

When the protein involved is a regulatory protein, these fluctuations in concentration from cell to cell cause dispersion in the time to complete regulated events in different cells, for example, different times to complete regulatory cascades. A particularly interesting case occurs when two independently produced regulatory proteins competitively control a developmental switch. The independent, stochastic temporal patterns of production of each regulatory protein can vary widely from cell to cell. In this case, the path choice from the competitively regulated switch would not be deterministic. Rather, the choice would be random with the probabilities of alternative choices dependent on the stochastic properties of the gene expression mechanisms and the design of the switch circuit. As a result an initially homogeneous cell population would partition into subpopulations following different pathways. The phenotypes on each path could be radically different. In many pathogenic organisms random variation of surface features assists in evasion of host defenses or otherwise enhances virulence (PUTTE and GOOSEN 1992 ; ROBERTSON 1992 ; FINLAY and FALKOW 1997 ; STRAUSS and FALKOW 1997 ). We suggest in this article that one source of the randomness expressed in the phenotype variations can be the random thermal fluctuations in the reaction rates of the chemical reactions comprising the regulatory circuit.

To examine this phenomenon, we analyze herein the effect of fluctuations in gene expression rates and other molecular-level fluctuations on lysis or lysogeny pathway selection statistics by phage -infected Escherichia coli cells. This path selection is made by the lysis-lysogeny decision circuit wherein a well-characterized competitive regulatory mechanism is central to the regulatory circuit that partitions the population between lytic and lysogenic outcomes.

	APPROACH

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

Numerous studies have shown that the fraction of -infected E. coli cells that become lysogenic is influenced by environmental parameters, especially the nutritional state of the cell and the ratio of phage particles to cells at the time of infection. [We follow KOURILSKY 1973 and call this ratio the average phage input (API). The API is distinguished from the number of post-infection phage particles in a specific cell which we call the multiplicity of infection (MOI)]. We investigate here (i) the molecular mechanisms that cause the lysis-lysogeny decision circuit (described below) to select randomly different pathways in different cells, and (ii) use of a stochastic kinetic model of the circuit to predict the fraction selecting each pathway. The predicted fraction of lysogens over a range of APIs is compared with the fraction assayed by KOURILSKY 1973 . The Kourilsky experiments were selected for comparison because his experimental method was designed to minimize error in the percent lysogenization measured, and the measurements were made over a two-decade API range.

Our approach to kinetic analysis of the lysis-lysogeny decision outcome is as follows: The cell-level kinetic model of the phage lysis-lysogeny decision circuit (details below) uses the stochastic formulation of chemical kinetics (GILLESPIE 1976 ), includes stochastic mechanisms of gene expression (MCADAMS and ARKIN 1997 ), and models promoter regulation using the statistical-thermodynamic approach described by SHEA and ACKERS 1985 . The set of coupled stochastic equations comprising the stochastic kinetic model is solved using a standard Monte Carlo algorithm (GILLESPIE 1976 ) for systems of coupled chemical equations to predict the fraction of cells that will commit to lysogeny for various MOI. Over a wide API range, infecting phage particles distribute randomly among the target cells to produce a Poisson distribution of phage particles per cell as predicted theoretically (ELLIS and DELBRUCK 1939 ). Thus, given the predicted probability of lysogeny in individual cells at different MOIs, the predicted percentage of lysogens in an experimental cell population can be computed as the Poisson-weighted average of cell level predictions. This estimate is shown to compare favorably with experimental results.

	STOCHASTIC KINETICS

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

When concentrations of the reacting species are low and reaction rates are slow, conventional deterministic chemical kinetics may not describe the development of systems of coupled reactions correctly (MCQUARRIE et al. 1964 ; ZHENG and ROSS 1991 ). Rather, for such chemical systems, one has to recognize that the individual chemical reaction steps occur discretely and are separated by time intervals of random length. Bacterial genetic regulatory mechanisms typically involve low intracellular concentrations of the reacting species and relatively slow reaction rates; the concentrations are low because the majority of regulatory molecules are produced in low quantities per cell (GUPTASARMA 1995 ) and most individual genes (and hence their regulatory regions where the regulatory molecules bind) are present in only one or two copies per cell. The rates of genetic reactions are frequently so slow that many minutes may separate individual transcripts successfully initiated (i.e., neglecting abortive initiations) from an activated promoter. The expected erratic pattern of protein production cited earlier results from production of random numbers of proteins in quick succession from these intermittently produced transcripts (MCADAMS and ARKIN 1997 ). [Note that the erratic time pattern, which we postulate for protein production from the operons comprising the decision circuit, will result so long as the mechanisms of transcript initiation and translation control have the following broad statistical characteristics: (i) the statistical distributions of intertranscript intervals and proteins per transcript are skewed and have long tails, and (ii) the mean intertranscript time interval is relatively long (MCADAMS and ARKIN 1997 ).]

Kinetics of conventional macroscopic coupled chemical reaction systems is modeled using systems of ordinary differential equations, and there is an implicit assumption of continuously varying chemical concentration and deterministic dynamics. Two critical characteristics of chemical systems compatible with these assumptions are: (i) that the number of molecules of each type in the reaction mix is large compared to thermal fluctuations in concentration, and (ii) for each type of reaction in the system, the number of reactions is large within each observation interval. For genetic circuits both of these presumptions are frequently invalid so that the deterministic approach to chemical kinetics breaks down.

In small, low-rate chemical systems it is necessary to pay attention to the fact that changes in chemical population levels really occur in integral numbers of molecules, and are occasioned by essentially random distinct reaction events. It has been shown that the time evolution of such a chemical system is a stochastic process of the Markov type (GILLESPIE 1992A , GILLESPIE 1992B ; VAN KAMPEN 1992 ) that is described by the chemical master equation (GILLESPIE 1992B ; VAN KAMPEN 1992 ). The solution to the chemical master equation can be calculated using a stochastic simulation algorithm devised by GILLESPIE 1976 , GILLESPIE 1977 , GILLESPIE 1992B . The Gillespie algorithm correctly accounts for the fluctuations that occur in a well-stirred chemically reacting system.

	PHAGE LAMBDA DECISION CIRCUIT

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

The regulatory mechanisms controlling the phage lysis or lysogeny decision are generally known (HERSKOWITZ and HAGEN 1980 ; FRIEDMAN and GOTTESMAN 1983 ; ECHOLS 1986 ; FRIEDMAN 1992 ; PTASHNE 1992 ; MCADAMS and SHAPIRO 1995 ). The mechanisms of the "switch" that locks the phage into one or the other of the alternate pathways are described by MEYER et al. 1980 , MEYER and PTASHNE 1980 , PTASHNE 1992 , and SHEA and ACKERS 1985 . It is established that the stability of the CII protein is an important determinant of the number of lysogens produced in an infected cell population; that the HflA and HflB proteins affect CII degradation; and that CII is stabilized somehow in the presence of the CIII protein (HOYT et al. 1982 ; RATTRAY et al. 1984 ). The integrated functions of the switch together with other coupled regulatory mechanisms that determine commitment to execution of the lysogenic and lytic pathways are described by MCADAMS and SHAPIRO 1995 . Immediately after lambda phage infection of a target E. coli population, the cell population partitions into lytic and lysogenic phenotypes following mutually exclusive regulatory pathways.

The core of the lysis-lysogeny decision circuit is the four-promoter, five-gene regulatory network shown in Figure 1A. The organization of the genetic elements of the decision circuit in the phage DNA is shown in Figure 1B. Reinforcement of the path commitment and initiation of the pathway-specific actions associated with the selected pathway are accomplished by other coupled genes not shown (HERSKOWITZ and HAGEN 1980 ; MCADAMS and SHAPIRO 1995 ). The circuit shown in Figure 1A has two key subsystems: (i) the P_RM-based switch that creates the circuit's bistability, and (ii) the Hfl proteolytic system, which integrates environmental signals into the circuit's behavior.

View larger version (34K): In this window In a new window Download PPT slide   Figure 1. The phage lysis-lysogeny decision circuit. (a) Bold horizontal lines indicate stretches of double-stranded DNA. Arrows over genes indicate direction of transcription. Dashed boxes enclose operator sites that comprise a promoter control complex. The three operator sites, OR1–3, of the "lambda switch" implement concentration-dependent logic controlling promoters PRM and PR. Cro and CI dimers bind to the three sites with different affinities and in opposite order to control the activation level of the PRM and PR promoters (PTASHNE 1992 ; SHEA and ACKERS 1985 ). The five boxes R1–R5 contain nongenetic protein reaction subsystems. In R1, R2, and R5, "deg" indicates degradation. When protein N is available, transcribing RNAPs can be antiterminated at the NUTR and NUTL sites; termination sites TR1 and TL1 are inoperative for antiterminated RNAPs. The CI dimer acts as either a repressor or activator of promoter PRM, depending on its concentration. See text for discussion of the proteases labeled as P1 and P2 in R3 and R4. (b) decision circuit DNA organization. Phage-encoded genetic elements of the decision circuit are located in a 5000 nucleotide region of the phage DNA. Genes are separated onto leftward and rightward transcribed strands as indicated by the arrows. Rightward extensions of the antiterminated PR transcript transcribe the O and P genes essential for phage genome replication and the Q gene that controls transcription of later genes on the lytic pathway. Leftward extension of the antiterminated PL transcript transcribes xis and int genes essential for phage chromosome integration and excision into and out of the host chromosome. Locations of four termination sites are indicated by TR1–2 and TL1–2.

The core of the bistable switch is the complex biochemistry of the P_R and P_RM promoters' operator regions, which share three overlapping operator sites (Figure 1A, O_R1, O_R2, O_R3), where Cro and CI dimers bind competitively and in sequence, but in opposite order (MAURER et al. 1980 ; MEYER et al. 1980 ; MEYER and PTASHNE 1980 ; PTASHNE 1992 ). Figure 2, a and b, shows contour maps of the activation level of P_R and P_RM, respectively, as a function of CI and Cro dimer concentration. The activation levels are calculated using the model and parameters in SHEA and ACKERS 1985 .

View larger version (83K): In this window In a new window Download PPT slide   Figure 2. Activation contours of the PR and PRM promoters versus Cro2 and CI2 concentrations. Contour intervals are 5% of maximum activation of the respective promoters. Activation levels are calculated using the method and parameters in SHEA and ACKERS 1985 . (a) PR promoter. Peak activation, immediately after infection, is ~0.013 open complexes/sec (OC/s) or an average of one OC every 1.3 min. With no CI2, 50% repression by ~170 molecules of Cro2; with no Cro2, 50% repression by ~36 molecules of CI2. Number of molecules corresponding to the molarity on the y-axis is indicated. (b) PRM promoter. Immediately after infection, PRM has a low basal level of activity. The maximum activation in the absence of Cro2 is ~0.0063 OC/s or an average of one OC about every 2.6 min. (One molecule in an E. coli cell is approximately a 1 nM concentration.) The positive slope of PRM activity versus CI2 concentration around 100 nM enables autoregulatory maintenance of CI concentration in the stable lysogenic cell. (c) "Trajectory" of one cell in concentration space from infection to initiation of lysogeny. Points plotted are at 1-sec intervals. The erratic path is a consequence of the randomness of the protein production and degradation reactions. (d) Postinfection trajectories of the respective lysogenic-fated and lysis-fated cell subpopulations. Bold lines are the mean trajectories. The downward slope of the trajectory of the lysogenic subpopulation arises principally because Cro2 concentration is decreased by dilution in the growing cell after repression of PR.

The lysis or lysogeny outcome in each cell is determined by the specific temporal pattern of CI and Cro accumulation in that cell after infection. Immediately after infection, there are no CI or Cro molecules in the cell so the regulatory circuit is in the state labeled "S" in the lower left corners of Figure 2, a and b. At that point, P_R is fully activated; P_RM has only a low basal activation, and promoter P_L is also activated. Transcription and translation of the phage DNA is accomplished by the host cell's machinery.

Cro and N proteins are produced from transcripts initiated at promoters P_R and P_L and both proteins begin to accumulate immediately after infection. Initially terminators T_R1 and T_L1 partially block RNA polymerase (RNAP) transcription: about 50% at T_R1 (FRIEDMAN and GOTTESMAN 1983 ) and 80% at T_L1 (DRAHOS and SZYBALSKI 1981 ). However, as the concentration of N increases, the N protein (with other molecules from the host cell) acts to antiterminate RNAP at NUT sites upstream from T_R1 and T_L1 (DAS 1992 ), and the average rate of transcription of downstream genes, including cII and cIII, increases.

The CI production rate is determined by the combined transcript initiation rates from P_RE and P_RM (Figure 1A). Since P_RM is essentially OFF until there is some CI in the cell (Figure 2B), there will be no accumulation of CI (and thus potential for lysogeny) unless P_RE-initiated transcripts produce enough CI to get the cell into the concentration state where P_RM is activated and P_R is repressed, that is, into the regime where Log[CI₂] > ~-6.8 and Log[Cro₂] < ~-7.2 (or roughly more than 145 CI dimers and less than 55 Cro dimers, respectively). This can occur only in those cells where, by chance, there is early, strong CII production that persists long enough to activate P_RE. If, however, more than about 55 Cro dimers accumulate first, lysogeny will be precluded. (A nominal cell volume of 1.4 x 10^-15 liters is used here to relate molecular concentration to molecule count.) In cells where, by chance, enough early CI production from P_RE transcripts occurs to repress P_R and activate P_RM, then CI concentration will tend to continue to increase automatically due to positive autoregulation leading to ever increasing repression of P_R and P_L by CI₂. Eventually CI₂ concentration will rise into the negative feedback region of the P_RM repression curve and stabilize by autoregulation at a concentration in a range of 140–200 dimers per cell (REICHARDT and KAISER 1971 ; LEVINE et al. 1979 ), and Cro will eventually disappear due to dilution and degradation.

Though CII is produced from the same P_R-initiated transcripts that encode Cro, CII accumulation is initially attenuated by termination of about 50% of the transcripts at T_R1 and by CII's relatively short half-life (~2 min). Thus, initially, only Cro accumulates in the infected cells. In the presence of CIII, degradation of CII is reduced (HOYT et al. 1982 ). In the lysis-lysogeny decision, the Hfl system integrates two environmentally dependent signals into the circuit function: (i) nutritional state of the cell, and (ii) the level of the phage population in the cell's surrounding vicinity. This article examines how the latter sensing mechanism works. Both sensing functions depend on active control of CII proteolysis (Figure 1A). At higher levels of nutrition, Hfl-related proteolytic activity is higher so that CII and CIII have shorter lifetimes (GRODZIKER et al. 1972 ; BELFORT and WULFF 1974 ). This tends to reduce the mean and peak CII concentration levels, and thus the probability of P_RE activation. So the probability of lysogeny is reduced in well-fed cells. Dependence of lysogeny on MOI arises because the concentration of the host-encoded Hfl system is independent of MOI, while the number of cII and cIII genes is proportional to MOI. Thus, cells with higher MOI have a higher probability of achieving CII concentration necessary to activate P_RE and kickstart CI production.

If a cell reaches a state where (i) the Cro feedback loop is established, (ii) P_RM and P_L are repressed, and (iii) CII concentration is low, there is a high probability that the cell will continue on the lytic path. On the other hand, if a cell reaches a state where (i) the CI feedback loop is established, and (ii) P_R and P_L are repressed, there is a high probability that the cell will continue to lysogeny (MCADAMS and SHAPIRO 1995 ).

	SOURCE DATA AND GENETIC CIRCUIT MODEL

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

Kourilsky's measurements of lysogeny versus API:
We use the experimental assays of percent lysogeny versus API in KOURILSKY 1973 to compare with predictions of the stochastic kinetic model. In Kourilsky's experiments, lambda phages were added at various API to exponentially growing E. coli cultures for an incubation time that ensured near 100% phage absorption. Kourilsky's measurements included O^- and P^- strains incapable of phage chromosome replication. Since the phage chromosome count does not increase, the postinfection distribution of the O^- or P^- phage particles among the target cells can be computed using the Poisson infection statistics model described below. Selection of the O^- and P^- mutants for modeling eliminates the need to model chromosome replication.

In Kourilsky's plots of log API versus the log of the percent cells lysogenized, the shape of the rate of lysogenization versus API curves was similar for starved and unstarved cells, but the starved curves were systematically shifted to a 50–100 times higher lysogenization rate with little effect on the qualitative dependence of lysogenization rate on the infection ratio [ Figure 2 in KOURILSKY 1973 ]. The starved cell results with 50–100 times higher rates of lysogeny are used for comparison since the number of simulation runs necessary to estimate the fraction, f, of lysogeny varies as 1/f.

Stochastic kinetic model:
The stochastic kinetic model used here to analyze operation of the lysis-lysogeny decision circuit includes the genetic mechanisms and the coupled protein dimerization and degradation reactions shown in Figure 1A. Genetic mechanisms are modeled using explicit, though approximate, reaction models of each submechanism and explicitly including features such as termination sites. Thus, promoter operator sites are modeled using the statistical-thermodynamic approach described by SHEA and ACKERS 1985 . (The stochastic version of the Shea and Ackers model of promoter kinetics is produced by calculating the instantaneous probability of each distinct transcriptionally active state of a promoter using the partition function, and then using this probability in calculating the reaction probabilities for the transcript initiation reactions.) Transcript elongation is modeled as a sequence of individual nucleotide steps. Translation control is modeled as described by MCADAMS and ARKIN 1997 . Assumptions and reaction models used for elements of the system are described below and listed in Table 1 and Table 2. The reaction models are represented as a set of coupled stochastic kinetic equations.

Analytical solution to such systems of stochastic reaction equations is only practical for simple reaction systems. However, numerical solutions can be computed for complex systems of coupled stochastic reactions using the Monte Carlo algorithm described by GILLESPIE 1977 . The Gillespie algorithm produces a stochastic realization of the temporal behavior of the system by calculating the probabilistic outcome of each discrete chemical event and the resulting changes in the number of each molecular species. In the application of the Gillespie algorithm to simulation of bacterial regulation, each simulation run provides a representative case of the sequence and timing of events and the regulatory outcome in an individual cell starting from specified conditions. Multiple runs with the same initial conditions (e.g., the same MOI) are used to estimate the probability that cells will enter lysogeny for these conditions. [About samples are required to estimate the probability, P, of a binary random event with 95% confidence where f_e is the desired maximum fractional error in P (FELLER 1968 ). In the present case, P is the probability of lysogeny in each cell.] Computations were performed on SGI workstations and parallel array supercomputers (200 node Cray T3D machine at Eglin AFB, and 400 node SP2 machine at Maui High Performance Computer Center). Additional details on the software are available on the website: www.lbl.gov/~aparkin/LambdaMod.html/. Criteria described below were used to categorize the outcome of each individual run as a lysogenic outcome or not.

To compare the percent lysogenization predicted by the simulation at various infection levels to experimental observations obtained by KOURILSKY 1973 , the Poisson-weighted average of the probability of lysogeny at each MOI is computed to estimate the expected percent lysogeny as a function of API. The theoretical Poisson probability that a given cell will be infected with exactly MOI = M phage when API = A is

(1)

,

where P(M,A) is the probability of a cell having MOI = M, at API = A (ELLIS and DELBRUCK 1939 ). The expected fraction of lysogens at a given API, F_lysogens, is then

(2)

,

where F(M) is the estimated probability of lysogeny for cells with various MOIs as estimated using the stochastic kinetic model.

Modeling Assumptions

1. Dispersion in cell generation times can be neglected. Thus all runs used a cell cycle time of 35 min, consistent with the cell cycle time reported by KOURILSKY 1973 . E. coli cell generation times are observed to be approximately normal distributed with standard deviation of about 22% of the mean (PLANK and HARVEY 1979 ). The neglect of dispersion in cell generation times is equivalent to assuming that any growth rate-related effects in faster growing cells roughly offset the effects in slower growing cells.

2. The volume of the cell grows approximately linearly from 1 x 10^-15 to 2 x 10^-15 liters. (The maximum difference in volume between linear and exponential cell growth models is <6% with negligible effect on simulation results.)

3. Host housekeeping molecules relevant to phage gene expression and phage protein degradation are constitutively expressed and regulated at constant concentration, which is the same in all cells. This implies, for example, that all enzymes required for metabolic pathways, etc., are expressed at levels consistent with a healthy bacterium and that cytoplasmic concentrations of proteases, RNAP, ribosomes, and metabolic substrates are maintained during the early postinfection period when the lysis-lysogeny decision is being resolved. Both RNA polymerase and ribosomes are present in the cell in relatively large numbers, however, the free polymerase and ribosome concentrations are thought to be a fraction of the total and to be buffered by exchange with units that are engaged in other reactions. Under these conditions, fluctuations in polymerase and ribosome concentrations would be relatively small.

4. Regulatory effects of host proteins such as integration host factor and RNase III on phage gene expression are assumed to be equivalent for all cells and constant over time. For example, the effect of integration host factor on P_L activity (GILADI et al. 1990 ) is assumed to be included in the kinetic parameters of the promoter and to be independent of phage MOI.

5. Effects such as macromolecular crowding or two-step binding to DNA or RNA that might affect reaction kinetics are assumed to be subsumed into the kinetic parameter characterizing the reaction.

6. Intermediate reactions (as with sigma-factors or other subelements) in assembly of the RNAP and ribosome complexes are not rate limiting. Instead, we assume either (i) that the host cell maintains an effective concentration of transcriptionally and translationally available concentrations of these molecular complexes that are in rapid exchange with their binding sites on the DNA or RNA, or (ii) that the component subunits are in rapid equilibrium with functionally active assemblies (SHEA and ACKERS 1985 ; PTASHNE and GANN 1997 ). The rate-limiting step in transcript initiation is assumed to be the closed- to open-complex isomerization reaction (MCCLURE 1980 ).

7. Phage gene expression is stochastic, consistent with the mechanisms described by MCADAMS and ARKIN 1997 .

8. An average of 10 proteins are produced per transcript for all genes (KEPES 1963 ; SHEA and ACKERS 1985 ). Transcript degradation rates and ribosome binding rates are chosen to produce that average yield.

9. In Kourilsky's experiment the E. coli cells were unsynchronized, hence they were presumably infected at random times in the cell cycle (KOURILSKY 1973 ). We assume all infections occur early enough in the cell cycle so that cell growth only affects operation of the decision logic by dilution effects on concentrations of phage-encoded molecules. The initial rates of phage protein production from P_R- and P_L-initiated transcripts in each host cell are independent of cell volume. However, for the same rate of protein production, the consequent rate of change in phage protein concentration is cell size dependent so that timing of subsequent events could be slowed somewhat for larger cell size at infection time. Most cells that are fated to become lysogens are committed by 10 to 15 min (causing, for example, the cessation of Cro₂ production shown in Figure 3C). Most infections early in the cell cycle are thus resolved before the next cell division. For infections occurring late in the cycle, if commitment has not occurred before division, the phage chromosomes and proteins at division are randomly shared between the daughter cells when the cell divides. Then the phage infection continues, but with lower MOI. For the O^- or P^- phage mutants, the average postdivision MOI is halved, since phage chromosome replication is not possible. Halving the MOI reduces the probability of lysogeny in the daughter cells. This suggests that neglecting cell division leads to some degree of overestimation of the probability of lysogens in the simulation.

   
   

   View larger version (58K): In this window In a new window Download PPT slide   Figure 3. The solid lines in (a) show the time course of the average intracellular Cro2 and CI2 concentrations at MOI 6. The shaded region indicates the ±1 range as estimated by determining the 16th and 84th percentile points in the population at each time. (b, c) show the same data, but for the two subpopulations with different phenotypic fates. The concentration profiles of the two regulatory dimers in each subpopulation are similar for the first few minutes, but diverge into a substantially different time pattern after about 7 min. Common experimental methods for assaying time evolution of protein concentrations would yield data equivalent to the average value curves in (a), masking the differences in the diverging subpopulations.

10. The target cells are infected effectively simultaneously so that no temporal infection effects or phage infection-dependent immunity occurs.

11. The cell is assumed to be a homogeneous, well-stirred medium so the concept of "protein concentration" is valid and spatial effects are averaged out. [E. coli signaling proteins have been shown to diffuse distances comparable to the cell dimensions in much less than a second (ISHIHARA et al. 1983 ).]

12. A cell becomes committed to lysogeny if there is (i) a sufficient time-integrated concentration of CII to activate P_RE, and (ii) [CI₂] > [Cro₂] at the end of the 35-min cell cycle. Activation of P_RE was defined as an average activation level of one open-complex per 2 min over a contiguous 4-min period. This level of CII production would also activate the other CII-dependent promoters, P_anti-Q and P_I, that function in execution of the lysogenic pathway (MCADAMS and SHAPIRO 1995 ). CI₂ concentration greater than that of Cro₂ at the end of the cell cycle is an additional indication that activation of P_RE occurred early enough and was productive enough to lock on the CI feedback loop.

Reaction Models
The following paragraphs describe the models used for the mechanisms of the lysis-lysogeny decision circuit. Reactions and parameters are listed in Table 1 Table 2 Table 3. Parameters in the kinetic model are derived from the sources cited in Table 1 Table 2 Table 3. Considerations underlying selection of the CII and CIII proteolysis reaction models and rate parameters are described below.

Phage gene expression: The genetic mechanisms associated with transcript initiation and translation control produce the largest component of the stochastic effects that lead to divergent phenotypes in the infection system. The following genetic reactions are modeled: operator/promoter binding, transcript initiation, transcription, initiation of translation, translation, and initiation of mRNA degradation. The transcription model includes mechanisms for the two RNAP termination sites, T_R1 and T_L1, and antitermination at the NUT_R and NUT_L sites (Figure 1A).

Operator/promoter binding and control of transcript initiation: SHEA and ACKERS 1985 describe a statistical mechanical/thermodynamic approach to modeling the P_R/P_RM promoter complex. We use the same approach to model P_L repression by CI₂ and Cro₂ and P_RE activation by CII using the parameters in Table 1. The instantaneous probability of each distinct occupancy state of a promoter is assumed to be determined by the partition function defined in accord with the Shea/Ackers formulation, and we use the probability in the stochastic formulation of kinetics as defined by GILLESPIE 1977 , GILLESPIE 1992B . A key assumption is that effector molecule binding reactions at a promoter occur much faster than the rate of transcript initiation at the promoter. With this rapid equilibrium assumption, the binding state of the promoter is modeled by randomly choosing the promoter state at each instant using the probabilities given by the partition function (SHEA and ACKERS 1985 ). If the promoter state selected is one from which RNAP can initiate transcription, then that transcript initiation reaction is included in the list of possible reactions for the next Monte Carlo calculation in accord with the Gillespie algorithm. Whenever the Monte Carlo calculation determines that a transcript is initiated from one of the promoters, a new transcribing RNAP is "initiated" on the corresponding DNA transcription object. The promoter activation functions in Figure 2 show the resulting average rates of transcript initiation as a function of CI₂ and Cro₂ concentrations. Occlusion of the promoter site by the footprint of a recently launched RNAP is included in the model.

Transcription: A transcript elongation model estimates the time delays between transcript initiation and arrival at the end of each coding region on the operon. This delay, plus the delay until an effective level of signaling molecules is accumulated, determines the timing of regulatory molecule concentrations that control regulatory networks. The movement of transcribing RNAP along the DNA is modeled as a sequence of independent one-nucleotide reaction steps. Each such reaction is assumed to be unidirectional, that is, RNAP movement is assumed to be strongly forward-biased. It is assumed that there is a single rate-determining reaction for each RNAP step and that each forward step has constant probability of occurring per unit time, leading to an exponential distribution of interstep times. The exponential is characterized with an average step-time parameter. The same average step time was used at each nucleotide position and for coding and noncoding regions, neglecting the differences between transcription rates for different nucleotides. (Transcription through termination and antitermination sites is described below.)

Termination: At termination sites, transcribing RNAPs slow down (i.e., the average interstep time parameter is larger) and there is a probability of transcript termination at the site. When the RNAP is antiterminated upstream of the terminator site, the termination site is then modeled as the "normal" DNA described above.

Antitermination: The reactions to assemble the antiterminated form of RNAP at NUT_L and NUT_R sites (Figure 1A) depend on N protein concentration. The antitermination reaction complex also involves Rho and at least four additional host factors: NusA, NusB, NusG, and S10 (WHALEN et al. 1988 ; MASON and GREENBLATT 1991 ; LI et al. 1992 ; DEVITO and DAS 1994 ). These factors are assumed to be constitutively expressed and present in the necessary concentrations so that the concentration of N is limiting. In vitro studies suggest that in the presence of the proper host factors, the fractional readthrough of a downstream terminator increases in proportion to the concentration of N until there is full antitermination at N concentrations between 50 and 100 nM (WHALEN et al. 1988 ; LI et al. 1992 ; DEVITO and DAS 1994 ). Under these conditions the antitermination process at the NUT site is modeled as a single-step reaction assumed to be a pseudo-first-order reaction with the rate chosen such that antitermination of RNAP is near 100% for N concentrations above about 75 nM. This parameter choice fits the experimental dependence of fractional readthrough of a downstream terminator as a function of N concentration (WHALEN et al. 1988 ; LI et al. 1992 ; DEVITO and DAS 1994 ).

Translation: Translation control is modeled as described by MCADAMS and ARKIN 1997 , based on the mechanism described by YARCHUK et al. 1992 . In that model, a competition between ribosome and RNase E binding determines the average number of proteins produced per transcript. In the stochastic kinetic model of the circuit, this ribosome-RNase E competitive binding reaction is treated as a stochastic chemical reaction. The temporary occlusion of the ribosome binding site after a successful ribosome binding event is modeled. Motion of a translating ribosome on a transcript is modeled similarly to the model of motion of a transcribing RNAP on DNA described above. If one ribosome by chance overtakes another in the model, the progression of the former is halted until the latter moves ahead. The average ribosome step time is selected to be shorter than the RNAP step-time parameter, producing ribosome queuing as is observed (KENNELL and RIEZMAN 1977 ; YARCHUK et al. 1992 ). As with the transcription model, the statistics of the interstep times are assumed to be described by the exponential probability function.

Phage protein dimerization and degradation: The principal reactions involving phage-encoded proteins in the decision circuit are identified in the boxes labeled R1 to R5 in Figure 1A. Reactions in R1 include degradation and dimerization of CI; R2 includes dimerization and degradation of Cro; R3 and R4 include competitive degradation of CII and CIII by the two host cell proteases (see below); and R5 includes degradation of N.

CI and Cro dimerization and degradation: Degradation of CI and Cro is modeled as occurring predominantly by proteolysis of the monomeric form, a common degradation mode for multimeric proteins (SHEA and ACKERS 1985 ; GOTTESMAN and MAURIZI 1992 ). The dimerization reactions are assumed to be characterized by fast forward and reverse reaction rates whose ratio, the dissociation constant, is 20 nM (SHEA and ACKERS 1985 ). The dynamics of decay of multimeric protein populations is not exponential in general. Rather, the dynamics and specifically, the measured lifetimes, depend strongly on the forward and reverse rates of the multimerization reactions and on the initial concentration of each subspecies. Reaction parameters for the stochastic kinetic model here were selected to match experimental measurements of lifetimes.

Degradation of CII and CIII: Two membrane-bound protein complexes, HflA (CHENG et al. 1988 ) and HflB (BANUETT et al. 1986 ; HERMAN et al. 1995 ; KIHARA et al. 1997 ), act to degrade CII. HflB is thought to degrade CIII (HERMAN et al. 1995 ). Initial studies identified HflA as a nonessential membrane-bound GTP-utilizing protease (CHENG et al. 1988 ; ZORICK and ECHOLS 1991 ; NOBLE et al. 1993 ). However, later experiments report that HflA is probably a regulator of HflB activity and that HflB is the major protease responsible for CII degradation. Additional evidence for proteolytic activity of HflB is provided by in vitro experiments with purified protein (HERMAN et al. 1995 ; KIHARA et al. 1997 ; SHOTLAND et al. 1997 ). The Hfl proteolytic system has other host protein targets in addition to CII (CHENG and ECHOLS 1987 ; HERMAN et al. 1993 ; HERMAN et al. 1995 ). Both HflA and HflB respond to host environmental signals. There is some evidence that HflA activity is directly or indirectly affected by the catabolite-activating protein/cAMP system, which has been shown to reduce proteolytic activity in response to carbon-source starvation (HOYT et al. 1982 ; BANUETT and HERSKOWITZ 1987 ). CIII protects CII from proteolysis (HOYT et al. 1982 ; RATTRAY et al. 1984 ; BANUETT et al. 1986 ) even in the absence of HflA and HflB activity, which implies the existence of yet another proteolytic pathway for CII degradation (KIHARA et al. 1997 ). In summary, although phenomenology of the proteolysis of CII and CIII is relatively well characterized, the exact mechanisms whereby HflA, HflB, and perhaps another unidentified protein control degradation of CII and protection of CII by CIII are not known.

The half-life of unprotected CII has been observed to be anywhere from 5 min (HOYT et al. 1982 ) to less than 30 sec depending on conditions (RATTRAY et al. 1984 ). The short half-life of CII and the relatively low concentrations of CII protein and HflA/B suggest that binding of CII to the Hfl proteins is tight and fast. Two alternative mechanisms have been hypothesized for the protection of CII by CIII (CHENG et al. 1988 ): (i) competitive binding of CII and CIII to the Hfl (and perhaps another) proteolytic complex, or (ii) direct CIII binding to CII to form a proteolysis-resistant complex. Available experimental data does not differentiate between the two alternatives. In order to select among candidate models, we investigated alternative reaction mechanisms seeking a reaction system that meets four constraints: (a) yields a 2-min half-life for CII in the range of initial concentrations spanning 50–100 nM (GOTTESMAN and GOTTESMAN 1981 ; CHENG et al. 1988 ), (b) produces an approximately 6-min half-life for CIII in the absence of CII (KORNITZER et al. 1991A ), (c) yields CIII protection of CII consistent with HOYT et al. 1982 and RATTRAY et al. 1984 , and (d) functions in the overall simulation model to produce the simulated percent lysogeny as a function of API consistent with KOURILSKY 1973 . Constraint (d) proved the most restrictive. The best satisfaction of the constraints (a) through (d) was obtained using a proteolytic system in which CII and CIII are competitive substrates for two independent proteases. The resulting reaction mechanism is shown in Figure 1A and rate parameters used are in Table 3. One protease (called P1) is more specific than the other but saturates at very low concentrations of CII; the other (called P2) is only slightly less specific, does not saturate, and has a much higher maximal activity. These proteases correspond to HflB and the putative second protease identified by KIHARA et al. 1997 .

Degradation of N: The half-life of the antitermination-controlling protein N is approximately 5 min. Degradation is by the Lon protease, and Lon is also thought to be responsible for degrading the Hfl proteins (GOTTESMAN and GOTTESMAN 1981 ). Reported data on N degradation fit a first-order decay curve (GOTTESMAN and GOTTESMAN 1981 ); thus N probably does not saturate the protease. Accordingly, N degradation is modeled as a first-order reaction.

Cell growth: The linear cell growth assumption was implemented as a constant probability of adding a small fixed volume increment each instant of time. Each run was started at an initial cell volume of 1 x 10^-15 liter and continued until the volume doubled to 2 x 10^-15 liter over 35 min of simulated cell time.

	RESULTS

TOP ABSTRACT APPROACH STOCHASTIC KINETICS PHAGE LAMBDA DECISION CIRCUIT SOURCE DATA AND GENETIC... RESULTS DISCUSSION LITERATURE CITED

Time course of pathway selection:
Figure 4 and Figure 5 show the temporal trajectory of the concentration of key protein molecules in one lytic and one lysogenic case selected from runs at MOI 6. (Figure 2C is based on the same two cases.) The two cases show the randomness in the intracellular regulatory protein concentration trajectories and the differences in the trajectories for the divergent developmental paths possible in two initially identical cells. Of the phage-encoded proteins shown in Figure 4 and Figure 5, Cro₂ and CII are expressed earliest in both the lytic and lysogenic cases. Cro₂ appeared within 1 min of infection (Figure 5B) and CII appeared within 2 min (Figure 4A). Protein expression in the two cases began to diverge after about 5 min. Both the lytic- and lysogeny-fated cases experienced a nearly equal burst of CII production at this time (Figure 4A), however, in the lysogeny-fated case, there was a simultaneous burst of CIII production (Figure 4B). So lysogeny resulted in this case because, by chance, the bursts of CII and CIII were both large and simultaneous so that CII degradation was slowed and it survived long enough to activate P_RE and kickstart CI production. Figure 4C shows that the CII/CIII proteases were strongly inhibited by the bursts of CIII production in the lysogenic case. CI₂ concentration (Figure 5A) in the lysogenic case began to grow at about 12 min just after CII concentration peaked. The growing CI₂ concentration repressed P_R and stopped Cro production. As a result, the Cro₂ concentration declined in the lysogeny-fated case after 12 min (Figure 5B). In contrast, in the lytic-fated case no CIII production occurred so the unprotected CII rapidly degraded and did not activate P_RE enough to start the CI expression feedback loop. Without expression of CI, Cro₂ production continued (Figure 5B) and lysis ensued.

View larger version (32K): In this window In a new window Download PPT slide   Figure 4. Time evolution of CII and CIII concentration for two cases at MOI 6 illustrating a lytic and a lysogenic (bold) outcome. The pattern of protein concentration growth is distinctive and different for every simulation run. Lysogeny requires early, higher CII concentration as in (a) so that promoter PRE is activated, and protein CI and its dimer begin to accumulate to turn on PRM and repress Cro production. The high CIII concentration in the lysogenic case in (b) protected CII from degradation. The percent proteolytic activity in (c) calculated by ((k10 · ([P1] + [P1 · CII])/(total P1)) + (k16 · ([P2] + [P2 · CII])/(total P2)))/(k10 + k16) indicates the percentage of the total protein activity available for the degradation of CII. (Figure 5 shows additional protein concentrations from the same simulation runs.)

View larger version (32K): In this window In a new window Download PPT slide   Figure 5. Time evolution of Cro and CI dimer concentrations for the same two simulation runs at MOI 6 as Figure 4. For the lysogenic case (bold), the high CII concentration after 6 min (Figure 4A) leads to the accumulation of CI2 (a) and cessation of Cro production (b). Dilution and degradation causes Cro2 concentration to decline thereafter. For the lytic case, in contrast, the initial burst of CII is not sustained (Figure 4A) so that PRE is not significantly activated and CI production is negligible (a). Cro2 growth begins immediately after infection (b) and, in lytic cases, continues building until it represses both PL and PRM thus ending the possibility of lysogeny.

Figure 2C shows the intracellular CI and Cro dimer concentration trajectory for the lysogenic-fated case at MOI 6 (identical data as Figure 5) superimposed on the P_R and P_RM promoter activation contours. The CI₂ repressor concentration began to autoregulate its own concentration 20 min after infection; thereafter, the CI₂ concentration remained constant and P_RM activation slowly increased as Cro₂ concentration was diluted in the growing cell (Figure 2C, arrow). The concentration trajectory for the lytic case in Figure 5 is not shown in Figure 2C to avoid confusing the figure. However, the oval on Figure 2C indicates the region where the concentrations stabilized at 12 min after infection.

Estimated statistics of concentration trajectories:
The Monte Carlo solution to the stochastic kinetic equations produces a database of representative time-dependent samples of the concentration trajectories as the infection progresses for each molecular species in the reaction system. Analysis of this database provides estimates of the statistical parameters of the infection progress in cells with corresponding initial conditions, e.g., a particular MOI. Figure 3A shows the estimated statistical distribution of the CI₂ and Cro₂ concentration trajectories for the subset of cells at MOI 6. (For all plots in Figure 3, the bold lines are the average concentration of the indicated species and the lighter lines are the ±1 range.) The lysis- and lysogeny-fated subsets shown in Figure 3B and Figure C, each experience a different pattern of Cro₂ and CI₂ concentration growth statistics, distinct from each other and from the combined statistics. Figure 2D shows the same average Cro₂ and CI₂ concentration trajectories for the lysogenic-fated and lytic-fated cases at MOI 6 superimposed on the P_R and P_RM promoter activation contours.

Lysogenic fraction: Kinetic model estimates compared to experiment:
The experimental lysogeny fraction data shown in Figure 6B for starved O^- () and P^- () mutants are from Figure 2 of KOURILSKY 1973 . For the higher API values in Figure 6B, Poisson statistics of infection (Equation 1) predicts that some cells would have quite high infection levels. At some infection level, the phage processes must become disruptive to the host cell processes so that the assumptions underlying the kinetic model become invalid. To address this possibility we solved for two sets of CII/CIII proteolysis parameters: first, for the best match to the experimental O^- and P^- data for all available API (Figure 6A, labeled "Full", symbol: ), and, second, for the best match considering only API values 6 (Figure 6A, labeled "Lower", symbol: ). Points in Figure 6A ( and ) reflect the estimated probability of lysogeny in individual infected cells vs. MOI from solution of the stochastic kinetic model equations for these two different choices of Hfl parameters. The vertically hatched area in Figure 6B indicates the range of difference between the resulting estimates of the fraction of lysogens. The hatched area in Figure 6A indicates the corresponding range of differences in the probability of lysogeny vs. MOI resulting from the different proteolytic model parameters. Both curves in Figure 6A show negligible lysogenization at MOI < 3 and a rapid increase in lysogeny for MOI > 3. Corresponding points in Figure 6B yield the solid lines bounding the hatched region. These estimates of the fraction of lysogens in an infected cell population versus API are calculated as the Poisson-weighted sum of points for different MOIs in Figure 6A for corresponding cases using Equation 2. The match with experiment is good at the critical low MOI values, but falls above observed values at high MOI. We attribute the overestimation of the percentage lysogeny at high API in Figure 6B predominantly to disruption of host cell processes at high infection levels.

View larger version (24K): In this window In a new window Download PPT slide   Figure 6. Comparison of results from stochastic kinetic model and experimental results from KOURILSKY 1973 for fraction of lysogens produced vs. MOI. (a) Simulation results for fraction of cells producing lysogens (, ). Curve labeled "Full" results from choice of proteolysis parameters to match the full experimental data set in (b); curve labeled "Lower" results from proteolysis parameters chosen for best match to experimental points at lower API values in (b). Vertically hatched area in (b) indicates the range of difference between the resulting estimates of the fraction of lysogens. Results with several "digital mutants" are shown in (a): O-T-, termination sites removed (x); O-N-, hence no-antitermination (); O-Coop-, noncooperative binding of CI dimers at OR1–3 (). (b) Solid lines bounding the hatched region are the predicted fraction of lysogens for the Full and Lower cases in (a) calculated by weighting the results shown in (a) by the theoretical Poisson statistical distribution of the number of phage per cell at each API. Experimental points for the fraction of lysogens for O- () and P- () strains. Experimental points are from Figure 2 in KOURILSKY 1973 for cells starved before infection. The three curves labeled Poisson-n in (b) show the hypothetical fraction of lysogens vs. API expected for a "threshold model" where all cells with MOI n are assumed to become lysogens. (c) Curves labeled O-N-/50 and O-T-/50 are predictions for "digital mutants." See text for explanation.

The three curves labeled Poisson-n show the hypothetical fraction of lysogens vs. API expected for a "threshold model," where all cells with MOI n are assumed to become lysogens. The solution of the stochastic kinetic model exhibits rapid onset of lysogeny for MOI > 2 (Figure 6A), representing an approximation to a threshold process in the decision circuit produced by the reinforcing effects of production from multiple promoters and earlier antitermination as MOI increases. At low APIs, both the experimental points and the stochastic model predictions for the O^- and P^- mutants lie between the idealized threshold model predictions for thresholds at MOI 3 and MOI 4.

Digital mutants:
Additional tests of the kinetic model by predictions of other experimental observations are needed; however, we are unaware of additional, independent measurements for similar strains and conditions. Accordingly, we include in Figure 6C testable predictions of rates of lysogeny for several "digital mutants" based on changes in the stochastic kinetic model reflecting several mutant cases. The curve labeled O^-N^- in Figure 6C reflects our prediction of the percent lysogeny for a digital mutant with the function of the N protein disabled in the kinetic model and all other parameters as for the curve labeled "Full". (We use the "O^-" notation to indicate replication deficient, i.e., either O^- or P^-, mutants.) This is the prediction of percent lysogeny from the kinetic model for a starved O^-N^- mutant; the curve labeled N^-/50 is the corresponding prediction for an unstarved O^-N^- mutant. [The "unstarved" estimate is derived by dividing the "starved" estimate by 50, consistent with the observation by KOURILSKY 1973 that starved cells systematically exhibited increased lysogeny by 50–100x.] The predicted level of lysogeny is reduced for the O^-N^- case because the transcribing polymerases are never antiterminated and production of CII and CIII is always lower than when N production leading to RNAP antitermination is possible.

The curve in Figure 6A labeled "O^-Coop^-" shows the predicted probability of lysogeny for a digital mutant where CI binding to the P_R–P_RM operator sites is made noncooperative in the kinetic model, reducing the effectiveness of positive autoregulation of P_RM. The predicted experimental fraction of lysogens (not shown) is close to the curves for the O^-N^- cases in Figure 6C. The dashed line labeled O^-T^- in Figure 6C is the estimate for another starved digital mutant with the T_R and T_L termination sites disabled (the line labeled O^-T^-/50 is for unstarved mutants). The reduced slope of lines for the O^-T^- mutant in Figure 6C is due to the predicted increase in the estimated probability of lysogeny for cells with MOI 1 and 2 for this mutant as shown in Figure 6A.