Exact results in the large system size limit for the dynamics of the Chemical Master Equation, a one dimensional chain of equations

# Exact results in the large system size limit for the dynamics of the Chemical Master Equation, a one dimensional chain of equations

## Abstract

We apply the Hamilton-Jacobi equation (HJE) formalism to solve the dynamics of the Chemical Master Equation (CME). We found exact analytical expressions (in large system-size limit) for the probability distribution, including explicit expression for the dynamics of variance of distribution. We also give the solution for some simple cases of the model with time-dependent rates. We derived the results of Van Kampen method from HJE approach using a special ansatz. Using the Van Kampen method, we give a system of ODE to define the variance in 2-d case. We performed numerics for the CME with stationary noise. We give analytical criteria for the disappearance of bi-stability in case of stationary noise in 1-d CME.

###### pacs:
02.50.Ey, 87.10.-e

## I Introduction

The statistical physics of a living cell requires a theory for chemical reactions with few molecules ko08 (),qi08 (). One of mathematical tools here is the Chemical Master Equation (CME) ga (); ka92 (), describing the dynamics of probability of having different (integer) number of molecules. The same CME equation could be applied to other areas of science as well, even in the financial market theory lu09 (). Recently in qi09 () has been considered a kinetic model fe01 () for the phosphorylation-dephosphorylation cycle in the cell, and the corresponding CME was investigated. The authors considered the model with the bi-stability phenomenon, and have been derived some results for the existence of bi-stability. Here we solve exactly the dynamics of CME, a very important topics according to qi09 (). The accurate (exact) solution of the CME dynamics is important for financial market modelling lu09 (). In this article we will derive exact dynamics for the CME.

The master equation is formulated as a system of linear differential equations for of having molecules, , where N is a large integer:

 dP(X,t)Ndt=R+(X−1N)P(X−1,t)+ R−(X+1N)P(X+1,t)−(R+(XN)+R−(XN))P(X,t) (1)

Here is the growth rate and is the degradation rate.

Actually we should modify the equation at the border:

 dP(0,t)Ndt=R−(1N)P(1,t)−R+(0)P(0,t) (2)

and

 dP(N,t)Ndt=R+(1−1N)P(1−1N,t)−R−(1)P(N,t) (3)

The large parameter N describes the system volume.

A close master equations have been considered in evolution theory in hu88 ()-sa08 ().

Let us introduce the coordinate and function :

 x=X/N,p(x,t;N)≡NP(X,t) (4)

Assuming that the probability distribution is a smooth function of ,

 P(X+1,t)−P(X±1,t)≪1 (5)

one gets the Fokker-Plank equation for the model by Eq.(1).

The investigation of CME via Fokker-Plank equations meet some problems ke (). An alternative approach is to assume that is not a smooth function of , but the function is, where

 p(x,t;N)=exp[Nu(x,t)] (6)

Thus the might be un-smooth in the limit of , but still have smooth u(x,t). Such an ansatz has been assumed for the statics in hu88 (), and for the dynamics in sa07 (); ka07 (), while considering the evolution models. This ansatz gives the solution of the dynamics and steady state with the accuracy , while the approximation of the master equation via Fokker-Plank equation, assuming the smoothness of , is certainly wrong. This topic we already discussed in sa08 (), then in as10 ().

In section II we calculate the dynamics of variance for population distribution, using HJE. In section III we solve the CME dynamics for the simple case of time dependent rates. In appendix we re-derive the variance of distribution using Van Kampen method, also used that method to calculate the variance for general 1-d multi-step CME, as well as give ODE to derive the variance in case of 2-d CME.

## Ii The master equations with constant rates

### ii.1 Hamilton-Jacobi equation for the Chemical Master Equation

Using ansatz by Eq. (6), the model equations (1) can be written as Hamilton-Jacobi equations for

 ∂u∂t+H(u′,x)=0, (7)

where ,

 H(u′,x)=R+(x)+R−(x)−R+(x)e−u′−R−(x)eu′ (8)

