Analytical study of an exclusive genetic switch

Following the stochastic dynamical processes described above the system of master equations that defines the model can be written as follows:

$$\begin{eqnarray} &&\fl \frac{\partial P_0}{\partial t}(N_1,N_2) = g[P_0(N_1-1,N_2)+P_0(N_1,N_2-1)-2P_0(N_1,N_2)]\nonumber \\ &&+\,d[(N_1+1)P_0(N_1+1,N_2)+(N_2+1)P_0(N_1,N_2+1)\nonumber \\ &&-\,(N_1+N_2)P_0(N_1,N_2)]-b(N_1+N_2)P_0(N_1,N_2)\nonumber \\ &&+\,u[P_1(N_1-1,N_2)+P_2(N_1,N_2-1)] \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\fl \frac{\partial P_1}{\partial t}(N_1,N_2) = g[P_1(N_1-1,N_2)-P_1(N_1,N_2)]+d[(N_1+1)P_1(N_1+1,N_2)\nonumber \\ &&+\,(N_2+1)P_1(N_1,N_2+1)-(N_1+N_2)P_1(N_1,N_2)]\nonumber \\ &&+\,b(N_1+1)P_0(N_1+1,N_2)-uP_1(N_1,N_2) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\fl \frac{\partial P_2}{\partial t}(N_1,N_2) =g[P_2(N_1,N_2-1)-P_2(N_1,N_2)]+d[(N_1+1)P_2(N_1+1,N_2)\nonumber \\ &&+\,(N_2+1)P_2(N_1,N_2+1)-(N_1+N_2)P_2(N_1,N_2)]\nonumber\\ &&+\, b(N_2+1)P_0(N_1,N_2+1)-uP_2(N_1,N_2) \end{eqnarray} \tag{ 3 }$$

where g is the generation rate of a protein (when the generation is not suppressed by the switch state), d is the degeneration rate of a single protein, b is the binding rate of a single protein and u is the unbinding rate of the bound protein. The whole problem is clearly symmetric with respect to the variables 1 and 2. Also note that, the degeneration term is the same for the three equations, since it does not depend on the value of S.

In order to illustrate the qualitative behaviour, we first present stochastic simulations of the Exclusive Switch. These simulations are performed with a Gillespie algorithm [25, 26], in which the reactions described in the model are given a certain probability, depending on the state of the system and the value of the parameters. These probabilities are used to determine stochastically which reaction is going to happen next, and when it will happen. The time that the system spends in a given state N₁, N₂, S is normalized to obtain the probability distributions. (An alternative approach would be to solve the master equation numerically by using e.g. the method of finite state projection [27].)

The reference values for the simulations are the typical values for bacteria such as Escherichia coli (E. coli) (in s⁻¹) [1, 7]

$$\begin{equation} g=0.05, \qquad d=0.005, \qquad b=0.1, \qquad u= 0.005 . \end{equation} \tag{ 4 }$$

In the system there is a clear symmetry between the two protein species, since they undergo the same microscopic reactions. However, at any given time the system is typically dominated by one of the proteins. The reason for this is that, once a protein, e.g. of type X₁, binds to the promoter site, proteins X₂ start disappearing, while the number of X₁ fluctuates around a steady value. That means that when the bound protein unbinds (as it will do eventually due to stochasticity of the system), it is much more probable for proteins X₁ to bind to the promoter site again, since there are more of them. At the same time it is more difficult for proteins X₂ to bind to the promoter site. However the X₂ will not disappear permanently from the system (there is no absorbing state), and will be produced again the moment the bound protein X₁ unbinds. Thus, with a small probability, proteins X₂ will be able to bind again to the promoter site, and become the dominant species, as proteins X₁ start to degenerate. This means that there are two long-lived symmetry-related states, in which one species is much more abundant than the other.

With regard to the probability distribution P(N₁, N₂), this bistability is translated into two peaks, concentrated around the axes, i.e., where one of the protein numbers is almost zero. This is illustrated in a contour plot in figure 1 for u = 0.05. However, this bistability depends strongly on the value of the parameters: the bigger the value of u, the more irrelevant is the switch state for the dynamics of the protein, and the less important is the bistability we have described. For example in figure 1, when g, b and d are kept constant and u is increased, the peaks move together and eventually merge at the value of u = 0.145 ± 0.005 (see table 1 for the value u at which the transition happens, for other values of g, b, d). Thus there is an apparent transition from a distribution with two symmetry-related peaks to a distribution with one symmetric peak. We wish to study the nature of this transition, i.e. is there an underlying transition at a finite value of u where the system changes from bistable behaviour to symmetric behaviour, or is the transition simply due to two peaks coming closer together and no longer being resolved?

Zoom In Zoom Out Reset image size

Although the whole probability distribution is always symmetric, distributions P₁ and P₂ will be a priori asymmetric, since they describe the probability of the number of proteins when a protein 1 or 2 are bound, which are not symmetric situations. Let us now define r_A and r_B as the probability masses of P₁ on either sides of the diagonal N₁ = N₂:

$$\begin{equation} \begin{array}{l} r_A =\dsty\sum _{N_1>N_2}P_1(N_1,N_2)+\dsty\frac{1}{2} \dsty\sum _{N_1=N_2}P_1(N_1,N_2) \\\ms r_B =\dsty\sum _{N_1<N_2}P_1(N_1,N_2)+\dsty\frac{1}{2} \dsty\sum _{N_1=N_2}P_1(N_1,N_2). \end{array} \end{equation} \tag{ 5 }$$

We can now study how r_A, r_B change with u, and whether there is a clear transition between the situation in which they are different, and the one in which they are equal to each other (if there is any). Figure 2 shows that these two quantities approach each other in a continuous way, and they are equal only when u → ∞, that is, when the only possible state of the switch is S = 0 and it no longer has any effect on the protein dynamics. Therefore, the probability distributions P₁, P₂ and hence r_A, r_B tend to zero, all the probability being concentrated in the distribution P₀(N₁, N₂), which is always symmetric.

Zoom In Zoom Out Reset image size

The marginal distribution P₁ appears to deform continuously into a symmetric distribution when u → ∞. We conclude that (at least for these parameter values) P₁(N₁, N₂) remains asymmetric for u < ∞, and as a consequence so does P₂, and there is no transition to a symmetric state in this marginal distribution. As figure 3 shows, P₁ has only one peak for different values of u, and even if it becomes smaller as u increases, it does not become symmetric at any point.

Zoom In Zoom Out Reset image size

We deduce that, even though P(N₁, N₂) appears to become unimodal at some finite value of u (see figure 1), there is no transition between symmetric and asymmetric regimes, since P₁ and P₂ always remain asymmetric. Thus, the bistability of the switch is always present, with the asymmetry in distributions P₁, P₂ decreasing as the switch state becomes less important, that is, as u increases.

We now consider the system in the limit u, b → 0 with the binding constant

$$\begin{equation} k= \frac{b}{u} \end{equation} \tag{ 6 }$$

held fixed. In this limit the system will equilibrate between binding/unbinding events. When the switch is in state 1 the number of type 2 proteins decays to zero and the number of type 1 is given by Poissonian statistics for the generation/degeneration processes:

$$\begin{equation} P_1(N_1,N_2) = r_1 \frac{1}{N_1!}\left(\frac{g}{d}\right)^{N_1} {\rm e}^{-g/d} \delta _{N_2,0} \end{equation} \tag{ 7 }$$

