DNA looping provides stability and robustness to the bacteriophage (cid:1) switch

The bistable gene regulatory switch controlling the transition from lysogeny to lysis in bacteriophage (cid:1) presents a unique challenge to quantitative modeling. Despite extensive characterization of this regulatory network, the origin of the extreme stability of the lysogenic state remains unclear. We have constructed a stochastic model for this switch. Using Forward Flux Sampling simulations, we show that this model predicts an extremely low rate of spontaneous prophage induction in a recA mutant, in agreement with experimental observations. In our model, the DNA loop formed by octamerization of CI bound to the O L and O R operator regions is crucial for stability, allowing the lysogenic state to remain stable even when a large fraction of the total CI is depleted by nonspeciﬁc binding to genomic DNA. DNA looping also ensures thattheswitchisrobusttomutationsintheorderofthe O R binding sites. Our results suggest that DNA looping can provide a mechanism to maintain a stable lysogenic state in the face of a range of challenges including noisy gene expression, nonspeciﬁc DNA binding, and operator site mutations.

. Schematic drawing of tripartite RND multidrug efflux system AcrAB-TolC of the Gram-negative bacterium Escherichia coli (courtesy of Andrea Eberle, ETH Zü rich, Switzerland). Suggestions on the stoichiometry of the adaptor AcrA to inner membrane RND component AcrB (or to outer membrane channel TolC) vary between 1 and 4. Edited by Peter H. von Hippel, Institute of Molecular Biology, Eugene, OR, and approved March 16, 2009 (received for review October 16, 2008) The bistable gene regulatory switch controlling the transition from lysogeny to lysis in bacteriophage presents a unique challenge to quantitative modeling. Despite extensive characterization of this regulatory network, the origin of the extreme stability of the lysogenic state remains unclear. We have constructed a stochastic model for this switch. Using Forward Flux Sampling simulations, we show that this model predicts an extremely low rate of spontaneous prophage induction in a recA mutant, in agreement with experimental observations. In our model, the DNA loop formed by octamerization of CI bound to the O L and OR operator regions is crucial for stability, allowing the lysogenic state to remain stable even when a large fraction of the total CI is depleted by nonspecific binding to genomic DNA. DNA looping also ensures that the switch is robust to mutations in the order of the O R binding sites. Our results suggest that DNA looping can provide a mechanism to maintain a stable lysogenic state in the face of a range of challenges including noisy gene expression, nonspecific DNA binding, and operator site mutations.
biochemical networks ͉ genetic switches ͉ rare events ͉ stochastic modeling ͉ nonspecific DNA binding T he bistable developmental switch controlling the transition from lysogeny to lysis in bacteriophage is one of the best characterized gene regulatory networks (1). In the lysogenic state, the phage genome is integrated into the chromosome of the Escherichia coli host cell and is essentially dormant because of expression of the cI repressor gene, the product of which represses cro and other genes (Fig. 1). A transition to the lytic switch state can occur in response to DNA damage (via UV irradiation), when CI molecules are cleaved by RecA. Transcription of the cro gene from P R then triggers a cascade of gene activation, leading to phage excision, replication, and cell lysis. In mutants where this cascade is blocked, a state with elevated Cro levels is stable for several cell generations. This is known as the anti-immune state (2). A simple and intuitive explanation has been presented for this bistability (1): in the lysogenic state, CI is dominant and cro is repressed; yet, once cro begins to be expressed, Cro dimers repress transcription of cI, making the transition to lysis inevitable. However, quantitative measurements have revealed a puzzle: the lysogenic state of the phage switch is both extremely stable (3)(4)(5) and robust to rewiring of its transcriptional regulatory interactions (3). These features have not yet been explained by mathematical models. Here, we present dynamical simulations of a stochastic mathematical model that reproduces this seemingly mysterious behavior. Our simulations provide evidence that DNA looping plays a key role in ensuring the stability and robustness of the switch. The molecular interactions controlling the transcription of cI and cro have been studied in great detail. CI and Cro bind as dimers to the operator sites O R 1, O R 2, and O R 3 ( Fig. 1), which control expression of the cI and cro genes from the P RM and P R promoters. Transcription from the unactivated P RM promoter is Ϸ10 times weaker than from P R ; however, when a CI dimer is bound at O R 2, transcription from P RM is enhanced to approximately the same level as from P R (i.e., the 2 promoters can compete with one another only when CI is bound at O R 2). CI dimers bind preferentially and cooperatively to the O R 1 and O R 2 sites, which overlap the P R promoter. When CI is bound to these 2 sites, cro (P R ) is repressed and cI (P RM ) is activated. Cro dimers bind preferentially to O R 3, which overlaps P RM , so that when this site is occupied, cI is repressed (1). An important additional component of the network architecture is a DNA loop that can form between the O R site and the left operator site O L , located 2,400 bp from O R , which also has 3 adjacent binding sites for CI dimers (O L 1, O L 2, and O L 3). This loop is mediated by octamerization between pairs of CI dimers bound at O R and O L (6). The role of this DNA looping interaction in the function of the phage switch remains a subject of debate (7)(8)(9)(10)(11).
Quantitative measurements have revealed several intriguing features of the phage switch. First, the lysogenic state is extremely stable (3)(4)(5), despite the stochastic nature of the underlying gene regulatory network. Stochastic fluctuations in gene regulation (''noise'') might be expected to cause spontaneous transitions from the lysogenic to lytic states, even in lysogens lacking RecA; yet the rate of these transitions is so low that it is almost unmeasurable. Second, recent measurements of P RM and P R promoter activity suggest that only a small fraction of the total CI in the cell is available for binding to O R (7,8,(12)(13)(14), whereas other measurements show that the total concentration of CI in the lysogen varies dramatically from cell to cell (15). Taken together, these results suggest that the stability of the lysogen is rather insensitive to the number of free intracellular CI molecules. Despite this stability, This article is a PNAS Direct Submission. 1 Present address: Division of Ecology and Evolutionary Biology, University of Glasgow, Glasgow G12 8QQ, United Kingdom. 2 To whom correspondence may be addressed. E-mail: tenwolde@amolf.nl or rallen2@ ph.ed.ac.uk.
This article contains supporting information online at www.pnas.org/cgi/content/full/ 0810399106/DCSupplemental. transition to lysis occurs readily in wild-type phage on UV irradiation, which leads to cleavage of CI by RecA. Finally, and remarkably, switch function is robust to changes in the gene network architecture itself: when the order of the 3 O R binding sites is altered so that O R 1 is replaced by O R 3 or vice versa, the network remains functional (3). Computer simulations should be an excellent tool for explaining this behavior. Although stochastic simulations have successfully been used to model the initial developmental choice between lysogeny and lysis for lambda-infected cells (16), modeling spontaneous switching of an already established lysogen has proved problematic. Despite the fact that a wealth of biochemical parameters is available for this network, no model has yet reproduced its extraordinary stability and functional robustness (4,12,(17)(18)(19)(20)(21). In this article, we present a model that takes into account the stochastic character of the chemical reactions and includes DNA looping and depletion of free CI and Cro by nonspecific binding to genomic DNA. Our model also explicitly describes the detailed dynamics of the binding of transcription factors to the promoters. Accurate computation of spontaneous switching rates for this large reaction set is achieved using the Forward Flux Sampling (FFS) rare event simulation method (22,23), in combination with temporal coarse graining of dimerization and nonspecific DNA binding reactions. Our simulations show that this stochastic model can reproduce the bistability of the switch, its robustness to operator site mutations, and the extreme stability of the lysogenic state, even in the presence of nonspecific DNA binding.
In this work, we study the effect on switch function of 2 key parameters: the strengths of the DNA looping and nonspecific binding interactions. We find that the DNA looping interaction plays a crucial role. In the absence of the looping interaction, a highly stable lysogenic state can be achieved, but this state is very sensitive to depletion of free CI and to operator site mutations. When looping is included in the model, the lysogen is insensitive to CI depletion and robust to rearrangement of the operator sites. We conclude that DNA looping may play an important role in allowing the phage switch to function reliably even under highly destabilizing conditions in the host cell.

The Model
Our model consists of a set of chemical reactions, simulated using the Gillespie algorithm (24). The components of the model are: dimerization of CI and Cro proteins, binding of CI and Cro dimers to specific DNA binding sites and O L 3, binding of RNA polymerase (RNAp) to promoters P RM , P R , and P L , transcription of cI and cro, translation of the corresponding mRNA transcripts, degradation of mRNA transcripts, and removal of CI and Cro monomers and dimers from the cell. Our model also includes formation of a DNA loop between O R and O L , mediated by a CI octamer, and nonspecific binding of CI and Cro dimers to genomic DNA. The key parameters that we vary are the strength of the nonspecific DNA binding interaction ⌬G NSB and the strength of the DNA looping interaction ⌬G loop . Other parameters are fixed using biochemical data as far as possible. The model parameters are discussed in brief here and are described in full in the supporting information (SI) Text.
Host Cell Parameters. We assume that the E. coli host cell is growing rapidly [doubling time, 34 min (3)], and has 3 copies of each of the O R and O L operators (4), in a cell volume of 2 m 3 (4). The concentration of free RNAp in the cytoplasm is taken to be 50 nM (25), but our conclusions are not sensitive to this parameter, as we demonstrate in the SI.
Operator Binding Dynamics. Equilibrium constants from the literature were used for CI (26,27) and Cro (28) for RNAp binding to P RM and P R (29) and for CI binding to the O L sites (8). Cro is assumed to bind to O L and O R sites identically. The total number of possible (unlooped) configurations of the O R and O L operators are, respectively, 40 and 36. Because we are performing dynamical simulations, we require rate constants k a and k d for association and dissociation. For all association rates, we used the diffusion-limited value k a ϭ 4D ϭ 0.314 m 3 s Ϫ1 (taking the diffusion constant D ϭ 5 m 2 s Ϫ1 and the molecular size ϭ 5 nm). The rate constant for dissociation, in s Ϫ1 , was then deduced from the equilibrium constant, using k a /k d ϭ (exp[Ϫ⌬G/RT])/ (6.023 ϫ 10 8 ) m 3 , where k a is in m 3 s Ϫ1 , ⌬G is in kcal/mol, RT ϭ 0.616 kcal/mol at 37°C, and 6.023 ϫ 10 8 m 3 is a volume conversion factor.
Protein and mRNA Production and Removal. We model transcription as a single reaction in which an mRNA molecule is produced when RNAp is bound to a promoter. P RM activity is enhanced when a CI dimer is bound at O R 2. Transcription rates are 0.014 s Ϫ1 for P R , 0.001 s Ϫ1 for unstimulated P RM , and 0.011 s Ϫ1 for stimulated P RM (29)(30)(31)(32). All mRNA transcripts are degraded with a half-life of 2 min. Translation and protein folding are combined into a single step. The model generates a statistical distribution for the number of proteins produced per transcript, which is governed by the balance between the translation and mRNA degradation rates. The average of this distribution (the ''burst size'') is 6 and 20 for CI and Cro, respectively. The burst size for Cro follows ref. 4. The value for CI is based on the observation that the CI ribosome binding site (RBS) is Ϸ7-fold weaker than that of LacZ (I. Dodd, personal communication), and that the LacZ burst size is 30-40 (33, 34). A somewhat weaker CI RBS was observed by Shean and Gottesman (35). Protein monomers and dimers are removed with rate constant ln(2)/T c , where the cell cycle time T c ϭ 34 min (3). We also include active degradation of Cro monomers with half-life of 42 min (4, 13).
Dimerization. Dimerization free energies are taken to be Ϫ11 kcal/mol for CI (36,37) and Ϫ8.7 kcal/mol for Cro (38). The association reaction is again assumed to be diffusion-limited. To increase the efficiency of our simulations, we coarse-grain the monomer-monomer association and dissociation reactions for both CI and Cro (39), as described in the SI. DNA Looping. When an O R and an O L operator each carry at least 2 adjacent CI dimers (at the 1-2 or 2-3 sites), these operators can associate to form a ''looped state,'' with rate k loop , which dissociates with rate k unloop . The total number of possible looped states is 49. Binding of CI dimers to the 2 nonoctamerized sites in a loop occurs with a cooperativity factor exp(Ϫ⌬G tet /RT), where ⌬G tet ϭ Ϫ3 kcal/mol (8). Because we assume fast growth of the host cell, our model contains 3 copies of the host genome and 3 copies of the phage switch. Because the loop is much longer than the persistence length of DNA, we assume that any O R -O L combination can form a loop. The strength and dynamics of the DNA loop in vivo are unknown. We therefore test the effects of DNA looping on the network behavior, by varying the ratio k loop /k unloop ϵ exp(Ϫ⌬G loop / RT). We generally assume a fixed value 62.1 s Ϫ1 for k loop (arising from considerations of polymer dynamics as discussed in the SI), but we find that only the ratio is important.
Nonspecific DNA Binding Dynamics. We model nonspecific DNA binding (8,14) by including in our reaction set association and dissociation of CI and Cro dimers to 10 7 genomic DNA sites, corresponding to 2-3 copies of the bacterial genome. The association rate k a is assumed to be diffusion-limited, and we assume identical nonspecific binding affinities for CI and Cro. Nonspecifically bound dimers are removed from the cell with rate constant ln(2)/T c . We do not model nonspecific binding of RNAp, because our value of 50 nM corresponds to the free RNAp concentration (25). To investigate the effects of nonspecific DNA binding on the model switch, we vary the parameter ⌬G NSB where k a /k d ϭ (exp[Ϫ⌬G NSB /RT])/(6.023 ϫ 10 8 ) m 3 . We assume that these reactions are fast compared with the other reactions in the network, and can be coarse-grained (39) as described in the SI.

Bistability
Our model represents a mutant phage in which the lytic pathway is nonfunctional, because we do not model the lytic genes downstream of cro (2). Such mutants show 2 stable states: the very stable lysogenic state, with high CI and low Cro levels, and the anti-immune state, with elevated Cro and little CI, which is less stable, but is nevertheless maintained for several generations (2). For our model to reproduce this bistability, simulations initiated in either of the lysogenic or anti-immune states should remain stable, only rarely making a spontaneous transition to the other state. Fig. 2 shows that this is indeed the case: in Fig. 2 A, a simulation initiated in the lysogenic state remains in that state, whereas in Fig. 2B, a simulation run with the same parameters but initiated with a high Cro concentration and little CI remains in the anti-immune state. Fig. 3 shows the range of parameter values for which our model shows bistability. Our simulations give steady-state concentrations of CI and Cro in good agreement with measured values for the lysogenic and anti-immune states. For the lysogenic state, we obtain Ϸ200-400 CI per cell, compared with a measured value of 220 (40). For the antiimmune state, Cro per cell ranges from Ϸ250-900 (depending on ⌬G NSB ), corresponding to concentrations of 200-750 nM, compared with a measured concentration of Ϸ400 nM (12). These values are presented as functions of ⌬G loop and ⌬G NSB in the SI (Fig. S1). In our model, the DNA looping interaction decreases the lysogenic CI concentration by as much as a factor of 2, in agreement with the observations of Dodd et al. (7,8). We have also simulated mutants without the O R 3 or the O L 3 binding sites, corresponding approximately to the OR3-r1 and OL3-4 mutants of refs. 7, 8, and 11. For these mutants, we find lysogenic CI concentrations 1.9 and 1.8 times the wild-type values (for ⌬G loop ϭ Ϫ3.7 kcal/mol and ⌬G NSB ϭ Ϫ2.8 kcal/mol), in reasonable agreement with the results of ref. 8 (factors of 2.8 and 3 respectively), even though our model does not include the up-regulation of P RM by the loop identified in ref. 11, or the effect on P RM of the OR3-r1 mutation. For the same parameters, the P R promoter is repressed by a factor of 1.6 in the antiimmune state, in good agreement with ref. 9.

Extreme Stability of the Lysogenic State
The spontaneous switching rate from lysogeny to lysis in recAϪ mutants is so low as to be almost undetectable, being Ͻ2 ϫ 10 Ϫ9 per cell per generation (3,4). Reproducing such low spontaneous switching rates while maintaining bistability provides an extreme challenge for a computational model (4). We quantified the stability of the lysogenic and anti-immune states for our model by using the FFS rare event simulation method (22,23). This method allows us to compute the rate of spontaneous switch flips, even though such flips would hardly ever be observed in a typical ''brute-force'' simulation run. FFS uses a series of interfaces (defined by an order parameter) between the initial (CI-rich) and final (Cro-rich) states to split a switching event into a number of more likely transitions between successive interfaces. Here, our order parameter was the difference between the total number of Cro and CI molecules. This method has been used for a simple model of a mutually repressing genetic switch (22,23). Details of the FFS method and its implementation are given in the SI. Our results are listed in Table 1. In the absence of DNA looping (Table 1, top 2 rows), the lysogen is stable for Ϸ10 9 generations only when nonspecific DNA binding is absent. When nonspecific binding is included in the model, the lysogen becomes much less stable: a relatively modest nonspecific binding strength ⌬G NSB ϭ Ϫ2.1 kcal/mol produces a lysogen that is only stable for Ϸ1,700 generations: approximately a million times less stable than the experimental observation. In contrast, when DNA looping is included in the model (Table 1, bottom 3 rows), extremely stable lysogens can be achieved for a wide range of DNA looping and nonspecific DNA binding strengths. These lysogenic states are in fact even more stable than those observed experimentally (3). One possible explanation might be the effects of passing DNA replication forks (not included in the model), which might be expected to destabilize looped DNA configurations and/or remove bound CI from the operator sites. The anti-immune state for our model is much less stable than the lysogenic state, in agreement with experimental observations (2,9).

DNA Looping Allows the Switch to Function Despite CI Depletion
It is believed that a large fraction of the transcription factors in a bacterial cell are unavailable for binding to their specific binding sites because they are nonspecifically bound to genomic DNA (41). Recent expression measurements for the P RM promoter have suggested that this is the case for CI (7,8,14) (see also Fig. S2). This nonspecific binding poses a severe challenge to the phage switch. Although both CI and Cro are depleted B A   3. Bistability of the model switch as a function of the nonspecific binding strength ⌬G NSB and the DNA looping strength ⌬Gloop. Blue squares represent parameter sets for which only the lysogenic state is stable (no stable anti-immune state); for red diamonds only the anti-immune state is stable (no stable lysogen), whereas green circles represent parameter combinations where the model switch is bistable. To determine bistability, independent simulations are run, starting in the lysogenic and anti-immune states. If the system remains in the lysogenic state for at least 1,000 generations, and in the anti-immune state for at least 15 generations, it is considered to be bistable. by nonspecific binding, the P RM promoter is intrinsically weak and requires activation by CI at O R 2 to compete effectively with P R . One would therefore expect nonspecific DNA binding by CI to drastically destabilize the lysogen, and to compromise the bistability of the switch. Table S1 shows how the concentration of free CI depends on the nonspecific DNA binding strength. To characterize the effects on switch function, we investigated the range of nonspecific binding strengths over which our model gave bistability for different values of the DNA looping parameter ⌬G loop . Our results are shown in Fig. 3. In the absence of DNA looping (Fig. 3, top row), bistability is indeed strongly compromised as the parameter ⌬G NSB is increased in magnitude (left to right). For ͉⌬G NSB ͉ Ͼϳ3 kcal/mol, the model is no longer bistable; the lysogenic state cannot be sustained, and only the anti-immune state is stable. However, when the DNA looping interaction is included in the model, bistability is maintained over a much wider range of ⌬G NSB values. In Fig. 3, the width of the green (bistable) region increases dramatically as the parameter ͉⌬G loop ͉ increases. Our model therefore suggests that one role of the DNA looping interaction may be to ensure that the switch continues to function even when the level of intracellular free CI is depleted by nonspecific DNA binding (7,8) or by cell-to-cell fluctuations (15). We note that this stability to CI fluctuations is not expected to prevent lysis from occurring on UV irradiation of the wild-type phage: even for a strong looping interaction, rapid degradation of CI by RecA will eventually lead to so little CI being present that the loop cannot be maintained, on which lysis will occur. To check this, we simulated a version of the model in which we fixed the total CI concentration, for a typical parameter set (⌬G loop ϭ Ϫ3.7 kcal/mol, ⌬G NSB ϭ Ϫ2.8 kcal/mol). When we artificially lower the CI concentration to Ϸ10% of the lysogenic steady-state level (30 CI molecules per cell), the system flips to the anti-immune state. This result is in good agreement with the observation of Bailone et al. (42) that prophage induction occurs at a CI concentration Ϸ10% of that in the lysogen.

DNA Looping Causes Robustness to Operator Mutations
In an important series of experiments, Little et al. (3) showed that the basic functions of the phage regulatory network are robust to changes in network architecture. Mutants O R (121) and O R (323), in which the order of the O R 1, O R 2, and O R 3 binding sites was altered compared with the wild-type O R (321), formed stable lysogens (although less stable than the wild type) that could be induced to enter the lytic pathway on UV irradiation (3). To our knowledge, this robustness has not been reproduced in computer models (4).
We tested whether our model was able to produce stable lysogens for the O R (121) and O R (323) mutants, which were created in our simulations by changing the operator binding site affinities. We neglect possible changes in the properties of P RM (43), whose DNA sequence is also affected by the Little et al. (3) substitutions. The range of parameters for which our model mutants are bistable is shown in Fig. 4. We find that the DNA looping interaction plays a key role in ensuring robustness. In the absence of DNA looping (Fig. 4, top row), our model produces a stable lysogen for O R (121) only if the nonspecific binding strength ⌬G NSB is below a critical value, and cannot sustain a stable lysogen for O R (323) at all. However, when DNA looping is included in the model (Fig. 4, lower rows), the O R (323) mutant can achieve a stable lysogenic state, and the O R (121) lysogen becomes able to tolerate stronger nonspecific DNA binding interactions. In fact, if the looping strength is too strong, our model predicts that the O R (121) mutant may lose the ability to sustain a stable anti-immune state. The steady-state concentrations of CI in the lysogen predicted for these mutants are in reasonable agreement with the observations of Little et al. (3). For a typical parameter set (⌬G loop ϭ Ϫ3.7 kcal/mol, ⌬G NSB ϭ Ϫ2.8 kcal/mol), we obtained steady-state lysogenic CI concentrations 38% and 81% of the wild-type values for O R (121) and O R (323), respectively, compared with 25-30% and 60-75% measured by  Table  S2). For all combinations of ⌬G loop and ⌬G NSB , the O R (323) lysogen is less stable than O R (121), which in turn is less stable than the wild type, although quantitatively the magnitude of this destabilization is greater than the factors of 10 observed by Little et al. (3). Our model allows us to make the tentative prediction that (in a recAϪ host) the order of stability of the anti-immune states will be O R (323) Ͼ wild-type Ͼ O R (121). These stabilities have not yet been measured (J. Little, personal communication).

Alternatives to Looping
Our results indicate that the DNA looping interaction allows the phage switch to maintain an extremely stable lysogenic state, which is robust to operator mutations, while retaining its essen- A cell generation time of 34 min is assumed. FFS calculations used 8 -85 interfaces, 1,000 -10,000 configurations at the first interface, and 500 -10,000 trials per interface, and were averaged over 10 -100 runs. tial bistability. Could this result have been achieved by evolution in any other way? Recent work by Babić and Little (44) has shown that a triple mutant lacking CI cooperativity but with increased P RM activity and enhanced binding of CI to O R 2, can maintain a stable lysogen, with increased CI levels compared with the wild-type, and can switch into the lytic state. It is very unlikely that this mutant can form a DNA loop (J. Little, personal communication). This suggests that increased lysogenic CI levels could substitute for cooperative interactions (including DNA looping). To test whether our model could reproduce these experiments, we modified our simulation parameters to represent the 2 triple mutants cI Y210N prm252 O R 2up and cI Y210N prmup-1 O R 2up. In both cases, we removed all cooperative interactions between CI dimers, including looping, and increased the affinity of CI for O R 2 by a factor of 5. For the prm252 mutation, the basal and stimulated P RM transcription rates were increased by factors of 6 and 2.5, respectively, and for the prmup-1 mutation these factors were 10 and 5 (44). With nonspecific binding strength ⌬G NSB ϭ Ϫ2.1 kcal/mol, both these ''mutants'' formed stable lysogens, with CI levels Ϸ3 and Ϸ6 times that of the ''wild-type.'' The mutant representing cI Y210N prmup-1 O R 2up also formed a stable lysogen for ⌬G NSB ϭ Ϫ2.8 kcal/mol, but cI Y210N prm252 O R 2up did not. These results are in reasonable agreement with those reported by Babić and Little (44). Furthermore, our results support the view that increased lysogenic CI levels could provide an alternative mechanism for maintaining a stable lysogenic state, although this mechanism has apparently not been selected by evolution.

Discussion
The phage switch represents an important test for computational modeling of gene regulatory networks. This network has been a paradigm in molecular biology for many years: its architecture and biochemical parameters have been extensively studied and its behavior has been thoroughly and quantitatively characterized. Nevertheless, computational models have failed to produce results in agreement with experimental observations. If modeling is to prove a useful tool in the analysis of larger and more complex gene regulatory networks, it is essential that it should produce convincing results for phage . The computer simulation model presented here reproduces both the bistability of the phage switch and its extremely low spontaneous flipping rate. Our model includes only known biochemical interactions and uses, as far as possible, biochemically measured parameters, and is qualitatively insensitive to parameter variation (Table S3 and Figs. S3 and S4). The FFS rare event sampling method allows us to quantify the stability of the lysogenic and anti-immune states, even though spontaneous switching rates would never be observed in a typical brute-force simulation run. An important difference between our model and previous theoretical studies is that we simulate the full transcription factor-DNA binding dynamics. These ''operator state fluctuations'' were eliminated in previous models using the physical-chemical quasi-equilibrium assumption of Shea and Ackers (4,18,29). Our previous work (45), and that of Walczak and Wolynes (46), has shown that operator state fluctuations can drastically affect the switching pathways and the switching rate for bistable genetic networks. Moreover, Vilar and Leibler (47) have shown that operator state fluctuations are coupled to DNA looping dynamics. We therefore model explicitly both the protein-DNA binding dynamics and the DNA looping dynamics. Our model also includes nonspecific DNA binding, as does ref. 4. Although DNA looping and nonspecific binding are modeled in a simplified manner, our model involves many reactions and is computationally expensive to simulate. This problem is alleviated by the FFS method, in combination with the coarsegraining of dimerization and nonspecific binding reactions. We note that a semianalytical approach, based on Large Deviation Theory, which has been applied successfully to a simplified model of the phage switch (10), may provide a promising alternative to direct simulations as performed here.
Our results suggest a key role for the DNA looping interaction in ensuring that the network retains its bistability even when exposed to perturbations. Our model requires DNA looping to achieve bistability in the presence of nonspecific DNA binding and/or operator state mutations. The P RM promoter is intrinsically weak and requires CI to be bound at O R 2 to compete effectively with P R . The DNA loop enhances CI occupancy at O R 2, providing a way to achieve sustained activation of the P RM promoter, even when the free intracellular CI concentration is very small. We also find, in agreement with Dodd et al. (7), that the presence of the loop increases autorepression of P RM by enhancing CI binding at O R 3. Looping thus increases the strength of both autoactivation via O R 2 and autorepression via O R 3; this reduces the fluctuations in CI, which enhances the stability of the switch (48). Our model does not include the loop-mediated increase in maximal P RM activity recently discovered by Anderson and Yang (11); we expect that including this in the model would only strengthen our conclusions.
Our simulations and the experiments of Babić and Little (44) show that a stable lysogen can be obtained in the absence of DNA looping and other cooperative interactions, by raising the level of CI. This observation is perhaps not so surprising, because the stability of genetic switches is predicted to depend exponentially on the expression levels of gene regulatory proteins (48,50). What is then the role of DNA looping? DNA looping not only reduces fluctuations in CI levels, as discussed above, but also reduces operator state fluctuations, which have been predicted to limit the stability of genetic switches (22,45,46). This suggests that the looping interaction allows a stable lysogen to be achieved with lower CI levels. One advantage of this may be a reduction in the energetic cost to the host cell of producing CI (49,51). Alternatively, the dynamical pathways for the transition to lysis may be more favorable for a switch with DNA looping.
Several predictions emerge from our simulations. First, we predict that a mutant phage that cannot form the DNA loop will form a lysogen with much reduced stability compared with the wild type, and may not be able to sustain a lysogen at all. We further predict that depleting free CI (for example, by introducing a large number of specific binding sites) will strongly destabilize the lysogen in a nonlooping mutant but will have a much less severe effect on the wild-type lysogen. Finally, our simulations allow us to tentatively suggest an order for the stability of the anti-immune states in a nonlysing version (for example, O R (121) N Ϫ O Ϫ ) of the Little mutants in a recAϪ host.
This work presents an encouraging picture of the potential for computational modeling to unravel the contribution of different biochemical mechanisms to observed biological behavior. By using a simplified representation of the core part of the phage switch, we are able to obtain behavior in qualitative and partial quantitative agreement with experimental results. Future work should include more sophisticated models for DNA looping and nonspecific binding, explicit modeling of cell growth, DNA replication and cell division, and detailed characterization of the switching pathways. Such work should make stochastic modeling, in combination with experiments, an important tool in unraveling the mechanisms behind the complex behavior of biochemical networks.