# Bifurcation Analysis for Phage Lambda with Binding Energy Uncertainty

Authors: Ning Xu; Xue Lei; Ping Ao; Jun Zhang
Journal: Computational Biology Journal
 (2014)

---

## Abstract

In a phageλ genetic switch model, bistable dynamical behavior can be destroyed due to the bifurcation caused by inappropriately chosen model parameters. Since the values of many parameters with biological
significance often cannot be accurately acquired, it is thus of fundamental importance to analyze how and to
which extent the system dynamics is influenced by model parameters, especially those parameters pertaining
to binding energies. In this paper, we apply a Jacobian method to investigate the relation between bifurcation
and parameter uncertainties on a phage λ OR model. By introducing bistable range as a measure of system robustness, we find that RNA polymerase binding energies have the minimum bistable ranges among all the
binding energies, which implies that the uncertainties on these parameters tend to demolish the bistability
more easily. Moreover, parameters describing mutual prohibition between proteins Cro and CI have finite
bistable ranges, whereas those describing self-prohibition have infinity bistable ranges. Hence, the former are
more sensitive to parameter uncertainties than the latter. These results help to understand the influence of
parameter uncertainties on the bifurcation and thus bistability.


---

## Body

## 1. Introduction

Bistability is a salient feature of phageλ genetic switch. Mathematical models with a number of parameters have been established to describe its bistable behavior. However, parameter variation in phage λ model may considerably change the system dynamics. Therefore, in the case when the model parameters are not precisely known, it is useful to study the influence of parameter uncertainty on system dynamics.Since 1980s, researchers have developed several models for phageλ genetic switch. Reference [1] proposed a quantitative model for phage λ according to the principles of statistical thermodynamics. The model by [2] included two variables, that is, concentrations of proteins CI and Cro. Reference [3] extended this model and pointed out the discrepancy between experimental and theoretical results. Moreover, models with four variables were built to include more ingredients such as mRNA concentration in phage λ [4]. Other useful models include Markov chain model [5], delayed reaction stochastic model [6], fuzzy logic model [7], decision making [8], and qualitative model [9]. Based on these models, significant research progresses on phage λ have been made [10–15].Most of these models contain parameters that cannot be precisely determined. For example, the model by [2] has more than 10 binding energy related parameters, each of which may vary up to ±0.06 kcal/mol and thus affects system dynamics. Nonspecific binding also introduces additional factors that influence the system models [16]. Reference [17] developed a model that has 13 binding energy related parameters with uncertainty ranges around 0.5 kcal/mol. Moreover, [18] pointed out that some parameters such as protein synthesis rates can be changed as surrounding environment changes. Reference [19] reported that single point mutation of DNA leads to ±2 kcal/mol change for binding free energy. Due to these parameter uncertainties, it is often desired to analyze the effect of these parameters on system dynamics.Researchers have been using sensitivity and robustness to study the influence of parameter uncertainty. Sensitivity and robustness are closely related because a robust system is not sensitive to parameter variations [20]. By using lysis frequency as a sensitivity measure, [21] simulated a mutation model introduced by [22] and concluded that spontaneous switching from lysogen to lysis depends sensitively on model parameters. By systematically perturbing binding energies with ±1 kcal/mol, [23] concluded that lysogenic activity at promoter PRM remains stable under this perturbation.In this paper we investigate bifurcation by using bistable range as a system robustness measure. In a typical case, phageλ genetic switch model has three equilibria. Two of them are stable corresponding to the lysogenic and lytic states, and the third one is a saddle point. When some parameters are perturbed, it is possible that one stable equilibrium and the saddle point coalesce, then both of them disappear, resulting in the loss of bistable behavior. Hence, bistable region can be defined as the parameter region in which bistable behavior exists [18]. The bistable region for a parameter is large if the system dynamics is not sensitive to that parameter.We are interested in analyzing bistable ranges for binding energy related parameters. Some other parameters such as protein degradation rate and protein synthesis rate with respect to gene transcription have already been proved that their variations may lead to the loss of bistability [3, 24]. Our calculations are based on a phage λ OR model [3]. Although DNA looping formed by cooperations between the left operators and the right operators affects the stability of the genetic switch [25], the OR model only includes the effects of the right operators, which are believed to play a critical role for lysogeny maintenance [26]. We will calculate bistable ranges for all binding energies, which are crucial in determining the affinity that proteins bind to phage λ genes. Furthermore, we define Z-shaped bifurcation and calculate Jacobian matrix to study robustness. We find out that different binding energy affects lysogeny stability differently: some of them may destroy bistability, whereas others have little influence on system dynamics. Therefore, we need to precisely determine the values of those sensitive parameters in the model because otherwise we may not get an effective bistable switch.


## 2. Phageλ Modeling