Eqs. (7),(8) has been derived in es09 (),as10 (), where has been investigated mainly the corresponding Hamilton equations to calculate the time of extinction from meta-stabile state. We are interested in the investigation of the whole distributions, using the traditional mathematical method of characteristics. To solve the HJE, we consider the Hamilton equation for x and corresponding momentum, getting the system for the characteristics me98 (); ev02 ()

 ˙x=Hv(x,v)=R+(x)e−v−R−(x)ev, ˙v=−Hx(x,v), ˙u=vHv(x,v)−H(x,v)=v˙x+q, (9)

subject to initial conditions: , Here

The respective solution to Eq. (II.1) in - space is called the characteristic of Eq. (7).

Our Hamiltonian is time independent. Then Eq. (7) and Eq. (II.1) result in

 ˙q=0 (10)

Along the characteristic the variable is constant, so is selected to parameterize these curves.

Consider the equation

 q=−(R+(x)+R−(x))+R+(x)e−v+R−(x)ev (11)

It has a solution for

 q≥0, (12)

if at some point

 R+(x)=R−(x), (13)

and we take .

Using Eq. (11), we transform the first equation in (II.1) into

 ˙x=±√(q+R++R−)2−4R+R− (14)

Consider the following initial distribution

 u0(x)=−a(x−x0)2. (15)

with large .

The maximum of the distribution corresponds to the point , therefore .

Thus for the maximum of distribution we should consider a characteristic with .

Integrating the Eq.(14), we derive:

 T=∫xx0dyR+(y)−R−(y) (16)

Such equation was derived in qi09 ().

We can define the dynamics of the full distribution.

Let us define the function:

 T(x,q)=∫xx0dy√(q+R+(y)+R−(y))2−4R+(y)R−(y) (17)

To calculate we first calculate q for the given x from the equation

 T(x,q)=t (18)

Eq.(18) defines an implicit function

 q=Q(x,t) (19)

### ii.2 The elasticity

It is important to calculate at the point of the maximum of distribution, the ”elasticity” qi09 ().

We calculate from Eq.(19):

 ∂q∂x=∂Q(x,t)∂x (20)

We can calculate the last derivative at fixed t from the expression:

 T(x,q)=const ∂T(x,q)∂x+∂T(x,q)∂qq′x=0 (21)

It is equivalent to calculate from Eq. (17) at the fixed t.

Thus we get the following equation for at the point of maximum ():

 1R+(x)−R−(x)−q′∫xx0dy(R+(y)+R−(y))[(R+(y)−R−(y))]3=0 (22)

From Eq.(11) we can obtain:

 q′x=−(R+(x)−R−(x))v′ (23)

Thus eventually we get:

 −1v′x=(R+(x)−R−(x))2∫xx0dy(R+(y)+R−(y))[(R+(y)−R−(y))]3 (24)

Eq. (24) is the main result of our work.

In Fig. 1 is given the comparison of analytical result with the numerics.

### ii.3 Probability distribution

Eq.(17) defines the for given , we can define also using Eq.(7).

To calculate at the given , let us consider the trajectory of points connecting that point with the starting point .

We have chosen the trajectory to have . We take just as a solution of the equation

 τ=T(x(τ),q) (25)

