Biophysical Journal 86:75-84 (2004)
© 2004 The Biophysical Society
Why the Lysogenic State of Phage 
Is So Stable: A Mathematical Modeling Approach
Moisés Santillán * and
Michael C. Mackey 
* Centre for Nonlinear Dynamics, McGill University, H3G 1Y6 Montreal, Quebec, Canada; and Departments of Physiology, Physics & Mathematics and Centre for Nonlinear Dynamics, McGill University, H3G 1Y6 Montreal, Quebec, Canada
Correspondence: Address reprint requests to Moisés Santillán, Depto. de Física, Esc. Sup. de Física y Matemáticas, Inst. Politécnico Nal. 07738 México D.F., México. Tel.: +52-55-57296000 ext. 55319; Fax: +52-55-57296000 ext. 55051; E-mail: moyo{at}esfm.ipn.mx.
 |
ABSTRACT
|
|---|
We develop a mathematical model of the phage 
lysis/lysogeny switch, taking into account recent experimental evidence demonstrating enhanced cooperativity between the left and right operator regions. Model parameters are estimated from available experimental data. The model is shown to have a single stable steady state for these estimated parameter values, and this steady state corresponds to the lysogenic state. When the CI degradation rate (
cI) is slightly increased from its normal value (cI 
0.0 min-1), two additional steady states appear (through a saddle-node bifurcation) in addition to the lysogenic state. One of these new steady states is stable and corresponds to the lytic state. The other steady state is an (unstable) saddle node. The coexistence these two globally stable steady states (the lytic and lysogenic states) is maintained with further increases of cI until cI 0.35 min-1, when the lysogenic steady state and the saddle node collide and vanish (through a reverse saddle node bifurcation) leaving only the lytic state surviving. These results allow us to understand the high degree of stability of the lysogenic state because, normally, it is the only steady state. Further implications of these results for the stability of the phage switch are discussed, as well as possible experimental tests of the model.
|
INTRODUCTION
|
|---|
Bacteriophage (or simply phage) is a virus capable of infecting Escherichia coli bacteria. After infection, the virus can follow either one of two different pathways: 1), the virus integrates its DNA into the host bacterial DNA and duplicates when the bacterium divides (this pathway is known as lysogeny); 2), the virus uses the bacterial molecular machinery to make many viral copies and leave (after killing the host bacterium) to infect other bacteria (the so-called lysis pathway). Once the virus is in the lysogenic state, it can shift to the lysis state under certain conditions, e.g., if the bacterial culture is irradiated with ultraviolet (UV) light. The molecular regulatory mechanism responsible for the lysogeny/lysis decision is known as the phage switch.
The two patterns of phage behavior, lytic and lysogenic, and the subtle ways in which it subverts its host, E. coli, have made it a paradigm for many biological pathways (Gottesman, 1999
). One of the most striking characteristics of the phage switch is the fact that the intrinsic loss rate of lysogeny is of the order of 10-7 per cell and generation (Aurell et al., 2002; Little et al., 1999; Rozanov et al., 1998). In contrast, the mutation rate in the portion of the genome involved in lysogeny is between 10-6 and 10-7 per generation (Aurell et al., 2002; Little et al., 1999). Thus, the lysogenic state is more stable than the genome itself.
The large volume of experimental data on the behavior of this system and the logical structure underlying the switch performance make it appealing for mathematical modeling. Several mathematical models have been proposed in the last two decades to explain the astonishingly high degree of stability of the phage lysogenic state. Nevertheless, to our knowledge, both the stability of the lysogenic state and the efficiency of the switch still lack a proper explanation. The aim of this paper is to offer an explanation of the lysogenic state stability based on a new modeling effort. Below, we briefly review some of the existing models of the phage switch before turning to an exposition of our results.
The first model incorporating quantitative biochemical information was proposed by Ackers et al. (1982). This is an equilibrium model that explains the existence of the lytic and lysogenic steady states, but not their relative stability. Later, Shea and Ackers (1985) improved their model making it fully dynamical. This model gave the expected qualitative behavior for stable maintenance of lysogeny, as well as for the induction of lysis. Reinitz and Vaisnys (1990) extended the dynamical model, and found a quantitative inconsistency between their experiments and the model predictions. The most extensive model was developed by McAdams and Shapiro (1995) who included all of the proteins involved in the lysis/lysogeny switch, as well as the genes that become active after the fate of the phage is determined. Arkin et al. (1998) published a model of the switch based on a stochastic representation of transcription, translation, and interaction between proteins. They accurately predicted the fraction of lysogens developed after infection. Further, their simulations clearly show how two identical cells in identical conditions, infected with the same number of phage, can still meet different fates. More recently, Aurell and Sneppen (2002) and Aurell et al. (2002) studied the stability of the lysogenic state using a stochastic mathematical model. From their results, they suggested that the current view of the phage switch is incomplete, given the difference they observed between the model predictions and the experimentally observed behavior of a mutant virus strain.
Here, we extend the Shea and Ackers (1985) model to account for recently discovered (Dodd et al., 2001; Ptashne and Gann, 2000) interactions (cooperativity) between regulatory molecules bound at two different operator regions. We pay special attention to the estimation of all of the model parameters from published experimental results. An analysis of the model steady states reveals that under normal conditions it only has one stable steady state, corresponding to lysogeny. When the degradation rate of the protein CI is slightly increased from its normal value (0.0 min-1), two additional steady states appear through a saddle node bifurcation. One of these new steady states is also stable and corresponds to lysis, whereas the other is an (unstable) saddle node. This bistable behavior is maintained with further increases of the CI degradation rate until it reaches 
0.35 min-1, when the saddle node and the lysogenic steady state collide and annihilate each other through a reverse saddle node bifurcation. For even larger degradation CI rates there is only one stable steady state, corresponding to lysis. The consequences of these results on the performance and stability of the phage switch are finally discussed.
|
THE PHAGE SWITCH
|
|---|
An excellent review of the molecular regulatory mechanisms in the phage switch is given by Ptashne (1986). A schematic representation of the switch performance in the lysogenic and lytic steady states is shown in Fig. 1. All of the switch regulatory processes take place in the right operator (OR), which is composed of three regions designated OR1, OR2, and OR3. The promoter PR completely overlaps OR1 and partially overlaps OR2. RNA polymerase enzymes that bind to the promoter PR initiate transcription of gene cro. Similarly, the promoter PRM completely overlaps OR3 and partially overlaps OR2. RNA polymerase enzymes bound to the promoter PRM initiate transcription of the cI gene. Dimers of cI product (denoted CI2) can bind to OR1, OR2, and OR3, in order of increasing affinity. Conversely, dimers of cro product (denoted Cro2) bind to OR3 with the highest affinity, then to OR2, and finally to OR1. When CI2's bind to adjacent OR locations, they do so cooperatively. Thus, the binding energy when there are CI2's bound to OR1 and OR2 or to OR2 and OR3 is smaller (larger in absolute value given that interaction and binding energies are negative) than the sum of the individual binding energies. Recently, Darling et al. (2000b) observed that there is also cooperativity between Cro2's bound to OR1 and OR2, to OR2 and OR3, and to OR1, OR2, and OR3. All of the individual-site binding and interaction energies are given in Darling et al. (2000b).
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 1 A schematic representation of the phage switch in the two stable states: lysogeny (top) and lysis (bottom).
|
|
In the lysogenic state (see Fig. 1), gene cI is on while gene cro is off. Monomers of cI product (denoted by CI) spontaneously combine to form dimers CI2. Similarly, monomers of cro product (denoted by Cro) spontaneously form dimers Cro2. Due to cooperativity, most of the OR1 and OR2 sites are occupied by CI2's in the normal lysogenic state. This has two effects: 1), promoter PR is repressed (because it is blocked by a CI2); and 2), the initiation of transcription at promoter PRM is enhanced. One CI2 bound to OR2 does not affect the probability that a RNA polymerase will bind promoter PRM and form a closed complex, but it increases the probability for the closed complex to isomerize into an open complex to start transcription. In other words, dimers CI2 repress the production of Cro and enhance the production of CI. Nevertheless, if the concentration of CI2 reaches very high values, the probability for CI2 to bind OR3 will be increased, which has the effect of repressing RNA polymerase binding to PRM. Thus CI2 regulates its own concentration by enhancing CI production if its concentration is not too high, and otherwise repressing transcription of gene cI.
If the CI2 concentration decreases, for instance by the cleavage of CI by RecA proteins (activated by UV light), the probability for OR1 and OR2 to be free from CI2 is increased. This, on its own, creates the possibility that a polymerase will bind PR and start transcription of gene cro and, in the long run, leads to an increasing Cro2 concentration. At a high enough Cro2 concentration, a Cro2 can bind OR3 and repress CI production, establishing the lytic state. In this state, gene cro is on while gene cI is off. When the concentration of Cro2 is too high, a Cro2 can bind to OR2 and even to OR1, repressing the production of Cro.
New experimental evidence (Dodd et al., 2001; Ptashne and Gann, 2000) suggests that there is cooperativity between CI2's bound to the OR and OL operators. As depicted in Fig. 2, the phage DNA folds in such a way that CI2's bound to OR and OL sites can interact due to their proximity. Ptashne and Gann (2000) speculate that this newly discovered cooperativity may help explain the incredibly high stability of the lysogenic state.
|
MODEL DEVELOPMENT
|
|---|
The complex formed by the operator OR and the promoters PR and PRM has 40 different binding states. Twenty-seven of these states arise when we consider that OR1, OR2, and OR3 can be either empty, bound by a CI2, or bound by a Cro2, giving 33 = 27 different possible combinations without considering the binding of RNA polymerase molecules to promoters PR and PRM. The fact that whenever OR3 is free (no matter what the 32 = 9 states of OR1 and OR2 are), promoter PRM can be bound by a RNA polymerase, accounts for nine more binding states. Similarly, three more states appear because whenever OR1 and OR2 are free (regardless of the state of OR3 (three combinations)), promoter PR can be bound by a RNA polymerase. Finally, the last binding state (to complete all 40 of them) is that in which all three OR1, OR2, and OR3 are free and both PR and PRM are occupied by RNA polymerase molecules.
The left operator of phage (OL) consists of three regions OL1, OL2, and OL3, that can be bound by CI2 and Cro2 molecules. There is only one promoter regulated by this operator, promoter PL. This promoter completely overlaps OL1 and partially overlaps OL2. Thus, the left operator-PL complex has 30 different binding states. Twenty-seven of them arise from the fact that each one of OL1, OL2, and OL3 has three binding states. They can be either free, bound by a CI2, or bound by a Cro2. The other three states arise because whenever OL2 and OL3 are free, a RNA polymerase molecule can bind PL, independent of the state of OL1.
The energy (Ei) of every one of the 40 x 30 = 1200 binding states of the complex formed by the right and left operators plus the promoters PR, PRM, and PL can be calculated given the cooperativity and individual binding energies of CI2, Cro2, and mRNAP. The equation to calculate these energies is:
 | (1) |