We first present phageλ gene regulatory networks and then introduce their mathematical model. Lysogeny and lysis are two possible states for phage λ. It is believed that wild-type phage λ has quite a stable lysogenic state since lysogeny loss rate is even lower than 10-8 per cell generation [27]. However, certain conditions such as the exposure to UV light [28, 29], the presence of mitomycin C [30], or the starvation of the host cells [31] can induce phage λ from lysogeny to lysis. Mathematical models have been built to successfully explain stability and high efficiency of genetic switch.The concentrations of two proteins CI and Cro can serve as an indication of which state phageλ is in, lysogeny or lysis. In lysogeny, the concentration of protein CI in the cell is significantly higher than that of protein Cro, whereas in lysis, protein Cro dominates the cell. The amount of CI and Cro thus implies which state phage λ is in.Figure1 illustrates the gene segment regulating cro and cI gene expression. Three binding sites OR1, OR2, and OR3 control the expression of two kinds of regulatory proteins. The PRM promoter covers the whole binding site OR3 and part of OR2, controlling cI gene expression. When RNA polymerase (RNAP) binds to the site OR3, the gene expression is on, resulting in the production of protein CI. The PR promoter covers all of OR1 and part of OR2, controlling cro expression. When RNAP binds to sites OR1 and OR2 simultaneously, the cro gene expression is on. The binding sites may be occupied by either Cro or CI protein. Usually Cro and CI first form their dimers so that they are able to bind to OR1, OR2, and OR3, prohibiting the gene expressions on these sites.Figure 1
Three binding sites on a DNA segment. The sitesOR1, OR2, and OR3 can be bound by CI dimer, Cro dimer, or RNA polymerase. This DNA segment consists of promoters PRM and PR. If the RNA polymerase binds to OR3, cI gene is expressed; if the RNA polymerase binds to OR1 and OR2, cro gene is expressed. The amount of Cro and CI in the cell directly determines which state phage λ is in.Several phageλ models have been established based on this regulating mechanism. In this paper we adopt the well-known two-dimensional OR model stated in [3] with updated parameter values given in Table 1. This model involves three OR gene binding sites, which are believed to be the most critical factors for the stability of genetic switch [33]. This two-dimensional model is in line with the general biological switch model for two-gene networks [34]. To keep it simple, we eliminate the effects from other activities in the cell, such as nonspecific binding [16, 17], stochastic gene expression [35], and intrinsic and extrinsic noises [36]. The model is described by the following system equations:
(1)dNrdt=A1fr-μ1Nr,dNcdt=A2fc-μ2Nc,
where
(2)fr=P28+P29+P30+P40,fc=∑i=3133‍Pi+k(∑j=3440‍Pj),(3)Ps=exp(-ΔGs/RT)Nr2i(s)Nc2j(s)NRNAPk(s)∑m=140exp(-ΔGm/RT)Nr2i(m)Nc2j(m)NRNAPk(m),hhhhhhhhhhhhhhhhllhhhhh(s=1,2,…,40).
The subscript s corresponds to 40 binding configurations shown in Table 2. The conventions between CI dimer (Nr2), Cro dimer (Nc2), and their monomers (Nr and Nc) are
(4)Nr2=12Nr+18Kr-18Kr2+8NrKr,Nc2=12Nc+18Kc-18Kc2+8NcKc.Table 1
Parameter values for our model.
Δ
G
1
c
=
-
12.0 kcal/mol†
Δ
G
2
c
=
-
10.8 kcal/mol†
Δ
G
3
c
=
-
13.4 kcal/mol†
Δ
G
1
r
=
-
12.5 kcal/mol‡
Δ
G
2
r
=
-
10.5 kcal/mol‡
Δ
G
3
r
=
-
9.5 kcal/mol‡
Δ
G
12
r
=
-
2.7 kcal/mol‡
Δ
G
23
r
=
-
2.9 kcal/mol‡
Δ
G
R
12
=
-
12.5 kcal/mol†
Δ
G
R
3
=
-
11.5 kcal/mol†
K
r
=
5.56
×
1
0
-
9 M†
K
c
=
3.26
×
1
0
-
7 M†
N
RNAP
=
3.0
×
1
0
-
8 M§
μ
1
=
2.0
×
1
0
-
2 min−1†
μ
2
=
3.6
×
1
0
-
2 min−1†
R
=
1.99 cal/(mol·K)
T
=
310 K
k
=
1
/
11
¶
†denotes the parameter values from [14], ‡from [32], §from [21], and ¶from [3].Table 2
40 Binding configurations for the three sites.
s
O
R
3
O
R
2
O
R
1
1
2
CI2
3
CI2
4
CI2
5
CI2↔
CI2
6
CI2
CI2
7
CI2↔
CI2
8
CI2↔
CI2↔
CI2
9
Cro2
10
Cro2
11
Cro2
12
Cro2
Cro2
13
Cro2
Cro2
14
Cro2
Cro2
15
Cro2
Cro2
Cro2
16
Cro2
CI2
17
Cro2
CI2
18
Cro2
Cro2
CI2
19
Cro2
CI2
20
CI2
Cro2
21
Cro2
CI2
Cro2
22
CI2
Cro2
23
CI2
Cro2
24
CI2
Cro2
Cro2
25
Cro2
CI2↔
CI2
26
CI2
Cro2
CI2
27
CI2↔
CI2
Cro2
28
RNAP
29
CI2
RNAP
30
Cro2
RNAP
31
RNAP
CI2
32
RNAP
CI2↔
CI2
33
RNAP
CI2
Cro2
34
RNAP
35
RNAP
CI2
36
RNAP
Cro2
37
RNAP
Cro2
38
RNAP
Cro2
CI2
39
RNAP
Cro2
Cro2
40
RNAP
RNAPFrom (1)–(4), Nr, Nc represent concentrations of proteins CI and Cro. The parameters A1 and A2 are synthesis rates of CI and Cro, respectively; μ1 and μ2 are the overall rates that decrease protein concentration, which are caused by normal protein degradation and cell growth [37]. The function fr (or fc) represents the overall probability that the PR (or PRM) promoter is occupied by RNAP.   Ps is the probability that sth configuration occurs, and ΔGs is the total binding energy needed for this configuration. The exponents i(s), j(s), and k(s) represent the numbers of CI dimers, Cro dimers, and RNAP at sth configuration, respectively. ΔGs is calculated by summing up the needed single binding energies listed in Table 1, including cooperation energy. For example, ΔG5=ΔG12r+ΔG1r+ΔG2r. The left right arrows in Table 2 indicate the cooperations. k is the coefficient for positive feedback. When CI dimer binds to OR2, positive feedback makes the transcription of cI gene more efficient, hence promotes its expression. Kr, Kc are dissociation constants for CI and Cro proteins, respectively. R represents the universal gas constant and T the absolute temperature.


