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.

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 [1015].

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 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 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.

Figure 1 illustrates the gene segment regulating cro and cI gene expression. Three binding sites , , and control the expression of two kinds of regulatory proteins. The promoter covers the whole binding site and part of , controlling cI gene expression. When RNA polymerase (RNAP) binds to the site , the gene expression is on, resulting in the production of protein CI. The promoter covers all of and part of , controlling cro expression. When RNAP binds to sites and 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 , , and , prohibiting the gene expressions on these sites.

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 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: where The subscript corresponds to 40 binding configurations shown in Table 2. The conventions between CI dimer (), Cro dimer (), and their monomers ( and ) are

From (1)–(4), , represent concentrations of proteins CI and Cro. The parameters and are synthesis rates of CI and Cro, respectively; and are the overall rates that decrease protein concentration, which are caused by normal protein degradation and cell growth [37]. The function (or ) represents the overall probability that the (or ) promoter is occupied by RNAP.    is the probability that th configuration occurs, and is the total binding energy needed for this configuration. The exponents , , and represent the numbers of CI dimers, Cro dimers, and RNAP at th configuration, respectively. is calculated by summing up the needed single binding energies listed in Table 1, including cooperation energy. For example, . The left right arrows in Table 2 indicate the cooperations. is the coefficient for positive feedback. When CI dimer binds to , positive feedback makes the transcription of cI gene more efficient, hence promotes its expression. , are dissociation constants for CI and Cro proteins, respectively. represents the universal gas constant and 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 and 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 . Set  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 is small in region A, for example,  M/min, there exists only one stable point and thus one single curve segment for both and . At the boundary of regions A and B, saddle node bifurcation occurs: increases to region B, for example, a typical value  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.

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, say  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 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 system at equilibrium, namely , has no zero or purely imaginary eigenvalues, then the phase trajectory of the system resembles that of the linearized system near if none of the eigenvalues of are on the axis.

For a given equilibrium , in most cases, we can determine its dynamical behavior qualitatively by calculating its Jacobian matrix . For a two-dimensional system, its Jacobian matrix can be calculated as

The following theorem gives the necessary conditions for the critical point that bifurcation occurs [38].

Theorem 3. For any nonlinear system , the system goes through a saddle node bifurcation at , only if the following relations are satisfied:

From (6), 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 , the system in (1) can be rewritten as From Theorems 2 and 3, for a given , the Jacobian matrix is evaluated at equilibrium points and . These two equilibrium points can be represented as a function of as shown in Figure 2.

By Theorem 3, 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  M/min and  M/min.

To determine the stability of an equilibrium point, we calculate the eigenvalue of the Jacobian matrix. Figure 3 plots the eigenvalues as a function of . The Jacobian matrix can be calculated analytically (see the Appendix for details). By solving the value of that makes the determinant of the Jacobian matrix vanish, we get the critical value of at which bifurcation occurs. When is small (less than  M/min in Figure 3), there are only two negative eigenvalues corresponding to one equilibrium point. This equilibrium is thus stable [38].

