Solution of the master equation of a bistable reaction system

Physica 128A (1984) 150-163 North-Holland, Amsterdam

SOLUTION OF THE MASTER EQUATION REACTION Christian

OF A BISTABLE

SYSTEM WISSEL

Fachbereich Physik, Philipps-Unioersitiit, D-3550 Marburg, Fed. Rep. Germany

Received 21 February 1983 Revised 11 April 1984

In this paper several dynamical properties of the stochastic Schlogl model are calculated. A semi-numerical method of solving the eigenvalue problem of the corresponding master equation enables us to determine all interesting quantities to a very high accuracy. The dependence of the metastable state on the various parameters is shown. The transient evolution in a quenching process is investigated. It is shown how the deterministic discontinuous transition is altered by the stochastics. Finally, the dynamical correlation function has been determined in the bistability region. The demonstrated method can generally be used for the solution of a one-dimensional single step master equation.

1. Introduction In recent years there analysis of a trimolecular

has been chemical

a growing reaction,

number

of papers3-*391~~) on the

*+2x23x; *2x. kz

Schlbgl’)

(1)

k4

was the first to investigate

this reaction

scheme

by deterministic

rate

equations. The steady states show a behaviour analogous to first order phase transitions: For certain parameter constellations a bistability is found. Because of the simplicity of these mode1 reactions the corresponding master equation has been used to investigate qualitatively the stochastic influence on discontinuous transitions in homogeneous systems. A general discussion on the validity of a master equation description for chemical reactions can be found in ref. 2. The stationary probability distribution 12-13)for the particle number of the molecules X can easily be found from a recursion relation. The study of the dynamics of the system is more difficult. Several approaches have been tried: approximations for the escape timelg), approximative solutions of the corresponding eigenvalue equation for the bistability region3) and near the critical pain?), the Mori-Zwanzig projection operator formalism for the dynamical 0378-4371/84/$03.00 @ Elsevier Science Publishers (North-Holland Physics Publishing Division)

B.V.

151

MASTER EQUATION OF A BISTABLE REACTION SYSTEM

correlation

function4),

computer

simulations’).

Some

exact

results

have

been

found, too’l). Finally, the analytic expression for the mean first passage time”‘0’2’ ) has been used to get some information on the slow part of the dynamics. There is also an abundant literature on the Fokker-Planck description of the same problem. In this paper a seminumerical method for the solution of the corresponding eigenvalue-problems is presented. This method is well conditioned, fast and very accurate. From that all dynamical quantities of interest could be determined. The approximations mentioned above are compared with this exact result. In section 2 the master equation corresponding to the reaction (1) is transformed so that a suitable numerical method for the solution of the eigenvalue problem can be applied. In section 3 it is discussed how the properties of the lowest eigenvalue leads to the existence of a metastability which is investigated in detail. In section 4 the influence of the stochastics on the discontinuous transition is determined. In section 5 the dynamical correlation function is represented.

2. Master equation The master

y

=

equation

describing

the reaction

(1) can be written

in the form

2 W,,Pmr(t). m=O

Here P,,,(t) is the conditional probability to find 12 particles of the molecule X the transition matrix W,,, can at the time t. Using standard techniques 14*1s722), be expressed by the birth and death rates

W n+l,n= b, = Kq[n/K(nlK

- l/K) + I/a] ,

Wn_l,n = d, = Kn/K[lla(n/K

- lIK)(n/K

- 2/K) + I] ,

(3)

W,,,=-b,-d,,. We have used dimensionless A,

a=

quantities4) (z,“‘, 2

The quantity

K is proportional

K= 3

(s)1’4_ 1

(4

2

to the size of the system

i.e. to the volume

152

C. WISSEL

whereas

the parameters

4 and a are independent

of it. The time t is measured

in units of Ki’. The stationary

p:+1

distribution

P*, is determined

by the recursion

relation’4~15~22)

P~b,ld,,, .

=

(5)

of it can be found This unique distribution is stable “,*‘). The discussion elsewhere.~~.‘*.‘3.‘9 ). For a > 3 and for values of 9 near 1 the deterministic description yields a bistability for the particle number N of the species X (see fig. 4). Here the master equation gives a stationary distribution Pz showing two peaks (see fig. 7). In order to find the time dependent matrix W,,, into a symmetric form A,, nm =