## 3. Bifurcation in Phageλ

The model used in this paper is two-dimensional nonlinear ordinary differential equations. To get all the equilibria, we set the left hand side of (1) to zero and solve the transcendental equations. It can be found that in general the model has three equilibria. However, the number of equilibria may vary due to bifurcation under parameter variation. In phage λ model, frequently occurring is the saddle node bifurcation, whose precise mathematical definition is given as follows [38].Definition 1 (saddle node bifurcation).
A saddle node bifurcation in dynamical system takes place when two equilibrium points coalesce and disappear.In (1), two parameters A1 and A2 are not specified in Table 1. They represent proportional constants relating the transcription initiation rate to the absolute rate of CI and Cro protein synthesis [3]. We use Figure 2 to illustrate Z-shaped bifurcation with the variation of A2. Set A1=8×10-9 M/min and divide the whole plane into three regions as shown in Figure 2, labeled by A, B, and C. Initially, when the value of A2 is small in region A, for example, A2≈5×10-9 M/min, there exists only one stable point and thus one single curve segment for both Nr and Nc. At the boundary of regions A and B, saddle node bifurcation occurs: A2 increases to region B, for example, a typical value A2≈5×10-8 M/min, and two new equilibria (one saddle and one stable equilibrium) emerge, resulting in three equilibria. Note that this is a saddle node bifurcation with the reversed direction, in comparison to Definition 1.Figure 2
Saddle node bifurcation with the variation of Cro synthesis rateA2. This is a schematic plot and not drawn to scale.At the boundary of regions B and C, saddle node bifurcation takes place again, leading to the vanishing of two equilibria (one saddle point and one stable equilibrium). In region C, sayA2≈5×10-7 M/min, only one stable equilibrium remains. We notice the motion of the saddle point. When the saddle point shows up at the boundary of regions A and B, the lysis stable equilibrium appears as well. As parameter A2 increases, the saddle point moves toward the lysogen stable point and finally merges with it at the boundary of regions B and C. Two successive bifurcations form a Z-shaped curve are shown in Figure 2.We have not yet specified the boundaries between these regions quantitatively. Next we solve the values of the boundaries by using the Jacobian method. Let us first introduce the following two theorems [38].Theorem 2 (Hartman-Grobman Theorem).
If the Jacobian of the nonlinear systemx˙=g(x) at equilibrium, namely J(x0), has no zero or purely imaginary eigenvalues, then the phase trajectory of the system x˙=g(x) resembles that of the linearized system z˙=J(x0)z near x0 if none of the eigenvalues of J(x0) are on the jω axis.For a given equilibriumx0, in most cases, we can determine its dynamical behavior qualitatively by calculating its Jacobian matrix J(x0). For a two-dimensional system, its Jacobian matrix can be calculated as
(5)J=[∂g1∂Nr∂g1∂Nc∂g2∂Nr∂g2∂Nc].The following theorem gives the necessary conditions for the critical point that bifurcation occurs [38].Theorem 3.
For any nonlinear systemz˙=g(z,μ), the system goes through a saddle node bifurcation at z=z†, μ=μ† only if the following relations are satisfied:
(6)g(z†,μ†)=0,(7)|∂g∂z(z†,μ†)|=0.From (6), z† is an equilibrium point. In a planar system, (7) is the determinant of system Jacobian. When the determinant of Jacobian matrix vanishes, at least one eigenvalue is zero. Jacobian is a characterization of planar system dynamics. For variation of parameter A2, the system in (1) can be rewritten as
(8)dNrdt=g1(Nc,Nr,A2),dNcdt=g2(Nc,Nr,A2).
From Theorems 2 and 3, for a given A2, the Jacobian matrix is evaluated at equilibrium points Nc† and Nr†. These two equilibrium points can be represented as a function of A2 as shown in Figure 2.By Theorem3, we know that the bifurcation appears when Jacobian has a zero eigenvalue. From Figure 3, we can obtain the approximate values of critical points by checking the points of the eigenvalue curves crossing zero line. These values are 5.6×10-9 M/min and 3.4×10-7 M/min.Figure 3
Jacobian eigenvalues at equilibria. The figure shows how the eigenvalues of Jacobian at equilibria change withA1 fixed at 8×10-9 M/min and A2 varying. The solid line (black) shows the eigenvalues of lysis equilibrium. The dashed line (blue) is for lysogeny equilibrium. The dot-dashed line (red) is for saddle.To determine the stability of an equilibrium point, we calculate the eigenvalue of the Jacobian matrix. Figure3 plots the eigenvalues as a function of A2. The Jacobian matrix can be calculated analytically (see the Appendix for details). By solving the value of A2 that makes the determinant of the Jacobian matrix vanish, we get the critical value of A2 at which bifurcation occurs. When A2 is small (less than 5.6×10-9 M/min in Figure 3), there are only two negative eigenvalues corresponding to one equilibrium point. This equilibrium is thus stable [38].AsA2 exceeds 5.6×10-9 M/min, we have the case when there are three equilibria and correspondingly six eigenvalues. Among these three equilibria, two are stable and the third is unstable (saddle point). This is because positive eigenvalue appears at saddle point (see dash-dot line in Figure 3). For a two-dimensional linear system, if the Jacobian has one positive and one negative eigenvalues at the equilibrium, the equilibrium is a saddle point. If A2 continues increasing, only one equilibrium point is left (A2 greater than 3.4×10-7 M/min in this case), and both of its eigenvalues are negative. The curves in Figure 3 have two intersections with zero eigenvalue, which indicates that the eigenvalues of Jacobian vanish twice, forming a Z-shaped bifurcation.


## 4. Bifurcation Analysis of Binding Energy Uncertainty