where r₁ is the probability of the switch being in state 1 (r₁ = ∑^∞_{N1, N2 = 0}P₁(N₁, N₂)). Similarly, in state 2 the number of type 2 is given by Poissonian statistics

$$\begin{equation} P_2(N_1,N_2) = r_2 \frac{1}{N_2!}\left(\frac{g}{d}\right)^{N_2} {\rm e}^{-g/d} \delta _{N_1,0} \end{equation} \tag{ 8 }$$

and in state 0 both N₁ and N₂ are given by Poissonian statistics

$$\begin{equation} P_0(N_1,N_2) = r_0 \frac{1}{N_1!N_2!}\left(\frac{g}{d}\right)^{N_1+N_2} {\rm e}^{-2g/d}. \end{equation} \tag{ 9 }$$

In order to fix the probabilities r₀, r₁, r₂ we consider the master equation for r₁ which can be obtained from summing (2) over the variables N₁, N₂:

$$\begin{equation} \dot{r_1} = -u r_1 + b \langle N_1\rangle_0 \end{equation} \tag{ 10 }$$

where the quantity 〈N₁〉₀ is the mean value of N₁ given the system is in switch state S = 0, multiplied by the probability of being in switch state S = 0, i.e., 〈N₁〉₀ = ∑^∞_{N1, N2 = 0}N₁P₁(N₁, N₂). In the stationary state we have 〈N₁〉₀ = r₀g/d and (10) becomes

$$\begin{equation} 0 = -u r_1 + \frac{bg}{d} r_0. \end{equation} \tag{ 11 }$$

A similar equation holds for r₂ and the normalisation of probability yields

$$\begin{equation} r_0 = \frac{1}{1+ 2kg/d}\qquad r_1=r_2 = \frac{kg/d}{1+ 2kg/d}. \end{equation} \tag{ 12 }$$

Thus we see that the system has three long-lived states: when S = 1 the system is dominated by type 1 proteins; when S = 2 the system is dominated by type 2 proteins, and when S = 0 the system is in a symmetric state.

When k → ∞ the symmetric state has zero weight (r₀ → 0) and the system is either in the S = 1 state or the S = 2 state. The switching time between the two states is expected to diverge in this limit therefore the system exhibits symmetry breaking. Note that the transition from the S = 1 state to the S = 2 state is through the symmetric S = 0 state.

In this limit the unbinding and binding events happen on a faster time scale than the growth and degradation of proteins. (This may be a physically relevant limit since it is often assumed that binding-unbinding processes are faster than other processes [9].) Therefore, the switch state becomes decoupled from the numbers of protein and the probabilities obey

$$\begin{equation} P_S(N_1,N_2) = P(N_1,N_2) r_S(N_1,N_2) \end{equation} \tag{ 13 }$$

where P(N₁, N₂) is the probability of N₁, N₂ proteins being in the system regardless of switch state, and r_S(N₁, N₂) is the probability for the switch to be in state S, given that the numbers of proteins are N₁, N₂. That is, for given number of proteins N₁ and N₂ (which include any bound protein) the switch probabilities equilibrate and obey

$$\begin{equation} u r_1(N_1,N_2) = b N_1 r_0(N_1,N_2) \qquad u r_2(N_1,N_2) = b N_2 r_0(N_1,N_2) \end{equation} \tag{ 14 }$$

which may be solved to obtain

$$\begin{equation} r_0 = \left[ 1+ k(N_1+N_2)\right]^{-1} \qquad r_1= k N_1 r_0 \qquad r_2= k N_2 r_0 \end{equation} \tag{ 15 }$$

where the binding constant k is given by (6).

A master equation for the evolution of P(N₁, N₂) on the slower timescale on which generation and degeneration events occur may then be written down:

$$\begin{eqnarray} \fl \frac{\partial P}{\partial t}(N_1,N_2) = & g[1-r_2(N_1-1,N_2) ]P(N_1-1,N_2) -d[N_1-r_1(N_1,N_2)]P(N_1,N_2) \nonumber\\ &-g[1-r_2(N_1,N_2)P(N_1,N_2) ] +d[N_1+1-r_1(N_1+1,N_2)]P(N_1+1,N_2)\nonumber \\ &+ g[1-r_1(N_1,N_2-1) ]P(N_1,N_2-1) -d[N_2-r_2(N_1,N_2)]P(N_1,N_2)\nonumber \\ &-g[1-r_1(N_1,N_2) ]P(N_1,N_2) +d[N_2+1-r_2(N_1,N_2+1)]P(N_1,N_2+1). \nonumber\\ \end{eqnarray} \tag{ 16 }$$

The terms with coefficient g in (16) represent generation of N_i when the switch is not in switch state i. The terms with coefficient d in (16) represent reduction of N_i with rate dN_i when the switch is not state i and reduction with rate d(N_i − 1) when the switch is in state i (this is because the bound protein does not degrade). The stationary state of (16) obeys detailed balance with respect to generation and degradation of N₁ and N₂ individually, thus

$$\begin{eqnarray} &&\fl P(N_1,N_2)= P(N_1-1,N_2)\frac{g}{d}\frac{[1-r_2(N_1-1,N_2) ]}{[N_1-r_1(N_1,N_2)]} \nonumber\\ &&\hphantom{\fl P(N_1,N_2)}= P(N_1-1,N_2)\frac{g}{d}\frac{[1+k(N_1-1)]}{[1+k(N_1+N_2-1)]} \frac{1}{N_1} \frac{[1+k(N_1+N_2)]}{[1+k(N_1+N_2-1)]} \end{eqnarray} \tag{ 17 }$$

with a similar equation holding for detailed balance in N₂. Then these equations may be iterated to obtain

$$\begin{eqnarray} &&\fl P(N_1,N_2) = \left(\frac{g}{d}\right)^{N_1} \frac{(1+k(N_1+N_2))}{(1+kN_2)} \frac{1}{N_1!} \prod _{n_1=0}^{N_1-1}\frac{(1+k n_1)}{1+k(N_2+ n_1)} P(0,N_2)\nonumber \\ &&\hphantom{\fl P(N_1,N_2) }= \left(\frac{g}{d}\right)^{N_1+N_2} \frac{(1+k(N_1+N_2))}{N_1! N_2!} \frac{\prod _{n_1=0}^{N_1-1}(1+k n_1) \prod _{n_2=0}^{N_2-1}(1+k n_2)}{\prod _{n=0}^{N_1+N_2-1}(1+k n)} P(0,0). \end{eqnarray} \tag{ 18 }$$

The constant P(0, 0) will be determined by normalisation of the sum of probabilities to one.

The expression (18) is a single-peaked distribution which is symmetric in N₁ and N₂. To see this one can identify the stationary points of P(N₁, N₂) through the conditions P(N₁ + 1, N₂) = P(N₁, N₂) and P(N₁, N₂ + 1) = P(N₁, N₂) which yield

$$\begin{equation} \begin{array}{l} \dsty\frac{g}{d} \dsty\frac{1}{(N_1+1)} \dsty\frac{(1+kN_1)}{(1+k(N_1+N_2))} \dsty\frac{(1+k(N_1+N_2+1))}{(1+k(N_1+N_2))} =1 \\\ms \dsty\frac{g}{d} \dsty\frac{1}{(N_2+1)} \dsty\frac{(1+kN_2)}{(1+k(N_1+N_2))} \dsty\frac{(1+k(N_1+N_2+1))}{(1+k(N_1+N_2))} =1. \end{array} \end{equation} \tag{ 19 }$$

