Bistability, Probability Transition Rate and First-Passage Time in an Autoactivating Positive-Feedback Loop

A hallmark of positive-feedback regulation is bistability, which gives rise to distinct cellular states with high and low expression levels, and that stochasticity in gene expression can cause random transitions between two states, yielding bimodal population distribution (Kaern et al., 2005, Nat Rev Genet 6: 451-464). In this paper, the probability transition rate and first-passage time in an autoactivating positive-feedback loop with bistability are investigated, where the gene expression is assumed to be disturbed by both additive and multiplicative external noises, the bimodality in the stochastic gene expression is due to the bistability, and the bistability determines that the potential of the Fokker-Planck equation has two potential wells. Our main goal is to illustrate how the probability transition rate and first-passage time are affected by the maximum transcriptional rate, the intensities of additive and multiplicative noises, and the correlation of additive and multiplicative noises. Our main results show that (i) the increase of the maximum transcription rate will be useful for maintaining a high gene expression level; (ii) the probability transition rate from one potential well to the other one will increase with the increase of the intensity of additive noise; (iii) the increase of multiplicative noise strength will increase the amount of probability in the left potential well; and (iv) positive (or negative) cross-correlation between additive and multiplicative noises will increase the amount of probability in the left (or right) potential well.


Introduction
Bistability arises within a wide range of biological systems from the bacteriophage l to cellular signal transduction pathways in mammalian cells [1,2]. As a fundamental behavior of biological system, bistability has been studied extensively through experiments, theoretical analysis and numerical simulations. Hasty et al. [3] considered a single network derived from bacteriophage l and constructed a two-parameter deterministic model describing the temporal evolution of the concentration of l repressor protein. They showed how additive and multiplicative external noise can be used to regulate gene expression. In the case with only additive noise, they demonstrated the utility of such control through the concentration of protein switch, whereby protein production is turned ''on'' and ''off'' by using short noise pulse. In the case with multiplicative noise, they showed that small deviations in the transcription rate can lead to large fluctuations in the production of protein. Combining theory and experiments, Isaacs et al. [4] investigated the dynamics of an isolated genetic module, an in vivo autoregulatory gene network. As predicted by their theoretical model, temperature-induced protein destabilization led to the existence of two expression states. The result of Isaacs et al. shows clearly the effects of varying the strength of feedback activation on population heterogeneity (see also [5]). Recently, Acar et al. [6] experimentally explored how switching affects population growth by using the galactose utilization network of Saccharomyces cerevisiae.
Kaern et al. [2] pointed out that a hallmark of positive-feedback regulation is bistability, which gives rise to distinct cellular states with high and low expression levels, and that stochasticity in gene expression can cause random transitions between the two states, yielding bimodal population distributions (see also [7][8][9]). So, for the positive-feedback regulation with bistability, a challenging question is how to determine probability transition rates between the two states, or how the random transitions between the two states are affected by the transcription rate and noise strength. In this paper, a simple theoretical model for an autoactivating positive-feedback loop is investigated. Our main goal is to provide a theoretical analysis for the probability transition rate and firstpassage time in the system with external noise. The paper is organized as follows. In section 2, the basic model and its bistability is presented. The analysis of the probability transition rate and first-passage time for the situations with additive and multiplicative external noises are given in section 3.