In this section, we will analyze bifurcation forΔG’s by using the Jacobian method. In Table 1, the symbol ΔG with subscripts such as ΔG1r and ΔG2c represents Gibbs free energy required for CI or Cro dimers to bind to the OR sites. Because of the difficulty in determining the precise values of the binding energy, it is useful to study how the equilibria change as these Gibbs free binding energies vary. In this section, we set synthesis rates A1 and A2 at appropriate values, A1=8×10-9 M/min and A2=5×10-8 M/min. To calculate the equilibria, we set the left hand side of (1) to zero and obtain the expression of Nr and Nc as functions of ΔG’s:
(9)Nr*=f1(ΔG1c,ΔG2c,…,ΔGR3),Nc*=f2(ΔG1c,ΔG2c,…,ΔGR3),
where Nr* and Nc* are the concentrations of proteins CI and Cro at the equilibrium point. Once (9) is obtained, we can determine the quantitative behavior around equilibria when changing ΔG’s. The relative sensitivities are given by
(10)SNr(ΔGs)=∂f1∂ΔGs,SNc(ΔGs)=∂f2∂ΔGs,
where the subscript s∈{1r,2r,…,R12,R3}. Since the system equilibria cannot be analytically solved, we may apply the numerical procedure to calculate (10) as a measure of parameter sensitivity. As to the study on system robustness, the Jacobian method can be used to study bifurcation introduced by parameter variation.Since there are ten differentΔGs, we divide them into four groups:(1)
ΔG1r, ΔG2r, and ΔG3r: binding free energies of CI;(2)
ΔG1c, ΔG2c, and ΔG3c: binding free energies of Cro;(3)
ΔG12r and ΔG23r: the cooperation between proteins CI’s if both OR1 and OR2, or OR2 and OR3 are bound by CI dimers;(4)
ΔGR12 and ΔGR3: RNAP binding energies. We then study each of these four groups.


### 4.1.ΔG1r, ΔG2r, and ΔG3r

We systematically perturb the values ofΔG1r, ΔG2r, and ΔG3r to investigate their impact on the system dynamics. Two proper robustness measures are the range of parameters in which lysogeny exists (stable range), and in which bistable dynamics exists (bistable range). In stable range, it is possible to keep phage λ dormant and integrate into the gene of the host. On the other hand, bistable range defines the parameter range to keep bistability and the maximal uncertainty tolerance. In this case, it implies that the system favors dynamics with two stable equilibria and a saddle. Further, if the stochastic force is strong enough, phage λ can switch from the lysogenic stable point to the lytic stable point even without the parameter changes.Figures4(a)–4(c) illustrate the impact of ΔG1r, ΔG2r, and ΔG3r variations. We set parameter values as in Table 1 and find out all the equilibria when varying one ΔGs. We find that the energy ΔG3r does not show Z-shaped bifurcation, while the other two parameters ΔG1r and ΔG2r show quite similar behaviors since Figures 4(a) and 4(b) are almost identically the same. The result is not surprising because ΔG1r and ΔG2r play similar roles in system dynamics. These two parameters determine the likelihood that CI dimer binds to sites OR1 and OR2. When either of the two sites is occupied by protein dimers, the expression of cro gene is prohibited, resulting in lysogenic state. Nonetheless, the parameter ΔG3r indicates the likelihood of CI dimer binding to the site OR3. If the absolute value of ΔG3r is too large, it is easy for CI dimer to bind to OR3; thus, cI gene expression is prohibited, resulting in lysis. If the absolute value of ΔG3r is small, the Z-shaped bifurcation vanishes. In this case, cI gene expression is not prohibited, but it has no influence on cro gene expression. Bistable behavior always exists when ΔG3r is less than its nominal value in Table 1. Without stopping or weakening the expression of cro gene, single lysogen stable point layout cannot be achieved just by increasing ΔG3r. Compared with other parameters, the system dynamics is more robust to ΔG3r variation.Movement of equilibria under parameter variation. Z-shaped bifurcation can be generated by variation ofΔG1r, ΔG2r, and ΔG3c. All vertical axes are plotted on a log scale.
(a)
Z-shaped bifurcation forΔG1r
(b)
Z-shaped bifurcation forΔG2r
(c)
Bifurcation forΔG3r
(d)
Bifurcation forΔG1c
(e)
Bifurcation forΔG2c
(f)
Z-shaped bifurcation forΔG3cTable3 lists the stable range, bistable range, and the length of bistable range for all ten ΔGs. It can be shown by Theorem 2 that one of the three equilibria is a saddle and the other two are stable points. Jacobian matrix is also used to determine the boundary of stable and bistable ranges according to Theorem 3. As long as ΔG1r is less than the nominal value at -11.33kcal/mol, or ΔG2r is less than -9.34kcal/mol, the system always has lysogenic state. They are critical bifurcation points. For ΔG3r, however, lysogenic state always exists if it is greater than -12.33kcal/mol. Consequently, for CI protein, a larger binding energy to OR1 and OR2 helps to keep the DNA of phage λ silent in the host. The system tends to be unstable if the values of binding energies for CI to OR1 and OR2 become small, thus forming Z-shaped bifurcation. Parameter variations may lead to Z-shaped bifurcation if and only if the bistable range is bounded.Table 3
Robustness measures of binding energies.
Variable
Stable Range (kcal/mol)
Bistable Range (kcal/mol)
Length (kcal/mol)
Δ
G
1
r
(
-
∞
,
-
11.33
)
(
-
15.49
,
-
11.33
)
4.16
Δ
G
2
r
(
-
∞
,
-
9.34
)
(
-
13.04
,
-
9.34
)
3.70
Δ
G
3
r
(
-
12.33
,
+
∞
)
(
-
12.33
,
+
∞
)
∞
Δ
G
1
c
(
-
∞
,
+
∞
)
(
-
15.14
,
+
∞
)
∞
Δ
G
2
c
(
-
∞
,
+
∞
)
(
-
15.12
,
+
∞
)
∞
Δ
G
3
c
(
-
15.67
,
+
∞
)
(
-
15.67
,
-
12.24
)
3.43
Δ
G
12
r
(
-
∞
,
-
1.55
)
(
-
6.03
,
-
1.55
)
4.48
Δ
G
23
r
(
-
7.71
,
∞
)
(
-
7.71
,
+
∞
)
∞
Δ
G
R
12
(
-
13.67
,
∞
)
(
-
13.67
,
-
10.86
)
2.81
Δ
G
R
3
(
-
∞
,
-
10.50
)
(
-
12.66
,
-
10.50
)
2.16When the stochastic noise is strong enough, the noise may induce the system to lytic state. Moreover, if free energy varies due to gene mutation, bistable behavior might be destroyed and only one stable equilibrium is left. This is what the mathematical model implies about the parameter influence on the system dynamics.