As exceeds  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 continues increasing, only one equilibrium point is left ( greater than  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 ’s by using the Jacobian method. In Table 1, the symbol with subscripts such as and 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 and at appropriate values,  M/min and  M/min. To calculate the equilibria, we set the left hand side of (1) to zero and obtain the expression of and as functions of ’s: where and 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 ’s. The relative sensitivities are given by where the subscript . 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 s, we divide them into four groups:(1), , and : binding free energies of CI;(2), , and : binding free energies of Cro;(3) and : the cooperation between proteins CI’s if both and , or and are bound by CI dimers;(4) and : RNAP binding energies. We then study each of these four groups.

4.1. , , and

We systematically perturb the values of , , and 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.

Figures 4(a)4(c) illustrate the impact of , , and variations. We set parameter values as in Table 1 and find out all the equilibria when varying one . We find that the energy does not show Z-shaped bifurcation, while the other two parameters and show quite similar behaviors since Figures 4(a) and 4(b) are almost identically the same. The result is not surprising because and play similar roles in system dynamics. These two parameters determine the likelihood that CI dimer binds to sites and . 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 indicates the likelihood of CI dimer binding to the site . If the absolute value of is too large, it is easy for CI dimer to bind to ; thus, cI gene expression is prohibited, resulting in lysis. If the absolute value of 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 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 . Compared with other parameters, the system dynamics is more robust to variation.

Table 3 lists the stable range, bistable range, and the length of bistable range for all ten s. 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 is less than the nominal value at kcal/mol, or is less than kcal/mol, the system always has lysogenic state. They are critical bifurcation points. For , however, lysogenic state always exists if it is greater than kcal/mol. Consequently, for CI protein, a larger binding energy to and 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 and become small, thus forming Z-shaped bifurcation. Parameter variations may lead to Z-shaped bifurcation if and only if the bistable range is bounded.

When 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. , , and

The bifurcation results of , , and are shown in Figures 4(d)4(f). The equilibria behaviors with respect to and variation are similar to . represents the likelihood that Cro dimer binds to site . If Cro dimer binds to site , the expression of cro is prohibited, resulting in only one single lysogenic stable point. When or 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 is similar to the influence of and . The bistable range of is about 3.43 kcal/mol, which is comparable to those of and .

The six parameters we have shown so far are the basic binding energies. The values of , , and decide the likelihood that negative feedback (i.e., self-prohibition) happens. It is known that the sites and are used to transcript and then translate protein Cro. The site contributes to the production of protein CI. If protein CI binds to , or Cro binds to and , 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 readily occurs when the amount of protein CI is too large. On the other hand, , , and represent the likelihood of positive feedback (i.e., mutual prohibition). For instance, with large , the binding between Cro and is favored, stopping the expression of cI gene. If the value of 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. and

The parameters and 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 and . The nominal values of and are close to zero. Since protein cooperation process is an exothermic reaction, and must be less than zero. We find out that the variation of leads to Z-shaped bifurcation whereas does not. The value of affects the expression of both cro and cI genes, whereas only promotes mutual prohibition from CI to Cro. with large value prohibits the binding of RNAP to both sites and , thus both cI and cro gene expressions are closed.

Moreover, we know the length of bifurcation range for is kcal/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 is from kcal/mol to kcal/mol (Table 3). In our model, if is equal to the center value at kcal/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 kcal/mol, in comparison to the in vitro value kcal/mol used in this work. The kcal/mol difference between these two values might be caused by the looping effect [40]. Hence, the looping effect makes the system more robust to from this point of view.

4.4. and

Figure 6 illustrates the influence of and 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 is 2.81 kcal/mol, and for is only kcal/mol, which is the minimum among all ten parameters.

Notice that and are the binding free energies of RNAP to promoter () and to promoter ( and ). 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 and than any other parameters. For example, the bistable ranges of and are around 4 kcal/mol, while those of and are less than 3 kcal/mol. Parameter uncertainties of and most easily destroy bistable behavior. The concentrations of CI and Cro at equilibrium point also change significantly under or uncertainty.

4.5. Discussion

From the above analysis, we find that the binding free energies to sites and have similar effects on system dynamics. As we have shown in Table 3, the results of , , , and are similar. This is mainly because of the similar status for binding sites and , both of which control the expression of Cro protein.

Moreover, the bistable ranges have infinite length for , , , and . 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 and 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 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 Table 3, 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 , , and indicate the likelihood of protein CI binding to sites and , or the likelihood of mutual prohibition occurring. Similarly, indicates the likelihood of protein Cro binding to site 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 and . Since and 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 and , 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 and . Moreover, parameters describing mutual prohibition can form Z-shaped bifurcation (i.e., , , , and ), whereas parameters describing self-prohibition (i.e., , , and ) 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.

Appendix

Calculation of Jacobian

From (1)–(5), we may represent system Jacobian by Since and are sums of , , as (2) shows, the key part to evaluate their partial derivatives in (A.1) is and . Since is a function with respect to Cro and CI dimers, we evaluate the partial derivatives by Chain Rule: where Then we have to deal with other terms in (A.2). Note that is a fractional function in (3), so we write its numerator and denominator as The partial derivatives for with respect to and are given by The partial derivative of can be obtained as The subscript corresponds to the 40 binding configurations. The partial derivative of is exactly the sum of 40 partial derivatives. From the above derivation, we may obtain the Jacobian matrix analytically. This is quite helpful when we calculate the eigenvalues of Jacobian matrix at equilibrium points. By referencing Table 2, we may count the exponents (i.e., , , and ) in (A.4) and list them in Table 4.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research has been supported by the China National 973 Project under Grant no. 2010CB529200, the Innovation Program of Shanghai Municipal Education Commission under Grant no. 11ZZ20, the Shanghai Pujiang Program under Grant no. 11PJ1405800, NSFC under Grant no. 61174086, and Project-sponsored by SRF for ROCS SEM.