where
The terms 
(X = RM, R, L) represent the binding energy of a RNAP molecule bound to site PX, the terms 
(X = R, L) represent the binding energy of molecule Y to site 
the terms 
represent the interaction energy between two Y molecules bound to sites 
and 
and the terms 
represent the interaction energy between three Cro2 molecules bound to OX1, OX2, and OX3. All of these individual binding and interaction energies have been measured elsewhere and are given in Appendix A. To our knowledge, the interaction energy between CI2's bound to operators OR and OL has not yet been measured. The present model assumes that whenever there is one CI2 bound to OR
and a second CI2 is bound to OL, they interact with an energy 
GRL whose value is also estimated in Appendix A.
Following the same technique used by Ackers et al. (l982) and Shea and Ackers (1985), under the assumption of thermodynamic equilibrium, the probability of every one of the 1200 binding states can be calculated from
 | (2) |
In Eq. 2, Pi and Ei are the probability and the energy of the i-th binding state, respectively. R is the ideal gas constant and T is the absolute temperature. Here we take T = 37°C so RT 0.617 kcal/mol (Ackers et al., 1982; Shea and Ackers, 1985). The partition function Z is given by
 | (3) |
Finally, [CI2], [Cro2], and [RNAP] are, respectively, the concentration of CI2, Cro2, and RNA polymerase molecules, whereas 
i, ßi, and i are, respectively, the number of CI2, Cro2, and RNA polymerase molecules bound to the operator-promoter complex in the i-th state.
Note that the thermodynamic equilibrium assumption does not mean that the probabilities remain constant in time. Rather, it means that the probabilities considered are those that would be attained by the real system if it would have enough time to reach the state of thermodynamic equilibrium, given the concentrations [CI2], [Cro2], and [RNAP]. If these concentrations are time dependent, the probabilities given under the thermodynamic equilibrium assumption will change accordingly. If the chemical processes involving association and dissociation of CI2 and Cro2, as well as RNA polymerase molecules to the right operator-promoter complex, take place rapidly (relative to the temporal changes of the concentrations [CI2], [Cro2], and [RNAP]), the thermodynamic equilibrium assumption is appropriate.
The probability fR for PR to be bound by a polymerase can be calculated by adding the probability of all compatible states, i.e., all the states in which a polymerase is bound to PR while OR1 and OR2 are free, independent of the states of OR3, PRM, OL1, OL2, OL3, and PL. Similarly, the probability 
for a polymerase to be bound to PRM with a CI2 bound to OR2 can be calculated by adding the probability of all the binding states in which OR3 is free, a polymerase is bound to PRM, and a CI2 is bound to OR2, no matter what the states of OR1, PR, OL1, OL2, OL3, and PL are. Finally, the probability 
for a polymerase to be bound to PRM without a CI2 bound to OR2 is calculated by adding the probability of all the states in which OR3 is free, a polymerase is bound to PRM, and OR2 is either free or bound by a Cro2 independent of the states of OR1, PR, OL1, OL2, OL3, and PL.
The rate at which RNA polymerase molecules start transcription of the genes cro and cI can be written as kcrofR[OR] and 
respectively (Ackers et al., 1982; Shea and Ackers, 1985). kcro is the transcription initiation rate constant of promoter PR, whereas 
and 
are the transcription initiation rate constants of promoter PRM when OR2 is bound by or free from a CI2 molecule, respectively. OR is the right operator concentration.
Let [McI] and [Mcro] be the concentration of cI and cro mRNA molecules, respectively. Let [CIT] and [CroT], respectively, denote the total concentration of CI and Cro monomers. These concentrations are calculated by adding the free monomer plus twice the dimer concentrations. Thus, from the considerations in the previous paragraphs, we propose the following equations to model the dynamic evolution of [McI], [Mcro], [CIT], and [CroT]:
 | (4) |
 | (5) |
 | (6) |