The solution of these equations is N₁ = N₂ = N/2 which means

$$\begin{equation} \frac{g}{d}\frac{1}{\left(\frac{N}{2}+1\right)} \frac{\left(1+k\frac{N}{2}\right)}{(1+kN)} \frac{(1+k(N+1))}{(1+kN)} =1 \end{equation} \tag{ 20 }$$

and when N is large,

$$\begin{equation} \frac{g}{d} \frac{\left(2+kN\right)}{N(1+kN)} \simeq 1. \end{equation} \tag{ 21 }$$

The remaining quadratic

$$\begin{equation} kN^2 + \left[1- k\frac{g}{d}\right]N - 2 \frac{g}{d} =0 \end{equation} \tag{ 22 }$$

has a single positive root. Thus, in the u, b → ∞ limit the system has reached a symmetric state with the probability distribution peaked at N₁ = N₂ = N/2 given by (22).

We start from the exact master equations (1)–(3) for the evolution of probabilities P_S(N₁, N₂). The zeroth moments of N_i are the probabilities r_S, i.e.

$$\begin{equation} r_S = \sum _{N_1=0,N_2=0}^{\infty } P_S(N_1,N_2)\qquad \mbox{for} \quad S= 0,1,2. \end{equation} \tag{ 23 }$$

We now define the first moments of 〈N_i〉_S of N_i as follows:

$$\begin{equation} \langle N_i\rangle_S = \sum _{N_1=0,N_2=0}^{\infty } N_i P_S(N_1,N_2)\qquad \mbox{for}\quad i=1,2\qquad S= 0,1,2 \end{equation} \tag{ 24 }$$

and the second moments

$$\begin{equation} \fl \langle N_iN_j\rangle_S= \sum _{N_1=0,N_2=0}^{\infty } N_i N_j P_S(N_1,N_2)\qquad \mbox{for}\quad i,j=1,2\qquad S= 0,1,2. \end{equation} \tag{ 25 }$$

Summing (1) gives

$$\begin{eqnarray} &&\fl \frac{\partial r_0}{\partial t} = -b(\langle N_1\rangle_0+ \langle N_2\rangle_0) +u[r_1+r_2] \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} &&\fl \frac{\partial \langle N_1\rangle_0}{\partial t} = gr_0 -d\langle N_1\rangle _0 -b[\langle N_1N_1\rangle _0+\langle N_1N_2\rangle _0] + u[\langle N_1\rangle_1+\langle N_1\rangle_2+r_1] . \end{eqnarray} \tag{ 27 }$$

Note that the physical meaning of e.g. 〈N_i〉_S is the probability of being in switch state S (r_S) multiplied by the mean number of type i, given that the switch is in state S.

Similarly, summing (2) gives

$$\begin{eqnarray} \frac{\partial r_1}{\partial t} &=& +b\langle N_1\rangle_0 -ur_1 \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} \frac{\partial \langle N_1\rangle_1}{\partial t} &=&gr_1 -d\langle N_1\rangle_1 +b\langle N_1(N_1-1)\rangle_0 - u\langle N_1\rangle_1 \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&\frac{\partial \langle N_2\rangle_1}{\partial t} = -d\langle N_2\rangle_1 +b\langle N_1N_2\rangle_0 - u\langle N_2\rangle_1. \end{eqnarray} \tag{ 30 }$$

We now invoke symmetry between switch state 1 and 2:

$$\begin{equation} \fl r_1=r_2 = (1-r_0)/2\quad \langle N_1\rangle _0 =\langle N_2\rangle_0 \quad \langle N_1\rangle_1 =\langle N_2\rangle_2 \quad \langle N_1\rangle_2 =\langle N_2\rangle_1. \end{equation} \tag{ 31 }$$

Then the exact steady-state versions of equations (26)–(30) read

$$\begin{eqnarray} &&r_1= \frac{b}{u}\langle N_1\rangle_0 \qquad r_0 = 1-2 r_1 \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} &&b[\langle N_1N_1\rangle_0+\langle N_1N_2\rangle_0] =gr_0 -d\langle N_1\rangle_0 + u[ \langle N_1\rangle_1+ \langle N_1\rangle_2+r_1] \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} &&b\langle N_1N_1\rangle_0 =-gr_1 +(d+u)\langle N_1\rangle_1 + b\langle N_1\rangle_0 \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} &&b\langle N_1N_2\rangle_0 =(d+u)\langle N_2\rangle_1. \end{eqnarray} \tag{ 35 }$$

Note that if we sum (33)–(35) we obtain the exact relation

$$\begin{equation} d( \langle N_1\rangle_0+\langle N_1\rangle_1+\langle N_1\rangle_2) =g(1-r_1)=g(1-r_2) \end{equation} \tag{ 36 }$$

which simply gives the overall birth/death balance for N₁. Also note that (34) and (35) give exact relations between the second moments and first moments. However, to actually evaluate these moments one would have to consider equations for higher moments, leading to a hierarchy of equations.

We now make a mean-field approximation that expresses second moments in terms of first moments:

$$\begin{eqnarray} \langle N_1N_1\rangle_0 &=& \frac{ \langle N_1\rangle_0\ \langle N_1\rangle_0}{r_0} +\langle N_1\rangle_0 \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} \langle N_1N_2\rangle_0 &=& \frac{ \langle N_1\rangle_0\ \langle N_2\rangle_0}{r_0}=\frac{ (\langle N_1\rangle_0)^2}{r_0}. \end{eqnarray} \tag{ 38 }$$

Note that a symmetry condition from (31) has been explicitly used in the last equation of (38). The first relation (37) comes from the assumption that N₁ has a Poisson distribution when the switch is in the 0 state. This means that the second moment is equal to the square of the mean plus the mean itself. However, there is an important factor r₀, which comes from the fact that 〈N₁N₁〉₀/r₀ is the mean square value of N₁ given that the switch is in the 0 state and 〈N₁〉₀/r₀ is the mean value of N₁ given that the switch is in the 0 state. The Poisson approximation is in fact exact in the limit where u, b tend to zero (see section 2.3).

The second relation (38) is a simple factorization scheme which ignores correlations between the values of N₁ and N₂ when the switch state is 0.

Using this approximation scheme (34) becomes

$$\begin{equation} \langle N_1\rangle_1 =\frac{b}{d+u} \frac{(\langle N_1\rangle_0)^2}{r_0} + \frac{gr_1}{d+u} =\frac{br_0}{d+u}\left[ \left(\frac{\langle N_1\rangle_0}{r_0}\right)^2 + \frac{g}{u} \left(\frac{\langle N_1\rangle_0}{r_0}\right) \right] \end{equation} \tag{ 39 }$$

and (35) becomes

$$\begin{equation} \langle N_1\rangle_2 = \frac{br_0}{d+u} \left(\frac{\langle N_1\rangle_0}{r_0}\right)^2. \end{equation} \tag{ 40 }$$

Using expressions (39) and (40) in (33), combined with the previous approximations (37) and (38) yields the following quadratic equation for 〈N₁〉₀/r₀:

$$\begin{equation} \frac{2bd}{d+u}\left(\frac{\langle N_1\rangle_0}{r_0}\right)^2+\left[d-\frac{bg}{d+u}\right] \left(\frac{\langle N_1\rangle_0}{r_0}\right) -g=0. \end{equation} \tag{ 41 }$$