### 4.2.ΔG1c, ΔG2c, and ΔG3c

The bifurcation results ofΔG1c, ΔG2c, and ΔG3c are shown in Figures 4(d)–4(f). The equilibria behaviors with respect to ΔG1c and ΔG2c variation are similar to ΔG3r. ΔG1c represents the likelihood that Cro dimer binds to site OR1. If Cro dimer binds to site OR1, the expression of cro is prohibited, resulting in only one single lysogenic stable point. When ΔG1c or ΔG2c is small, the prohibition effect is weak, but this is not enough to drive the system dynamics to only one single lysis stable point. Hence, bistable layout always exists if these two parameters are small. The influence of ΔG3c is similar to the influence of ΔG1r and ΔG2r. The bistable range of ΔG3c is about 3.43 kcal/mol, which is comparable to those of ΔG1r and ΔG2r.The six parameters we have shown so far are the basic binding energies. The values ofΔG3r, ΔG1c, and ΔG2c decide the likelihood that negative feedback (i.e., self-prohibition) happens. It is known that the sites OR1 and OR2 are used to transcript and then translate protein Cro. The site OR3 contributes to the production of protein CI. If protein CI binds to OR3, or Cro binds to OR1 and OR2, they can stop their own transcription. Thus, we call the phenomenon negative feedback. Figure 4 shows that such negative feedback cannot take the system to Z-shaped bifurcation. The negative feedback occurs when the amount of Cro or CI is too large. For instance, the binding between CI dimer and OR3 readily occurs when the amount of protein CI is too large. On the other hand, ΔG1r, ΔG2r, and ΔG3c represent the likelihood of positive feedback (i.e., mutual prohibition). For instance, with large ΔG3c, the binding between Cro and OR3 is favored, stopping the expression of cI gene. If the value of ΔG3r is small, and the negative feedback is weakened, but it is not enough to make the system dynamics change from single lysis equilibrium to single lysogenic equilibrium. Parameters associated with positive feedback can show Z-shaped bifurcation, whereas parameters associated with negative feedback cannot. Hence, the system is more robust to negative feedback parameters than to positive feedback parameters.


### 4.3.ΔG12r and ΔG23r

The parametersΔG12r and ΔG23r determine the energy difference caused by the cooperation of adjacent binding CI dimers. Although [32] discovered that adjacent Cro dimers also have cooperations, we ignore this cooperation effect since it has a small cooperation energy and focus on CI cooperation in our analysis. Figure 5 shows the bifurcation results for the parameters ΔG12r and ΔG23r. The nominal values of ΔG12r and ΔG23r are close to zero. Since protein cooperation process is an exothermic reaction, ΔG12r and ΔG23r must be less than zero. We find out that the variation of ΔG12r leads to Z-shaped bifurcation whereas ΔG23r does not. The value of ΔG23r affects the expression of both cro and cI genes, whereas ΔG12r only promotes mutual prohibition from CI to Cro. ΔG23r with large value prohibits the binding of RNAP to both sites OR2 and OR3, thus both cI and cro gene expressions are closed.Z-shaped bifurcation generated by variation ofG12r and G23r. The variation of ΔG12r also leads to Z-shaped bifurcation, while ΔG23r is not sensitive enough to create Z-shaped bifurcation.
(a)
Z-shaped bifurcation forΔG12r
(b)
Bifurcation forΔG23rMoreover, we know the length of bifurcation range forΔG12r is 4.48kcal/mol, which is the maximal value among all the finite values in Table 3. Compared with the influence from other parameters, the system is relatively more robust to the cooperation energy variations. The bistable range for ΔG12r is from -6.03kcal/mol to -1.55kcal/mol (Table 3). In our model, if ΔG12r is equal to the center value at -3.79kcal/mol, the system is most robust to parameter uncertainty because of maximum tolerance of uncertainty. This result is consistent with the calculation by [39], where they calculated biological behaviors in phage λ and found that cooperation energy is -3.7kcal/mol, in comparison to the in vitro value -2.7kcal/mol used in this work. The 1.0kcal/mol difference between these two values might be caused by the looping effect [40]. Hence, the looping effect makes the system more robust to ΔG12r from this point of view.


### 4.4.ΔGR12 and ΔGR3