p;-ll2p

cp

nm

The transformed

$& dt

p*1/2 m

solution of (2) we bring the transition 15). To this end we introduce (9)

7

distribution

(P”,~ fulfills a master

equation

c A”,%?n.

=

(11)

r=O

Now the transition

matrix

A,,=-b,-d,, A

nn+l

=

A,+1,

=

@,4,+l)1’2

is tridiagonal and symmetric. eigenvalue equation

c

A,wZ” = r=O The eigenvalues

(12)

This

makes

it easy to solve

the corresponding

h,w(,“). A, are real and because

A,SO.

Since for n > 0 the birth and death coincide’6*21). Therefore the eigenvectors

(13)

of the stability

of P*, it is”) (14)

rates do not vanish no eigenvalues w,” ( ) form a complete orthonormal set.

MASTER

We choose

Then

EQUATION

OF A BISTABLE

REACTION

SYSTEM

153

the numbering

we find

(16)

A,=O, w(o) =

p*w .

n

(17)

”

As we have a tridiagonal symmetric matrix we can determine the other eigenvalues by a very easy and well conditioned numerical method, the socalled bisection method”). The accuracy of the result can easily be estimated. In our case its absolute value is of the order of magnitude lo-“. The eigenvectors w(ny)have been determined by the inverse iteration method18). A cut-off for the particle number it has been used in order to get a vector space of finite dimension. It has been checked that in this way the eigenvalue equation (13) has been solved with a relative error of less than lo-‘. By this method as many eigenvalues and eigenvectors have been calculated as are necessary for a reliable determination of the interesting quantities. In no case more than 20 eigenvalues are necessary and often the first two or three are sufficient. The conditional distribution P,,,(t) can be expressed by

= 2

p,,,(t)

pE1’2w(,y)erJW(,Y)pG-1’2.

(18)

3. Metastability In

fig.

1 the

parameters

lowest

eigenvalues

are

shown

as functions

of

q for

the

a = 7 and K = 20. In the range

q3 = 0.74 < q < 1.35 = q4, where the deterministic description shows the bistability the decrease remarkable, there A, is lower than the rest of the eigenvalues by several of magnitude. For q = 0.89 we find A, = -1.451 x 10m4; A,= -5.150 A, = -7.180 X 10-l; A, = -1.046. Janssen’$) approximative formula value of A, which deviates by a factor of about 1.6 in the middle of the (q3, q4). Near the boundaries the deviation increases considerably. The result

of the inverse

system

size expansion3*4*6)

of A, is orders x 10-l; gives a interval

154

C. WISSEL

Fig. 1. Eigenvalues Al and A2 for a = 7 and metastability are Al = 0.826 and A2 = 1.110.

K = 20. The

boundaries

of the

range

of the

Al - exp(- const K) for larger values of K is confirmed by our results. This asymptotic behaviour is reached already for K = 30 (compare refs. 3 and 4). The asymptotic dependence on the parameter a (see fig. 6) is exponential, too. The enormous difference between A, and A, has the following consequences: There is a time interval where the terms v = 2,3, . . . in (18) are vanishingly small already and where e-*lr equals 1 to a good accuracy. Thus in this time interval p(m) n

we get a time-indpendent =

p*

+ n

which describes

distribution

p*l/2w(1)w(1)p*-1/2 n nmm

a metastable

(19) state.

It is reached

in a characteristic

time of the

order of magnitude A;‘. The time AT’ to reach the completely stable state P*, is greater by orders of magnitudes. Depending on the parameter values (see above) it can become so long that it may go beyond the realizable duration of an experiment. occupied