One must take the positive root of this quadratic which yields

$$\begin{equation} \fl \frac{\langle N_1\rangle_0}{r_0} = \frac{1}{4bd}[ bg -d(d+u) + ( (d(d+u)-bg)^2 + 8bdg(d+u))^{1/2}]. \end{equation} \tag{ 42 }$$

Then using (32) one obtains

$$\begin{eqnarray} r_0 &=& \left[ 1 + \frac{2b}{u}\frac{\langle N_1\rangle_0}{r_0}\right]^{-1} \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} \langle N_1\rangle_0 &=& \left(\frac{\langle N_1\rangle_0}{r_0}\right)\left[ 1 + \frac{2b}{u}\frac{\langle N_1\rangle_0}{r_0}\right]^{-1}. \end{eqnarray} \tag{ 44 }$$

One may check the limits of section 2 from the quadratic (41). In the limit b, u → 0 with k = b/u one obtains 〈N₁〉₀/r₀ = g/d and $r_0 = \big[1 + \frac{2kg}{d}\big]^{-1}$ in agreement with section 2.3 where it is shown that N₁ follows a Poisson distribution with mean g/d when the switch is in state S = 1 or S = 0.

In the limit b, u → ∞ with k = b/u fixed the quadratic (41) reduces to

$$\begin{equation} \left(\frac{\langle N_1\rangle_0}{r_0}\right)^2 2kd + \frac{\langle N_1\rangle_0}{r_0} (d- kg) -g =0. \end{equation} \tag{ 45 }$$

This quadratic for $\frac{\langle N_1\rangle_0}{r_0}$ is the same as the quadratic (22) for the value of N = 2N₁ that maximises P(N₁, N₂) in the exact solution of section 2.4.

The mean field theory we have developed can be compared with the simulations by studying the zeroth and first order moments of the probability distributions, i.e. r₀, 〈N₁〉₀, 〈N₁〉₁, 〈N₁〉₂. (The probabilities r₁ and r₂ may automatically be obtained from r₀, since $r_1=r_2=\frac{1-r_0}{2}$.)

Different values of the mentioned quantities are given in table 2, where the parameters of the models have different values. The reference for this table is the set of E. coli values, and only the parameters that are changed are written. For the simulations performed, r₀ is always in good agreement with the mean field theory approximation, and so are r₁ and r₂. Also, 〈N₁〉₀ is in quite good agreement, too.

Table 2. Results from different quantities from simulations and mean field theory approach. The values of r₀ and 〈N₁〉₀ predicted by the mean field theory are always in good agreement with the simulations. 〈N₁〉₁ and 〈N₁〉₂ predictions are quite far from the simulation results, except on the limit u, b → 0, where our approximation is exact.

	E. coli	MFT	u = 0.05	MFT
r0	(4.9186 ± 0.0007) × 10−3	4.92705 × 10−3	$\hspace{6pt}(4.5263\pm 0.0008)\times 10^{-2}$	4.54635 × 10−2
〈N1〉0	(2.489 ± 0.005) × 10−2	2.48768 × 10−2	$\hspace{6pt}0.23862\pm 0.00014$	0.238634
〈N1〉1	4.942 ± 0.008	3.74372	$\hspace{6pt}4.113\pm 0.003$	2.71128
〈N1〉2	(5.908 ± 0.007) × 10−2	1.25604	$\hspace{6pt}0.8734\pm 0.0005$	2.2774
	u = 50, b = 1000	MFT	$\hspace{6pt} u=5\times 10^{-8}, b=10^{-6}$	MFT
r0	(4.941 ± 0.004) × 10−3	4.95073 × 10−3	$\hspace{6pt}(2.48 \pm 0.03)\times 10^{-3}$	2.49872 × 10−3
〈N1〉0	(2.55 ± 0.06) × 10−2	2.48762 × 10−2	$\hspace{6pt}(2.477\pm 0.024)\times 10^{-2}$	2.49375 × 10−2
〈N1〉1	9.89 ± 0.11	2.50019	$\hspace{6pt}4.92\pm 0.05$	4.98751
〈N1〉2	0.227 ± 0.013	2.49969	$\hspace{6pt}(5.012\pm 0.007)\times 10^{-5}$	4.97754 × 10−5

〈N₁〉₁ and 〈N₁〉₂ are different in the mean field theory, which is an improvement over simpler approximations where they have the same value. However, the values of 〈N₁〉₁ and 〈N₁〉₂ are rather different from the simulations, and this comes from ignoring higher correlations of the numbers of proteins and the state of the switch. Only when u, b → 0 are the values in close agreement with the simulation values as expected in this limit where the mean field theory is exact.

Mean field theory can, in principle, be improved by considering higher order moments and correlations. However, the algebra soon gets quite complicated.

We begin by considering the formal solution of the master equation system (1)–(3). To transform the system of equations into a system of partial differential equations, we take the generating function of the different probability distributions:

$$\begin{equation} K_{S}(z_1,z_2)=\sum ^{\infty }_{N_1=0}\sum ^{\infty }_{N_2=0}z_1^{N_1}z_2^{N_2}P_{S}(N_1,N_2), \end{equation} \tag{ 46 }$$

where S = 0, 1, 2. In this way we obtain a system of linear partial differential equations with non-constant coefficients:

$$\begin{eqnarray} &&\fl g(z_1+z_2-2)K_0+ [d-(d+b)z_1] \frac{\partial K_0}{\partial z_1} +[d-(d+b)z_2]\frac{\partial K_0}{\partial z_2} +uz_1K_1+uz_2K_2 = 0 \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} &&\fl [g(z_1-1)-u]K_1+d(1-z_1)\frac{\partial K_1}{\partial z_1} +d(1-z_2)\frac{\partial K_1}{\partial z_2} +b\frac{\partial K_0}{\partial z_1} = 0 \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} &&\fl [g(z_2-1)-u]K_2+d(1-z_1)\frac{\partial K_2}{\partial z_1}+d(1-z_2) \frac{\partial K_2}{\partial z_2}+b\frac{\partial K_0}{\partial z_2} =0, \end{eqnarray} \tag{ 49 }$$

where the right-hand side terms have been set to 0, since the stationary probabilities P_S(N₁, N₂) are the main quantities to be determined in this paper.

The second and the third equations of the system work in a completely analogous way, because of the symmetry of species 1 and 2, so it will be enough to deal with the first two equations. We also note the symmetries K₂(z₁, z₂) = K₁(z₂, z₁) and K₀(z₁, z₂) = K₀(z₂, z₁).

In appendix A we give a formal solution to the system (47)–(49). However, in practice it is not clear how to actually compute e.g. probability distributions from this solution. In order to do this we develop instead a perturbative approach.

In this section we develop a perturbative approach to the problem of finding the exact stationary state. To do so we require a suitable small parameter of the model which we choose to be u, the unbinding parameter. In the u → 0, the exact solution is simple: if, for example, one protein of type 1 is bound, the proteins of this kind will obey the usual Poisson distribution regulated by death and birth terms, while the number of proteins of type 2 will just decay to 0. This limit is the starting point for a perturbative solution, wherein the probability distribution will be expanded in a power series of u:

$$\begin{equation} P_S = \sum _{n=0}^{\infty } u^n P_S^{(n)}. \end{equation} \tag{ 50 }$$