and
 | (7) |
The meaning of all parameters appearing for the first time in these equations is as follows: M is the rate of mRNA degradation; cro is the rate of CroT degradation; cI is the rate of CIT degradation; 
M is time required after transcription initiation to have a mRNA ready to be bound by a ribosome; cI is the time it takes to translate a monomer CI; cro is the time it takes to translate a monomer Cro; 
cI is the rate of cI translation initiation; and cro is the rate of cro translation initiation. The notation X means X 
X(t - ).
The dimer concentrations ([CI2] and [Cro2]) can be calculated in terms of the total monomer concentrations ([CIT] and [CroT]) if we make a quasi-steady-state assumption for the dimerization reactions (see Appendix B):
 | (8) |
and
 | (9) |
where 
and 
respectively, denote the CI and Cro dimerization dissociation constants.
The values of all of the parameters in Eq. 4–9 are estimated in Appendix A and tabulated in Table 1.
|
NUMBER AND LOCATION OF THE STEADY STATES
|
|---|
The model given by Eq. 4–9 reaches a steady state whenever d[McI]/dt = 0, d[Mcro]/dt = 0, d[CIT]/dt = 0, and d[CroT]/dt = 0. Given that, in the steady state, all of the delayed variables attain their steady-state value, the time delays M, cI, and cro may be ignored. It can be shown after a little algebra that the CIT and CroT steady-state values satisfy the following equations:
 | (10) |