In this

case

Pz has no physical

relevance.

Its position

is

by P’,“’ 19V21).

In order to find the properties of the metastable state Pim’ we have shown ,+,‘,“P;‘n and ,,,(,l)p;-1” . fj gs. 2 and 3, respectively. m From fig. 3 it is obvious that the form of the peaks in fig. 2 are the same as in PE. Thus the metastable Pim’ has almost the same shape as P*, with the same form of the peaks3720). Only the relative heights of the peaks differ. They are determined by the initial condition i.e. by the factor w(1)P*-1’2. In fig.“h &me consequences for the mean values are shown. Let us start at small values of q and increase q so slowly, that the system reaches the corresponding metastable state. As the starting point m of (19) always is near

MASTER

EQUATION

OF A BISTABLE

.o{~~~-,7--~ n1

60

n2

Fig.2. Eigenvector

t-q

120

~(nl)P:“~for a = 7, K = 20 and 9 = 0.89. The position of the extrema of PE are nr, nr and

n3.

REACTION

SYSTEM

155

_,,j, i, , , _ nl

n2

60

n3 120 n

Fig. 3. Eigenvector w$‘)P:-“~for a = 7, K = 20 and q = 0.89. The position of the extrema of P: are ni, 112and ns.

Fig. 4. Deterministic equilibrium values (full curve) of the particle curve stable, intermediate one unstable). The stable curve coincides metastable distribution PLrn)except for the upper curve for low values the mean value of the stationary distribution Pz. The boundaries of are q, and q2.At qothe duration of the metastability is maximal and same statistical weight (a = 7, K = 20).

number N (upper and lower with the mean values of the of q. The dotted curve shows the range of the metastability the two peaks of P: have the

to the lower peak n, (see fig. 7) the statistical weight of the upper peak remains vanishingly small. The corresponding mean value of the particle number follows almost exactly the lower full curve for the deterministic equilibrium value until at higher values of q the transition to the upper state takes place (see section 4). The corresponding decrease of q delivers the analogous result for the upper full curve but with a small deviation of the mean value at small values of q. If the alteration of q is performed extremely slowly the corresponding stationary state P: is reached which gives the dashed curve in fig. 4 for the mean value. It follows from fig. 1 that for low and high values of q we find the same order

156

C. WISSEL

I

I

-

1

5

3

a

7

Fig. 5. Boundary values of q. Between q3 and 44 deterministic there is metastability. At 90 the duration of the metastability the same statistical weight (a = 7, K = 20).

bistability is maximal

exists. Between 41 and q2 and the peaks of PE have

IOOO-

500-

7

5 Fig. 6. Exponential

increase

of the maximal

duration

a of the metastability

(K = 20).

of magnitude for A, and A,. Here no metastability exists. The boundaries of the range of metastability are q1 and q2 which are given in fig. 4, too. The maximal duration

t,-

t, of the metastability

is achieved

at qo. If it reaches

high values

it

equals A;‘. Therefore t2- t1 increases exponentially with K and a (fig. 6). In addition q,, labels the position of a sort of Maxwell construction3). That means that here the statistical weight of both of the peaks of Pz are equal. In the limit K + 00, fig. 4 remains almost unchanged. Only the dashed curve approaches a step function at q,,. In fig. 5 the various q-values discussed above are plotted versus a. Only for a 3 5.2 where the bistability is pronounced a metastability is found. Next we discuss a change of the parameters which is equivalent to a quenching process in a thermal system with a first order phase transition. We start with the stationary distribution Pz (single peak in fig. 7) at a = 2, q = 1.2 and K = 20, i.e. in the range where the system shows a normal stability

MASTER

nl 10

EQUATION

“2 40

OF A BISTABLE

90 \“j

REACTION

SYSTEM

157

n