Figure6 illustrates the influence of ΔGR3 and ΔGR12 variations. From Figure 6, we find that Z-shaped bifurcation appears for both parameters. Note that the appearance of Z is different from the previous parameters as discussed above. As in Table 3, the length of bistable range for ΔGR12 is 2.81 kcal/mol, and for ΔGR3 is only 2.16kcal/mol, which is the minimum among all ten parameters.Z-shaped bifurcation generated by variation ofGR12 and GR3. Both of them are able to generate Z-shaped bifurcation. Their bistable ranges are obviously smaller than any other parameters. The appearance of “Z” shows steep slope.
(a)
Z-shaped bifurcation forΔGR12
(b)
Z-shaped bifurcation forΔGR3Notice thatΔGR3 and ΔGR12 are the binding free energies of RNAP to PRM promoter (OR3) and to PR promoter (OR1 and OR2). Because RNAP is necessary for gene expression, the result from RNAP could be the most influential. Table 3 summarizes all the sensitivity measures for the above ten parameters. We observe that the last two parameters have the shortest length of bistable ranges, which implies that bistable behavior of phage λ is more sensitive to ΔGR12 and ΔGR3 than any other parameters. For example, the bistable ranges of ΔG1r and ΔG2r are around 4 kcal/mol, while those of ΔGR12 and ΔGR3 are less than 3 kcal/mol. Parameter uncertainties of ΔGR12 and ΔGR3 most easily destroy bistable behavior. The concentrations of CI and Cro at equilibrium point also change significantly under ΔGR12 or ΔGR3 uncertainty.


### 4.5. Discussion

From the above analysis, we find that the binding free energies to sitesOR1 and OR2 have similar effects on system dynamics. As we have shown in Table 3, the results of ΔG1r, ΔG2r, ΔG1c, and ΔG2c are similar. This is mainly because of the similar status for binding sites OR1 and OR2, both of which control the expression of Cro protein.Moreover, the bistable ranges have infinite length forΔG3r, ΔG1c, ΔG2c, and ΔG23r. All of the above four parameters do not have Z-shaped bifurcation. When the values of those parameters become close to zero, bistable behavior always exists. The parameters ΔG1c and ΔG2c can be used to maintain lysogenic state when they are less than a critical value because their stable ranges do not have any limits. For example, if we introduce gene mutation to make ΔG1c lower than −15.14 kcal/mol, bistable behavior disappears and only one lysogenic state is left. Consequently, the state of phage λ stabilizes at lysogeny.As to the other six parameters with limited bistable ranges in Table3, their variations lead to Z-shaped bifurcations. We are interested in why the original system seems more sensitive to these six parameters. The free energies ΔG1r, ΔG2r, and ΔG12r indicate the likelihood of protein CI binding to sites OR1 and OR2, or the likelihood of mutual prohibition occurring. Similarly, ΔG3c indicates the likelihood of protein Cro binding to site OR3 and also results in mutual prohibition. Our result has shown that, compared with self-prohibition, the mutual prohibition is more influential on system dynamics. The binding energy uncertainty related to mutual prohibition may lead to the loss of bistable behavior.Note that the system dynamics of this model is even more sensitive to the parametersΔGR12 and ΔGR3. Since ΔGR12 and ΔGR3 have the minimum bistable ranges, the uncertainties of these two parameters are more likely to lead to the loss of bistable behavior. This is consistent with the conclusion by [23] that perturbation of RNAP binding energy significantly changes promoter activity, thus affecting gene expression. The binding free energies of RNAP to cro and cI promoters play a key role for the generation of proteins.In building mathematical models for phageλ genetic switch, it is crucial to determine the values for those parameters that are more sensitive to parameter uncertainties such as ΔGR12 and ΔGR3, because otherwise the system dynamics may deviate from normal bistable behavior. Furthermore, our results provide guidelines for the experimentalists in laboratory on how to choose those parameters that can alter mutations more easily.


## 4.1.ΔG1r, ΔG2r, and ΔG3r