and
 | (11) |
where
 | (12) |
and
 | (13) |
The equation 
([CIT],[CroT],cI) = 0 determines a family of curves in the ([CIT],[CroT]) space (one for every value of cI), whereas the equation 
([CIT],[CroT]) = 0 determines one single curve in the same space. For a given value of cI, the points where the curves = 0 and = 0 intersect correspond to the steady states of the system.
The curve = 0 is shown in Fig. 3 in the [CroT] vs. [CIT] phase space, along with = 0 curves for various values of cI. For the normal lysogenic value cI 0.0 min-1 (see Table 1), both curves intersect at only one point implying that for this value of cI there is only one steady state. This steady state is stable (this point is discussed below in more detail) and corresponds to the lysogenic steady state. Nevertheless, this behavior is at the verge of a bifurcation because a slight increment of cI makes the curve = 0 tangentially intersect the curve = 0. At this point, two more steady states appear via a saddle node bifurcation. One of these steady states is stable and corresponds to the lytic state (see below for more detail), whereas the other is a saddle node. With further increases in cI, the lytic steady state separates from the saddle node, and the saddle node approaches the lysogenic steady state. When cI reaches a value of 0.35 min-1, the saddle node and the lysogenic state collide and annihilate through a reverse saddle node bifurcation.
|
NUMERICAL TEST OF THE STEADY STATES' STABILITY PROPERTIES
|
|---|
Consider the system of time-delay differential equations given by Eq. 4–7. In Table 1, it is clear that the time delay M is one order of magnitude smaller than cI and two orders of magnitude smaller than cro. This observation allows us to slightly simplify the model by ignoring the time delay M wherever it appears. After this simplification, the system of equations was numerically solved by means of the fourth-order Runge-Kutta algorithm (adapted to account for the time delays), implemented in FORTRAN.
Projections of the phase space trajectories onto the ([CIT],[CroT]) space are shown in Fig. 4 for various values of cI. The curves in this figure exemplify a large amount of numerical experiments that we performed to test the steady states' stability properties. These properties can be described as follows: for cI = 0.0 min-1 the system has one single globally stable steady state corresponding to lysogeny (Fig. 4 A). A slight increment of cI beyond this value produces two more steady states via a saddle node bifurcation. One of these is stable and corresponds to lysis whereas the other is an unstable saddle node (Fig. 4 B). This bistable behavior is maintained with further increments of cI until about cI 0.35 min-1, when the saddle node and the lysogenic steady state collide annihilating each other by means of a reverse saddle node bifurcation. For cI's larger than 0.35 min-1 the system has one single globally stable lytic steady state (Fig. 4 C).
We numerically solved the model equations by changing the value of parameters other than cI to test the robustness of our results. We found that of all the parameters, the model is particularly sensitive to changes in GRL. Indeed, for absolute values of GRL's smaller than 3.1 kcal/mol the system shows a bistable behavior even for cI = 0.0 min-1. When the absolute value of GRL is increased beyond 3.1 kcal/mol, the system has one single steady state (corresponding to lysogeny) for cI = 0.0 min-1, and the values of cI at which the bifurcations take place increase. This behavior is pictured as a bifurcation diagram in Fig. 5. The above described results indicate that cooperativity between the right and left operator regions plays a very important role in the stability of the lysogeny steady state. Moreover, they also highlight the necessity of detailed experiments to clarify the nature of the interaction between CI dimers bound to the left and right operator regions.
Some of the estimated parameters were obtained by averaging previously reported values. In some cases, the relative errors are higher than 30%. This is in particular the case with 
and 
To test how this variability affects the system dynamic behavior we redraw the plots of Fig. 3 (which illustrate the number and location of the system's steady states for different values of cI) using the extreme values of the PRM transcription initiation rates and 
These plots are presented in Fig. 6. There, notice that with the lowest and values, the system shows bistability even for cI = 0 min-1. Moreover, the cI value at which the lysogenous steady state and the saddle node collide and annihilate decreases. A comparison with the bifurcation diagram of Fig. 5 reveals that the only effect of decreasing the values of and is to slightly displace the bifurcation curves of Fig. 5 to the left along the GRL axis. Similarly, increasing the and values has the single effect of slightly displacing the bifurcation curves to the right. Thus, we conclude that changes of parameters and even larger than 30% do not alter the overall qualitative dynamic behavior of the system.
|
CONCLUSIONS
|
|---|
Our results provide a possible answer to the question of why the lysogenic steady state is so stable. Namely, from our model and analysis the phage switch has only one steady state under normal conditions. Moreover, based on our numerical results the lysogeny steady state appears to be globally stable. This is in agreement with the fact that the lysogenic state is even more stable than the genome itself. Only that fraction of mutations that typically increase the CI degradation rate can cause some lysogens to become lytic. Conversely, if after being irradiated with UV light the activated RecA proteins cleave CI monomers at a rate >0.35 min-1, all of the infected bacteria will go into lysis. This will happen because lysis will be the only available steady state, which is globally stable under this condition. Thus, by taking into account the cooperativity between CI2's bound at the OR and OL operators, it is possible to explain the astonishingly high stability of the phage lysogenic state, as well as the almost perfect efficiency of the switch.
The estimated energy of interaction between CI2's bound to the left and right operators (GRL), relies on an assumed interaction that has not been experimentally demonstrated. Nevertheless, we believe our model is robust in this respect based on the fact that the value we estimate for GRL is similar to the measured interaction energy between CI2's bound to adjacent OR sites (see Appendix A). Furthermore, our assumption can be tested experimentally.
In this work, we predict bistable behavior (between the lysogenic and lytic states) for CI degradation rates (cI) larger than 0.0 but smaller than 0.35 min-1. To test this prediction, mutant phage strains in which CI is not as stable as it is in the wild type would be useful. If for such mutant strains cI is large enough to induce bistable behavior, a rate of lysogen loss of the same order of magnitude as that predicted by Arkin et al. (1998) should be observed, instead of the (orders of magnitude smaller) wild-type loss rate. Other possible tests of the model would involve mutant E. coli strains, in which the activated RecA proteins fail to cleave CI at a rate high enough to induce one single lytic stable steady state. Then, the system would still be in the bistable region, and some of the lysogens would always survive after having been irradiated with UV light.
The dynamic influence of the thermodynamic equilibrium approximation, on which this and other models rely, should be analyzed. In our particular case, this approximation is equivalent to a quasi-steady-state assumption for the chemical reactions by means of which CI2, Cro2, and RNA polymerase molecules bind to, and detach from, the promoter-operator regions. Its validity relies upon the assumption that those chemical reactions are rapid, relative to the transcription and translation processes. In other words, the feasibility of a thermodynamic equilibrium (or a quasi-steady-state) assumption is a question of having phenomena with orders-of-magnitude different characteristic times. In our opinion, although it is not plausible that a phenomenological model based on a thermodynamic equilibrium approximation will generally work for living biological systems, it seems to work well here because the required difference in characteristic times is fulfilled.
The model here presented does not allow us to make quantitative predictions concerning the lysogen loss and survival rates of wild-type and mutant strains, under different experimental conditions. To make them, the fluctuations inherent to molecular systems should be accounted for in a stochastic model.
Finally, we call attention to the fact that the model is highly sensitive to changes in the parameter GRL. This indicates that the interaction between CI dimers bound to the left and right operator regions plays an important role in the stability of the lysogeny state. However, it also highlights the necessity of detailed experiments concerning this interaction to validate our results.
|
APPENDIX A: PARAMETER ESTIMATION
|
|---|
The CI2 and Cro2 binding energies to OR1, OR2, and OR3, the interaction energy between CI2's or Cro2's bound to adjacent sites, and the RNA polymerase binding energy to PR and PRM have been measured elsewhere. The most recent, and presumably most accurate, values are reported by Darling et al. (2000b) and tabulated in Table 2.
The CI2 and Cro2 single-site binding energies at OL are reported by Aurell and Sneppen (2002). However, these energies were measured before the cooperativity between Cro2's bound to adjacent sites was discovered. Cro2 single-site binding energies at OR reported before and after cooperativity was discovered are very different. This is not the case with CI2 (Darling et al., 2000b). Moreover, Cro2 single-site binding energies at OR, that also ignore cooperativity, are very similar to those of OL (Aurell et al., 2002). From this, we take the CI2 single-site binding energies reported by Aurell and Sneppen (2002), but for the CI2 cooperativity and Cro2 single-site binding and cooperativity energies at OL we assume that they are equal to those of OR. Giladi et al. (1990) measured the association constant of the polymerase-PL closed complex formation. They report a value of KB 8.94 x 107 mol-1. From this and taking RT 0.617 kcal/mol, we can estimate the RNA polymerase-PL binding energy to be 
kcal/mol. These energies are tabulated in Table 3.
The value of the probabilities fR, 
and can be calculated for every value of [CI2] and [Cro2] from the energy values tabulated in Tables 2 and 3, as well as from the energy GRL and the RNA polymerase concentration [RNAP] (the last two parameters are estimated below). We also need estimates for the right operator concentration [OR], the growth rate µ, the degradation rates M, cI, and cro, the transcription initiation rates 
and kcro, and for the translation initiation rates cI and cro.
Growth rate
Reinitz and Vaisnys (1990) and Little et al. (1999) performed experiments on the shift of the phage switch from the lysogeny to the lytic state. In these experiments, they employ a bacterial generation time of 34 min. For the purpose of the present work, we will consider the same generation time, which corresponds to a growth rate of µ 2.0 x 10-2 min-1.
E. coli volume
E. coli are rodlike bacteria 3–5 µm long and 0.5 µm in diameter, so they have a volume in the range from 6.0 x 10-16 L to 9.8 x 10-16 L. We considered a mean volume of 8.0 x 10-16 L.
RNA polymerase concentration [RNAP]
According to Bremer and Dennis (1996), there are 1,500 RNA active polymerase molecules per cell in E. coli bacterial cultures growing at the rate µ as estimated above. This leads to a concentration [RNAP] 3.0 µM.
Right operator concentration [OR]
According to Bremer and Dennis (1996), there are 2.5 genome equivalents per average E. coli cell at the growth rate determined by µ. Assuming one operator OR per genome equivalent, the right operator concentration can be estimated as [OR] 5.0 x 10-3 µM.
Average lysogenic CIT concentration
According to Ptashne (1986), there are between 200 and 350 CI molecules per cell (in monomer units) in lysogenic bacteria. Here, we take the mean number 275, which corresponds to a concentration of 
PR transcription initiation rate kcro
Hawley et al. (1985) and Fong et al. (l994) report measurements of the closed-to-open complex isomerization rate at the PR operator kcro. The values they report are respectively: kcro 5.0 x 10-2 s-1 and kcro 4.12 x 10-2 s-1. We take the mean value kcro 4.6 x 10-2 s-1 2.76 min-1.
PRM transcription initiation rates 
and 
Hwang et al. (1988) and Li et al. (1997) report measurements of the closed-to-open complex isomerization rate at the PRM operator. According to Hwang et al. (1988), this rate in the absence of CI2 repressor is around 
and 
in the presence of CI2 repressor. Li et al. (1997) report the following values: 
and 
We consider here the mean values: 
and 
mRNA degradation rate M
Court et al. (1980) performed experiments to measure the half-life of PL transcripts. They report 1/2 6 min, and assert that the half-lives of the PR and PRM transcripts are similar. From this, the mRNA degradation rate can be estimated to be M = ln 2/1/2 0.12 min-1.
CIT degradation rate cI
Following Hawley et al. (l985), Reinitz and Vaisnys (1990), Aurell et al. (2002), and Aurell and Sneppen (2002), we ignore the degradation rate of CIT because it is very small compared with other rates in the model, e.g., the growth rate µ. Therefore, we take cI = 0.
CroT degradation rate cro
Pakula et al. (1986) measured a Cro half-life in vivo of 2600 s 43.3 min. From this, cro 1.6 x 10-2 min-1.
cI translation initiation rate cI
Under normal conditions, the system of time-delay differential equations (Eq. 4–7) model has one single stable steady state corresponding to lysogeny. In this lysogenic state, PR is very tightly repressed whereas PRM is repressed only 
20% (Dodd et al., 2001; Meyer et al., 1980). Because 
the maximum activity of the promoter PRM corresponds to 
and 
and is given by 
Thus, the normalized activity of the promoter PRM in the lysogenic state is 
with which the dynamic equation for [McI] becomes
By solving for the cI translation initiation rate cI from this last equation and the dynamic equation for [Mcro], and taking into account that the concentration of all chemical species is constant in a steady state (d[McI]/dt = 0 and d[CIT]/dt = 0), we obtain
cro translation initiation rate cro
Aurell et al. (2002) estimate that, when translated, a cro transcript produces 51% of the polypeptides an ideal lacZ transcript produces. From this and taking into account that ribosomes load onto lacZ transcripts at 3.2 s = 5.33 x 10-2 min intervals (Kennell and Riezman, 1977), that the lacZ transcript half-life is 2 min (Leive and Kollin, 1967), and that the cro transcript half-life is of the order of 6 min (Court et al., 1980), the cro translation initiation rate can be estimated as cro (0.51 x 2 min)/(6 min x 5.33 x 10-2 min) 3.2 min-1.
The time delays due to cI and cro translation cI and cro
CI and Cro are proteins that are 267 and 66 amino-acid-long, respectively. This means that gene cI is 711 basepairs long, whereas cro gene is 198 basepairs long. From this and taking into account that, according to Bremer and Dennis (1996), the mRNA chain elongation rate is 50 nucleotide/s, the time it takes for genes cI and cro to be translated is cI 0.24 min and cro 6.6 x 10-2 min, respectively.
The time delay due to transcription M
Once a RNA polymerase has transcribed a mRNA chain long enough for a ribosome to bind to it, translation can start. According to Draper (1996), efficient mRNA's can initiate translation every 3 s. From this and the fact that the mRNA chain elongation rate is of the order of 50 nucleotide/s (Bremer and Dennis, 1996), <150 nucleotides are required for a ribosome to bind a mRNA and start translation. Furthermore, the DNA chain elongation rate is at least 490 nucleotide/s (Bremer and Dennis, 1996). Thus it takes <0.31 s after transcription initiation to have a mRNA ready for translation initiation. i. e., M 5.1 x 10-3 min.
The CI and Cro dimerization dissociation constants 
and 
Burz et al. (1994) measured the CI dimerization dissociation constant obtaining 
The Cro dimerization dissociation constant was measured by Darling et al. (2000a) to give 
Interaction energy between CI2 dimers bound to OR and OL, GRL
According to Meyer et al. (1980) and Dodd et al. (2001), PR is very tightly repressed whereas PRM is repressed only 20% in the lysogenic state. From this, the PRM transcription initiation rate 
should be 80% of the maximum transcription initiation rate That is, 
After substitution of the steady state [CIT] 0.55 µM and [CroT] 0.0 µM values this becomes a nonlinear algebraic equation of GRL which, when solved with the aid of Octave's algorithm "fsolve" gives:
|
APPENDIX B: DIMERIZATION KINETICS
|
|---|
Consider the following chemical reaction:
In this equation, a2 represents a dimer of chemical species a, and k+ and k- are the forward and backward reaction rates, respectively.
In a state of chemical equilibrium the following relation is satisfied,
where KD = k-/k+ is the so-called dissociation constant and [x] represents the concentration of species x. Let [aT] be the total monomer concentration:
The last two equations constitute a complete set of equations for the variables [a] and [a2]. By solving for [a2] we obtain the following expression for the dimer concentration in terms of aT and KD,
|
ACKNOWLEDGEMENTS
|
|---|
We wish to thank Olivier Landry and Miryam Elouneg-Jamroz for their collaboration in this project.
This work was supported by COFAA-IPN (México), EDI-IPN (México), CONACyT (México), MITACS (Canada), the Natural Sciences and Engineering Research Council (NSERC grant OGP-0036920, Canada), and Le Fonds pour la Formation de Chercheurs et l'Aide à la Recherche (FCAR grant 98ER1057, Québec).
|
FOOTNOTES
|
|---|
Moisés Santillán is currently an invited researcher at Departamento de Matemáticas, Centro de Investigación y Estudios Avanzados del IPN.
Submitted on May 12, 2003;
accepted for publication September 16, 2003.
|
REFERENCES
|
|---|
Ackers, G. K., A. D. Johnson, and M. A. Sea. 1982. Quantitative model for gene regulation by phage repressor. Proc. Natl. Acad. Sci. USA. 79:1129–1133.[Abstract/Free Full Text]
Arkin, A., J. Ross, and H. H. McAdams. 1998. Stochastic kinetic analysis of developmental pathway bifurcation in phage -infected Escherichia coli cells. Genetics. 149:1633–1648.[Abstract/Free Full Text]
Aurell, E., S. Brown, J. Johnson, and K. Sneppen. 2002. Stability puzzles in phage . Phys. Rev. E65:051914.
Aurell, E., and K. Sneppen. 2002. Epigenetics as a first exit problem. Phys. Rev. Lett. 88:048101.[Medline]
Bremer, H., and P. P. Dennis. 1996. Modulation of chemical composition and other parameters of the cell by growth rate. In Escherichia coli and Salmonella thyphymurium: Cellular and Molecular Biology, Vol. 2. F. C. Neidhart, R. Curtis, J. L. Ingraham, E. C. C. Lin, K. B. Low, B. Magasanik, W. S. Reznikoff, M. Riley, M. Schaechter, and H. E. Umbarger, editors. American Society for Microbiology, Washington, DC. 1553–1569.
Burz, D., D. Becket, N. Benson, and G. K. Ackers. 1994. Self-assembly of bacteriophage CI repressor: effects of single-site mutations on the monomer-dimer equilibrium. Biochemistry. 33:8399–8405.[Medline]
Court, D., B. Crombrugghe, S. Adhya, and M. Gottesman. 1980. Bacteriophage lambda Hin function II. Enhanced stability of lambda messenger RNA. J. Mol. Biol. 138:731–743.[Medline]
Darling, P. J., J. M. Holt, and G. K. Ackers. 2000a. Couple energetics of cro repressor self-assembly and site-specific DNA operator binding. I: Analysis of cro dimerization from nanomolar to micromolar concentrations. Biochemistry. 39:11500–11507.[Medline]
Darling, P. J., J. M. Holt, and G. K. Ackers. 2000b. Coupled energetics of cro repressor self-assembly and site-specific DNA operator binding. II: Cooperative interactions of cro dimers. J. Mol. Biol. 302:625–638.[Medline]
Dodd, I. B., A. J. Perkins, D. Tsemitsidis, and B. Egan. 2001. Octamerization of CI repressor is needed for repression of PRM and efficient switching from lysogeny. Genes &. Development. 15:3013–3022.[Abstract/Free Full Text]
Draper, D. E. 1996. Translational initiation. In Escherichia coli and Salmonella thyphymurium:
TOP
ABSTRACT
INTRODUCTION
and 
i = 1,2,3), interaction energies between adjacent CI2's or Cro2's (
and 
i < j) (in this case the dimer-operator binding energies are not considered), and RNA polymerase (RNAP) binding energies to PR and PRM (
and 
)
and 
i = 1,2,3), interaction energies between adjacent CI2's or Cro2's (
and 
i < j) (in this case the dimer-operator binding energies are not considered), and energies of RNA polymerase (RNAP) binding to PL (
)