Owing to the symmetry of the system P₂(N₁, N₂) = P₁(N₂, N₁) we need only consider P₀ and P₁.

Writing out the expansion explicitly we have

$$\begin{equation} \begin{array}{l} P_1 = P_1^{(0)}+ uP_1^{(1)} \dots \\\ms P_0 =uP_0^{(1)}+ u^2 P_0^{(2)} \ldots. \end{array} \end{equation} \tag{ 51 }$$

Note that the constant term P⁽⁰⁾₀ = 0 since P₀ = 0 in the limit of no unbinding.

This approach also makes sense when the typical E. coli values for the parameters are considered [7]:

$$\begin{equation} g=0.05 \qquad d=0.005 \qquad b=0.1 \qquad u=0.005 \end{equation} \tag{ 52 }$$

u is, along with d, the smallest of the parameters. An expansion in 1/b (appendix B) could also be developed, but we find the expansion in u more convenient.

In the zeroth order of the u expansion of the stationary master equation (2) (with left-hand side set to zero) we find

$$\begin{eqnarray} &&\fl 0=g\big[P_1^0(N_1-1,N_2)-P_1^0(N_1,N_2)\big]+d\big[(N_1+1)P_1^0(N_1+1,N_2)\nonumber \\ &&+\,(N_2+1)P_1^0(N_1,N_2+1)-(N_1+N_2)P_1^0(N_1,N_2)]. \end{eqnarray} \tag{ 53 }$$

We define the generating function

$$\begin{equation} K_1^{(0)}(z_1,z_2) = \sum ^{\infty }_{N_1=0}\sum ^{\infty }_{N_2=0} z_1^{N_1}z_2^{N_2}P_{1}^{(0)}(N_1,N_2) \end{equation} \tag{ 54 }$$

which obeys

$$\begin{equation} 0=g(z_1-1)K^{(0)}_1+ d(1-z_1)\frac{\partial K^{(0)}_1}{\partial z_1} + d(1-z_2)\frac{\partial K^{(0)}_1}{\partial z_2} \end{equation} \tag{ 55 }$$

the solution of which is independent of z₂:

$$\begin{equation} K^{(0)}_1 = c_1\exp \left(\frac{g}{d}z_1\right), \end{equation} \tag{ 56 }$$

where c₁ is a constant to be determined. Analogously $K_2^{(0)}(z)=c_2\exp \big(\frac{g}{d}z_2\big)$. Applying the normalization condition K⁽⁰⁾₁(1) + K⁽⁰⁾₂(1) = 1 and the symmetry consideration c₁ = c₂ leads to

$$\begin{equation} K_1^{(0)}=\frac{1}{2}\exp \left(-\frac{g}{d}\right)\ \exp \left(\frac{g}{d}z_1\right). \end{equation} \tag{ 57 }$$

Expanding as a power series in z₁, z₂ yields

$$\begin{equation} P_1^{(0)}(N_1,N_2)=\frac{1}{2}\frac{\exp \big({-}\frac{g}{d}\big)}{N_1!}\left(\frac{g}{d}\right)^{N_1}\delta _{N_2,0}. \end{equation} \tag{ 58 }$$

This is a Poisson distribution for N₁ with mean g/d, with N₂ fixed to be zero. The normalisation factor 1/2 is so that P⁽⁰⁾(N₁, N₂) = P⁽⁰⁾₁(N₁, N₂) + P⁽⁰⁾₂(N₁, N₂) is normalised to unity.

We substitute the expansion (50) into the stationary master system (1)–(3) with time derivatives set equal to zero. Arranging orders of u the equations may be written as

$$\begin{equation} \mathcal {L}_S P_S^{(n)}(N_1,N_2) = -f_S^{(n)}(N_1,N_2) \end{equation} \tag{ 59 }$$

for S = 0, 1 where the action of the linear operators $\mathcal {L}_S$ is

$$\begin{eqnarray} &&\fl \mathcal {L}_0 P_0^{(n)}(N_1,N_2) =g\big[P_0^{(n)}(N_1-1,N_2)+P_0^{(n)}(N_1,N_2-1)-2P_0^{(n)}(N_1,N_2)\big]\nonumber \\ &&+\,d\big[(N_1+1)P_0^{(n)}(N_1+1,N_2)+(N_2+1)P_0^{(n)}(N_1,N_2+1) \nonumber \\ &&-\,(N_1+N_2)P_0^{(n)}(N_1,N_2)\big]-b(N_1+N_2)P_0^{(n)}(N_1,N_2) \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} &&\fl \mathcal {L}_1 P_1^{(n)}(N_1,N_2) =g\big[P_1^{(n)}(N_1-1,N_2)-P_1^{(n)}(N_1,N_2)\big] \nonumber \\ &&+\,d\big[(N_1+1)P_1^{(n)}(N_1+1,N_2)+(N_2+1)P_1^{(n)}(N_1,N_2+1) \nonumber \\ &&-\,(N_1+N_2)P_1^{(n)}(N_1,N_2)\big] \end{eqnarray} \tag{ 61 }$$

and the inhomogenous terms are

$$\begin{eqnarray} &&f_0^{(n)}(N_1,N_2)= P_1^{(n-1)}(N_1-1,N_2)+P_2^{(n-1)}(N_1,N_2-1) \end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray} &&f_1^{(n)}(N_1,N_2)=-P_1^{(n-1)}+b(N_1+1)P_0^{(n)}(N_1+1,N_2) . \end{eqnarray} \tag{ 63 }$$

As noted above, we can first determine the zeroth order P⁽⁰⁾₁(N₁, N₂) and P⁽⁰⁾₂(N₁, N₂) = P⁽⁰⁾₁(N₂, N₁) as functions of the parameters of the model. Then, owing to the form of the equation (62), f⁽¹⁾₀ is determined. This allows us to solve for P⁽¹⁾₀ which in turn determines f⁽¹⁾₁ and allows us to solve for P⁽¹⁾₁. Continuing in this fashion the rest of the probabilities will be found following the structure:

$$\begin{equation} P_1^{(0)}\rightarrow P_0^{(1)} \rightarrow P_1^{(1)} \rightarrow P_0^{(2)} \rightarrow P_1^{(2)} \cdots \end{equation} \tag{ 64 }$$

In general, P⁽ⁿ⁾₀ will be found before P⁽ⁿ⁾₁.

Having laid out the general perturbation scheme and established the zeroth order (n = 0) solution, we now outline how equations (59) can be solved. We note that for each switch state S the same linear operator $\mathcal {L}_S$ appears at all orders n. This means that the homogeneous parts of the equations are independent of the order (with only the inhomogeneous term on the right-hand side varying between orders) and that they only have to be solved once.

Let us define a Green function Q₀(N₁, N₂|N⁰₁, N⁰₂) for the operator $\mathcal {L}_0$ through

$$\begin{equation} \mathcal {L}_0 Q_0\big(N_1,N_2\big|N^0_1,N^0_2\big) = -\delta _{N_1,N_1^0}\delta _{N_2,N_2^0} \end{equation} \tag{ 65 }$$

so that the solution of (59) may be written

$$\begin{equation} P_{0}^{(n)}(N_1,N_2)= \sum _{N_1^0,N_2^0}f_0^{(n)}\big(N_1^0,N_2^0\big)Q_0\big(N_1,N_2\big|N_1^0,N_2^0\big). \end{equation} \tag{ 66 }$$