We systematically perturb the values ofΔG1r, ΔG2r, and ΔG3r to investigate their impact on the system dynamics. Two proper robustness measures are the range of parameters in which lysogeny exists (stable range), and in which bistable dynamics exists (bistable range). In stable range, it is possible to keep phage λ dormant and integrate into the gene of the host. On the other hand, bistable range defines the parameter range to keep bistability and the maximal uncertainty tolerance. In this case, it implies that the system favors dynamics with two stable equilibria and a saddle. Further, if the stochastic force is strong enough, phage λ can switch from the lysogenic stable point to the lytic stable point even without the parameter changes.Figures4(a)–4(c) illustrate the impact of ΔG1r, ΔG2r, and ΔG3r variations. We set parameter values as in Table 1 and find out all the equilibria when varying one ΔGs. We find that the energy ΔG3r does not show Z-shaped bifurcation, while the other two parameters ΔG1r and ΔG2r show quite similar behaviors since Figures 4(a) and 4(b) are almost identically the same. The result is not surprising because ΔG1r and ΔG2r play similar roles in system dynamics. These two parameters determine the likelihood that CI dimer binds to sites OR1 and OR2. When either of the two sites is occupied by protein dimers, the expression of cro gene is prohibited, resulting in lysogenic state. Nonetheless, the parameter ΔG3r indicates the likelihood of CI dimer binding to the site OR3. If the absolute value of ΔG3r is too large, it is easy for CI dimer to bind to OR3; thus, cI gene expression is prohibited, resulting in lysis. If the absolute value of ΔG3r is small, the Z-shaped bifurcation vanishes. In this case, cI gene expression is not prohibited, but it has no influence on cro gene expression. Bistable behavior always exists when ΔG3r is less than its nominal value in Table 1. Without stopping or weakening the expression of cro gene, single lysogen stable point layout cannot be achieved just by increasing ΔG3r. Compared with other parameters, the system dynamics is more robust to ΔG3r variation.Movement of equilibria under parameter variation. Z-shaped bifurcation can be generated by variation ofΔG1r, ΔG2r, and ΔG3c. All vertical axes are plotted on a log scale.
(a)
Z-shaped bifurcation forΔG1r
(b)
Z-shaped bifurcation forΔG2r
(c)
Bifurcation forΔG3r
(d)
Bifurcation forΔG1c
(e)
Bifurcation forΔG2c
(f)
Z-shaped bifurcation forΔG3cTable3 lists the stable range, bistable range, and the length of bistable range for all ten ΔGs. It can be shown by Theorem 2 that one of the three equilibria is a saddle and the other two are stable points. Jacobian matrix is also used to determine the boundary of stable and bistable ranges according to Theorem 3. As long as ΔG1r is less than the nominal value at -11.33kcal/mol, or ΔG2r is less than -9.34kcal/mol, the system always has lysogenic state. They are critical bifurcation points. For ΔG3r, however, lysogenic state always exists if it is greater than -12.33kcal/mol. Consequently, for CI protein, a larger binding energy to OR1 and OR2 helps to keep the DNA of phage λ silent in the host. The system tends to be unstable if the values of binding energies for CI to OR1 and OR2 become small, thus forming Z-shaped bifurcation. Parameter variations may lead to Z-shaped bifurcation if and only if the bistable range is bounded.Table 3
Robustness measures of binding energies.
Variable
Stable Range (kcal/mol)
Bistable Range (kcal/mol)
Length (kcal/mol)
Δ
G
1
r
(
-
∞
,
-
11.33
)
(
-
15.49
,
-
11.33
)
4.16
Δ
G
2
r
(
-
∞
,
-
9.34
)
(
-
13.04
,
-
9.34
)
3.70
Δ
G
3
r
(
-
12.33
,
+
∞
)
(
-
12.33
,
+
∞
)
∞
Δ
G
1
c
(
-
∞
,
+
∞
)
(
-
15.14
,
+
∞
)
∞
Δ
G
2
c
(
-
∞
,
+
∞
)
(
-
15.12
,
+
∞
)
∞
Δ
G
3
c
(
-
15.67
,
+
∞
)
(
-
15.67
,
-
12.24
)
3.43
Δ
G
12
r
(
-
∞
,
-
1.55
)
(
-
6.03
,
-
1.55
)
4.48
Δ
G
23
r
(
-
7.71
,
∞
)
(
-
7.71
,
+
∞
)
∞
Δ
G
R
12
(
-
13.67
,
∞
)
(
-
13.67
,
-
10.86
)
2.81
Δ
G
R
3
(
-
∞
,
-
10.50
)
(
-
12.66
,
-
10.50
)
2.16When the stochastic noise is strong enough, the noise may induce the system to lytic state. Moreover, if free energy varies due to gene mutation, bistable behavior might be destroyed and only one stable equilibrium is left. This is what the mathematical model implies about the parameter influence on the system dynamics.


## 4.2.ΔG1c, ΔG2c, and ΔG3c

The bifurcation results ofΔG1c, ΔG2c, and ΔG3c are shown in Figures 4(d)–4(f). The equilibria behaviors with respect to ΔG1c and ΔG2c variation are similar to ΔG3r. ΔG1c represents the likelihood that Cro dimer binds to site OR1. If Cro dimer binds to site OR1, the expression of cro is prohibited, resulting in only one single lysogenic stable point. When ΔG1c or ΔG2c is small, the prohibition effect is weak, but this is not enough to drive the system dynamics to only one single lysis stable point. Hence, bistable layout always exists if these two parameters are small. The influence of ΔG3c is similar to the influence of ΔG1r and ΔG2r. The bistable range of ΔG3c is about 3.43 kcal/mol, which is comparable to those of ΔG1r and ΔG2r.The six parameters we have shown so far are the basic binding energies. The values ofΔG3r, ΔG1c, and ΔG2c decide the likelihood that negative feedback (i.e., self-prohibition) happens. It is known that the sites OR1 and OR2 are used to transcript and then translate protein Cro. The site OR3 contributes to the production of protein CI. If protein CI binds to OR3, or Cro binds to OR1 and OR2, they can stop their own transcription. Thus, we call the phenomenon negative feedback. Figure 4 shows that such negative feedback cannot take the system to Z-shaped bifurcation. The negative feedback occurs when the amount of Cro or CI is too large. For instance, the binding between CI dimer and OR3 readily occurs when the amount of protein CI is too large. On the other hand, ΔG1r, ΔG2r, and ΔG3c represent the likelihood of positive feedback (i.e., mutual prohibition). For instance, with large ΔG3c, the binding between Cro and OR3 is favored, stopping the expression of cI gene. If the value of ΔG3r is small, and the negative feedback is weakened, but it is not enough to make the system dynamics change from single lysis equilibrium to single lysogenic equilibrium. Parameters associated with positive feedback can show Z-shaped bifurcation, whereas parameters associated with negative feedback cannot. Hence, the system is more robust to negative feedback parameters than to positive feedback parameters.


## 4.3.ΔG12r and ΔG23r

The parametersΔG12r and ΔG23r determine the energy difference caused by the cooperation of adjacent binding CI dimers. Although [32] discovered that adjacent Cro dimers also have cooperations, we ignore this cooperation effect since it has a small cooperation energy and focus on CI cooperation in our analysis. Figure 5 shows the bifurcation results for the parameters ΔG12r and ΔG23r. The nominal values of ΔG12r and ΔG23r are close to zero. Since protein cooperation process is an exothermic reaction, ΔG12r and ΔG23r must be less than zero. We find out that the variation of ΔG12r leads to Z-shaped bifurcation whereas ΔG23r does not. The value of ΔG23r affects the expression of both cro and cI genes, whereas ΔG12r only promotes mutual prohibition from CI to Cro. ΔG23r with large value prohibits the binding of RNAP to both sites OR2 and OR3, thus both cI and cro gene expressions are closed.Z-shaped bifurcation generated by variation ofG12r and G23r. The variation of ΔG12r also leads to Z-shaped bifurcation, while ΔG23r is not sensitive enough to create Z-shaped bifurcation.
(a)
Z-shaped bifurcation forΔG12r
(b)
Bifurcation forΔG23rMoreover, we know the length of bifurcation range forΔG12r is 4.48kcal/mol, which is the maximal value among all the finite values in Table 3. Compared with the influence from other parameters, the system is relatively more robust to the cooperation energy variations. The bistable range for ΔG12r is from -6.03kcal/mol to -1.55kcal/mol (Table 3). In our model, if ΔG12r is equal to the center value at -3.79kcal/mol, the system is most robust to parameter uncertainty because of maximum tolerance of uncertainty. This result is consistent with the calculation by [39], where they calculated biological behaviors in phage λ and found that cooperation energy is -3.7kcal/mol, in comparison to the in vitro value -2.7kcal/mol used in this work. The 1.0kcal/mol difference between these two values might be caused by the looping effect [40]. Hence, the looping effect makes the system more robust to ΔG12r from this point of view.


