Cent. Eur. J. Phys. • 12(9) • 2014 • 615-621 DOI: 10.2478/s11534-014-0497-0

VERS ITA

Central European Journal of Physics

Self-regulating genes.

Exact steady state solution by using Poisson representation

Research Article

István P. Sugár1*, István Simon2t

1 Department of Neurology and Center for Translational Systems Biology, Ichan School of Medicine at Mount Sinai, New York, NY 10029

2 Institute of Enzymology, Research Center for Natural Sciences, Hungarian Academy of Sciences, Budapest, Hungary

Received 20 December 2013; accepted 28 May 2014

Abstract: Systems biology studies the structure and behavior of complex gene regulatory networks. One of its aims

is to develop a quantitative understanding of the modular components that constitute such networks. The self-regulating gene is a type of auto regulatory genetic modules which appears in over 40% of known transcription factors in E. coli. In this work, using the technique of Poisson Representation, we are able to provide exact steady state solutions for this feedback model. By using the methods of synthetic biology (P.E.M. Purnickand Weiss, R., Nature Reviews, Molecular Cell Biology, 2009, 10: 410-422) one can build the system itself from modules like this.

PACS C200B): 87.10.-e; 87.16.-b

Keywords: genetics • master equations • Poisson representation • probability distribution

© Versitasp. zo.o.

1. Introduction

An Important aim of systems biology Is to dissect gene regulatory networks to modular components [1]. In the last decade, by using graph theoretical methods, important functional modules were found [2-4]. These extracted, relatively simple modules can be studied theoretically. Once the modules are well characterized they can be coupled

^E-mail: Lstvan.sugar@mssm.edu (Corresponding author) ^E-mail: simon@enzLm.hu

to Increase network complexity and to better understand cellular behavior.

Because of the low reactant number in the cell the fluctuations should be included into the theoretical description of a module, i.e. a stochastic theory is needed. The stochastic processes can be either simulated by using Gillespie's [5] method or by analytically solving the respective master equations.

Among the identified modules, feedback loops are present in many cellular networks [2, 4, 6, 7] and self-regulating genes are simple examples of these modules. So far two models of self-regulating genes have been developed and studied theoretically: Hornos' model [8] and Iyer-Biswas'

Springer

model [9]. These are two-state models of a gene which produces a protein that regulates its own activity. In Iyer-Biswas' model the gene is either in 'off' state (no protein production at all) or in 'on' state (proteins are produced). In the positive (negative) feedback model, the amount of protein produced proportionally increases the propensity of the gene to be in the 'on (off)' states. In Hornos' model the gene is in a state when the operator site is free, while it is in p state when the protein is bound to the operator site. The rates for protein production, ga and gp, are different for the free and bound state. ga > gp and ga < gp in the case of a self-repressing and self-activating gene, respectively.

In this article we consider a slightly modified version of the self-regulating gene model originally developed by Hornos et al. [8]. We also correct the master equations of Hornos' model and calculate the steady state solution by using the method of Poisson Representation [10, 11]. Note that the same corrections of the master equations have been made by Grima et al. [12], and solved the equations by using the generating function method1.

2. Model

In our self-regulating gene model a single gene produces a protein that represses or enhances its own activity. The six elementary reactions where a and ¡3 denote the protein unbound and protein bound states of the gene, respectively are:

a + P—^ß

ß —-—> a + P

protein binding (1)

protein dissociation (2)

ß--—> a degradation of bound protein

a ———> a + P unbound gene produces protein

ß-> ß + P bound gene produces protein

P ———> 0 degradation of unbound protein

1 The errors in Hornos' master equations were noticed by us independently from Grima et al. We were not aware of Grima's work that came out online at July 20th, 2012, until the referee of J. Phys. Rev. E drew our attention to Grima's work on October 24, 2012

where P is the unbound protein, h is the bimolecular rate of protein binding to the gene, f is the unimolecular rate of protein release from the gene. The rates for protein production ga and gp are different for protein free and protein bound states of the gene. ga > gp and ga < gp in the case of self-repression and self-activation, respectively. k is the degradation rate of an unbound protein. The above reactions and notations are similar to the ones in Hornos' model [8] except Eq. 3 describes the degradation of a bound protein with a rate of kt. In Hornos' model this reaction is taken into consideration only when there is no unbound protein present and the protein degradation rate is kt = k. In all subsequent version of Hornos' model [13-16] the degradation rate of the bound protein is taken to be zero. In our model, for the sake of generality, any value of kt is acceptable. This is the case because the degradation rate may depend on the actual gene-protein interaction.

2.1. Master equations

The probability of finding n proteins (unbound) in the system at time t can be calculated by solving the master equations. The master equations of our model at n > 0 are:

dPa (n,t) dt

= fPp(n - 1,t) + kbPp(n, t)+ gaPa(n - 1 , t) + k[n + 1]Pa(n + 1 , t) - [hn+ga+ kn]Pa(n, t)

dPß (n,t)

= k[n +1]Pß(n + 1,t) + gßPß(n - 1,t)

+ h[n+ 1]Pa(n + 1, t) - [f + kb+ gß + nk]Pß(n, t)

where Pa(n, t) and Pp(n, t) are the Individual probabilities that the gene is unbound and bound, respectively, while immersed in a solution containing n unbound proteins at time t. In addition it is physically plausible to make the following additional restrictions:

Pa(n, t) = Pß(n, t) = 0 at n < 0

Since we intend to get the steady state solutions of the above master equations, initial conditions are not needed. In the Supplemental Material (Part 1) detailed comparisons are made between our master equations, Eqs. 7-8, and the master equations utilized in the works of Hornos and his coworkers [8, 13-16]. We point out differences

between the two sets of equations that are either because of the extra reaction in our model (Eq. 3) or because of the erroneous terms in Hornos' equations.

2.2. Equations of generating functions

In order to get the steady state solutions of the master equations, we rewrite Eqs. 7 - 9 in the form of generating equations, where the generating functions are defined as:

G¡(s, t) = ¿ s"Pi(n,t)

and for generating function :

dC, (s't)= KddC^+ g,sC, {s,t) + hdC* M)

- (f + kb + g,)C,(s,t) - ks

ds dC,(s, t) ds

2.3. Poisson representation at steady state

Let us solve Eqs. 7, 8 at steady state by using the method of Poisson representation [10]. The method assumes the existence of pa (A) and p$ (A) functions that yield the steady state solutions, Pa(n, to) and P$(n, to), by

where i = a or $.

We obtain the following equation for generating function Ga:

Pi(n, to) = J dAp¡(A)e

dCa(t' t) = fsCe (s, t) + kbCe (s, t) + gasCa (s, t)

hdCa(s,t) dCa(s, t) + k----(h + k)s-—--gaCa(s, t)

where i = a or $.

After substituting the above forms of the probability functions into the generating functions in Eqs. 11, 12, we get the following equations for the pa(A) and p$(A) functions:

0= [fpe (B) + {g a - hB - kB^pa (B) eB<s-1> - fp, (A) + {ga - hA - kA^pa (A)

3A(s-1)

\a\\(„ (dP$ r, dPa + Ld(APa)+u d(APa) ,, I dAj fp, - f— + kbP, - ga^Y + k dA + ^^--hAPa

pA(s-1)

0 = [g, - kB]pe(B)eB<s-1> - [g, - kÀ\pe(A)eA<s-1> + j d^k

d¥ - g>dS+hAp. - fp, -

aA(s-1)

The detailed derivation and the solutions of Eqs. 14, 15 are given in the Supplemental Material (Pa rt 2).

mental Material (Part 2)). The functions of the expansion coefficients, pa(A) and (A) are:

3. Results

The steady state solution of the master equations, Eqs. 7 -8, is based on an expansion of the probability distribution in Poisson distributions (see Eq. 13). The functions of the expansion coefficients, pa(A) and (A), are the solutions of Eqs. 14,15. With the proper selection of the integration boundaries, A and B, the first two terms of each equation, Eqs. 14, 15, become zero (see derivation in the Supple-

Pa (A) = D

f + g, - kA kA + hA - ga

" D ' f + g, - kA'

k + h A - g- Ia k+h J

h g,\F-L\i ga \

Ia - k\ \A-^h\

„,n nil g,\F-Lh g a \C „EA p,(A) = D\A - k \ \A - ^ \ e

Table 1. Integration Boundaries

Conditions Integration boundaries Comments

Sß < ga k < k+h F - L> 0 and G - 1 > 0 II H« b — H+h pa ■ Pß < 0 if '+«> X > 0 if f+gß<A > 0

gß > ga k > k+h A — ga A — k+h B — k Pa ■ Pß > 0 if f+gß>X < 0 if f+ß<A > 0

F - L + G < 0 F - L > 0 and G - 1 < 0 a — 9ß B — TO Self contradicting conditions*

F - L < 0 and G - 1 > 0 A — g» A — k+h B — TO Self contradicting conditions**

All other conditions Poisson representation does not exist

* From conditions F — L + G < 0 and G — 1 < 0 follows F — L < —1, and this contradicts with condition F — L > 0. ** From conditions F — L + G < 0 and F — L < 0 follows G < 0, and this contradicts with condition G — 1 > 0.

where D is the integration constant and

(k + h)

k 2(k+h) ts - f] h

ga - f + gß k + h k

ga k+h

ga __ gß

(k - kb - f)

By means of a normalization condition we can calculate the integration constant, D (Supplemental Material (Part 3))

Table 1 lists the integration boundaries, A and B, at which the first two terms in both Eq. 14 and Eq. 15 become zero. After substituting Eqs. 16, 17 into Eq. 13 one can calculate the steady state solutions of the master equations (Eqs. 7, 8). The integration boundaries in Eq. 13 are given in Table 1. The integrals can be obtained in closed forms when either > g = 0or ^ > g = 0and thus one can get Pa(n, to) and (n, to) in closed form too (Supplemental aterial (Part 4)).

In general the integral in Eq. 13 has been calculated numerically by Romberg's method [17]. In Supplemental Material (Part 5) we point out that the solution obtained from the numerical integration is consistent with the master equations (Eqs. 7, 8).

The subfigures of Figure 1 show the calculated steady state probability distributions at twelve different parameter sets. All but one of these parameter sets are characteristic to self-restricting genes, i.e. where ga > . Namely, the parameters at subfigure 1D refer to a self-activating gene.

Subfigures 1A-D belong to increasing values of the g-parameter from 0 to 200, while the other parameters are fixed. With increasing g one expects an increasing number of proteins produced by the bound gene, while with increasing number of protein the overall probability of bound gene state increases too. These qualitative expectations are supported by the calculated values (see figure legends) and the probability distributions in Subfigures 1A-D.

Subfigures 1E-H belong to increasing values of the g parameter from 0.2 to 20, while the other parameters are fixed. With increasing g one expects an increasing probability of the unbound gene state. Since ga > g = 0, the increasing probability of unbound gene state results in an increasing number of proteins produced by the unbound gene. Again these qualitative expectations are supported by the calculated overall probabilities of the unbound gene state, Pa(= 1 — ), and probability distributions in Subfigures 1E-H.

Finally, subfigures 1I-L belong to increasing values of the J parameter from 0.015 to 0.9, while the other parameters are fixed. With increasing j one expects an increasing overall probability of the bound gene state. Since ga > Sgr = 0, the increasing probability of the bound gene state results in a decreasing number of proteins produced by the unbound gene. These qualitative expectations are supported by the calculated overall probabilities of the bound gene, , and probability distributions in Subfigures 1I-L.

The subfigures of Figure 2 show the calculated steady state probability distributions at twelve different parameter sets. All these parameter sets are characteristic to self-activating genes, i.e. where ga < . We have similar qualitative expectations regarding the change of the probability distributions with the change of parameters

Figure 1. Probability distributions in the case of self-restricting genes.

Calculated steady state probability distributions at twelve different parameter sets. Each parameter set refers to self-restricting genes except the parameter set at subfigure 1D. Blue: P(n) = Pa(n, to), and red: P(n) = P$(n, to) where n is the number of unbound proteins. At each subfigure the inset shows the respective p(A) = pa(A) + p$(A) function. The total probability of the bound state of the gene, i.e.:

P$= Y. P$(n, to), and the model parameter values are given in the legends to the subfigures.

A) P$=0.728; k =0.9,k =0.2,k =80, =0,% =0.1;

B) P$=0.779; k =0.9,k =0.2,k =80,=20,% =0.1;

C) P$=0.853; k =0.9,k =0.2,-a =80,=60,% =0.1;

D) P$=0.944; k =0.9,k =0.2,-a =80,=200,% =0.1;

E) P$=0.854; k =0.9,k =0.2, k =80,=0,% =0.4;

F) P$=0.801; k =0.9,k = 1,T =80,=0,% =0.4;

G) Pp=0.66; k =0.9, k =5, k =80,

H) Pp=0.456; k =0.9, k =20, k

I) Pp=0.435; k =0.9,k =0.2,if =80, J) Pp=0.522; k =0.9,k =0.2,k =80.

L) Pp=0.899; k =0.9,k =0.2, k =80

=0, k =0.4;

'i =0, k=0.4;

i =0, k =0.015;

7i =0, k =0.025;

0.635; k =0.9,k =0.2,k =80,k =0,k =0.05;

=0, k =0.9.

as in the case of the subfigures of Figure 1. However, there is one important difference between the parameter sets listed in Figure 1 and Figure 2. In Figure 2 (except for subfigures 2B-D ) both the k and -p parameters are small. Consequently one can expect smaller aver-

age protein numbers, {n} = ^ n[Pa(n, to) + Pp(n, to)].

This expectation is supported by the calculated probability distributions. In the subfigures of Figure 2 (except for subfigures 2B-D ) the average protein number is about 5, while in the case of subfigures to Figure 1 it is, in most

Figure 2. Probability distributions in the case of self-activating genes.

Calculated steady state probability distributions at twelve different parameter sets each referring to self-activating genes. Blue: P(n) = Pa(n, to), and red: P(n) = Pp(n, to) where n is the number of unbound proteins.

At each subfigure the inset shows the respective p(A) = pa(A) + pp(A) function. The total probability of the bound state of the gene, i.e.: Pp = Y. Pp(n, to), and the model parameter values are given in the legends to the subfigures. The model parameter values in the

subfigures are:

A) Pp=0.213; ^ =0.66,{ =9.3,if =4.6,^ =4.8,£ =0.6;

B) Pp=0.311; j =0.66,j =9.3,j =4.6,M = 16,£ =0.6;

C) Pp=0.542; j =0.66, j =9.3, j =4.6,j =36,£ =0.6;

D) Pp=0.742; j =0.66, j =9.3, j =4.6,j =66,£ =0.6;

E) Pp=0.592; j =0.66, j = 1.3, j =4.6, j- =6,£ =0.6;

F) Pp=0.518; j =0.66, j =2,j- =4.6, j- =6,£ =0.6

G) Pp=0.222; j =0.66, j =9.3,j- =4.6, j- =6,£ =0.6

H) Pp=0.0835; j =0.66, j =30, j =4.6, j =6, £ =0.6

I) Pp=0.0232; j =0.66, j =9.3,j =4.6, j- =6,£ =0.05 J) Pp=0.0856; j =0.66, j =9.3,^ =4.6, j =6,£ =0.2 K) Pp=0.276; j =0.66, j =9.3,j =4.6, j =6,£ =0.8 L) Pp=0.749; j =0.66, j =9.3,j =4.6, j =6, £ =6

of the cases, much higher. Finally, it is important to note that at model parameters satisfying the conditions in Table 1 the calculated probability distributions have been always unimodal ones.

4. Conclusions

Systems biology studies the structure and behavior of complex gene regulatory networks. One of its aims is to develop a quantitative understanding of the modular com-

ponents that constitute such networks. The self-regulating gene is a type of auto regulatory genetic modules which appears in over 40% of known transcription factors in E. coli. In this work using the technique of Poisson Representation, we are able to provide exact steady state solutions for this feedback model. The model is an extended and corrected version of Hornos' model [8, 13-16]. The respective master equations are extended and corrected versions of Hornos' [8] and Ramos' equations [1416]. The corrections are presented in Table S1 (Supplemental Material) on a clear and condensed fashion. The steady state solutions of our equations, i.e. the probability distributions, are obtained in integral forms (Eq. 13). In general these integrals have been calculated numerically, while they have been obtained in closed forms when either k > it = 0 or ^ > ^ = 0. The behavior of the obtained probability distributions has been investigated at many parameter sets for both self-repressing and self-activating genes. Finally it is important to note that the steady state solution of the master equations provided by the Poisson Representation is valid at a limited parameter space (see Table 1), while the general steady state solution exists at the entire parameter space of the self-regulating gene model [12].

References

Acknowledgements

I.P.S. acknowledges the support by contract NIH/NIAID HHSN272201000054C. I.S. acknowledges the support from the Hungarian Science Research Fund (OTKA-NK 100482).

L. H. Hartwell et al., Nature 402, C47 (1999) S. S. Shen-Orr et al., Nature Genetics 31, 64 (2002) T. I. Lee et al., Science 298, 799 (2002) A. Ma'ayan et al., Science 309, 1078 (2005) D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977) F. J. Isaacs et al., Proceedings of the National Academy of Sciences of the United States of America 100, 7714 (2003)

D. Thieffry et al., Bioessays 20, 433 (1998)

J. E. M. Hornos et al., Physical Review E 72, 051907

(2005)

S. Iyer-Biswas, C. Jayaprakash, arXiv preprint arXiv:1110.2804 (2011)

C. W. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer, Berlin, 2004), Series in Synergetics

C. W. Gardiner, S. Chaturvedi, Journal of Statistical Physics 17, 429 (1977)

R. Grima, D. R. Schmidt, T. J. Newman, J. Chem. Phys. 137, 035104 (2012)

G . C. P. Innocentini, J. E. M. Hornos, J. Math. Biol. 55, 413 (2007)

A. F. Ramos, J. E. M. Hornos, Phys. Rev. Lett. 99, 108103 (2007)

A. F. Ramos, G. C. P. Innocentini, J. E. M. Hornos, Phys. Rev. E 83, 062902 (2011) A. F. Ramos et al., Iet Systems Biology 4, 311 (2010) W. H. Press et al., Numerical Recipes (Cambridge Press, New York, 1989).

Cent. Eur. J. Phys. • 12(9) • 2014 • 622-627 DOI: 10.2478/s11534-014-0497-0

VERS ITA

Central European Journal of Physics

SUPPLEMENTAL MATERIAL TO

Self-regulating genes. Exact steady state solution by

using Poisson Representation

Supplemental Material

Istvan P. Sugar1*, Istvan Simon2t

1 Department of Neurology and Center for Translational Systems Biology, Ichan School of Medicine at Mount Sinai, New York, NY 10029

2 Institute of Enzymology, Research Center for Natural Sciences, Hungarian Academy of Sciences, Budapest, Hungary

1. Comparison of the master equations

Here we compare the master equations In the main text, Eqs. 7-8, with the master equations in the papers of Hornos and his coworkers [1-5]. Our self-regulating gene model, given by Eqs. 1-6, is equivalent with the models discussed in the above papers when the degradation rate of the bound protein is zero, i.e. k = 0. In spite of the similarity of the model, the master equations in the above papers are different from Eqs. 7-8.

By using the notations of Eqs. 7-8 the master equations in the papers of Ramos et al. [3-5] are:

t] = fPß(n, t) + gaPa(n - 1 , t)

+ k[n + 1]Pa(n + 1, t) - [hn + ga + kn]Pa(n, t)

*E-mail: istvan.sugar@mssm.edu (Corresponding author) ^E-mail: simon@enzim.hu

dP^nt't] = k[n + 1]Pe(n + 1,f) + gePe(n- 1,f)

+ hnPa(n,t) - [f + ge + nkP(n,t) (S2)

It is important to mention that in Ramos' equations, Eqs. S1 - S2, n is defined as the number of 'free' proteins, i.e. unbound proteins.

ote that the 1st term on the right hand side and the 3rd term on the right hand side of Ramos' first and second equations, respectively are different from the respective terms in Eqs. 7-8. In Ramos' first equation (Eq. S1) it is not taken into consideration that the protein dissociation process increases the number of unbound proteins by 1, while in Ramos' second equation (Eq. S2) it is not taken into consideration that the protein binding process decreases the number of unbound proteins by 1. Ramos et al. [4] obtained the exact time dependent solutions of Eqs. S1, S2 by using the Heun functions. However, because of the above mentioned differences between our and Ramos' master equations we are unable to get the exact time dependent solutions of our master equations (Eqs. 7 - 8).

The master equations in Hornos' paper [1] look similar to Ramos' equations [3-5] but n is defined as the total number

Springer

Table S1. Summary of differences between the terms of the master equations in this, Homos' and Ramos' paper

This paper Hornos' paper (ref. [1]) Ramos' papers (refs. [3-5]) Comment

n = number of unbound proteins

1si master Eq. term 1 fPß(n - 1,t) fPß(n, t) k =0

2nd master Eq. term 3 h[n + 1]P„ (n + 1,t) hnPa(n, t) k =0

n = total number of proteins

2nd master Eq. term 1 knPß (n + 1,t) k[n + 1]Pß (n + 1,t) k=0

2nd master Eq. term 4 -[f + gß + k(n - 1)]Pß(n, t) -[f + gß + kn]Pß(n, t) k=0

1si master Eq. term 2 kPß (n + 1,t) missing kb - 1 k ~ 1

2nd master Eq. term 1 knPß (n + 1,t) k[n + 1]Pß (n + 1,t) kb - 1 k ~ 1

of proteins (bound and unbound). For clarity let n be the total number of proteins. Thus, by using our notations (in Eqs. 7 - 8), Hornos' equations are:

dPa^"' t] = fPp(n, t) + gaPa(n - 1 , t)

+ k[n + 1]Pa(n+ 1 , t) - [hn+ga+ kn]Pa(n, t)

dPf t) = k[n + 1]Pp(n + 1,t) + gP(n - 1,t)

+ hnPa(n, t) - [f + gp + nk]Pp(n, t) (S4)

In order to compare Hornos' equations with our master equations, let us rewrite Eqs. 7 - 8 by using ~n instead of n:

^MnJl = fPp (n, t) + kbPPj (n + 1 , t) + gaPa (n - 1 , t) + k[n + 1]Pa(n+ 1 , t) - [hn+ga+ kn]Pa(n, t)

first equation (Eq. S5). However, in the second set of equations (Eq. S4 and Eq. S6) there are two differences: one in the 1st term and one in the last term of the right hand side. This is because Hornos does not take into consideration that in the ¡3 state out of n proteins, only the unbound proteins, i.e. ~n - 1 proteins, can degrade. Finally, let us compare Hornos' equations (Eqs. S3 - S4) with our rewritten master equations (Eqs. S5 - S6), when k = 1, i.e. the degradation rate of the bound and unbound proteins are similar. We consider both conditions: k = 0 and k = 1 because it is not clear in Hornos' paper which condition was really applied. In this case our rewritten master equations are:

dPadnt,t] = fPp (n, t) + kPp (n + 1 , t) + gaPa (n - 1 , t)

+ k[n + 1]Pa(n + 1, t) - [hn + ga + kn]Pa(n, t)

dPdnt,t) = knPp (n + 1, t) + gpPp (n - 1, t)

+ hnPa(n, t) - [f + ge + kn]Pp(n, t) (S8)

^d: t) = knPp(n + 1,t) + gpPp(n - 1, t)

+ hnPa(n, t) - [f + gp+k(n - kb]Pp(n, t)

Let us compare now Hornos' equations (Eqs. S3 - S4) with our rewritten master equations (Eqs. S5 - S6) when k = 0. Hornos' first equation (Eq. S3) agrees with our

In Eq. S7 the 2nd term on the right hand side is missing in Hornos' 1st equation (Eq. S3), while in Eq. S8 the 1st term on the right hand side is different from the respective term in Hornos' equation (Eq. S4).

The following comparative Table S1 summarizes the above described differences between the terms of the master equations in our, Hornos' (Ref. [1]) and Ramos' (Refs. [3-5]) works on a condensed and clear fashion.

2. Solving the master equations by using Poisson Representation

Let us solve Eqs. 7, 8 at steady state by using the method of Poisson representation [6]. We assume the existence of pa(A) and pß(A) functions that yield Pa(n, to) and Pß(n, to) by

dApi(A)e-A-]

where i = a or ¡3.

After substituting the above forms of the probability functions into the generating functions in Eqs. 11, 12 we get the following equations for the pa(A) and (A) functions

[fpß(B) + [ga - hB - kB}pa(B) eB(s-1) - fpß(A) + [ga - hA - k^pa(A)

aA(s-1)

+ I dA| fpß - fd-§+ kbpß - gad-^ + kd(^ + h- hApa

j dA^f

aA(s-1)

0 = [gß - kB\pß(B)eB<s-1> - [gß - kA^pß(A)eA<s-1> + J dA^k^dr1 - gßdA + hApa - fpß - kbpß

d(APß) dpß

?A(s-1) (S11)

Note that Eqs. S10 and S11 are identical with Eqs. 14 and 15 in the main text, respectively. Deriving Eqs. S10, S11 the following relationships were utilized:

j dApt(A)t

,A(s-1)

eA(s-1)

to r A" f

G = C(s, TO) = ^s" dApi(A)e-A— = dAp(A)

"=0 A "' A

sG = (s - 1) J dAPieA(s-1) + Ci = -j dA^e^-1^ j dApe^ + pi(B)eB(s-1) - p(A)e

AA B B

dCi = (s - 1) J dA(piA)eA(s-1^ + JdA(pA)

AA B B

= - j dAd(dp)eAI-s-1)+ f dA(piA)eA{s-1^ pi(B)BeB(s-1) - p(A)Ae

= J dA(piA)eA<s-1>

,A(s-1)

eA(s-1)

where in the 2nd and 3rd equations we integrated by parts.

2.1. Eliminating terms at integration boundaries

In Eqs. S10, S11 the first two terms can be eliminated by properly choosing the integration boundaries A and B.

That is 0 =

fp, (B) + [ga-hB-k^pa (B) fp, M + {ga-hA-k^Pa (A)

0= [gp- kB]pp(B)eB^ - [g, - kA]pp(A),

After integrating Eq. S16 we get the steady-state solution for Pt:

p, (A) = D \A-f \

where D is the integration constant. By substituting ,A(s-i) Eq. S17 into Eq. S14 we get the steady-state solution

for Pa:

Pa (A) = Pt (A) = Pa (B) = Pt (B) = 0 (S12)

2.2. Determining Pa(A) and Pt(A) functions

By determining the p„ (A) and Pt (A) functions one can get the values of A and B where Eqs. S12 are valid. In order to calculate these functions let us take the sum of Eq. S10 and Eq. S11:

eA(s-i)

0= [dA{^A [(kX-f-9t )Pt (A)+(kA+hA-g„ )p„ (A)]

A (S13)

From Eq. S13 a relationship between p„ (A) and Pt (A) follows

[kA-f- g^Pt(A) + [kA + hA- g„]P„(A) = C (S14)

where, as a consequence of Eqs. S12, the constant, C, is equal to zero.

By properly choosing the integration boundaries the first two terms becomes zero in Eq. S11 and we get the following differential equation:

k^A^m - g^jdiA) + hAP„(A) - fPt(A) - kbPt(A) = 0

After substituting Eq. S14 into Eq. S15 we get:

dPl P,

k-kb-f

---+---

k(A-f) k(k + h)(A-k ) (A-gz)

A(kA -f-g, )

,r F-L E +--^ +

a g, a__ga-

A- k A k+h

(k + h) _fhg,

k 2(k+hn -g+h-f) h

k + h (k - kb - f)

ga - f + g, k + h k

ga k+h ga - g, k+h k

Pa (A) = D

f + gp- kA kA + hA - ga

\A-g-,\

\\ A -

" D ' 'f + g, - kA'

k + h aL A k+h J

\A-UrL\A- g a

k+h (S18)

Note that Eqs. S17 and S18 are identical with Eqs 17 and 16 in the main text, respectively.

2.3. Determining integration boundaries

From Eqs. S17, S18 one can get the integration boundaries A and B where Eqs. S12 are fulfilled. These integration boundaries are listed in the main text (Table 1).

3. Normalization condition

By means of the following normalization condition we can calculate the integration constant, D:

1 = X] Pa("■ ™) + Yi ("■

n=0 ™ r B

™ rB in ™ rB

(S15) =y dApa (A)e-A- + y dApp (A)

„-nJA n! „„ J A

f dApa (A)

J dApp (A)

f dA[pa (A)+ pp (A)

= D dA\A -

A --MlSeEA k + hi

kA - f - gp kA + hA - ga

4. Closed form solutions

When either gak > gpk = 0 or gpk > gak = 0 one can get Pa(n, <xi) and Pp(n, <xi) in closed form.

According to the Integral representations of the confluent hypergeometrlc function [7]

J dAAJ-1(b-A)v-1epA = r^r^l bJ+v-1M(j,j + v, bp)

where M is the Kummer's confluent hypergeometrlc function, and Eq. S20 holds when j > 0 and v > 0. When gak > gpk = 0 the exact steady-state solutions are:

Pp (n, to) = D fdAA-L(b - A)GeEAAne-A = DbG+n+1-' ,f(C + + 1 - L „ M(1 + n - L, 2 + n + G - L, Eb - b) p n! J r(n + 1)r(G + n + 2 - L)

Pa (n, TO) =

0 Df k + h Dk k + h

-kUA-L(b - A)G-1+ i^hA-L+1(b - A)G-1

r(n + 1 - L)f(G) l> + 1)|"(G + n + 1 - L) g„+- f(n + 2 - L)f(G)

bG+n-L b

f(n + 1)f(G + n + 2 - L)

M(1 + n - L, 1 + n + G - L, Eb - b) M(2 + n - L, 2 + n + G - L, Eb - b)

where b = -rh.

When gpk > gak = 0 the exact steady-state solutions are:

D [ n! J

Pp(n, to) = ^ I dAAG(b - A)f-LeEAAne-A

G+n+1 + F-L r(G + n + 1)r(1 + F - L) l~(n + 1)f(G + n + 2 + F - L)

M(1 G n, 2 n G F - L, Eb - b)

Pa(n, to) =

(k + h)n

! J dAAne-AeEA

(f + gp - kA)AG-1(b - A)f

D(T±^bG+n+F-L r(1+ F - L)r(G + n) M(G + n, 1+ n + G + F - L, Eb - b) k + h r(n + 1)r(G + n + 1+ F - L) v '

Dk i G+n+1+f-L 11 + F - L)r(G + n + 1)

f(n + 1)f(G + n + 2 + F - L)

M(1 G n, 2 n G F - L, Eb - b)

where b = ^.

5. Consistency check between the solution and the master equations

In the case of steady state and at n = 0 the master equations (Eqs. 7, 8) are:

Pa(0, to) = 0.00021673, Pp(0, to) = 0.018694 and Pp(1, to) = 0.020516. After substituting the above parameter values and solutions into Eq. S27 the numerical value of the left hand side, 0.011593559, agrees with that of the right hand side, 0.011566947, up to 4 decimals.

0 = kbPp(0, to) + kPa(1, TO) - gaPa(0, to) (S25)

0 = kPp(1, to) + hPa(1, to) - [f + kb + gp]Pp(0, to)

After eliminating Pa(1, to) from these equations we get

Pa(0, to) = kb kf kkb kgp k2 Pp(1, to) Pp(0, to) ga hga hga hga hga Pp(0, to)

The master equations are consistent with the solution (given by Eq. 13) if the above equality holds after substituting the values of Pa(0, to), Pp(0, to) and Pp(1, to). At the following parameter set: h = 0.1, k = 0.9, L = 0.2, k = 80, k = 0 the solution provides:

References

[1] J. E. M. Hornos et al., Phys. Rev. E 72, 051907 (2005)

[2] G. C. P. Innocentini, J. E. M. Hornos, J. Math. Biol. 55, 413 (2007)

[3] A. F. Ramos, J. E. M. Hornos, Phys. Rev. Lett. 99, 108103 (2007)

[4] A. F. Ramos, G. C. P. Innocentini, J. E. M. Hornos, Phys. Rev. E 83, 062902 (2011)

[5] A. F. Ramos et al., let Systems Biology 4, 311 (2010)

[6] C. W. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer, Berlin, 2004), Series in Synergetics.

[7] M. Abramowitz, I. A. Sregun, Handbook of Mathematical Functions (National Bureau of Standards, 1972).