Fig. 7. Particle number distribution P,(t) for a “quenching” process. At I = 0 we start with the stationary distribution P: for a = 2, 4 = 1.2 and K = 20. The other full curves show P,(t) at a = 7, 9 = 0.89 and K = 20 for t = 2; 2.5; 3.5; 5; 7.5. The heights of the peaks increase in the course of time until the metastable state (highest peak at nj) is reached. The dotted curve shows the stationary distribution P:. The positions of the extrema are nl, n2 and n3.

.6-

.5-

Fig. 8. The probability to find the system in the lower state around nl for a “quenching” process: the initial distribution is Pz at a = 2, 4 = 1.2 and K = 20. The evolution of the distribution takes place at a = 7, 9 = 0.89 and K = 20. Between t, and h the metastable phase appears. At t3 the stationary state P: is reached.

behaviour (see fig. 5). Then we “quench” the system by an instantaneous switching to the values a = 7, q = 0.89 and K = 20. Using 20 eigenvalues we can calculate the probability distribution for all values of the time t except for the extremely small ones. In fig. 7 the full curve with the highest upper peak

C. WISSEL

158

distribution P,cm)which is reached at the time t, = 7. For t > t, = 210 the system leaves this state and the stationary state PE (dotted line) is reached at t, = 2420. The area under the first peak gives the probability a,

shows the metastable

that the system

is in the lower state.

In fig. 8 the values

of a, are plotted

versus

time for the quenching process. Bearing in mind the suppressed zero point of a, and the logarithmic time scale one can see the metastable phase between t, and t2.

4. Discontinuous

transition

One of the most interesting problems is the question how the deterministic transition behaviour is changed by the stochastics. When in the deterministic description the value of q is slowly increased the value N of the particle number follows the lower curve in fig. 4. Just at the value q4 the discontinuous transition into the upper state takes place. For the stochastic description first of all one has to define how to increase the value of q in the course of time. A possibility is to hold the value of q fixed for a certain time interval At. One has to investigate the probability W(q) that in this time interval the transition takes place. If this is not the case the value of q is enlarged by Aq and the transition probability for this new value of q is investigated. This step-wise increase of q (see fig. 9) can also be used as an approximation of a linear increase when the step length is sufficiently diminished. For the calculation of the transition probability W(q) one needs the first passage time distribution F,,(t). It is F,,(t) dt the probability that the system will reach the particle number u the first time in the time interval (t, t + dt) presupposed that it starts with the particle number m at the time t = 0. It

/

-At-i

Fig. 9. Linear

and step wise increase

of q for the investigation

of the transition

behaviour.

MASTER

EQUATION

OF A BISTABLE

REACTION

SYSTEM

1.59

fulfills the backward equationt4.*‘)

(21) with the boundary conditions

(22)

F,,,(r) = 0 7 F,,,(O) =

LA-,.m .

(23)

We introduce the absorptive boundary condition d,=O.

(24)

Since we are asking for the first passage only the further fate of the system is of no interest once the state u is reached. Therefore the first passage results are not affected14) by the insertion (24). From the definition of F,,,(t) it is obvious that

where p,,,n is the solution of the master equation (2) with the additional absorption condition (24). The validity of (25) is simply checked by the fact2*,14) that p,,,n fulfills also the backward equation (21) for all n < U. The initial condition P,_,,,(O) = CYS~_~,~ yields (23). The absorption condition (24) gives P,_,,, = 0 leading to (22). From (25) and (18) we get

F,,,(t) = b,_, c PE!{2G,1”1, eiJG~)P~-‘“, v=o

(26)

where the tilde indicates that the absorption condition (24) is used in the eigenvalue equation (13). There we have to drop (16) and (17) as there is no stationary distribution in this absorptive case. Nevertheless, we can use the transformations (9) and (10) together with (5). Then P: is just a formal expression used for the symmetrization without the interpretation of a stationary state. A simple check of the formula (26) is the calculation of the first and the second moment of the distribution F,,(t). We have compared these results with the analytic formulas for these moments~‘0~‘4’“). The relative deviation is less than 10e4 in all treated cases.