We define a generating function

$$\begin{equation} K_0\big(z_1,z_2\big|N_1^0,N_2^0\big) = \sum ^{\infty }_{N_1=0}\sum ^{\infty }_{N_2=0} z_1^{N_1}z_2^{N_2}Q_{0}\big(N_1,N_2\big|N_1^0,N_2^0\big) \end{equation} \tag{ 67 }$$

the equation for which is obtained by summing (65)

$$\begin{equation} a(z_1)\frac{\partial K_0}{\partial z_1}+a(z_2)\frac{\partial K_0}{\partial z_2}+ g(z_1+z_2-2)K_0 = -z_1^{N_1^0}z_2^{N_2^0} \end{equation} \tag{ 68 }$$

with a(z_i) = d − (d + b)z_i.

In order to solve (68) we use the method of characteristics (see e.g. [28]). The characteristic equations are

$$\begin{eqnarray} &&\frac{{\rm d}z_1}{{\rm d}s} = a(z_1) \end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray} &&\frac{{\rm d}z_2}{{\rm d}s} = a(z_2) \end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray} &&\frac{{\rm d}K_0}{{\rm d}s} = -g(z_1+z_2-2)K_0 -z_1^{N_1^0}z_2^{N_2^0} \end{eqnarray} \tag{ 71 }$$

where s is a time-like parameter and (69), (70) define characteristic curves along which the partial differential equation (68) reduces to the ordinary differential equation (71). Solving (69), (70) yields the curves

$$\begin{equation} z_1 = \frac{d}{d+b} + A v \qquad z_2 = \frac{d}{d+b} + B v \end{equation} \tag{ 72 }$$

where

$$\begin{equation} v = {\rm e}^{-(d+b)s} \end{equation} \tag{ 73 }$$

and A, B are two constants to be fixed.

Equation (71) can be rewritten as an ordinary differential equation in v as

$$\begin{eqnarray} &&\fl \frac{{\rm d}}{{\rm d}v} \big[ K_0(v) v^{\frac{2gb}{(d+b)^2}} {\rm e}^{-\frac{g(A+B)}{(d+b)}v} \big] =\frac{1}{d+b} v^{\frac{2gb}{(d+b)^2}-1} {\rm e}^{-\frac{g(A+B)}{(d+b)}v} \left(\frac{d}{d+b} + A v\right)^{N_1^0} \left(\frac{d}{d+b} + B v\right)^{N_2^0}. \nonumber\\ \end{eqnarray} \tag{ 74 }$$

To integrate (74) we choose an end-point of the integration as v = 1 and set this to correspond to an arbitrary point in the z₁–z₂ plane. This fixes the two constants A, B as

$$\begin{equation} A= z_1 - \frac{d}{d+b} \qquad B= z_2 - \frac{d}{d+b}. \end{equation} \tag{ 75 }$$

Thus we obtain

$$\begin{eqnarray} &&\fl K_0(z_1,z_2) =\frac{1}{d+b} \int _{0}^1 {\rm d}v v^{\frac{2gb}{(d+b)^2}-1} \exp \left\lbrace -\frac{g}{d+b}\left( \frac{2d}{d+b} -z_1-z_2\right)(1-v) \right\rbrace \nonumber \\ &&\times \left(\frac{d}{d+b}(1-v) + vz_1\right)^{N_1^0} \left(\frac{d}{d+b}(1-v) + v z_2\right)^{N_2^0} . \end{eqnarray} \tag{ 76 }$$

We now expand as a power series in z₁, z₂ to obtain Q₀(N₁, N₂|N⁰₁, N⁰₂) from (67)

$$\begin{eqnarray} &&\fl Q_0\big(N_1,N_2\big|N_1^0,N_2^0\big) = \frac{{\rm e}^{ -\frac{2gd}{(d+b)^2}}}{d+b} \sum _{p=0}^{N_1} \frac{1}{p!}\left( \frac{g}{d+b}\right)^p \left(\begin{array}{@{}c@{}} N_1^0 \\ N_1-p\end{array} \right) \left(\frac{d}{d+b}\right)^{N_1^0-N_1+p} \nonumber\\ &&\times\, \sum _{r=0}^{N_2} \frac{1}{r!}\left( \frac{g}{d+b}\right)^r {N_2^0 \choose N_2-r} \left(\frac{d}{d+b}\right)^{N_2^0-N_2+r} \nonumber \\ &&\times\, I\left(\frac{2gb}{(d+b)^2}{+}N_1{-}p{+}N_2{-}r, N_1^0{+}N_2^0{-}N_1{-}N_2{+}2p{+}2r{+}1; \frac{2gd}{(d+b)^2} \right)\nonumber\\ \end{eqnarray} \tag{ 77 }$$

where I(α, β; x) is defined as the integral

$$\begin{equation} I(\alpha ,\beta ;x) =\int _0^1 {\rm d}v\, v^{\alpha -1}(1-v)^{\beta -1} {\mathrm{e}}^{x v} = \frac{\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}\ {}_1F_1(\alpha ,\alpha +\beta ;x) , \end{equation} \tag{ 78 }$$

and Γ(z) and ₁F₁(a, b; z) are the usual Euler Gamma function and confluent hypergeometric function respectively.

If we define a Green function for the operator $\mathcal {L}_1$ through

$$\begin{equation} \mathcal {L}_1 Q_1\big(N_1,N_2\big|N^0_1,N^0_2\big) =-\delta _{N_1,N_1^0}\delta _{N_2,N_2^0}, \end{equation} \tag{ 79 }$$

this equation only has a solution when N₂ ≠ 0. To see this we note that summing the left-hand side of (79) over all N₁, N₂ yields zero, i.e. the operator conserves probability. Therefore, one cannot solve (79) or indeed (59) for an arbitrary right-hand side; one requires that the sum of the right-hand side of (59) over all N₁, N₂ yields zero. However, since the null space of the operator $\mathcal {L}_1$ is concentrated on N₂ = 0 (i.e. the stationary state in (58) is proportional to $\delta _{N_2,0}$) we can find a solution of (79) for N₂ > 0.

It is simplest to proceed by considering the solution P₁ for an arbitrary right-hand side −h(N₁, N₂)

$$\begin{equation} \mathcal {L}_1 P_1(N_1,N_2) = -h (N_1,N_2) \end{equation} \tag{ 80 }$$

that satisfies

$$\begin{equation} \sum _{N_1=0}^{\infty }\sum _{N_2=0}^{\infty }h(N_1,N_2) =0. \end{equation} \tag{ 81 }$$

We define generating functions

$$\begin{eqnarray} K_1(z_1,z_2) &=& \sum ^{\infty }_{N_1=0}\sum ^{\infty }_{N_2=0} z_1^{N_1}z_2^{N_2}P_1(N_1,N_2) \end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray} H(z_1,z_2) &=& \sum ^{\infty }_{N_1=0}\sum ^{\infty }_{N_2=0} z_1^{N_1}z_2^{N_2}h(N_1,N_2) \end{eqnarray} \tag{ 83 }$$

then summing (80) yields

$$\begin{equation} d(1-z_1)\frac{\partial K_0}{\partial z_1} +d(1-z_2)\frac{\partial K_0}{\partial z_2}+ g(z_1-1)K_0 = -H(z_1,z_2). \end{equation} \tag{ 84 }$$

In order to solve (84) we again use the method of characteristics. The characteristic equations are this time

$$\begin{eqnarray} &&\frac{{\rm d}z_1}{{\rm d}s} = d(1-z_1) \end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray} &&\frac{{\rm d}z_2}{{\rm d}s} = d(1-z_2) \end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray} &&\frac{{\rm d}K_1}{{\rm d}s} = -g(z_1-1)K_1 - H(z_1,z_2) . \end{eqnarray} \tag{ 87 }$$

Solving the first two equations yields

$$\begin{equation} z_1 = 1 + A v \qquad z_2 = 1 + B v \end{equation} \tag{ 88 }$$

where now

$$\begin{equation} v = {\rm e}^{-ds} \end{equation} \tag{ 89 }$$

and A, B are constants to be fixed. The final equation becomes

$$\begin{equation} \frac{{\rm d}}{{\rm d}v}\big[ K_1\, {\rm e}^{-A\frac{g}{d}v}\big] = \frac{H(z_1,z_2)}{d v}\, {\rm e}^{-A\frac{g}{d}v}. \end{equation} \tag{ 90 }$$

We choose the integration to be from v = 0 to v = 1 where v = 0 corresponds to z₁ = z₂ = 1 and v = 1 corresponds to an arbitrary point in the z₁–z₂ plane which implies A = z₁ − 1 and B = z₂ − 1. We then obtain the solution of (90)

$$\begin{eqnarray} &&\fl K_1(z_1,z_2) = K_1(1,1)\, {\rm e}^{(z_1-1)\frac{g}{d}} +\frac{1}{d} \int _0^1 {\rm d}v \frac{{\rm e}^{(z_1-1)\frac{g}{d}(1-v)}}{v} H(1+(z_1-1)v, 1+(z_2-1)v). \nonumber\\ \end{eqnarray} \tag{ 91 }$$

Expanding as a power series in z₁, z₂ implies

$$\begin{eqnarray} &&\fl P_1(N_1,N_2) = K(1,1) {\rm e}^{-g/d} \frac{(g/d)^{N_1}}{N_1!} \delta _{N_2,0} \nonumber \\ &&+\, {\rm e}^{-g/d} \sum _{p=0}^\infty \sum _{q=0}^\infty \frac{h(p,q)}{d} {q \choose N_2} \sum _{r=0}^{N_1} \left(\frac{g}{d}\right)^{r} \frac{1}{r!} {p \choose N_1-r} \nonumber \\ &&\times\, I(N_1+N_2-r,p+q-N_1-N_2+2r+1; g/d). \end{eqnarray} \tag{ 92 }$$

As discussed above, for N₂ > 0 we can define the Green function (79) by means of which the solution of (80) may be written

$$\begin{equation} P_1(N_1,N_2) = \sum _{N_1^0,N_2^0} Q_1\big(N_1,N_2\big|N_1^0,N_2^0\big)h\big(N_1^0,N_2^0\big). \end{equation} \tag{ 93 }$$

Then we can read off the Green function from (92) as

$$\begin{eqnarray} &&\fl Q_1\big(N_1,N_2\big|N_1^0,N_2^0\big) = {\rm e}^{-g/d} \frac{1}{d}{N_2^0 \choose N_2} \sum _{r=0}^{N_1} \left(\frac{g}{d}\right)^{r} \frac{1}{r!} {N_1^0 \choose N_1-r} \nonumber\\ &&\times\, I\big(N_1+N_2-r,N_1^0+N_2^0-N_1-N_2+2r+1; g/d\big). \end{eqnarray} \tag{ 94 }$$

For N₂ = 0 the solution (92) reads

$$\begin{eqnarray} &&\fl P_1(N_1,0) = K(1,1)\, {\rm e}^{-g/d} \frac{(g/d)^{N_1}}{N_1!} \nonumber \\ &&+\, {\rm e}^{-g/d} \sum _{p=0}^\infty \sum _{q=0}^\infty \frac{h(p,q)}{d} \sum _{r=0}^{N_1} \left(\frac{g}{d}\right)^{r} \frac{1}{r!} {p \choose N_1-r} \nonumber \\ &&\times\, I(N_1-r,p+q-N_1+2r+1; g/d). \end{eqnarray} \tag{ 95 }$$

Some care is required with the r = N₁ term of the sum in (95) since the integral I(0, β; x) does not converge. However, the property (81) of the function h implies that the coefficient of the offending integral is zero. To see this one can write the r = N₁ term of (95) as

$$\begin{eqnarray} &&\fl {\rm e}^{-g/d} \sum _{p=0}^\infty \sum _{q=0}^\infty \frac{h(p,q)}{d} \left(\frac{g}{d}\right)^{N_1} \frac{1}{N_1!} I(0,p{+}q{+}N_1{+}1;g/d) \nonumber \\ &&= \sum _{p=0}^\infty \sum _{q=0}^\infty \frac{h(p,q)}{d} \left(\frac{g}{d}\right)^{N_1} \frac{1}{N_1!} \int _0^1 {\rm d}v\, {\rm e}^{-\frac{g}{d}(1-v)}v^{-1}(1-v)^{p+q+N_1} \nonumber \\ &&= \sum _{p=0}^\infty \sum _{q=0}^\infty \frac{h(p,q)}{d} \left(\frac{g}{d}\right)^{N_1} \frac{1}{N_1!} \sum _{s=1}^{p+q} (-1)^s\, {\rm e}^{-\frac{g}{d}} {p+q \choose s} I(s,N_1+1; g/d). \end{eqnarray} \tag{ 96 }$$

In the final equality, the binomial expansion of (1 − v)^{p + q} has been used with the term s = 0 not present since its coefficient vanishes due to (81). All the integrals in (96) then converge.

The first-order contribution to the stationary probability P⁽¹⁾₀(N₁, N₂) is given by:

$$\begin{equation} P_0^{(1)}(N_1,N_2) =\sum _{N_1^0,N_2^0}f_0^{(1)}\big(N_1^0,N_2^0\big)Q_0\big(N_1,N_2\big|N_1^0,N_2^0\big) \end{equation} \tag{ 97 }$$

with

$$\begin{equation} \fl f_0^{(1)}\big(N_1^0,N_2^0\big) =\frac{1}{2}\exp \left(-\frac{g}{d}\right) \left[ \left(\frac{g}{d}\right)^{N_1^0-1}\frac{\delta _{N_2^0,0}}{(N_1^0-1)!}+\left(\frac{g}{d}\right)^{N_2^0-1}\frac{\delta _{N_1^0,0}}{(N_2^0-1)!}\right]. \end{equation} \tag{ 98 }$$

Although Q₀(N₁, N₂|N⁰₁, N⁰₂) is expressed in terms of known integrals in (77), it is advisable to go back to the explicit expression of the integrals (78) to evaluate the sums appearing in (97) more easily. In that way, P⁽¹⁾₀(N₁, N₂) can be written as:

$$\begin{eqnarray} &&\fl P_0^{(1)}(N_1,N_2) =\frac{1}{2}\frac{{\rm e}^{-\frac{g}{d} -\frac{2gd}{(d+b)^2}}}{d+b} \sum _{N_1^0=0}^{\infty }\left(\frac{g}{d}\right)^{N_1^0-1}\frac{1}{(N_1^0-1)!}\nonumber\\ &&\times\,\sum _{p=0}^{N_1}\frac{1}{p!}\left(\frac{g}{d+b}\right)^{p}{N_1^0 \choose N_1-p}\left(\frac{d}{d+b}\right)^{N_1^0-N_1+p}\frac{1}{N_2!}\left(\frac{g}{d+b}\right)^{N_2}\nonumber \\ &&\times\, \int _0^1 {\rm d}v\, v^{\frac{2gb}{(d+b)^2}+N_1-p-1}(1-v)^{N_1^0-N_1+N_2+2p}\,{\rm e}^{\frac{2gd}{(d+b)^2}v}+ {\rm symm} \end{eqnarray} \tag{ 99 }$$

where the label symm refers to the fact that there will be another term equal to the written one, apart from a switch in the variables N⁰₁ and N⁰₂.

In appendix C it is shown how the expression may be simplified to the result

$$\begin{eqnarray} &&\fl P_0^{(1)}(N_1,N_2) = \frac{1}{2}\exp \left(\frac{g(b-d)}{(d+b)^2}-\frac{g}{d}\right)\left(\frac{g}{d}\right)^{-1}\frac{1}{d+b}\frac{1}{N_2!}\left(\frac{g}{d+b}\right)^{N_2}\sum _{m=0}^{N_1}\left(\frac{g}{d+b}\right)^{N_1-m}\nonumber \\ &&\times\, \left(\frac{g}{d}\right)^m \frac{1}{(N_1\,{-}\,m)!m!}\Bigg [\frac{g}{d\,{+}\,b}I\left(\frac{2gb}{(d\,{+}\,b)^2}\,{+}\,m,N_2\,{+}\,N_1\,{-}\,m\,{+}\,2;\frac{g(d\,{-}\,b)}{(d\,{+}\,b)^2}\right) \nonumber \\ &&+\,mI\left(\frac{2gb}{(d+b)^2}+m,N_2+N_1-m+1;\frac{g(d-b)}{(d+b)^2}\right)\Bigg ]+{\rm symm}. \end{eqnarray} \tag{ 100 }$$

Once we have P⁽¹⁾₀, we can plug this result into the P⁽¹⁾₁ equation which becomes for N₂ > 0

$$\begin{equation} P_1^{(1)}(N_1,N_2) =\sum _{N_1^0,N_2^0}f_1^{(1)}\big(N_1^0,N_2^0\big)Q_1\big(N_1,N_2\big|N_1^0,N_2^0\big) \end{equation} \tag{ 101 }$$

where f⁽¹⁾₁(N₁, N₂) = −P⁽⁰⁾₁(N₁, N₂) + b(N₁ + 1)P⁽¹⁾₀(N₁ + 1, N₂).

We consider separately the two terms of f⁽¹⁾₁. The first is

$$\begin{equation} {-}P_1^{(0)}\big(N_1^0,N_2^0\big)=-\frac{1}{2}\frac{\exp \big({-}\frac{g}{d}\big)}{N_1^0!}\left(\frac{g}{d}\right)^{N_1^0}\delta _{N_2^0,0}. \end{equation} \tag{ 102 }$$

Since the calculation with the Green function (see section 4.6) is only for N₂ ≠ 0 this term does not contribute and the only contribution comes from the second term in f⁽¹⁾₁ involving P⁽¹⁾₀. The resulting expression, obtained by substituting h = b(N₁ + 1)P⁽¹⁾₀(N₁ + 1, N₂) in (92) with P⁽¹⁾₀ given by (100), is

$$\begin{eqnarray} &&\fl P_1^{(1)}(N_1,N_2) =\frac{b}{2g} \exp \left(\frac{g(b-d)}{(d+b)^2}-\frac{2g}{d}\right)\sum _{N_1^0=0}^{\infty }\sum _{N_2^0=0}^{\infty }(N_1^0+1) \nonumber\\ &&\times\,\Bigg \lbrace \frac{1}{(d+b)}\frac{1}{N_2^0!}\left(\frac{g}{d+b}\right)^{N_2^0} \sum _{m=0}^{N_1^0+1}\left(\frac{g}{d+b}\right)^{N_1^0+1-m}\left(\frac{g}{d}\right)^m\frac{1}{\big(N_1^0+1-m\big)!m!} \nonumber \\ &&\times\,\Bigg [\frac{g}{d+b}I\left(\frac{2gb}{(d+b)^2}+m,N_2^0+N_1^0-m+3;\frac{g(d-b)}{(d+b)^2}\right)\nonumber \\ &&+\,mI\left(\frac{2gb}{(d+b)^2}+m,N_2^0+N_1^0-m+2;\frac{g(d-b)}{(d+b)^2}\right)\Bigg ]\nonumber\\ &&+\,{\rm symm}\big(N_1^0+1,N_2^0\big)\Bigg \rbrace {N_2^0 \choose N_2} \sum _{r=0}^{N_1} \left(\frac{g}{d}\right)^{r} \frac{1}{r!} {N_1^0 \choose N_1-r}\nonumber\\ &&\times\,I\big(N_1+N_2-r,N_1^0+N_2^0-N_1-N_2+2r+1; g/d\big), \end{eqnarray} \tag{ 103 }$$

where the symmetric term in this case corresponds to the term inside the curly brackets, exchanging the places of N⁰₁ + 1 and N⁰₂ coming from the symmetric term in P⁰₀. For N₂ = 0 we obtain the result shown in (95) and cannot be simplified further.

Equations (95), (100) and (103) are the main results of this section and give closed form expressions for the first-order contributions to the stationary probabilities. The normalization constant K(1, 1) appearing in (95) is obtained from the condition:

$$\begin{equation} \sum _{N_1=0}^{\infty }\sum _{N_2=0}^{\infty } \big[P_0^{(1)}(N_1,N_2)+P_1^{(1)}(N_1,N_2)+P_2^{(1)}(N_1,N_2)\big]=0. \end{equation} \tag{ 104 }$$

Once the values of the probabilities P⁽¹⁾_i(N₁, N₂) have been obtained numerically, they are multiplied by the unbinding parameter u and added to the zeroth order probabilities. In that way, the probability distributions P_i(N₁, N₂) can be computed and plotted up to the first order of the expansion. Figure 4 shows the probability distributions for different values. The agreement is visually good in the first two examples: typical E. coli values (4) and another interesting case, with smaller g. The order of magnitude is well reproduced in almost every point of the probability distribution. However, the approximation does not work accurately if the value of the unbinding parameter u is of the order or bigger than the other parameters.

Zoom In Zoom Out Reset image size

Although the first two examples of figure 4 seem visually in good agreement with the simulations, the best way to check this is to plot different slices of the probability distribution, that is the probability distribution of N₂ where N₁ is held constant (figure 5). We choose the values of N₁ to correspond to the slices with large probability mass, i.e. N₁ = 0, 1, 2 for E. coli values. Figure 5 shows that along the slice with greatest probability mass (N₁ = 0) the analytical and simulation plots show good agreement. For N₁ > 0, there is reasonable agreement for the E. coli values (figure 5(b)), whereas for larger u the agreement is not so good (figure 5(d)).

Zoom In Zoom Out Reset image size

Since the proposed method is general, calculations can be performed up to any necessary order to get better results. For example, in the case of E. coli in axes N₁ = 1, N₂ = 1 is enough to compute the second order to have very good results (figure 6). In general, the method can be iterated as many times as required.

Zoom In Zoom Out Reset image size