Basic model
It is well known that the simplest circuit motif able to exhibit multiple stable states is the autoactivating positive-feedback loop [5,[9][10][11], in which a single gene encodes a protein (activator), and the activator monomers bind into dimers that subsequently bind to the upstream regulatory site of the gene, activating production of the activator monomers (see Figure 1). For example, the autoactivation of CI protein by the P RM promoter of phage l [4]. The autoactivating positive-feedback loop is expected to exhibit bistability for the protein synthesis level, i.e., a higher level and a low level of protein concentration [2,9]. Let x(t) and y(t) denote the concentrations of mRNA and activator protein at time t, respectively. Then, in general the macroscopic rate equation for x(t) and y(t) can be expressed as where F (y) represents the mRNA transcription rate, which is defined as a function of activator protein concentration with dF (y)=dyw0, the parameter K denotes the translation rate, and the parameters c R and c P are the degradation rates of mRNA and protein, respectively, with c R &c P , i.e., the concentration of mRNA is a fast variable compared with the concentration of activator protein [12,13]. For the mRNA transcription rate F (y), we take it as a Hill-type function where k max is the maximum transcription rate, a H the Hill coefficient where we take a H~2 , k H the Hill constant, and k f the basal transcription rate with k f %k max [14]. In biology, these parameters mean that: (i) k max represents the mRNA transcription rate when the activator protein concentration is large enough; (ii) a H~2 implies that the activator binding processes are considered comparatively rapid and close to equilibrium, so the concentration of activator homodimer is proportional to the square of activator monomer concentration [14]; (iii) k H is the dissociation constant of activator dimer from the regulatory site; and (iv) k f is the mRNA transcription rate when the activator protein concentration is very low [15].
Notice that x(t) can be considered to be a fast variable compared with y(t) since c R &c P . Then, Eq. 1 can be reduced as i.e., the fast variable can be assumed to be at an effective equilibrium, whereas the slow variable is responsible for the dynamics of the system [5,16,17]. In mathematics, one of the most important properties of Eq. 3 is its bistability, i.e., Eq. 3 has at most three fixed points, denoted by y s 1 , y u and y s 2 with y s 1 vy u vy s 2 . For the stability of y s 1 , y s 2 and y u , it is easy to see that both y s 1 and y s 2 are locally asymptotically stable and y u is unstable since (K=c R )dF (y s i )=dyvc P for i~1,2 and (K=c R )dF (y u )=dywc P (see also [15]).
In biology, we are more interested in how the dynamic properties of Eq. 3 is affected by the maximum transcription rate k max (see also [15]). The relationship between the bistability and k max is plotted in Figure 2 (i.e. log y vs. k max for dy=dt~0) where, following Smolen et al. [15], the parameters are taken as and k H~1 0 (see also [16]) (in this paper, we keep these parameters to be fixed). Clearly, for the stability of Eq. 3, the parameter k max has two bifurcation values, denoted by k' max and k'' max , respectively, with k' max vk'' max (where k' max &6:11 h {1 and k'' max &25:10 h {1 ), i.e., the bistability exists if k max is in the interval k' max vk max vk'' max . For the situation with k max vk' max (or k max wk'' max ), the system will be monostable, i.e., the system has only one fixed point y s1 (or y s2 ) if k max vk' max (or if k max wk'' max ). On the other hand, if k max~k ' max (or k max~k '' max ), then the system will have two fixed points y s1 and y u (or y u and y s2 ), i.e., if k max exactly equal its bifurcation value, then Eq. 3 will have two fixed points, in which one is stable and the other semi-stable.
Similar to Hasty et al. [3], we here consider also how the dynamics of protein concentration is affected by the additive and multiplicative external noises. As discussed above, if the bistability exists, then in the absence of noise, the system state will evolve identically to one of the two fixed points. The presence of a noise source will lead to the fluctuation of the system state. Hasty et al. [3] pointed out that an additive noise source alters the and k H~1 0 (see also the main text). The bistability exists if k' max vk max vk'' max , where k' max &6:11 h {1 and k'' max &25:10 h {1 ). If k max vk' max (or k max wk'' max ), then the system will be monostable (see also Ref. [14]). doi:10.1371/journal.pone.0017104.g002 ''background'' protein production. This means that we need to consider the effect of a randomly varying external field on the biochemical reactions. For example, To et al. [18] provided an experiment evident to show that the change of the external noise can induce bimodality in positive transcriptional feedback loops without bistability. On the other hand, Hasty et al. [3] also pointed out that although transcription is represented by a single biochemical reaction, it is actually a complex sequence of reactions, and it is natural to assume that this part of the gene regulatory sequence is likely to be affected by fluctuations of many internal or external parameters. This implies that the transcription rate can be also considered to be a random variable.
In our model, the additive noise, denoted by j(t), alters the ''background'' protein production, and is defined as a white noise with vj(t)w~0 and vj(t)j(t')w~2D A d(t{t') where D A measures the level of additive noise strength. The multiplicative noise alters the transcription rate. We vary the transcription rate by allowing the parameter k max to vary stochastically, i.e., let k max~kmax zg(t), where g(t) is also a white noise with vg(t)w~0 and vg(t)g(t')w~2D M d(t{t') where D M measures the level of multiplicative noise strength. Here a natural question is whether the additive and multiplicative noises are statistically correlated. However we could imagine the correlation arising from the feedback regulation, i.e. the transcription rate (affected by noise) is chemically coupled to the protein concentration (also affected by noise). We define that the crosscorrelation of j(t) and where m is the cross-correlation intensity [19]. In fact, for the crosscorrelation between the additive noise j(t) and multiplicative noise g(t), we have no experimental evidence that indicate that the parameter m should be positive, or negative. We also noticed that in a previous model developed by Hasty et al. [3], the effect of the crosscorrelation between the additive and multiplicative noises on the stochastic gene expression is ignored. Thus, for the effect of m we will only provide some theoretical possibilities.
According to the above definitions about j(t) and g(t), the Langevin equation corresponding to Eq. 3 is given by where Let w(y,t) denote the probability density distribution that the concentration of activator protein exactly equals y at time t. Then, from Risken [20], the Fokker-Planck equation of w(y,t) corresponding to Eq. 4 can be given by The stationary distribution is w(y)~Ce {U FP (y) where C is the normalized constant and is called the potential of the Fokker-Planck equation.

Additive noise
In this subsection, we consider only the effect of additive noise on the bistability, and assume that D M = 0 (i.e. we here ignore the multiplicative noise). Clearly, for D M~0 , Eq. 6 can be reduced to and Eq. 7 can be rewritten as , stationary distribution w(y) has two peaks corresponding to the two stable points y s 1 and y s 2 , respectively. This also implies that the stable point y s i (i~1,2) must correspond to the local minimum of the potential U FP (y), and the unstable point y u to the local maximum of U FP (y). In physics, the local minimum of U FP (y) corresponding to y s 1 (or y s 2 ) is also called the potential well, and the local maximum of U FP (y) corresponding to y u the potential barrier [20]. Thus, for convenience we call the potential well corresponding to y s 1 (or y s 2 ) the left (or right) well. The depth of the right well is defined as d z~UFP (y u ){ U FP (y s2 ), and, similarly, the depth of the left well is d {~UFP (y u ){U FP (y s1 ). The depth of the potential well varies as the function of k max where we take the parameters K, k f , c R , c P and k H to be fixed. The bistable potential U FP (y) is plotted in Figure 3A for . This strongly implies that the depth of left well decreases with the increase of k max but the depth of right well increases with the increase of k max . The relationship between the depth of potential well and k max is plotted in Figure 3B, i.e. the system state should be more easily attracted by the left well with the increase of the maximum transcription rate. It is also easy to see that there must exist a k max value such that d {~dz (in Figure 3B, The effects of D A and k max on the stationary distribution are plotted in Figure 4. We show that the decrease of D A will increase the total probability in the left well if  Figure 4A, 4B and 4C, respectively. We noticed that this result has been used to explain how the external noise can be used to control the level of protein synthesis [2,3]. For example, Hasty et al. investigated the autoregulation of bacteriophage l repressor expression network which contains three operator sites known as OR1, OR2, and OR3 and constructed a two-parameter deterministic model describing the temporal evolution of the concentration of l repressor protein [3]. They showed how the bistable regime is enhanced with the addition of the first operator site in the promoter region through comparing two models with only the last two operator sites and the full operator regions, respectively. They also showed how external noise can be used to regulate expression through adding external additive noise or multiplicative noise using the stochastic simulations. For the case with additive noise, they demonstrated the utility of such control through the successful switch of the concentration of a protein, whereby protein production is turned ''on'' and ''off'' by using short noise pulses, where the short noise pulse means that the noise of the system is rapidly increased to a high level in a short period of time, and after this pulse, the noise is returned to its original value. When the system state is near the stable point y si (i~1,2), the steady-state statistics of the system can be given by vyw~y si and vy 2 w{vyw 2~{ D A =g'(y si ) (the proof is given in Text S1) (see also [21]). This result shows clearly that when the system state is near the stable point y si , the intensity of stochastic fluctuations around y si is proportional to the noise strength D A but inversely proportional to {g'(y si ). On the other hand, for the bistable potential function U FP (y), a more challenging question is how the system state jumps from one potential well to the other (i.e. the state is switches from one well to the other) because of the external noise.
Notice that the probability exchange between two potential wells occurs only near the potential barrier (i.e., at the unstable point y u ). Thus, the increase (or decrease) of the amount of probability in one potential well must result in the decrease (or increase) of the amount of probability in the other. Let P z (t) denote the total probability in the right well at time t, i.e., P z (t)~Ð ? yu w(y,t) dy, and, similarly, P { (t) the total probability in the left well at time t, i.e., P { (t)~Ð yu 0 w(y,t) dy (where we must have P z (t)zP { (t)~1). Then the master equations of P z (t) and P { (t) can be given by where R z is the probability transition rate from the right well to the left well, and R { the probability transition rate from the left well to the right well, which are given by respectively, (the mathematical derivation is given in Text S1) (see also [22]). Obviously, Eq. 10 has two eigenvalues, one is l 1~0 , and the other l 2~{ (R z zR { ). The first one is the eigenvalue of the Fokker-Plack equation corresponding to the stationary distribution, and the latter is the lowest non-vanishing eigenvalue where 1=l 2 measures the largest time scale of the probability transition between two potential wells. From Eq. 11, we have that: 1. For both R z and R { , we have LR z =LD A w0 and LR { =LD A w0. This means that the increase of D A will promote the probability exchange between two wells (see Figure 5A). 2. The transition rate will decrease with the increase of the well depth, i.e., LR z =Ld z v0 and ). Hence, from Figure 3B, we have also that R { will increase with the increase of k max but R z will decrease with the increase of k max . In biology, this means that the increase of k max will promote the probability transfer from the left well to the right well, or the increase of k max will be useful for maintaining a high gene expression level (see also Figure 4).

For convenience, we use ratio
to measure the relative intensity of probability exchange between two wells, i.e., if R z =R { w1, then the probability is more easily This result shows clearly that the increase of D A will promote probability exchange from the deep well to the shallow well. The relationship between k max , D A and R z =R { is plotted in Figure 5B. Clearly, the ratio R z =R { is an increasing function of In physics, the first-passage time is defined as the time at which the stochastic variable first leaves a given domain [20]. In general, the first-passage time can be used to measure the robustness of the system steady-state. For our model, let t z (t { ) denote the firstpassage time at which the system state first leaves the right (left) well across the potential barrier. Under the weak noise (i.e., D A %1), the expectations of t z and t { can be approximated as respectively, and the variances of t z and t { , denoted by V t z and V t{ , are respectively (the mathematical derivations of Eqs 13 and 14 are given in Text S1). This strongly implies that under the weak noise, both t z and t { should approximately obey the exponential distribution. Statistically, all transition events (called also the escape events) in a given direction (for example, from the left well to the right well, or from the right well to the left well) can be roughly considered to be independent of each other with a given average rate, i.e. the number of the transition events in a given time interval should be a Poisson process. Thus, the distribution of the firstpassage time should be an exponential distribution, and its scale parameter is the inverse of the probability transition rate [23]. The numerical simulation results for the statistical properties of t z and t { are given in Table 1 (The simulation algorithm is given in Text S1). It is easy to see that for different D A values (where we take k max~6 :6 h {1 ), the average of t z (t { ) (i.e. the mean first-passage time) and its standard deviation are almost same. The Monte Carlo simulations also show that under the weak noise, the relation t z =t {~( R z =R { ) {1 is true (see Table 2). All of these simulation results exactly match the theoretical predictions.

Multiplicative noise
For D M =0, we first consider the situation with m~0, i.e., j(t) and g(t) are independent of each other. According to this definition, we have that   Lw(y,t) Lt~{ and that (see Eqs. 6 and 7). In this situation, the effect of D M on U FP (y) is plotted in Figure 6A, where the parameters D A and k max are taken as D A~0 :1 and k max~6 :5861 h {1 . In subsection 3.1, we have shown that for D M~0 , both right and left wells have the same depth if k max~6 :5861 h {1 . We will select the special value in the following analysis. We can find that the depth of the right well will decrease with the increase of D M , but the change in the depth of the left well is very small. The stationary distribution and the Monte Carlo simulation corresponding to Figure 6A are plotted in Figure 6B- (the mathematical derivation is given in Text S1) (see also [20,22]). The Monte Carlo simulation results for the statistical properties of t z and t { are given in Table 3, in which we can find that not only both t z and t { should obey the exponential distribution but also t z is more sensitive for D M than t { . This implies that the increase of the multiplicative noise strength will be useful for maintaining the protein concentration at the low level. Secondly, for m=0, we are interested in how the stochastic dynamics of the system is affected by the correlation between j(t) and g(t). The effect of m on U FP (y) is plotted in Figure 7A where D A~0 :1, D M~0 :1 and k max~6 :5861 h {1 , in which we can also find that the depth of the right well is more sensitive for the change of m than the depth of the left well, and that if m is negative, then the depth of the right well will increase with the increase of m j j (i.e. absolute value of m), and, conversely, if m is positive, the depth of the right well will decrease with the increase of m. The effect of m on the stationary distribution and the Monte Carlo simulation corresponding to Figure 7A are plotted in Figure 7B-7D. It shows clearly that the positive correlation (mw0) will increase the amount of probability in the left well, and, conversely, the negative correlation (mv0) will increase the amount of probability in the right well.
For the effect of m on the first-passage time, the Monte Carlo simulation results are showed in Figure 8, in which both t z and t { will decrease with the increase of m, and the change rate of t z is obviously larger than that of t { . We can also notice that the ratio t z =t { will decrease with the increase of m (see Figure 8B). This also implies that the positive (or negative) correlation between j(t) and g(t) will promote the probability transition from the right (left) well to the left (right) well (see also Figure 7B-7D). Clearly, the dependence of the ratio t z =t { on m reflects how the stationary distribution is influenced by the cross-correlation between additive and multiplicative noises, or, theoretically, the cross-correlation between j(t) and g(t) should be also used to control the protein synthesis.

Discussion
In this paper, the probability transition rate and first-passage time in an autoactivating positive-feedback loop with bistability are investigated. In our model, similar to Hasty et al. [3], the gene expression is assumed to be disturbed by both additive and multiplicative external noises, and the bimodality in the stochastic gene expression is due to the bistatility. The bistability of the deterministic dynamics Eq. 3 implies that the potential of the Fokker-Planck equation Eq. 6 has two potential wells, which correspond to the two stable points of Eq. 3, respectively, and that the stationary solution (i.e. stationary distribution) of the Fokker-Planck equation is a bimodal distribution. For our main goal, we are interested in how the system state jumps from one potential well to the other because of the external noise. In subsection 3.1, for the situation with only additive noise, our main results show that (i) both probability transition rates R z and R { will increase with the increase of D A ; (ii) R { will increase with the increase of k max but R z will decrease with the increase of k max , i.e., the increase of k max will be useful for maintaining a high gene expression level; (iii) the ratio R z =R { measures the relative intensity of probability exchange between two potential wells, and there is a critical value of k max (which is 6:5861 h {1 in our case) such that the ratio R z =R { is an increasing function (or a decreasing function) of D A if k max is larger (or smaller) than the critical value; and (iv) for both firstpassage times t z and t { , if D A %1 (i.e. the additive noise is weak), then they obey the exponential distribution with expectations vt z w&1=R z and vt { w&1=R { , respectively. In subsection 3.2, for the situation with both additive and multiplicative noises, we show that (i) for m~0 (i.e. the additive noise and multiplicative noise are independent of each other), the increase of D M will be useful for maintaining the protein concentration at the low level; and (ii) for m=0, the positive (or negative) correlation between additive and multiplicative noises will promote the probability transition from the right (left) well to the left (right) well, i.e., the positive correlation will increase the amount of probability in the left potential well, and, conversely, the negative correlation will increase the amount of probability in the right potential well.
However, our results may provide some theoretical intuitions for real gene expression. For example, Hasty et al. [3] investigated the autoregulation of l repressor expression in the lysis-lysogeny pathway in the l virus, and they showed why the additive and multiplicative noises can be used to regulate expression (or to control the protein production). Our results further show how the probability transition between two stable states (i.e. two levels of protein synthesis) is affected, and why the cross-correlation between the additive and multiplicative external noises can be also used to regulate expression. Finally, we would like to say that the further analysis incorporating more realistic dynamics should be carried out in future studies.

Supporting Information
Text S1 Steady-state statistics, Probability transition rate, Firstpassage time, Method for stochastic simulation. (PDF)