160

C. WISSEL

We like to emphasize

that the mean first passage

AT’ as one may expect. distribution small

decays

region

We have found

exponentially

at small values

with

time does not coincide

with

that for q = q. the first passage

time

a decay

of t (compare

time

i-’

ref. 8). Thus

disregarding

a very

A;’ is the mean

first

passage time*l) which exceeds the time ATi to reach the stationary state PE by a factor of about 2. This surprising result is elucidated by the fact that for the exponentially

decaying

first passage

time distribution

in more than

60 percent

of all cases the transition is performed in a shorter time than the mean value Ai’. The inverse mean first passage time A, is of the same order of magnitude as A, for q = q0 only. We are interested in the transition from the lower state at the lower peak (see fig. 7) to the upper one. As we have found that F,,,(t) does not remarkably change when u varies in the range of the upper peak we can take u = n3 in the following. We start our investigation with a low value of q. We take the corresponding metastable distributionPkm)(q) with one lower peak only as our initial distribution. The probability W(q) that the transition to the upper state takes place in the time interval (0, At) is

W(q) =

cf n

~,Jt~ W’,“‘(q)dt.

(27)

0

At the end of the time interval (0, At) the probability distribution p”(At, q) with the presupposition that the value It3 has not been reached can be calculated analogously to (18) but with the absorption condition (24). p”(At, q) is the initial distribution for the next time interval (At, 2At). There the transition probability

W(q

is

+ &I=

c n

I

R,&

4 +“d!&Q

q)dt.

(28)

0

We continue in this introduce the probability T(q)Aq

= W(q).

time intervals. Finally, manner in the following distribution T(q) for the transition per q-interval

we by (29)

In fig. 10 this distribution T(q) is shown for different increases of q. In all shown cases the step lengths Aq and At are so small that a further diminishing of them does not change the probability distribution T(q). That means we have practically the infinitesimal limit, i.e. the linear increase of q (see fig. 9).

MASTER

EQUATION

OF A BISTABLE

REACTION

SYSTEM

161

Fig. 10. Probability distribution 7’(q) that in the interval (q, q + dq).the transition into the upper state takes place (n = 7, K = 20). The value of q is increased as shown in fig. 9. From the left to the right it is Aq = 0.01, At = 1000;Aq = 0.01, At = 100; Aq = 0.01, At = 40.

In the deterministic description one would have a S-peak at q4 for the transition probability distribution T(q). The stochastics broadens T(q) and shifts it to lower values of q. The broadening is the larger the quicker the increase of q is. The slower q increases the more the distribution T(q) is shifted to lower values of q. Thus T(q) depends strongly on the way q is increased.

5. Correlation

function

is often used to describe

of

physical system4). It is defined by C(t)

(N(t)N(O)) - (Ivy .

Because of (18) it is C(t) = 2

v=l

e*p’ C nPE’“w%)

[n

2.

1

(30)

In this expression we sum over the first 15 eigenvalues only. In order to see that this is a reasonable approximation we consider the case t = 0. There C(0) equals the variance which can be calculated with the help of (5). The comparision with the truncated sum (30) gives a relative deviation of less than 10Y3 in all investigated cases. As the error by the truncation decreases for increasing time t the approximation is satisfactory. The resulting correlation function has been fitted by C(t) = C(O)[(l - (Y)e-‘lt + ff e-“*‘I .

(31)

162

C. WISSEL

Fig. 11. Fitted relaxation frequencies A, and A, and weighting factor for a = 5 and K = 10. The boundaries of the range of bistability are the full curves it is -Al = h, and -A,= AZ.

In fig. 11 the fitted

values

of the relaxation

(Y of the correlation

function

q3= 0.857 and q4= 1.167. For

frequencies

A,, A2 and

of the

weighting

factor (Y are shown for the case a = 5 and K = 10. In the range q ~0.9 the fit shows that the first two terms of the sum (29) alone give sufficiently well the whole correlation function. The remaining terms give a vanishing contribution only. Therefore we have -A, = A, and -A2 = A, in this case. In the range q > 0.9 the fit gives -A,= A,. But A2 differs from all eigenvalues. Several terms of the sum (29) contribute. Eq. (31) gives an approximative description only with an effective relaxation frequency A,. For high values of q the weighting factor (Y approaches 1. This means that the slow part with -A, = A, does not contribute to the correlation function in this range. The comparison with the correlation function calculated by the MoriZwanzig projection operator formalism4) gives a satisfactory agreement of the relaxation frequencies A, and A, for small values of q only. For large values of q the agreement is very rough only. In the range of metastability and near its boundaries q1 and q2 a pronounced disagreement occurs. The projection operator formalism produces too small values of LYand too high values of Ar. This discrepancy grows for increasing a. In the case of a well pronounced bistability the Mori-Zwanzig formalism generally fails whereas Janssen’s approximation3) is good just here. But near the boundary of the metastability at q, and q2 where the first order transition takes place both of these approximations fail.

6. Conclusion In this paper it has been demonstrated that one-dimensional master equations of the birth and death type can be transformed so that the corresponding

MASTER

EQUATION

OF A BISTABLE

REACTION

SYSTEM

163

eigenvalue problem can be solved by a rather simple numerical method. In this way the complete time dependence of all quantities which are related to the master equation can precisely be determined. This is of special interest for birth and death processes in theoretical biology and in chemical reactions. The solution of the eigenvalue problem for the stochastic Schliigl model has been used to calculate the dynamical properties of the “first order type phase transition” and to compare them with existing approximative results. The dependence of the properties of the metastable state on various parameters have been discussed. The transient evolution in a “quenching” process has been investigated. The influence of the stochastics on the discontinuous transition has been determined. The dynamical correlation function has been calculated, too. Unfortunately we see no way to extend this method to problems of higher dimension.

Acknowledgement

The author is grateful to B. Sonneborn-Schmick

for helpful discussions.

References 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15) 16) 17) 18) 19) 20) 21) 22)

F. Schlbgl, Z. Phys. 253 (1972) 147. D.T. Gillespie, J. Stat. Phys. 16 (1977) 311. H.K. Janssen, Z. Phys. 270 (1974) 67. S. Grofimann and R. Schranner, Z. Phys. B30 (1979) 325. E. Ebeling and L. Schimansky-Geier, Physica 98A (1979) 587. I. Oppenheim, K.E. Shuler and G.H. Weiss, Physica 88A (1977) 191. D.T. Gillespie, Physica 95A (1979) 69. D.T. Gillespie, Physica 1OlA (1980) 535. V. Seshadi, B.J. West and K. Lindenberg, J. Chem. Phys. 72 (1980) 1145. B.J. West, K. Lindenberg and V. Seshadi, J. Chem. Phys. 72 (1980) 1151. D.T. Gillespie, J. Chem. Phys. 76 (1982) 3762. I. Matheson, D.F. Walls and C.W. Gardiner, J. Stat. Phys. 12 (1975) 21. G. Nicolis and J.W. Turner, Physica 89A (1977) 326. N.S. Gael and N. Richter-Dyn, Stochastic Models in Biology (Academic Press, New York, 1974). H. Haken, Synergetics. An Introduction (Springer, Berlin, 1977). J.H. Wilkinson, The Algebraic Eigenvalue Problem (Clarendon, Oxford, 1965) chap. 5, sec. 37. Ref. 16, chap. 5, sec. 39. Ref. 16, chap. 5, sec. 53. N. van Kampen, Suppl. Progr. Theor. Phys. 64 (1978) 389. K. Matsuo, J. Stat. Phys. 16 (1977) 169. S. Dambrine and M. Moreau, Physica 106A (1981) 559, 574. C.W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1983).