## 4.4.ΔGR12 and ΔGR3

Figure6 illustrates the influence of ΔGR3 and ΔGR12 variations. From Figure 6, we find that Z-shaped bifurcation appears for both parameters. Note that the appearance of Z is different from the previous parameters as discussed above. As in Table 3, the length of bistable range for ΔGR12 is 2.81 kcal/mol, and for ΔGR3 is only 2.16kcal/mol, which is the minimum among all ten parameters.Z-shaped bifurcation generated by variation ofGR12 and GR3. Both of them are able to generate Z-shaped bifurcation. Their bistable ranges are obviously smaller than any other parameters. The appearance of “Z” shows steep slope.
(a)
Z-shaped bifurcation forΔGR12
(b)
Z-shaped bifurcation forΔGR3Notice thatΔGR3 and ΔGR12 are the binding free energies of RNAP to PRM promoter (OR3) and to PR promoter (OR1 and OR2). Because RNAP is necessary for gene expression, the result from RNAP could be the most influential. Table 3 summarizes all the sensitivity measures for the above ten parameters. We observe that the last two parameters have the shortest length of bistable ranges, which implies that bistable behavior of phage λ is more sensitive to ΔGR12 and ΔGR3 than any other parameters. For example, the bistable ranges of ΔG1r and ΔG2r are around 4 kcal/mol, while those of ΔGR12 and ΔGR3 are less than 3 kcal/mol. Parameter uncertainties of ΔGR12 and ΔGR3 most easily destroy bistable behavior. The concentrations of CI and Cro at equilibrium point also change significantly under ΔGR12 or ΔGR3 uncertainty.


## 4.5. Discussion

From the above analysis, we find that the binding free energies to sitesOR1 and OR2 have similar effects on system dynamics. As we have shown in Table 3, the results of ΔG1r, ΔG2r, ΔG1c, and ΔG2c are similar. This is mainly because of the similar status for binding sites OR1 and OR2, both of which control the expression of Cro protein.Moreover, the bistable ranges have infinite length forΔG3r, ΔG1c, ΔG2c, and ΔG23r. All of the above four parameters do not have Z-shaped bifurcation. When the values of those parameters become close to zero, bistable behavior always exists. The parameters ΔG1c and ΔG2c can be used to maintain lysogenic state when they are less than a critical value because their stable ranges do not have any limits. For example, if we introduce gene mutation to make ΔG1c lower than −15.14 kcal/mol, bistable behavior disappears and only one lysogenic state is left. Consequently, the state of phage λ stabilizes at lysogeny.As to the other six parameters with limited bistable ranges in Table3, their variations lead to Z-shaped bifurcations. We are interested in why the original system seems more sensitive to these six parameters. The free energies ΔG1r, ΔG2r, and ΔG12r indicate the likelihood of protein CI binding to sites OR1 and OR2, or the likelihood of mutual prohibition occurring. Similarly, ΔG3c indicates the likelihood of protein Cro binding to site OR3 and also results in mutual prohibition. Our result has shown that, compared with self-prohibition, the mutual prohibition is more influential on system dynamics. The binding energy uncertainty related to mutual prohibition may lead to the loss of bistable behavior.Note that the system dynamics of this model is even more sensitive to the parametersΔGR12 and ΔGR3. Since ΔGR12 and ΔGR3 have the minimum bistable ranges, the uncertainties of these two parameters are more likely to lead to the loss of bistable behavior. This is consistent with the conclusion by [23] that perturbation of RNAP binding energy significantly changes promoter activity, thus affecting gene expression. The binding free energies of RNAP to cro and cI promoters play a key role for the generation of proteins.In building mathematical models for phageλ genetic switch, it is crucial to determine the values for those parameters that are more sensitive to parameter uncertainties such as ΔGR12 and ΔGR3, because otherwise the system dynamics may deviate from normal bistable behavior. Furthermore, our results provide guidelines for the experimentalists in laboratory on how to choose those parameters that can alter mutations more easily.


## 5. Conclusion

In this paper we studied the robustness of a phageλ model by investigating its bifurcation behavior induced by parameter uncertainties. We introduced bistable range as a measure of robustness and then applied the Jacobian method to calculate it. It can be concluded that binding energy uncertainties may destroy bistable behavior, especially for those sensitive parameters such as ΔGR12 and ΔGR3. Moreover, parameters describing mutual prohibition can form Z-shaped bifurcation (i.e., ΔG1r, ΔG2r, ΔG3c, and ΔG12r), whereas parameters describing self-prohibition (i.e., ΔG3r, ΔG1c, and ΔG2c) cannot. Hence, system dynamics is more sensitive to mutual prohibition parameters.There also exist other factors affecting the stability and robustness of a phageλ model. For example, DNA looping can promote CI expression in lysogenic state [14] and thus help phage λ keep stable and robust [25]. Our future work includes exploring the bifurcation and system robustness with respect to these new factors.



---

*Source: 465216-2014-02-03.xml*