At any point of our trajectory , we can calculate , while is constant. Eq.(11) gives

 R−e2v−(q+R++R−)ev+R+=0, v(y)=−ln2R−+ln[((q+R+(y)+R−(y)± √(q+R+(y)+R−(y))2−4R+(y)R−(y))] (26)

where we denoted .

We derive the solution of the original Eq. (2), integrating the equation along our trajectory (the characteristics connecting the points and ):

 u(x,t)=u(x,0)+∫xx0dyv(y)+qt= ∫xx0dy[−ln2R−(y)+ln[((q+R+(y)+R−(y)± √(q+R+(y)+R−(y))2−4R+(y)R−(y))]]+qt (27)

Having the expression we can calculate .

### ii.4 The restricted meaning of probability distributions in master equation

In case of evolution models [9-13], we have master equation similar to Eqs.(1)-(3), only the negative term has other coefficient than , and therefore there is no a balance condition

Contrary to the case of master equations in evolution models hu88 ()-sa08 () where all the initial distributions have a meaning, now there are some restrictions. We should clarify well the meaning of the probabilities . At every moments of time the system has only ONE value of X, and just gives such probabilities. We should solve the system (1),(2),(3) for the given initial value and for other , other initial distributions have not meaning.

Another difference is connected with the stable point solutions. Eq. (13) gives the steady state solution. If that equation has two stable solutions and the probability of these two positions is the same, i.e.

 ∫x2x1dxlogR+(x)R−(x)=0 (28)

we again should accurately interpret the HJE results qi09 (). For the rates given by Eq.(27) at , the transition is at . One can use the HJE method to calculate the mean period of time that solution, initially located at one stationary point, will move to the other stationary point qi09a (). Then it again should return back, as every moment the system could exists only with one value of X.

In case of evolution model, the system goes to the equilibrium state instead of oscillation between two stable solutions.

### ii.5 The dynamics for the stationary but random rates.

Consider now the case when the rates are smooth functions of plus some random noises. The noise in the rates is well confirmed experimentally xi10 (). We took the case of rates from qi09 (),

 R+(k1,x)=(1−x)(0.5+k1x2), R−(k2,x)=x(k2+0.01x2) (29)

where now and are random variables,

 k1=K1exp(aξ1(x)−a2/2)x2), k2=K2exp(aξ2(x)−a2/2) (30)

have normal distributions. We consider the model with , and performed numerics for different values of parameters . For , at there is a phase transition: instead of two steady state solution we get one steady state .

We can analytically estimate the transition point (from one stabile point to bi-stability), considering the behavior of

 U(x)=|k1,k2 (31)

We found that the function changes its behavior with the level of the noise. At it has only three roots , while for there is one root.

## Iii Master equation when the rates vary with time

### iii.1 The simplest solvable case

Consider the case when birth and death rate coefficients in CME change with the time as and , while they are the same for all the . We have the following Hamiltonian:

 −H=g(t)e−v+f(t)ev−g(t)−f(t),v=∂u∂x (32)

We get for the characteristics the following system of equations:

 dpdt=−∂H∂x=0 (33)

Thus the is constants:

 v=v0, (34)

For the coordinate we have a simple equation:

 dxdt=∂H∂v=g(t)e−v0−f(t)ev0 (35)

Thus we have a simple solution

 x−x0=e−v0∫t0g(τ)−ev0∫t0f(τ) (36)

We find as a function of from the solution of the last equation.

 u(x,t)=u(x0,0)+v0(x−x0)+ ∫t0[g(τ)e−v0+f(τ)ev0−g(τ)−f(τ)]dτ (37)

For the maximum of distribution we have , therefore we get an equation:

 x−x0=∫t0[g(τ)−f(τ)]dτ (38)

It is interesting to find the variance of distribution. Differentiating Eq. (36), we get, putting at the maximum point:

 −v′(x,t)=1∫t0[g(τ)+f(τ)]dτ (39)

Thus the variance of distribution decreases with the time.

### iii.2 Rates as linear functions of the number of molecules

Consider now the Hamiltonian

 −H=(ax+g(t))e−v+(bx+f(t))ev−(a+b)x−f(t)−g(t) (40)

We have

 dvdt=ae−v+bev−(a+b) (41)

We can solve this case:

 ∫vv0dvae−v+bev−(a+b)=t

We can define the function from the latter equation

 V(v0,t)=Log[aeat−aebt−beat+v0% +aebt+v0aeat−bebt−beat+v0+bebt+v0] (43)

We have also

 dVdv0(t)=−beat+v0+aebt+v0aeat−aebt−beat+v0+aebt+v0− −beat+v0+bebt+v0aeat−bebt−beat+v0+bebt+v0 (44)

Now we solve the equation for x:

 dxdt=(b−a)x+f(t)eV−e−Vg(t) (45)

Its solution gives

 x−x0=e(b−a)t∫t0[f(τ)eV(v0,τ)−g(τ)e−V(v0,τ)]e−(b−a)τdτ (46)

Eqs.(42),(46) together define the trajectory of characteristics.

To get the solution for the maximum, we put on the left hand side of Eq. (42). In this way we get a simple explicit equation for as a function of time t.

 aeat−aebt−beat+v0+aebt+v0aeat−bebt−beat+v0+bebt+v0=1 (47)

Putting that solution into Eq.(46), we find the trajectory of the maximum of distribution.

To calculate we should differentiate the expression in Eq.(46) via x. Using the relation , we derive

 1v′=e(b−a)t× ∫t0[f(τ)eV(v0,τ)+g(τ)e−V(v0,τ)]dVdv0(τ)e−(b−a)τdτ (48)

where is defined by Eq.(46) as a function of t.

## Iv CME in multi-dimensional space

Consider now the CME in multi-dimensional space, when we have d-dimensional .

 dP(→X,t)Ndt= K∑nl=−KR→n(→X−→nN)P(→X−→n,t)−∑→nR→n(→XN)P(→X,t) (49)

We get HJE for the function with following Hamiltonian:

 ∂H∂t+H(→x,→p)=0, −H(→x,→p)=K∑nl=−K,1≤l≤dR→n(→x)[exp[−→n→p]−1] (50)

Let us investigate the motion of the maximum of distribution. We denote the maximum of distribution at moment of time t, and assume the following ansatz for the near the maximum of distribution:

 u(→x,t)=−12 lt;→x−y(t)|V|x−y(t) gt; (51)

Putting this ansatz into

 u′′xlt+H′xl(→x,→p)+∑mH′→pm(→x,→p)dpmdxl=0 (52)

Consider the point , and use the identity:

 dpmdxl=−Vml, (53)

We get:

 ∑mVlmdyldt−(K∑nh=−KRn1...nd∑mnm)Vlm=0 (54)

We get the following ODE for the dynamics of the maximum of distribution:

 dyl(t)dt=∑mQmym(t)), Qm=∑1≤h≤d∑−K≤nh≤KRn1..ndnm (55)

