Genetics, Vol. 149, 1633-1648, August 1998, Copyright © 1998

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

Adam Arkin1,a, John Rossb, and Harley H. McAdamsa
a Department of Developmental Biology, Stanford University, Stanford, California 94305
b Department of Chemistry, Stanford University, Stanford, California 94305

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 {lambda} 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 Down, 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 Down, GILLESPIE 1992B Down). 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 Down). ("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 Down; ROBERTSON 1992 Down; FINLAY and FALKOW 1997 Down; STRAUSS and FALKOW 1997 Down). 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 {lambda}-infected Escherichia coli cells. This path selection is made by the {lambda} 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 {lambda}-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 Down 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 Down. 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 {lambda} lysis-lysogeny decision circuit (details below) uses the stochastic formulation of chemical kinetics (GILLESPIE 1976 Down), includes stochastic mechanisms of gene expression (MCADAMS and ARKIN 1997 Down), and models promoter regulation using the statistical-thermodynamic approach described by SHEA and ACKERS 1985 Down. The set of coupled stochastic equations comprising the stochastic kinetic model is solved using a standard Monte Carlo algorithm (GILLESPIE 1976 Down) 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 Down). 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 Down; ZHENG and ROSS 1991 Down). 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 Down) 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 Down). [Note that the erratic time pattern, which we postulate for protein production from the operons comprising the {lambda} 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 Down).]

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 Down, GILLESPIE 1992B Down; VAN KAMPEN 1992 Down) that is described by the chemical master equation (GILLESPIE 1992B Down; VAN KAMPEN 1992 Down). The solution to the chemical master equation can be calculated using a stochastic simulation algorithm devised by GILLESPIE 1976 Down, GILLESPIE 1977 Down, GILLESPIE 1992B Down. 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 {lambda} phage lysis or lysogeny decision are generally known (HERSKOWITZ and HAGEN 1980 Down; FRIEDMAN and GOTTESMAN 1983 Down; ECHOLS 1986 Down; FRIEDMAN 1992 Down; PTASHNE 1992 Down; MCADAMS and SHAPIRO 1995 Down). 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 Down, MEYER and PTASHNE 1980 Down, PTASHNE 1992 Down, and SHEA and ACKERS 1985 Down. 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 Down; RATTRAY et al. 1984 Down). 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 Down. 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 Down; MCADAMS and SHAPIRO 1995 Down). The circuit shown in Figure 1A has two key subsystems: (i) the PRM-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 {lambda} 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 Down; SHEA and ACKERS 1985 Down). 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) {lambda} 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 PR and PRM promoters' operator regions, which share three overlapping operator sites (Figure 1A, OR1, OR2, OR3), where Cro and CI dimers bind competitively and in sequence, but in opposite order (MAURER et al. 1980 Down; MEYER et al. 1980 Down; MEYER and PTASHNE 1980 Down; PTASHNE 1992 Down). Figure 2, a and b, shows contour maps of the activation level of PR and PRM, 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 Down.



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 Down. (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, PR is fully activated; PRM has only a low basal activation, and promoter PL 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 PR and PL and both proteins begin to accumulate immediately after infection. Initially terminators TR1 and TL1 partially block RNA polymerase (RNAP) transcription: about 50% at TR1 (FRIEDMAN and GOTTESMAN 1983 Down) and 80% at TL1 (DRAHOS and SZYBALSKI 1981 Down). 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 TR1 and TL1 (DAS 1992 Down), 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 PRE and PRM (Figure 1A). Since PRM 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 PRE-initiated transcripts produce enough CI to get the cell into the concentration state where PRM is activated and PR is repressed, that is, into the regime where Log[CI2] > ~-6.8 and Log[Cro2] < ~-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 PRE. 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 PRE transcripts occurs to repress PR and activate PRM, then CI concentration will tend to continue to increase automatically due to positive autoregulation leading to ever increasing repression of PR and PL by CI2. Eventually CI2 concentration will rise into the negative feedback region of the PRM repression curve and stabilize by autoregulation at a concentration in a range of 140–200 dimers per cell (REICHARDT and KAISER 1971 Down; LEVINE et al. 1979 Down), and Cro will eventually disappear due to dilution and degradation.

Though CII is produced from the same PR-initiated transcripts that encode Cro, CII accumulation is initially attenuated by termination of about 50% of the transcripts at TR1 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 Down). In the {lambda} 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 Down; BELFORT and WULFF 1974 Down). This tends to reduce the mean and peak CII concentration levels, and thus the probability of PRE 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 PRE and kickstart CI production.

If a cell reaches a state where (i) the Cro feedback loop is established, (ii) PRM and PL 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) PR and PL are repressed, there is a high probability that the cell will continue to lysogeny (MCADAMS and SHAPIRO 1995 Down).


*  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 Down 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 Down]. 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 {lambda} 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 Down. (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 Down. 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.


 
View this table:
In this window
In a new window

 
Table 1. Parameters for promoters PRE and PL


 
View this table:
In this window
In a new window

 
Table 2. Parameters for transcription and translation reactions

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 Down. 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 fe is the desired maximum fractional error in P (FELLER 1968 Down). 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 Down, 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 Down). The expected fraction of lysogens at a given API, Flysogens, 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 Down. 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 Down). 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 {lambda} gene expression are assumed to be equivalent for all cells and constant over time. For example, the effect of integration host factor on PL activity (GILADI et al. 1990 Down) 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 Down; PTASHNE and GANN 1997 Down). The rate-limiting step in transcript initiation is assumed to be the closed- to open-complex isomerization reaction (MCCLURE 1980 Down).

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

  8. An average of 10 proteins are produced per transcript for all genes (KEPES 1963 Down; SHEA and ACKERS 1985 Down). 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 Down). 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 PR- and PL-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 Cro2 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 {sigma} 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 Down).]

  12. A cell becomes committed to lysogeny if there is (i) a sufficient time-integrated concentration of CII to activate PRE, and (ii) [CI2] > [Cro2] at the end of the 35-min cell cycle. Activation of PRE 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 {lambda} promoters, Panti-Q and PI, that function in execution of the lysogenic pathway (MCADAMS and SHAPIRO 1995 Down). CI2 concentration greater than that of Cro2 at the end of the cell cycle is an additional indication that activation of PRE 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.


 
View this table:
In this window
In a new window

 
Table 3. Parameters for housekeeping and nongenetic reactions

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 {lambda} 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, TR1 and TL1, and antitermination at the NUTR and NUTL sites (Figure 1A).

Operator/promoter binding and control of transcript initiation: SHEA and ACKERS 1985 Down describe a statistical mechanical/thermodynamic approach to modeling the PR/PRM promoter complex. We use the same approach to model PL repression by CI2 and Cro2 and PRE 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 Down, GILLESPIE 1992B Down. 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 Down). 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 CI2 and Cro2 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 NUTL and NUTR sites (Figure 1A) depend on {lambda} 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 Down; MASON and GREENBLATT 1991 Down; LI et al. 1992 Down; DEVITO and DAS 1994 Down). 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 Down; LI et al. 1992 Down; DEVITO and DAS 1994 Down). 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 Down; LI et al. 1992 Down; DEVITO and DAS 1994 Down).

Translation: Translation control is modeled as described by MCADAMS and ARKIN 1997 Down, based on the mechanism described by YARCHUK et al. 1992 Down. 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 {lambda} 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 Down; YARCHUK et al. 1992 Down). 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 Down; GOTTESMAN and MAURIZI 1992 Down). 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 Down). 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 Down) and HflB (BANUETT et al. 1986 Down; HERMAN et al. 1995 Down; KIHARA et al. 1997 Down), act to degrade CII. HflB is thought to degrade CIII (HERMAN et al. 1995 Down). Initial studies identified HflA as a nonessential membrane-bound GTP-utilizing protease (CHENG et al. 1988 Down; ZORICK and ECHOLS 1991 Down; NOBLE et al. 1993 Down). 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 Down; KIHARA et al. 1997 Down; SHOTLAND et al. 1997 Down). The Hfl proteolytic system has other host protein targets in addition to CII (CHENG and ECHOLS 1987 Down; HERMAN et al. 1993 Down; HERMAN et al. 1995 Down). 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 Down; BANUETT and HERSKOWITZ 1987 Down). CIII protects CII from proteolysis (HOYT et al. 1982 Down; RATTRAY et al. 1984 Down; BANUETT et al. 1986 Down) 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 Down). 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 Down) to less than 30 sec depending on conditions (RATTRAY et al. 1984 Down). 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 Down): (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 Down; CHENG et al. 1988 Down), (b) produces an approximately 6-min half-life for CIII in the absence of CII (KORNITZER et al. 1991A Down), (c) yields CIII protection of CII consistent with HOYT et al. 1982 Down and RATTRAY et al. 1984 Down, and (d) functions in the overall simulation model to produce the simulated percent lysogeny as a function of API consistent with KOURILSKY 1973 Down. 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 Down.

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 Down). Reported data on N degradation fit a first-order decay curve (GOTTESMAN and GOTTESMAN 1981 Down); 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, Cro2 and CII are expressed earliest in both the lytic and lysogenic cases. Cro2 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 PRE 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. CI2 concentration (Figure 5A) in the lysogenic case began to grow at about 12 min just after CII concentration peaked. The growing CI2 concentration repressed PR and stopped Cro production. As a result, the Cro2 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 PRE enough to start the CI expression feedback loop. Without expression of CI, Cro2 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 PR and PRM promoter activation contours. The CI2 repressor concentration began to autoregulate its own concentration 20 min after infection; thereafter, the CI2 concentration remained constant and PRM activation slowly increased as Cro2 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 CI2 and Cro2 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{sigma} range.) The lysis- and lysogeny-fated subsets shown in Figure 3B and Figure C, each experience a different pattern of Cro2 and CI2 concentration growth statistics, distinct from each other and from the combined statistics. Figure 2D shows the same average Cro2 and CI2 concentration trajectories for the lysogenic-fated and lytic-fated cases at MOI 6 superimposed on the PR and PRM promoter activation contours.

Lysogenic fraction: Kinetic model estimates compared to experiment:
The experimental lysogeny fraction data shown in Figure 6B for starved O- ({square}) and P- ({circ}) mutants are from Figure 2 of KOURILSKY 1973 Down. 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: {blacksquare}), and, second, for the best match considering only API values <=6 (Figure 6A, labeled "Lower", symbol: {blacktriangledown}). Points in Figure 6A ({blacksquare} and {blacktriangledown}) 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 Down for fraction of lysogens produced vs. MOI. (a) Simulation results for fraction of cells producing lysogens ({blacksquare}, {blacktriangledown}). 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 ({bullet}); O-Coop-, noncooperative binding of CI dimers at OR1–3 ({diams}). (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- ({circ}) and P- ({square}) strains. Experimental points are from Figure 2 in KOURILSKY 1973 Down 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 Down 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 PRPRM operator sites is made noncooperative in the kinetic model, reducing the effectiveness of positive autoregulation of PRM. 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 TR and TL 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.


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

Stochastic gene expression and competitive genetic regulatory mechanisms:
In preceding sections a stochastic kinetic model of the {lambda} lysis-lysogeny genetic regulatory circuit is used to estimate the dynamical behavior of the circuit, including effects of random patterns of gene expression. The random developmental path choice between the lysogenic or lytic path in individual cells was shown to result from the inevitable fluctuations in the temporal pattern of protein concentration growth caused by the molecular-level thermal fluctuations in rates of rate-determining reactions within gene expression mechanisms. The resulting differences in concentration between the regulatory proteins controlling the bistable switching elements of the decision circuit led to different path selections in different cells. The estimated variation with API of the fraction of a phage {lambda}-infected cell population that become lysogenic was shown to be consistent with experimental observations.

This analysis indicates how molecular level thermal fluctuations can be exploited by the regulatory circuit designs of developmental switches to produce different phenotypic outcomes. Such regulatory mechanisms will produce diverse phenotypes even in clonal cell populations maintained in the most homogeneous laboratory environments. In such systems environmental signals can act on the parameters of the regulatory circuit to bias the probabilities of path choice under different conditions.

Even in nondifferentiated cell populations, the concentration of regulatory proteins within individual cells will vary widely from the average concentration measured by laboratory procedures. In situations where there are divergent phenotype subpopulations as analyzed herein, the concentration trajectories of each subpopulation can differ radically from each other and from the average measured for the full population (Figure 3). For these situations, experimental methods (such as various cell sorting techniques) that profile the distribution of individual cell parameters are necessary. In any case, the function of regulatory circuits that determine cell fates is determined by protein concentrations within each cell, and the temporal pattern of these intracellular concentration trajectories is completely different from any averaged measurement (compare Figure 3 with Figure 4 and Figure 5).

Role of termination sites in the {lambda} circuit:
The simulation results with hypothetical "digital-mutations" affecting termination effectiveness of the TR1 and TL1 termination sites demonstrate that both removal of these sites and removal of the antitermination effects of the N protein have major effects on the level of lysogeny, but in opposing directions (Figure 6A). We conjecture that the function of these phage-encoded features in the phage regulatory circuit design may be to adjust the level of the phage lysogenic response to an "optimum" range for phage survival.

Other stochastic switching mechanisms:
Although the specific analysis reported here deals with regulation of the phage {lambda} infection, regulatory circuits based on bistable genetic regulatory mechanisms are used in many organisms to produce subpopulations of distinct phenotypes by random phenotype switching. Table 4 shows a small sample of well-known cases of bistable regulatory mechanisms in regulatory circuits that produce stochastic phenotype outcomes. Many examples are found in pathogenic organisms. Conventional deterministic kinetics does not model statistics of regulatory systems that produce probabilistic outcomes. A stochastic kinetic analysis as used in this paper for the {lambda} decision circuit can be used to predict statistics of regulatory outcomes for these stochastically regulated systems. Such stochastic kinetic analyses may also permit improved exploitation of information in the statistics of phenotypic outcomes.


 
View this table:
In this window
In a new window

 
Table 4. Examples of other cases of bistable regulatory mechanisms producing stochastic phenotypic outcomes

For bistable switching phenomena that are integral to the operation of the cell, some complications encountered in analysis of the {lambda} decision circuit would not be present. One example is the uncertainty in phage gene dosage and of timing of infection in the cell cycle due to the random nature of the phage infection process. Also, phage infection is disruptive of "normal" cell processes, whereas the cell's integral processes would presumably be less so.

Switching dynamics:
We expect that many aspects of the {lambda} switch dynamics will be found in other multi-potent regulatory switches. Examples include: (i) transient, low-level expression of key regulatory proteins (as with CII, CIII, and N), (ii) stochastic progress toward commitment as concentrations of the controlling proteins in each cell change from moment to moment and the probability of each potential outcome also changes until eventually the cell's fate is determined, and (iii) design features that adjust timing of reactions as with the delay of CII and CIII concentration growth produced by termination sites TL1 and TR1.

Stochastic gene expression and regulatory determinism:
If random processes of gene expression tend to make the pattern of protein production inherently erratic, how do cells achieve the regulatory determinism necessary for most functions? One possibility is that the overall regulatory architecture can suppress deleterious effects of molecular-level reaction fluctuations. For example, consider a case where expression of a protein is controlled by a rate-limiting mechanism acting at the level of translation. For both the conventional and the stochastic formulation of chemical kinetics, the kinetics of the rate-limiting step(s) will dominate the overall kinetics of a series of cascaded chemical reactions. So, stochastic processes affecting transcription control, posttranscriptional editing, or message transport may be screened by the final translation control mechanism. Another possibility is that the stochastic pattern of signal protein production may only cause uncertainty in timing of regulatory events, not uncertainty in outcome. Within broad limits the duration of many cellular functions may be less important to proper cellular function than the proper sequencing of events. For example, cells halt at various checkpoints until conditions (e.g., restoration of essential nutrients, completion of precursor cellular events) for further progress are satisfied (HARTWELL and WEINERT 1989 Down; KAUFMANN and PAULES 1996 Down; WELLS 1996 Down). In this case, the indeterminism relates to whether the cell will progress or not progress along a developmental path at any instant, rather than to a choice among alternate stable pathways. So, the regulatory decision logic is: "HALT until CONDITIONS are met then PROCEED," where "CONDITIONS" are sensed environmental or cellular signals. The result is dispersion across the cell population in the rate of progression along prescribed pathways rather than dispersion in outcome.


*  FOOTNOTES

1 Present address: Division of Physical Biosciences, Calvin Labs 144, Lawrence Berkeley National Laboratory, Berkeley, CA 94720. Back


*  ACKNOWLEDGMENTS

This work was supported by Office of Naval Research Grant N00014-96-1-0564. A.A. was partially supported by National Science Foundation grant CHE9109301. The work was also supported in part by a grant of computer time from the Department of Defense High Performance Computing Centers at Eglin Air Force Base and Maui. We thank LUCY SHAPIRO for valuable comments on the manuscript.

Manuscript received March 5, 1998; Accepted for publication April 30, 1998.


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

ADHYA, S. and M. GOTTESMAN, 1982  Promoter occlusion: transcription through a promoter may inhibit its activity. Cell 29:939-944[Medline].

BANUETT, F. and I. HERSKOWITZ, 1987  Identification of polypeptides encoded by an Escherichia coli locus (hflA) that governs the lysis-lysogeny decision of bacteriophage {lambda}. J. Bacteriol. 169:4076-4085[Abstract/Free Full Text].

BANUETT, F., M. A. HOYT, L. MCFARLANE, H. ECHOLS, and I. HERSKOWITZ, 1986  hflB, a new Escherichia coli locus regulating lysogeny and the level of bacteriophage lambda CII protein. J. Mol. Biol. 187:213-224[Medline].

BELFORT, M. and D. L. WULFF, 1974  The roles of lambda cIII gene and the Escherichia coli catabolite gene activation system in the establishment of lysogeny by bacteriophage lambda. Proc. Natl. Acad. Sci. USA 71:779-782[Abstract/Free Full Text].

BURZ, D. S., D. BECKETT, N. BENSON, and G. K. ACKERS, 1994  Self-assembly of bacteriophage {lambda} cI repressor: effects of single site mutations on the monomer-dimer equilibrium. Biochemistry 33:8399-8405[Medline].

CHENG, H. H. and H. ECHOLS, 1987  A class of Escherichia coli proteins controlled by the hflA locus. J. Mol. Biol. 196:737-740[Medline].

CHENG, H. H., P. J. MUHLRAD, M. A. HOYT, and H. ECHOLS, 1988  Cleavage of the CII protein of phage lambda by purified HflA protease: control of the switch between lysis and lysogeny. Proc. Natl. Acad. Sci. USA 85:7882-7886[Abstract/Free Full Text].

DAMBLY-CHAUDIERE, C., M. GOTTESMAN, C. DEBOUK, and S. ADYHA, 1983  Regulation of the pR operon of bacteriophage lambda. J. Mol. Appl. Genet. 2:45-56[Medline].

DAS, A., 1992  How the phage lambda N gene product suppresses transcription termination: communication of RNA polymerase with regulatory proteins mediated by signals in nascent RNA. J. Bacteriol. 174:6711-6716[Free Full Text].

DEVITO, J. and A. DAS, 1994  Control of transcription processivity in phage {lambda}: Nus factors strengthen the termination-resistant state of RNA polymerase induced by N antiterminator. Proc. Natl. Acad. Sci. USA 91:8660-8664[Abstract/Free Full Text].

DRAHOS, D. and W. SZYBALSKI, 1981  Antitermination and termination functions of the cloned nutL, N, and tL1 modules of coliphage lambda. Gene 16:261-274[Medline].

ECHOLS, H., 1986  Multiple DNA-protein interactions governing high-precision DNA transactions. Science 233:1050-1056[Abstract/Free Full Text].

ELLIS, E. L. and M. DELBRUCK, 1939  The growth of the bacteriophage. J. Gen. Physiol. 22:365-384[Abstract/Free Full Text].

FELLER, W., 1968 An Introduction to Probability Theory. John Wiley & Sons, New York.

FINLAY, B. B. and S. FALKOW, 1997  Common themes in microbial pathogenicity revisited. Microbiol. Mol. Biol. Rev. 61:136-169[Abstract].

FRIEDMAN, D. I., 1992  Interaction between bacteriophage lambda and its Escherichia coli host. Curr. Opin. Genet. Dev. 2:727-738[Medline].

FRIEDMAN, D., and M. GOTTESMAN, 1983 Lytic mode of lambda development, pp. 21–52 in Lambda II, edited by R. HENDRIX, J. W. ROBERTS, F. W. STAHL and R. A. WEISBERG. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY.

GILADI, H., M. GOTTESMAN, and A. B. OPPENHEIM, 1990  Integration host factor stimulates the phage lambda PL promoter. J. Mol. Biol. 213:109-121[Medline].

GILLESPIE, D. T., 1976  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22:403-434.

GILLESPIE, D. T., 1977  Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81(25):2340-2361.

GILLESPIE, D. T., 1992a Markov Processes: An Introduction for Physical Scientists. Academic Press, San Diego.

GILLESPIE, D. T., 1992b  A rigorous derivation of the chemical master equation. Phys. A 188:404-425.

GOTTA, S. L., O. L. MILLER, JR., and S. L. FRENCH, 1991  Ribosomal RNA transcription rate in Escherichia coli. J. Bacteriol. 173:6647-6649[Abstract/Free Full Text].

GOTTESMAN, S. and M. GOTTESMAN, 1981  Protein degradation in E. coli: the lon mutation and bacteriophage lambda N and CII protein stability. Cell 24:225-233[Medline].

GOTTESMAN, S. and M. R. MAURIZI, 1992  Regulation by proteolysis: energy-dependent proteases and their targets. Microbiol. Rev. 56:592-621[Abstract/Free Full Text].

GRODZIKER, T., R. R. ARDITTI, and H. EISEN, 1972  Establishment of repression by lambdoid phage in catabolic activating protein and adenylate cyclase mutants of Escherichia coli.. Proc. Natl. Acad. Sci. USA 69:366-370[Abstract/Free Full Text].

GUPTASARMA, P., 1995  Does replication-induced transcription regulate synthesis of the myriad low copy number proteins of Escherichia coli? BioEssays 17:987-997[Medline].

HARTWELL, L. H. and T. A. WEINERT, 1989  Checkpoints: controls that ensure the order of cell cycle events. Science 246:629-634[Abstract/Free Full Text].

HERMAN, C., T. OGURA, T. TOMOYASU, S. HIRAGA, and Y. AKIYAMA et al., 1993  Cell growth and lambda phage development controlled by the same Escherichia coli gene, ftsH/hflB.. Proc. Natl. Acad. Sci. USA 90:10861-10865[Abstract/Free Full Text].

HERMAN, C., D. THEVENET, R. D'ARL, and P. BOULOC, 1995  Degradation of {sigma}32, the heat shock regulator in Escherichia coli, is governed by HflB. Proc. Natl. Acad. Sci. USA 92:3516-3520[Abstract/Free Full Text].

HERSKOWITZ, I. and D. HAGEN, 1980  The lysis-lysogeny decision of phage lambda: explicit programming and responsiveness. Annu. Rev. Genet. 14:399-445[Medline].

HOYT, M. A., D. M. KNIGHT, A. DAS, H. I. MILLER, and H. ECHOLS, 1982  Control of phage lambda development by stability and synthesis of CII protein: role of the viral cIII and host hflA, himA and himD genes. Cell 31:565-573[Medline].

ISHIHARA, A., J. E. SEGALL, S. M. BLOCK, and H. C. BERG, 1983  Coordination of flagella on filamentous cells of Escherichia coli.. J. Bacteriol. 155:228-237[Abstract/Free Full Text].

KAUFMANN, W. K. and R. S. PAULES, 1996  DNA damage and cell cycle checkpoints. FASEB J. 10:238-247[Abstract].

KENNELL, D. and H. RIEZMAN, 1977  Transcription and translation initiation frequencies of the Escherichia coli lac operon. J. Mol. Biol. 114:1-21[Medline].

KEPES, A., 1963  Kinetics of induced enzyme synthesis: determination of the mean life of the Galactosidase-specific messenger RNA. Biochem. Biophys. Acta 76:293-309[Medline].

KIHARA, A., Y. AKIYAMA, and K. ITO, 1997  Host regulation of lysogenic decision in bacteriophage lambda: transmembrane modulation of FtsH (HflB), the cII degrading protease, by HflKC (HflA). Proc. Natl. Acad. Sci. USA 94:5544-5549[Abstract/Free Full Text].

KORNBERG, A., and T. A. BAKER, 1992 DNA Replication. W. H. Freeman, New York.

KORNITZER, D., S. ALTUVIA, and A. B. OPPENHEIM, 1991a  The activity of the CIII regulator of lambdoid bacteriophages resides within a 24-amino acid protein domain. Proc. Natl. Acad. Sci. USA 88:5217-5221[Abstract/Free Full Text].

KORNITZER, D., S. ALTUVIA, and A. B. OPPENHEIM, 1991b  Genetic analysis of the cIII gene of bacteriophage HK022. J. Bacteriol. 173:810-815[Abstract/Free Full Text].

KOURILSKY, P., 1973  Lysogenization by bacteriophage lambda: I. Multiple infection and the lysogenic response. Mol. Gen. Genet. 122:183-195[Medline].

LEVINE, A., A. BAILONE, and R. DEVORET, 1979  Cellular levels of the prophage lambda and 434 repressors. J. Mol. Biol. 131:655-661[Medline].

LI, J., R. HORWITZ, S. MCCRACKEN, and J. GREENBLATT, 1992  NusG, a new Escherichia coli elongation factor involved in transcriptional antitermination by the N protein of phage {lambda}. J. Biol. Chem. 267:6012-6019[Abstract/Free Full Text].

MARRS, C. F., W. W. RUEHL, G. K. SCHOOLNIK, and S. FALKOW, 1988  Pilin-gene phase variation of Moraxella bovis is caused by an inversion of the pilin genes. J. Bacteriol. 170:3032-3039[Abstract/Free Full Text].

MASON, S. W. and J. GREENBLATT, 1991  Assembly of transcription elongation complexes containing the N protein of phage {lambda} and the Escherichia coli elongation factors NusA, NusB, NusG and S10. Genes Dev. 5:1504-1512[Abstract/Free Full Text].

MAURER, R., B. MEYER, and M. PTASHNE, 1980  Gene regulation at the right operator (OR) bacteriophage lambda. I. OR3 and autogenous negative control by repressor. J. Mol. Biol. 139:147-161[Medline].

MCADAMS, H. and A. ARKIN, 1997  Stochastic mechanisms in gene expression. Proc. Natl. Acad. Sci. USA 94:814-819[Abstract/Free Full Text].

MCADAMS, H. H. and L. SHAPIRO, 1995  Circuit simulation of genetic networks. Science 269:650-656[Abstract/Free Full Text].

MCCLURE, W. R., 1980  Rate-limiting steps in RNA chain initiation. Proc. Natl. Acad. Sci. USA 77:5634-5638[Abstract/Free Full Text].

MCCLURE, W. R., 1983 A biochemical analysis of the effect of RNA polymerase concentration on the in vivo control of RNA chain initiation frequency, pp. 207–217 in Biochemistry of Metabolic Processes, Proceedings of the 12th Steenbock Symposium, edited by D. LENNON, F. STRATMAN and R. ZAHLMAN. Elsevier Scientific, New York.

MCQUARRIE, D. A., C. J. JACHIMOWSKI, and M. E. RUSSELL, 1964  Kinetics of small systems II. J. Chem. Phys. 40:2914-2921.

MEYER, B. J. and M. PTASHNE, 1980  Gene regulation at the right operator (OR) of bacteriophage lambda. III. Lambda repressor directly activates gene transcription. J. Mol. Biol. 139:195-205[Medline].

MEYER, B. J., R. MAURER, and M. PTASHNE, 1980  Gene regulation at the right operator (OR) of bacteriophage lambda. II. OR1, OR2, and OR3: their roles in mediating the effects of repressor and cro. J. Mol. Biol. 139:163-194[Medline].

NOBLE, J. A., M. A. INNIS, E. V. KOONIN, K. E. RUDD, and F. BANUETT et al., 1993  The Escherichia coli hfla locus encodes a putative gtp-binding protein and two membrane proteins, one of which contains a protease-like domain. Proc. Natl. Acad. Sci. USA 90:10866-10870[Abstract/Free Full Text].

PLANK, L. D. and J. D. HARVEY, 1979  Generation time statistics of Escherichia coli B measured by synchronous culture techniques. J. Gen. Microbiol. 115:69-77[Abstract/Free Full Text].

PTASHNE, M., 1992 A Genetic Switch: Phage {lambda} and Higher Organisms. Cell Press and Blackwell Scientific Publications, Cambridge, MA.

PTASHNE, M. and A. GANN, 1997  Transcriptional activation by recruitment. Nature 386:569-577[Medline].

PUTTE, V. D. P. and N. GOOSEN, 1992  DNA inversions in phages and bacteria. Trends Genet. 8:457-462[Medline].

RATTRAY, A., S. ALTUVIA, G. MAHAJNA, A. B. OPPENHEIM, and M. GOTTESMAN, 1984  Control of bacteriophage lambda CII activity by bacteriophage and host functions. J. Bacteriol. 159:238-242[Abstract/Free Full Text].

REICHARDT, L. and A. D. KAISER, 1971  Control of {lambda} repressor synthesis. Proc. Natl. Acad. Sci. USA 68:2185-2189[Abstract/Free Full Text].

REINITZ, J. and J. R. VAISNYS, 1990  Theoretical and experimental analysis of the phage lambda genetic switch implies missing levels of cooperativity. J. Theor. Biol. 145:295-318[Medline].

ROBERTSON, B. D., 1992  Genetic variation in pathogenic bacteria. Trends Genet. 8:422-427[Medline].

SAUER, R. T., 1979 Molecular Characterization of the Lambda Repressor and Its Gene CI. Harvard University Press, Cambridge, MA.

SHEA, M. A. and G. K. ACKERS, 1985  The OR control system of bacteriophage lambda: a physical-chemical model for gene regulation. J. Mol. Biol. 181:211-230[Medline].

SHIH, M.-C. and G. N. GUSSIN, 1983  Differential effects of mutations on discrete steps in transcription initiation at the lambda PRE promotor. Cell 34:941-949[Medline].

SHIH, M.-C. and G. N. GUSSIN, 1984  Kinetic analysis of mutations affecting the CII activation site at the PRE promoter of bacteriophage {lambda}. Proc. Natl. Acad. Sci. USA 81:6432-6436[Abstract/Free Full Text].

SHOTLAND, Y., S. KOBY, D. TEFF, N. MANSUR, and D. A. OREN et al., 1997  Proteolysis of the phage {lambda} CII regulatory protein by Ftsh (HflB) of Escherichia coli.. Mol. Microbiol. 24:1303-1310[Medline].

SORENSEN, M. A. and S. PEDERSEN, 1991  Absolute in vivo translation rates of individual codons in Escherichia coli. J. Mol. Biol. 222:265-280[Medline].

STRAUSS, E. J. and S. FALKOW, 1997  Microbial pathogenesis: genomics and beyond. Science 276:707-712[Abstract/Free Full Text].

VAN KAMPEN, N. G., 1992 Stochastic Processes in Physics and Chemistry. North-Holland, Amsterdam.

VOGEL, U. and K. F. JENSEN, 1994  The RNA chain elongation rate in Escherichia coli depends on the growth rate. J. Bacteriol. 176:2807-2813[Abstract/Free Full Text].

WELLS, W. A. E., 1996  The spindle-assembly checkpoint: aiming for a perfect mitosis every time. Trends Cell Biol. 6:228-234.

WHALEN, W., B. GHOSH, and A. DAS, 1988  NusA protein is necessary and sufficient in vitro for phage {lambda} N gene product to suppress a {rho}-independent terminator placed downstream of nutL.. Proc. Natl. Acad. Sci. USA 85:2494-2498[Abstract/Free Full Text].

WOUDE, M. V. D., B. BRAATEN, and D. LOW, 1996  Epigenetic phase variation of the pap operon in Escherichia coli.. Trends Microbiol. 4:5-9[Medline].

YARCHUK, O., N. JACQUES, J. GUILLEREZ, and M. DREYFUS, 1992  Interdependence of translation, transcription and mRNA degradation in the lacZ gene. J. Mol. Biol. 226:581-596[Medline].

ZHENG, Q. and J. ROSS, 1991  Comparison of deterministic and stochastic kinetics for nonlinear systems. J. Chem. Phys. 94:3644-3648.

ZORICK, T. S. and H. ECHOLS, 1991  Membrane localization of HflA regulatory protease of Escherichia coli by immunoelectron microscopy. J. Bacteriol. 173:6307-6310[Abstract/Free Full Text].




This article has been cited by other articles:


Home page
Proc. Natl. Acad. Sci. USAHome page
M. J. Morelli, P. R. ten Wolde, and R. J. Allen
DNA looping provides stability and robustness to the bacteriophage {lambda} switch
PNAS, May 19, 2009; 106(20): 8101 - 8106.
[Abstract] [Full Text] [PDF]


Home page
J R Soc InterfaceHome page
M. Vellela and H. Qian
Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlogl model revisited
J R Soc Interface, February 9, 2009; (2009) rsif.2008.0476v1.
[Abstract] [Full Text] [PDF]


Home page
Brief BioinformHome page
J. Pahle
Biochemical simulations: stochastic, approximate stochastic and hybrid approaches
Brief Bioinform, January 16, 2009; (2009) bbn050v1.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
A. Prasad, J. Zikherman, J. Das, J. P. Roose, A. Weiss, and A. K. Chakraborty
Origin of the sharp boundary that discriminates positive and negative selection of thymocytes
PNAS, January 13, 2009; 106(2): 528 - 533.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
F. St-Pierre and D. Endy
Determination of cell fate selection during phage lambda infection
PNAS, December 30, 2008; 105(52): 20705 - 20710.
[Abstract] [Full Text] [PDF]


Home page
ScienceHome page
P. J. Choi, L. Cai, K. Frieda, and X. S. Xie
A Stochastic Single-Molecule Event Triggers Phenotype Switching of a Bacterial Cell
Science, October 17, 2008; 322(5900): 442 - 446.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
M. Turcotte, J. Garcia-Ojalvo, and G. M. Suel
A genetic timer through noise-induced stabilization of an unstable state
PNAS, October 14, 2008; 105(41): 15732 - 15737.
[Abstract] [Full Text] [PDF]


Home page
J R Soc InterfaceHome page
J. J Tyson, R. Albert, A. Goldbeter, P. Ruoff, and J. Sible
Biological switches and clocks
J R Soc Interface, August 6, 2008; 5(Suppl_1): S1 - S8.
[Abstract] [Full Text] [PDF]


Home page
J R Soc InterfaceHome page
D. Gonze, M. Jacquet, and A. Goldbeter
Stochastic modelling of nucleocytoplasmic oscillations of the transcription factor Msn2 in yeast
J R Soc Interface, August 6, 2008; 5(Suppl_1): S95 - S109.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
A. Sanchez and J. Kondev
Transcriptional control of noise in gene expression
PNAS, April 1, 2008; 105(13): 5081 - 5086.
[Abstract] [Full Text] [PDF]


Home page
DevelopmentHome page
A. C. Miller, H. Seymour, C. King, and T. G. Herman
Loss of seven-up from Drosophila R1/R6 photoreceptors reveals a stochastic fate choice that is normally biased by Notch
Development, February 15, 2008; 135(4): 707 - 715.
[Abstract] [Full Text] [PDF]


Home page
BioinformaticsHome page
T. Fournier, J.P. Gabriel, C. Mazza, J. Pasquier, J.L. Galbete, and N. Mermod
Steady-state expression of self-regulated genes
Bioinformatics, December 1, 2007; 23(23): 3185 - 3192.
[Abstract] [Full Text] [PDF]


Home page
Phil Trans R Soc BHome page
F. Mavelli and K. Ruiz-Mirazo
Stochastic simulations of minimal self-reproducing cellular systems
Phil Trans R Soc B, October 29, 2007; 362(1486): 1789 - 1802.
[Abstract] [Full Text] [PDF]


Home page
J. Biol. Chem.Home page
Y.-K. Oh, B. O. Palsson, S. M. Park, C. H. Schilling, and R. Mahadevan
Genome-scale Reconstruction of Metabolic Network in Bacillus subtilis Based on High-throughput Phenotyping and Gene Essentiality Data
J. Biol. Chem., September 28, 2007; 282(39): 28791 - 28799.
[Abstract] [Full Text] [PDF]


Home page
Brief BioinformHome page
B. S. Srinivasan, N. H. Shah, J. A. Flannick, E. Abeliuk, A. F. Novak, and S. Batzoglou
Current progress in network research: toward reference networks for key model organisms
Brief Bioinform, September 1, 2007; 8(5): 318 - 332.
[Abstract] [Full Text] [PDF]


Home page
Nucleic Acids ResHome page
J. Hu, S. C. Sealfon, F. Hayot, C. Jayaprakash, M. Kumar, A. C. Pendleton, A. Ganee, A. Fernandez-Sesma, T. M. Moran, and J. G. Wetmur
Chromosome-specific and noisy IFNB1 transcription in individual virus-infected human primary dendritic cells
Nucleic Acids Res., August 2, 2007; (2007) gkm557v1.
[Abstract] [Full Text] [PDF]


Home page
BioinformaticsHome page
M. Dobrzynski, J. V. Rodriguez, J. A. Kaandorp, and J. G. Blom
Computational methods for diffusion-influenced biochemical reactions
Bioinformatics, August 1, 2007; 23(15): 1969 - 1977.
[Abstract] [Full Text] [PDF]


Home page
Nucleic Acids ResHome page
V. B. Teif
General transfer matrix formalism to calculate DNA-protein-drug binding in gene regulation: application to OR operator of phage {lambda}
Nucleic Acids Res., June 28, 2007; 35(11): e80 - e80.
[Abstract] [Full Text] [PDF]


Home page
Brief BioinformHome page
D. J. Wilkinson
Bayesian methods in bioinformatics and computational systems biology
Brief Bioinform, April 12, 2007; (2007) bbm007v1.
[Abstract] [Full Text] [PDF]


Home page
BioinformaticsHome page
T. Tian, S. Xu, J. Gao, and K. Burrage
Simulated maximum likelihood method for estimating kinetic rates in gene expression
Bioinformatics, January 1, 2007; 23(1): 84 - 91.
[Abstract] [Full Text] [PDF]


Home page
Sci SignalHome page
M. S. Samoilov, G. Price, and A. P. Arkin
From Fluctuations to Phenotypes: The Physiology of Noise
Sci. Signal., December 19, 2006; 2006(366): re17 - re17.
[Abstract] [Full Text] [PDF]


Home page
Brief BioinformHome page
K. Burrage, L. Hood, and M. A. Ragan
Advanced computing for systems biology
Brief Bioinform, December 1, 2006; 7(4): 390 - 398.
[Abstract] [Full Text] [PDF]


Home page
J R Soc InterfaceHome page
F. J Doyle III and J. Stelling
Systems interface biology
J R Soc Interface, October 22, 2006; 3(10): 603 - 616.
[Abstract] [Full Text] [PDF]


Home page
BioinformaticsHome page
J. V. Rodriguez, J. A. Kaandorp, M. Dobrzynski, and J. G. Blom
Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (PTS) pathway in Escherichia coli
Bioinformatics, August 1, 2006; 22(15): 1895 - 1901.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
T. Tian and K. Burrage
Stochastic models for regulatory networks of the genetic toggle switch
PNAS, May 30, 2006; 103(22): 8372 - 8377.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
J. T. Mettetal, D. Muzzey, J. M. Pedraza, E. M. Ozbudak, and A. van Oudenaarden
Predicting stochastic gene expression dynamics in single cells
PNAS, May 9, 2006; 103(19): 7304 - 7309.
[Abstract] [Full Text] [PDF]


Home page
Phil Trans R Soc BHome page
S Ramsey, A Ozinsky, A Clark, K.D Smith, P de Atauri, V Thorsson, D Orrell, and H Bolouri
Transcriptional noise and cellular heterogeneity in mammalian macrophages
Phil Trans R Soc B, March 29, 2006; 361(1467): 495 - 506.
[Abstract] [Full Text] [PDF]


Home page
J. Virol.Home page
T. Subramanian, S. Vijayalingam, and G. Chinnadurai
Genetic Identification of Adenovirus Type 5 Genes That Influence Viral Spread
J. Virol., February 15, 2006; 80(4): 2000 - 2012.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
D. Bratsun, D. Volfson, L. S. Tsimring, and J. Hasty
Delay-induced stochastic oscillations in gene regulation
PNAS, October 11, 2005; 102(41): 14593 - 14598.
[Abstract] [Full Text] [PDF]


Home page
ScienceHome page
J. M. Raser and E. K. O'Shea
Noise in Gene Expression: Origins, Consequences, and Control
Science, September 23, 2005; 309(5743): 2010 - 2013.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
R. J. Johnston Jr., S. Chang, J. F. Etchberger, C. O. Ortiz, and O. Hobert
From the Cover: MicroRNAs acting in a double-negative feedback loop to control a neuronal cell fate decision
PNAS, August 30, 2005; 102(35): 12449 - 12454.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
N. E. Buchler, U. Gerland, and T. Hwa
Nonlinear protein degradation and the function of genetic circuits
PNAS, July 5, 2005; 102(27): 9559 - 9564.
[Abstract] [Full Text] [PDF]


Home page
BioinformaticsHome page
A. Chatterjee, K. Mayawala, J. S. Edwards, and D. G. Vlachos
Time accelerated Monte Carlo simulations of biological networks using the binomial {tau}-leap method
Bioinformatics, May 1, 2005; 21(9): 2136 - 2137.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
O. Kobiler, A. Rokney, N. Friedman, D. L. Court, J. Stavans, and A. B. Oppenheim
Quantitative kinetic analysis of the bacteriophage {lambda} genetic network
PNAS, March 22, 2005; 102(12): 4470 - 4475.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
S. Hooshangi, S. Thiberge, and R. Weiss
Ultrasensitivity and noise propagation in a synthetic transcriptional cascade
PNAS, March 8, 2005; 102(10): 3581 - 3586.
[Abstract] [Full Text] [PDF]


Home page
GeneticsHome page
C. A. Voigt, D. M. Wolf, and A. P. Arkin
The Bacillus subtilis sin Operon: An Evolvable Network Motif
Genetics, March 1, 2005; 169(3): 1187 - 1202.
[Abstract] [Full Text] [PDF]


Home page
J. Bacteriol.Home page
M. B. Yarmolinsky
Bacteriophage P1 in Retrospect and in Prospect
J. Bacteriol., November 1, 2004; 186(21): 7025 - 7028.
[Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
G. Ashkenasy, R. Jagasia, M. Yadav, and M. R. Ghadiri
Design of a directed molecular network
PNAS, July 27, 2004; 101(30): 10872 - 10877.
[Abstract] [Full Text] [PDF]


Home page
ScienceHome page
J. M. Raser and E. K. O'Shea
Control of Stochasticity in Eukaryotic Gene Expression
Science, June 18, 2004; 304(5678): 1811 - 1814.
[Abstract] [Full Text] [PDF]


Home page
SIMULATIONHome page
J. M. McCollum, G. D. Peterson, C. D. Cox, and M. L. Simpson
Accelerating Gene Regulatory Network Modeling Using Grid-Based Simulation
SIMULATION, May 1, 2004; 80(4-5): 231 - 241.
[Abstract] [PDF]


Home page
GeneticsHome page
M. Thattai and A. van Oudenaarden
Stochastic Gene Expression in Fluctuating Environments
Genetics, May 1, 2004; 167(1): 523 - 530.
[Abstract] [Full Text] [PDF]


Home page
GeneticsHome page
J. Peccoud, K. V. Velden, D. Podlich, C. Winkler, L. Arthur, and M. Cooper
The Selective Values of Alleles in a Molecular Network Model Are Context Dependent
Genetics, April 1, 2004; 166(4): 1715 - 1725.
[Abstract] [Full Text] [PDF]


Home page
J. Bacteriol.Home page
G. Bertani
Lysogeny at Mid-Twentieth Century: P1, P2, and Other Experimental Systems
J. Bacteriol., February 1, 2004; 186(3): 595 - 600.
[Full Text] [PDF]


Home page
Genes Dev.Home page
I. B. Dodd, K. E. Shearwin, A. J. Perkins, T. Burr, A. Hochschild, and J. B. Egan
Cooperativity in long-range gene regulation by the {lambda} CI repressor
Genes & Dev., February 1, 2004; 18(3): 344 - 354.
[Abstract] [Full Text] [PDF]


Home page
Genome ResHome page
P. M. Kim and B. Tidor
Limitations of Quantitative Gene Regulation Models: A Case Study
Genome Res., November 1, 2003; 13(11): 2391 - 2395.
[Abstract] [Full Text] [PDF]


Home page
Genome ResHome page
D. E. Zak, G. E. Gonye, J. S. Schwaber, and F. J. Doyle III
Importance of Input Perturbations and Stochastic Gene Expression in the Reverse Engineering of Genetic Regulatory Networks: Insights From an Identifiability Analysis of an In Silico Network
Genome Res., November 1, 2003; 13(11): 2396 - 2405.
[Abstract] [Full Text] [PDF]


Home page
GeneticsHome page
M. Louis, L. Holm, L. Sanchez, and M. Kaufman
A Theoretical Model for the Regulation of Sex-lethal, a Gene That Controls Sex Determination and Dosage Compensation in Drosophila melanogaster
Genetics, November 1, 2003; 165(3): 1355 - 1384.
[Abstract] [Full Text] [PDF]


Home page
J. Bacteriol.Home page
V. Sentchilo, R. Ravatn, C. Werlen, A. J. B. Zehnder, and J. R. van der Meer
Unusual Integrase Gene Expression on the clc Genomic Island in Pseudomonas sp. Strain B13
J. Bacteriol., August 1, 2003; 185(15): 4530 - 4538.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
F. J. Isaacs, J. Hasty, C. R. Cantor, and J. J. Collins
Prediction and measurement of an autoregulatory genetic module
PNAS, June 24, 2003; 100(13): 7714 - 7719.
[Abstract] [Full Text] [PDF]


Home page
Plant Physiol.Home page
K. E. Thum, D. E. Shasha, L. V. Lejay, and G. M. Coruzzi
Light- and Carbon-Signaling Pathways. Modeling Circuits of Interactions
Plant Physiology, June 1, 2003; 132(2): 440 - 452.
[Abstract] [Full Text] [PDF]


Home page
J. Bacteriol.Home page
J. L. Reed and B. O. Palsson
Thirteen Years of Building Constraint-Based In Silico Models of Escherichia coli
J. Bacteriol., May 1, 2003; 185(9): 2692 - 2699.
[Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
D. Battogtokh, D. K. Asch, M. E. Case, J. Arnold, and H.-B. Schuttler
An ensemble method for identifying regulatory circuits with special reference to the qa gene cluster of Neurospora crassa
PNAS, December 24, 2002; 99(26): 16904 - 16909.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
P. S. Swain, M. B. Elowitz, and E. D. Siggia
Intrinsic and extrinsic contributions to stochasticity in gene expression
PNAS, October 1, 2002; 99(20): 12795 - 12800.
[Abstract] [Full Text] [PDF]


Home page
Nucleic Acids ResHome page
B. A. Sokhansanj, G. R. Rodrigue, J. P. Fitch, and D. M. W. III
A quantitative model of human DNA base excision repair. I. mechanistic insights
Nucleic Acids Res., April 15, 2002; 30(8): 1817 - 1825.
[Abstract] [Full Text] [PDF]


Home page
J. Bacteriol.Home page
L. You, P. F. Suthers, and J. Yin
Effects of Escherichia coli Physiology on Growth of Phage T7 In Vivo and In Silico
J. Bacteriol., April 1, 2002; 184(7): 1888 - 1894.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
D. McMillen, N. Kopell, J. Hasty, and J. J. Collins
Synchronizing genetic relaxation oscillators by intercell signaling
PNAS, January 22, 2002; 99(2): 679 - 684.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
M. Thattai and A. van Oudenaarden
Intrinsic noise in gene regulatory networks
PNAS, June 28, 2001; (2001) 151588598.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
J. Paulsson, O. G. Berg, and M. Ehrenberg
Stochastic focusing: Fluctuation-enhanced sensitivity of intracellular regulation
PNAS, June 13, 2000; (2000) 110057697.
[Abstract] [Full Text]


Home page
ScienceHome page
C. Sander
Genomic Medicine and the Future of Health Care
Science, March 17, 2000; 287(5460): 1977 - 1978.
[Abstract] [Full Text]


Home page
J. Biol. Chem.Home page
A. M. Kierzek, J. Zaim, and P. Zielenkiewicz
The Effect of Transcription and Translation Initiation Frequencies on the Stochastic Fluctuations in Prokaryotic Gene Expression
J. Biol. Chem., March 9, 2001; 276(11): 8165 - 8172.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
D. Gonze, J. Halloy, and A. Goldbeter
Robustness of circadian rhythms with respect to molecular noise
PNAS, January 22, 2002; 99(2): 673 - 678.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
D. Endy, L. You, J. Yin, and I. J. Molineux
Computation, prediction, and experimental tests of fitness for bacteriophage T7 mutants with permuted genomes
PNAS, May 9, 2000; 97(10): 5375 - 5380.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
J. Paulsson, O. G. Berg, and M. Ehrenberg
Stochastic focusing: Fluctuation-enhanced sensitivity of intracellular regulation
PNAS, June 20, 2000; 97(13): 7148 - 7153.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
M. Thattai and A. van Oudenaarden
Intrinsic noise in gene regulatory networks
PNAS, July 17, 2001; 98(15): 8614 - 8619.
[Abstract] [Full Text] [PDF]