## V Derivation of elasticity in 1-d case

Let us get system of ODE for the . We consider the multi-step version of CME in 1-d.

 dP(X,t)Ndt= K∑n=−KRn(X−nN)P(X−n,t)−∑nRn(XN)P(X,t) (56)

Now we have the following Hamilton-Jacobi equation:

 −H(u′,x)=K∑n=−KRn(x)[e−nu′−1] (57)

Now we have to consider the higher terms in expansion of near the maximum.

 u=−V(x−y(t))2/2−T(x−y(t))3/6 (58)

Eq.(55) gives

 dy(t)dt=y(t))b b=∑−K≤n≤KRn(x)n (59)

We also denote:

 a(x)=∑−K≤n≤KRn(x)n2 (60)

therefore

 ˙y=b

With the ansatz by Eq.(56) we have:

 u′′′xxt=−˙V+Tb, u′′′xtt=2˙Vb−Tb2+V˙˙y= 2˙Vb−Tb2+Vbb′ (61)

From the other hand we have differentiating the right hand side of Eq.(50):

 −u′′′xxt=H′′ppV2−2H′′xpV−TH′p=aV2−2b′V−bT, −u′′′xtt=−H′′xpV˙y+H′′ppbV2−H′p˙V+H′pTb2= −bb′V+abV2−b˙V+Tb2 (62)

We derived Eqs.(58),(59) putting .

Then

 −˙V+2Tb=aV2−2b′V, 3˙Vb−2Tb2=aV2b−2b′V

Removing T we get

 2bddtV=−4bb′V+2baV2 (63)

or

 ddt1V=b′1V−2a (64)

Which gives Eq.(24) for 1-step CMW case.

For the multi-step CME we have the following expression for the elasticity:

 −1v′x=(b(x))2∫xx0dya(y)[b(y)]3 (65)

## Vi Conclusion

In conclusion, we calculated exact probability dynamics Eq.(27) for the Chemical Master Distribution using Hamilton Jacobi equation method. Using the methods of characteristcs for solution of HJE, we gave explicit expression for the variance of distribution Eq.(24). The latter is important both for chemical qi08 () and financial applications. We derived the variance also directly from HJE, using a special ansatz. The latter is equivalent to Van Kampen method ka92 (). Using Van Kampen method, in Appendix we give an exact expression for the variance for general 1-d model, as well as a system of ODE to define the variance in 2-d case. Both HJE and Van Kampen methods gave identical results for the variance and the maximum point dynamics of CME: HJE gives directly differential equation for the vasrianc, while the Van Kampen method originally gives a differential equation for the probability distributions, Eq.(A.5), which later proved the differential equation for the variance. HJE gives also exact steady state distribution, and is more adequate for investigation of meta-stabile points es09 () as10 (). Using HJE, we derived exact dynamics Eqs.(44)-(46) for the case of 1-d CME when the rates are linear functions of the coordinate (number of molecules) plus some time dependent functions.

We performed some numerics in case of static noise, also give an analytical criteria for the level of the noise when the bi-stability disappears. Our choice of potential for the averaging Eq.(31) is rather arbitrary (contrary to all other results in this article, which are derived rigorously and are exact). A further investigation of the problem is necessary. Perhaps it is possible to investigate the more realistic case of non-stationary noise xi10 ().

DBS thanks DARPA Prophecy Program and Academia Sinica for financial support.

## Appendix A The calculation of variance dynamics using the Van Kampen method

### a.1 1-d multi-step models.

Consider the CME defined trough the system of equations:

 dP(X,t)Ndt= K∑\lx@utf@l=−KRl(X+lN)P(X+l,t)−∑lRl(XN)P(X,t) (A.1)

Following to Van Kampen, we assume that near the maximum

 X=Nϕ(t)+√Nξ(t), P(X,t)=Π(ξ,t) (A.2)

For the maximum point we will derive an equation (4) later.

We consider terms with different scaling by in the equation:

 dΠ(ξ,t)dt−√Nϕ′(t)∂Π∂ξ= NK∑\lx@utf@l=−KRl(ϕ(t)+ξ√N−1N)P(X−l,t) −N∑lRl(ϕ(t)+ξ√N)P(X,t) (A.3)

Collecting together the terms, we get

 dϕ(t)dt=K∑l=−KlRl(ϕ(t)) (A.4)

The terms give an equation

 dΠ(t)dt=−b′∂∂ξ(ξ)Π+a2∂2∂ξ2Π b=∑llRl,a=∑ll2Rl (A.5)

We derive the following equations for :

 d<ξ>dt=b′(t)<ξi> (A.6)

The initial condition gives .

We get for the variance

 ddt<ξ(t)2>=a(ψ(t))+2b′(ψ(t))<ξi(t)>)2 (A.7)

We can consider ODE via instead of t. Then Eq. (A.4) gives , then:

 bddψ<ξ2>=a(ψ)+2b′(ψ)<ξ>)2 (A.8)

Eventually we get Eq.(24), with .

### a.2 2-d case

We will apply the Van Kampen method, giving ODE to calculate the variance of distribution in 2-d case.

Consider the following system of equations for :

 dP(X,Y,t)dt=−(R1+(XN,YN)+R1−(XN),XN) +R2+(XN,YN)+R2−(XN),XN))P(X,Y,t)+ R1+(X−1N,YN)P(X−1,Y,t)+ R1−(X+1N,Y−1N)P(X+1,Y,t) R2+(XN,Y−1N)P(X,Y−1,t)+ R2−(XN,Y+1N)P(X,Y+1,t) (A.9)

Following to ka92 (), we introduce the fluctuating variables and replace by ,see ka92 ():

 X=Nϕ1(t)+√Nξ1(t), Y=Nψ2(t)+√Nξ2(t), P(X,Y,t)=Π(ξ1,ξ2,t) (A.10)

where give the solution in case of infinite :

 dψ1(t)dt=b1(ψ1(t),ψ2(t)) dψ2(t)dt=b2(ψ1(t),ψ2(t)) b1(ψ1,ψ2)=R1+(ψ1,ψ2)−R1−(ψ1,ψ2), b2(ψ1,ψ2)=R2+(ψ1,ψ2)−R2−(ψ1,ψ2) (A.11)

We can solve Eq.(A.11) and calculate .

Following to the methods of ka92 (), we derive an equation:

 dΠ(t)dt= −∑α(R′α+−R′α−)∂∂ξα(ξαΠ(ξ1,ξ2,t))+ ∑αRα++Rα−2∂2∂2ξαΠ(ξ1,ξ2,t) ≡−∑α∂bα∂ξα∂∂ξαΠ+12∑αaα∂2∂ξ2αΠ (A.12)

where we denoted .

We derive the following equations for :

 d<ξi>dt=b′i(t)<ξi> (A.13)

The initial condition gives .

We get for the variance

 ddt<ξi(t)2>= ai(ψ1(t),ψ2(t))+2b′i(ψ1(t),ψ2(t))<ξi(t)>)2 (A.14)

We can calculate the variance numerically as function of , using the solution of Eq.(A.11).

### References

1. R. Phillips, J. Kondev, and J. Theriot, Physical Biology of the Cell (Garland Science, New York, 2008).
2. D. A. Beard, H. Qian, Chemical Biophysics: Quantitative Analysis of Cellular Systems. Cambridge Texts in Biomedical Engineering. Cambridge Univ. Press, New York (2008).
3. C. W. Gardiner, Handbook of stochastic methods. Berlin, Germany: Springer (1985).
4. Van Kampen N G,Stochastic processes in Physics and Chemistry,(North-Holland,Elsevier Science).
5. T. Lux, Journal of Economic Behavior &amp; Organization 72 638(2009).
6. H. Ge and H. Qian, for a Bistable Biochemical System, PRL 103, 148103 (2009).
7. J. E. Ferrell and W. Xiong, Chaos 11, 227 (2001).
8. J. Keizer, Proc. Natl. Acad. Sci. U.S.A. 75, 3023 (1978)
9. G. Hu, Physica A 136,607 (1986).
10. H. Woodcock and P. G. Higgs, J. Theor. Biol. 179, 61 (1996).
11. E. Baake and H. Wagner, Genet. Res. 78, 93 (2001).
12. D. B. Saakian, Journal of Stat. Physics, 128,781(2007).
13. K. Sato and K. Kaneko, Phys. Rev. E 75, 061909 (2007).
14. D. B. Saakian, O. Rozanova, A. Akmetzhanov, Phys.Rev. E 78, 041908(2008).
15. R. and S. J. Chapman, Euro. Jnl of Applied Mathematics, 16, p. 427(2005).
16. C. Escudero and A. Kamenev, Phys.Rev. E 79, 041149 (2009).
17. M. Vellela and H. Qian, J. R. Soc. Interface (2008)
18. M. Assaf and B. Meerson, Phys.Rev. E 81, 021116 (2010).
19. L. C. Evans, Partial Differential Equations AMS (2002).
20. A. Melikyan, Generalized Characteristics of First Order PDEs, Birkhäuser, Boston (1998).
21. X. S. Xie,pp. 435-448, Single Molecule Spectroscopy in Chemistry, Physics and Biology Nobel Symposium, Edited by Astrid Graslund, et al., Springer Series in Chemical Physics, Springer-Verlag Berlin Heidelberg (2010).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters