BULLETIN OF MATHEMATICALBIOLOGY VOLUME 38, 1976

STOCHASTIC EQUILIBRIA FOR MULTIVARIATE REACTION SYSTEMS

9 JoHn P. MVLLOOLY Health Services Research Center, Kaiser Foundation Hospitals, 4610 S.E. Belmont Street, Portland, Oregon 97215, U.S.A.

The equilibrium distribution for a general Qth-order multivariate reaction system is studied. The state transition i n t e n s i t y m a t r i x is developed a n d examples are given for small n u m b e r s of reaction components. A closed-form expression for the equilibrium distribution for systems which are symmetric with respect to the order of component reactions is presented. Numerical examples for three component systems are discussed.

1. Introduction. A summary of various stochastic approaches and application to chemical kinetics has been given by McQuarrie (1967). Stochastic kinetic models have also been used to describe the dynamics of small experimental biological populations (Bartlett, 1971). The equilibrium dynamics of multicomponent systems containing small numbers of particles is also of interest in microbiology and ecology. In this paper we are concerned with the equilibrium behavior of closed multicomponent systems. Chemical terminology is used in this discussion. The terms particle and reaction generalize to system unit, biological or physical, and transformation process respectively. Consider a closed reaction system consisting of v reaction species, X~, i = 1, 2 . . . . . c, which undergo transformation reactions of the type rtX~ rjX 1, i ~ j. Letting Y~ denote a reaction unit of the ith species, i.e. r~ particles of Xt, the set of allowable reactions may be written Y~ ~- Yj, i r j. In the development that follows it is assumed that the number of particles in the system is equivalent to an integer number of reaction units. Denote by N~(t) c 2Vt(t ) = /~, the constant the number of reaction units of type i at time t; ~'~4~1 total number of reaction units of the c types. The reaction rate for the trans597

598

J O H N P. MULLOOLY

formation of a reaction unit of t y p e i into a reaction unit of t y p e j, Yt --> YI, when there are Nl units of t y p e i present in the system is denoted byf~l(N d. In general, the reaction rate fij(Ni) is considered to be a non-negative, nondecreasing function of N t . The results presented in this paper concern Qthorder reo ctions which are defined byftj(Nt) = fq. N~Q,~. The reaction rate when N~ = 1 is given b y the rate constant fij, a non-negative constant independent of time. The order of the reaction, Qij, a non-negative constant independent of time, specifies the power function dependence of the reaction rate on Nt. The stochastic model t h a t we formulate for this reaction system is a finite continuous time Markov process. Such a formulation requires t h a t we enumerate the t o t a l i t y of m u t u a l l y exclusive states t h a t the system m a y occupy and t h a t the state transition intensities be specified in terms of the reaction rates f~j(Nl). States of the system at time t are the allowable values of the vector variable (Nl(t), N2(t) . . . . . Nc(t)), where ~ = 1 N i ( t ) = N. The n u m b e r of possible states is given b y the number of combinations of (N + c - 1) items chosen (c - 1 ) at a time, (Nc+Cll), which we denote by S. Denoting the j t h state of the system b y the vector S(J) = (n~i), n(~), 9 9 ", n(J)~c~, we index the states b y the following inductive definition : (i) S(1) = (n~ 1) = 0, /t~1) = 0 , . . . ,

n~l_)1 = 0, n~ 1) = N ) .

(ii) if n~j) => 1, then n7 +1) = n7 ), n~1+1) = n~J), 9 9 ", n(J+l)c-2 = n ~ 2 , o(i+1),Oc_1 ~-~ ao(l) ,re_ 1 .*. ~ 1, n~J+1) = n(l) _ 1

--~1) . > (iii) if n~/_)Z = 0, 1 = 0,1 . . . . . /c, for some k > 0, and '%-k-1 = 1, then n ~ + l ) = nT), n~l+l) = n~j) n 0c -'k+- 1l )

_-- 0, ~(i+1) _-- 0, '~c-k

""

'

"'

~(]+1)

" " "' " ~ c - k - 3

=

~ (J)

~c-k-3'

n~/_+1) = O, n~ 3')

~--

_- n(1)

~ (]+1) "~c-k-2

~(J) 're-k-l--

c-k-2*

a_ 1 '

1"

The intensity of the transition of the system from state S(~) to state S(J), i # j, is denoted by aij. Since the reaction rate f~l(Ni) depends on the number of reaction units of t y p e i, and not on time t, the state transition intensities are also independent of t. The probability t h a t the reaction system is in state j at time t + At given t h a t it is in state i at time t is, for small At, P[S(J) at time t+ ht ] S(~) at time t] = a~j. At+O(At), i # j. The event of two or more state transitions occurring in At has a probability of order (At) 2 or higher and is assumed zero if At is taken sufficiently small. State transition intensities a~i are defined in terms of the reaction rates fij(N d in the n e x t section. I t is known t h a t the state probabilities, P~(t) = probability t h a t the reaction system is in state i at time t, are related to the state transition intensities by the following set of S linear differential equations

STOCHASTIC EQUILIBRIA FOR MULTIVARIATE REACTION SYSTEMS

599

with constant coefficients: d s d--t Pl(t) = ~ Pt(t).a~j, i = 1, 2 . . . . . S, where alj = - ~ aj~. In order to solve this sot of differential equations for finite t > O, the initial conditions PI(0) must be specified for j = 1, 2 . . . . . S. Letting P(t) be the vector of state probabilities at time t, and denoting the S • S matrix of state transition intensities by A t ( N ) , the set of differential equations may be expressed in matrix form dR(t) -

dt

(1)

P(t)Ac(N).

This Markov process has a limiting distribution as t --> oo, which is independent of the initial distribution; the limiting distribution defines the stochastic equilibrium state of the model reaction system and is the subject of the discussion that follows. At equilibrium (d/dt)P(t) = 0 and the vector of equilibrium probabilities is given by the solution of P A c ( N ) = O, subject to the constraint that ~ = 1 PI = 1. The reader is referred to Chapter 6 of Chiang (1968) for a general discussion of Markov processes. 2. The transition Intensity Matrix. When c = 2 the transition intensities (aij) for the N + 1 states S(J) are given by 1

-f~l(-x) f12(1)

0

2

f~(N) -- [ f l 2 ( 1 ) + f 2 1 ( N

3

4

o

o

-- 1 ) ] f 2 1 ( N - 1)

f12(2) - [f~2(2) + f 2 1 ( N - 2)]f21(N - 2)

(2)

A2(N) = .

N-2 N--1

N-1

f 1 2 ( N - 2) - [f12(iv - 2) + f2~(2)]

N

0

N+I

0

N f~(2)

f ~ 2 ( ; v - 1) - [ f ~ 2 ( ; v - 1) +f~1(1)]

0

N+I 0 f21(1)

f12(.~) --f12(s

In the discussion that follows space is conserved by omitting the diagonal terms from transition intensity matrices.

600

JOHN P. MULLOOLY

Development of the state transition intensity matrix, A3(N), is aided by partitioning the totality of (N + 2)(N + 1)/2 states into N + 1 sections of states according to the number of reaction units of type 1. Section 1 states, having 0 units of type 1, are S(I), S(2). . . . , S(lV+l). Section 2 states, having 1 unit of type 1, are S(N+2), S(~r . . . . . S(2N+1). Section ( N + I ) consists of the single state S(s). In general, the /cth-section of states consists of states S(Lk), S(~+1),..., S(Uk), where Lk = 1 + (/c - 1 ) ( 2 N + 4 - k ) / 2 , U~ = k. ( 2 N + 3 -/c)/2, for k = 1, 2 . . . . . N + I . Since multiple transformations within infinitesimal At have probability zero, state transitions are possible only between states in contiguous sections. The general form of A3(N) for any set of transformation rates {fr and any N is seen to be: I 1

2

]-

,~ . 9 N N+I ,V+2N+3. 9 2N 2Nr 2N+2 ~N+3. . 3N-I 5N f 32 (N)

f ( )_\,

~

f3~ (N)

~

f~(l')'\

~

,

2,\\\ .-,.\\ ,,...\C\\\ ;v J ~ \ , "-f_(I)l~ ""\

f "~'(1)

N+'I (~

,.~r

"'~(N)

I L.P

0

L~_ E Lk.,+1 "" &~l-I Uk_ I L N L~+I'" Uk-I U/r LIf

L~I

f(k-I ) f(k- I ) r3 12

\

\ \\\\

\ \\%

\\\

Uk-I

e,(;-2

f(N-K+I) 32

I

. \N, N

\\x."

s-,

is

r,jN)

I

(N-~+J) f3E \

(3)

%\

N\\

x\\

L,~, Lk.rl'" Uk.,-I U~+,+I

"~NN

\\\

f(1)

N\ N

x"""x"~ \x-.. ~i )

Observe from Aa(N) and A2(N) that the submatrix for transitions within the kth-section of As(N) is identical to the matrix A u ( N - k + 1) with subscripts 2 and 3 replacing 1 and 2, respectively. In general, the submatrix for transitions within the kth-sections of Ac(N) is identical to the matrix A(c_x)(N-k+l) with subscripts 2, 3 . . . . . c replacing 1, 2 . . . . , c - l , respectively. Hence, denoting the transition intensity matrix for transitions between sections k and k' b y A~g')(N), we may write At(N) in the partitioned form (4). Each section of states m a y be further partitioned into subsections depending on the number of units of type 2 present. This nested classification m a y be further partitioned according to the number of type 3 units present, and so on. We now give a detailed discussion of the construction of A4(N). From expressions (4) and (3) we see that the construction of A4(N) is completed b y speci( ~ - 1 ) , k = 2, 3 , . . ", N + I . fication of A~(~+I)(N), k = 1, 2,. 9 ., N, and ~~4(N)

STOCHASTIC E Q U I L I B R I A F O R M U L T I V A R I A T E R E A C T I O N SYSTEMS

A(~_O (N)

601

A c~2!(tV)

(~)

A(=_O ( N - I )

,4(~_o (N-2) ,~4)(N

Ar (N)=

e,

~

.

The kth section of A4(N)is partitioned into subsections 1, 2 . . . . , ( N - - k + 2 ) having N~ = 0, 1 . . . . , N - k + l, respectively; the jtb-subsection of the kthsection containing ( N - k - j + 3) states. The composition of the states contained in the jth-subsection of the kth. section are : N1

N2

N3

k-1 k-1

j-1 j-]

0 ]

k-I

j-1

N-k-j+2

N4

N-k-j+2 N-k-j+l 0

where k = l, 2 . . . . , N + l , j = l, 2 . . . . , N - k + 2 . Since only a single transition may occur in infinitesimal At, the allowable transitions from section k to section (k+ 1) are seen to be from subsections j and (j + 1 ) of section k to subsectionj of section (k + 1). The submatrix of A~(~+I)(N) giving the intensities for transitions between the states of jth-subsection of the kth-section and the states of the jth-subsection of the (k + 1)th-section is given by {5). The only other possible transitions from states in the kth-section of A4(N) to states in the (k+ 1)th-section are from the ( j + 1)th-subseetion of section k to the jth-subsection of the (k+l)th-seetion. The diagonal submatrix (6) displays these transition intensities.

602

JOHN P. MULLOOLY

g

jth-subsection 2

3

1 f (N-k-j+2)

0

0

0

0

0

41

"5

2 fsl(I)

.8_

of ( k + l ) t D section

I

i

3

0

"4jf(N-k-J+l) f(2)

(N-k-j+l)(N-k-j+2)

..

f(N-k -j)

0

(5)

0

g I

---.,/V- k-i+ 2 N-k-j+3

~ " ~ /~sf-k-J+l)x~ f( t ) 4 0

0

f(N-k-j'l-2)

-Sl

I2

]l"h- subsecfion r

I

I ,ac

of (k, I)t-h-section

2

3

., ( N - k - j + l )

(N-k-j+2)

f(]) 21

"6

(J)

@ ~

0

f 21'~.,,,

(6)

I

"7 N-k-]+l

C

.,... f(j) 21

N-k-j+2

f (21j )

This concludes the definition of A~ (~+1)(N) and we now discuss the specification of A~(k-1)(N), # = 2, 3 , . . . , N + I , which will complete the definition of

A 4(/Y). The following submatrix of A4~(k-z)(N) gives the intensities for transitions between the states of the jth-subsection of the #th-section and the states of the jth-subsection of the (k-1)th-section.

I

g

j t h - subsection of (k-I)th section 2 3 (N-k-j+3) (N-k-j+4)

f(k-D 14

.

f(k-I)

.

0

i5

f

0 "6

.g_

==

4

k\

k',, %,,, \ kk \

g

s N+j-2 N-k-j+3

(7)

k% k\ k %\ -,, ,,. \\--(k -I ) %\-.(#-1 )

0 r,4 0

qs

(k-I)

f14

0 (k-I)

fl3

STOCHASTIC E Q U I L I B R I A FOR MULTIVARIATE REACTION SYSTEMS

603

The other possible transitions from the kth-section to the (/c-1)th-section are from the jth-subsection of section k to the (j + 1)th-subseetion of section (/c- 1). These transition intensities are seen to be:

(j+ I)th-subsectionof (k-I)l"h-section I

2

5

. . (N~-j+2)

(N-k-j+3)

f(k-I} 12.

flk

-1)

C (8)

I

f.

N-k-j+2

C

xxx'-f r I} f~-t)

N-k-j+3

This concludes the specification of A~(~-I)(N), k = 2, 3 . . . . . N + I , which completes the definition of A4(N). Expression (4) may then be employed for the definition of A 5 (N), although specification of the off-diagonal rectangular submatriees A5k(k+l)(N) and A~5(k-1)(N) is more complex than for the c = 4 case.

In a later section a closed form equilibrium solution is presented for systems which are order symmetric, i.e. Qil = Qi, which allows us to study the model behavior of c-component systems without having to generate the transition intensity matrices Ac(N). In the Appendix, A3(N) is used for the analytic verification of the closed form solution. Direct numerical solution employing A4(N) and the closed form solution gave identical results for four-component order symmetric systems. 3. Solution of PAc(N) = O. Substituting P s = 1 - ~ j ~ l Pt into P A t ( N ) = 0, we obtain S inhomogeneous linear equations in ( S - 1 ) unknowns. However, using the fact that ~ t aij = 0, i = 1, 2 . . . . . S, it may be shown that the Sthequation is redundant, being equal to the sum of equations 1 through ( S - 1). Hence the equilibrium distribution is given by the solution of the following set of ( S - 1) equations in ( S - 1) unknowns:

(a12 -as~)

(a~2 -as2)

9

(a(s-1)2 -as2)

P2

(a(s-1)(s-1) - as(s-l))

P -1

9

[ (al(s-1) - as(s_1)) (a2(s-t) - as(s-l)) 9

[aSl

]

604

JOHN P. MULLOOLY

The IBM scientific subroutine SIMQ was used for the numerical solution of this set of simultaneous linear equations. The utility of this method of solution is limited to very small N (say < 20). When N = 100, there are in excess of five thousand states for a three-component system so that an array larger than 25 million is required for the transition intensity matrix. For these dimensions, the matrix inversion required for the solution of (9) is clearly infeasible. Whereas a value of N in the range of hundreds or thousands is small for the size of a biological population and perhaps nondetectable as a chemical system, it is mathem~otically very large in the contest of the discrete state, continuous time, birth-death process under investigation. A closed-form solution for the equilibrium distribution would enable us to study the system behavior for values of N encountered in applications. We have been partially successful in achieving this objective and present in Section 5 a closed form solution for systems which satisfy a certain symmetry condition on the orders.

4. The Deterministic Analogue. Treating the number of reaction units of type i, :Vt(t), as a non-random real-valued function of t, the equilibrium composition is given b y the solution of the following set of algebraic equations: dN 1 dt

c

Xc f,J NQtl .

,=1

. ,=1

tcj

l~j

.

.

.

o, j.

. 1,

2,

c,

subject to the constraint that Z~=I Ng = N. Putting Nc = N-)~gc__-~ Nj, it is seen that the deterministic equilibrium solution is found by solving the following set of ( c - l ) equations in ( c - 1 ) unknowns:

c-1

[

~ftg.NQ'J+fcg. N I=1 tcJ

c2 N,

]Qoj o-1 - ~fl,.NO"-flo.NQ'c=

1=1

O,

/=1

l#J

j = 1,2 . . . . , c - 1 .

(10)

When the kinetics are first-order, i.e. Qtj = 1 for all i # j, i, j = 1, 2 . . . . . c, the set of linear equations has an explicit closed form solution. Before discussing the deterministic first-order solution for general c, the particular cases of c = 3 and 4 will be given. A schematic with nodes corresponding to reaction species, and connecting arrows indicating reactions will assist comprehension of the general formulation. Defining

f l = f21 "fs1 +f32 "f21 +f23"f31,

fz = f12 "f32 +f13 "f32 +f31 "f12,

= f23-f 3 +fx2.f 3 +f

'f13,

STOCHASTIC EQUILIBRIA FOR MULTIVARIATE REACTION SYSTEMS

605

Q

f~2

the deterministic equilibrium solution is seen to be N1 = N . F 1 , N2 = N . F 2 , N3 = N . F 3 ,

where F~ - ~ , -

i = 1,2,3.

l=t

In the general case of c types of reaction species, f~ is given by the sum of the products of ( c - 1 ) reaction rates corresponding to paths of ( c - 1 ) segments which lead to node i and have no closed loops. Details of the c = 4 case further illustrate which paths are to be included in the general definition of composite rate constants fi-

f43

f32

fL~

Here, there are twelve configurations of length 3 having no closed loops; eight of which are without branches, while four have single branch points.

606

JOHN P. MULLOOLY

The sum of products of rate constants over paths leading to component 3 are seen to be

f3 = I23.114 .f43 +f21-f14.143 +A1 .f12 ":23 +A3 "f12-A3 +f14 .f42 "123+f41.113 "123+/43 ./21 "f13+f12 ./24 "f43 -]-/14 "f24 "f43 "}-f41"f21 "f13 +f12 "f42 "f23 -]-/13 "123"f43The remaining composite rate constants are defined in a similar manner. With this general definition o f f , , the deterministic solution for a first-order system having c reaction species and a total of N reaction units is given by N~ = N.F~,

whereF~ = -f~ ,

i = 1,2,...,c.

(11)

t=1 In the non-linear case, i.e. some Q~j r 1, the set of simultaneous algebraic equations (10) must be solved numerically. The numerical results presented below were obtained b y iterative numerical solution using a Newton-Raphson procedure. For first-order stochastic reaction systems, all reaction rates are first order, QIj = 1, and the reaction units are mutually independent, so that the equilibrium distribution is multinominal, 5. Closed F o r m Solutions for Stochastic Equilibria.

Pr[N1 =i1,N2 =i2,...,Nc=ic]

=

-1 "-2 - . . F; '

i l I i 2 ! . . .ic!

(12)

0 < i I < N=, c ~I 9 = N , where F j is the function of rate constants defined in ~t=1 (11). This result is 9 consequence of the probability mechanism defining the multinomial distribution, i.e. N independent assignments to c mutually exclusive and exhaustive categories. Note that the mean value of the stochastic process is equal to the deterministic solution for first-order systems. It is also to be noted that upon specializing the transition intensity matrix A t ( N ) to the first-order case, the multinomial probabilities (12) may be shown to satisfy (9). In a recent paper, Mullooly (1975) discussed the c = 2 case, with Q12 general and Q21 = 1, and established the result that the equilibrium probability distribution is given b y

Pr IN1 = ~V-j, N2 = j] oc

[(AT_j) !]Q,~.j!

The argument used there can easily be extended to the ease where both Q12

STOCHASTIC E Q U I L I B R I A F O R M U L T I V A R I A T E R E A C T I O N SYSTEMS

607

and Q21 are general non-negative constants, Pr [N1 = N - j ,

= j]

(13)

We now address the question of the generalizability of this result to higher dimensions. Note that each reaction species contributed a factorial term to the denominator of (13) and that the power by which a factorial term is raised is the order of the reaction which acts to decrease the level of that reaction species. The above observation suggests that (13) may be directly extended to a reaction system with c types of species for which Qij = Qi, j = 1, 2 . . . . . c, j r i. This symmetry condition requires that there be a single order for transformations of species i; it does not, however, require that the reaction rates be equal, but allows for differences in the rate constants flj, j = 1, 2 . . . c, j r i. For the c = 3 symmetric system, Qij = QI, it may be verified analytically that the probability function 1 ~2 9 9 9Fc Pr [N1 = il, N2 = i2 . . . . . N c = ic] =

( i l !)Q'(Q!)Q~ . . . (ic!)Q~

Z Z . . . Z.

F

'.C ...

( i l !)Q'(ie!)Q~ . . . (ic!)Qo c

0 < ij < N,

~

ij = N

il+Q+

...

(14)

+ic = N

J=l

satisfies (9). Details of this verification are given in the Appendix. Furthermore, identical numerical results were obtained by this expression and the numerical solution of (9) for c = 4 over a wide range of values of {f~j}, (Q~}, for N = 5, 10. Although these results do not constitute a rigorous mathematical proof for (14) for general c, they do provide strong evidence for its validity. In evaluating (14), for even modest values of N, special care is required to avoid overflow of terms involving products of factorials. Hence, it is useful to develop an approximate expression which is calculated more conveniently than (14). Expression (14) may be written a~ : H(il,

Pr [N1 = il, Nu = i2 . . . . . N c = ic] = il+iu+

i2 . . . . .

. . . ZH(il,

i2 .....

...

= N

c

where H(il ....

ic) =

c

~j

YI I-[ ( i j - / + I ) Q j 3=1 l=1

ic)

+ic

ic)'

608

J O H N P. M U L L O O L Y

Taking the logarithm of H, we have: e

lnH=

c

Zis'lnF

s-

1=1

~j

Z Qs Zln(i l - k + l ) , t=1

k=l

~j

which, when the sum ~ = l In (iI - k + I ) is approximated by the integral

~[+ll~ In (it - Y + 1) dy = - 1/2 (ln 1/2) - i 1 + (iS+ 1/2) In (iS+ 1/2), /2

gives

H(il, i2 . . . . . ic) - E X P

1/2 (ln 1/2)

is" (Qs+ln FS)

Qs "EXP 1

Qs.(is+ll2).ln(is+ll2)]. -

(15)

Dropping the constant term, the joint probability function of (Nz, N2 . . . . . Nc) m a y be approximated by Pr [N1 = il . . . . . Nc = ic] EXP

i1" (Qs + In FS) - ~ Qs" (is+ 1/2) In (is + 1/2)] 1

N

N-ll

2..-

I1=0 12=0

ij.(Qj+lnFs)-

EXP l ~ - 1=0

Oj.(is+l/2 ) In (is+1/2) 1

where c-1

ic = N -

~ i i.

(16)

1

The form

EXP[~ is.(Qs+In F s ) j=l

~ Qs.(is+ l/2) ln (ii+ l/2)] J=l

is convenient for avoiding overflow problems. The second term is always negative while the sign of the first depends on the values of Q, and Ft. For particular values of Qj and F i , the range of values of the bracketed expression may be investigated. I f very large positive values are attained which would cause exponential overflow, an appropriate constant can be subtracted from the bracketed term and the resultant values compared against a large negative number, say -100. If the exponential argument is less than -100, the unnormalized solution can be put equal to zero, thus avoiding exponential underflow.

STOCHASTIC E Q U I L I B R I A FOR MULTIVARIATE REACTION SYSTEMS

609

The marginal probability function of N~ is found b y computing the sum

~1=o ~=o

~-1=o~

=o ~o-=o

-~1 Oj.(i~+l/2) l n ( i j + l / 2 ) ]

fori/~ = 0 , 1 , . . . , N ,

where the upper limits of summation are

ut = N - i s - i l - i 2 ut = N - i s - J 1 -

... -Q,

... -is-l-iS+l-

1 I.

%~1

~O~F3

N

Q

%AM1

-27.9 -30.4

428.2 417.1

10

1

-27.9 -30.4

1.5

-20.1

40.6

- 21.4

4 1.5

419.5 4 13.3

0.5

-41.7 - 42.8

-3.9 - 6.9

446.0 + 21.8

0

- 65.5 - 57.4

- 22.5 - 28.7

+ 88.3 + 22.3

1

- 27.9 - 30.4

0 0

+ 28.2 + 17.1

1.5

- 19.3 - 20.5

+ 0.5 - 1.4

+ 18.8 + 12.7

0.5

-50.0 - 54.4

-5.1 - 9.6

455.1 + 23.5

0

- 96.2 - 60.9

- 89.4 - 43.8

+ 185.5 + 2.5

100

~oAMu 0 0

%Adl//a 428.2 417.1

0.01

9

"9

9

(a) 0.1

9

)9.99

x"

X

9 xA

99.9

A

A x

99 A

x

A (g) ( [ [ )

9

(b)(i i) 9

(C) ( i i ) x

9 xA

il

&

90

Z~

9 x

A

ox~

5C - -

50

O•

9

• •

9C _' X~

g

9 x

i9

Z~

i

• 9s

99.E

3.1

99.9s

F i g u r e 1A.

'.I [ I 0,2 4

[ I l/,ll I I I 6s Io"o,z 4 6

t I ,,f[ I 1 I I 3.01 8 Io"ot2 4 6 s Jo

C u m u l a t i v e m a r g i n a l d i s t r i b u t i o n s o f N I ( . ) , N2( x ), N 3 ( A ) on normal probability paper

N = 10. (a) (ii): Qt = Q2 = Q8 = 1; F 1 = 0.240, F 2 = 0.333, F s = 0.427. (b) (ii): Qt = 1.25, Q~ = 1, Q8 = 0.75; F 1 = 0.240, F 2 = 0.333, F a = 0.427. (c) (ii): Q1 = 0.75, Q~ = 1, Q3 = 1.25; F 1 = 0.240, F 2 = 0.333, F 8 = 0.427. 0,01

o

(b) (f)

O.J--

(d) ( i i )

(ii)

99 9s

x

(e) ( i i )

x -- 9 9 9

A _

9

x

x

x

o k x

• x

x 9

--

x

9

&.

-- 90

x

x

• 9

x

A

x

--99

A

x

x

A

A

x

50--

•

9

~

L~ Zk

--50

o L

•

x

~ 0•

90--

A

--10 x

•

oL

99

99.9 --

99990 Figure lB.

i

oT'o' i

,o,o'"T

i

,to Do,

C u m u l a t i v e m a r g i n a l d i s t r i b u t i o n s of N I ( . ), N2( x ), N 3 ( A ) on normal probability paper

N = 10. (f) (ii): Q1 = Q2 = Q3 = 0; F 1 = 0.240, F ~ = 0.333, F 3 = 0.427. (d) (ii): Q1 = Q2 = Q8 = 1.5; F 1 = 0.240, F 2 = 0.333, F 8 = 0.427. (e) (ii): QI = Q2 = Q8 = 0.5; F 1 = 0.240, F 2 = 0.333, F 3 = 0.427.

o.o, (~)--a~@

A..'-'- )9.99

l---x

A

O. !

9

)9.9

9

X L&.

)9 =) (I.i)

(b)(li)

9

x

A

x

i

"~ 90

• x

50--

50 • X

90--

--10

A •

99--

LX x

99.9

--Lt

--

9e.99LIZ~l

] I I I ,t I Xl ~I

10205040506070

F i g u r e 1C.

Jr 0

I0

I lZ~l ] I I

3 0 4 0 5 0 6 0 7 0 8 0 90

Cumulative marginal distributions of Nt(.), n o r m a l p r o b a b i l i t y paper

o.ol I00

N 2( x ), N s ( A ) o n

N = 100. (a) (ii): Q1 = Q2 = Q3 = 1 ; F t = 0.240, F 2 = 0.333, F 8 = 0.427. (b) (ii): Q1 = 1.25, Q~ = 1, Qs = 0.75; F 1 = 0.240, F 2 = 0.333, F 3 = 0.427.

o/7

0.01

--

~

99.99

-

9 ~

99,9

)9

( f ) ClU)

)0

i~,2

zx

ix 5O

k

zX

•

50

Lx zx

9C

A A

9s

z~ /,,

99,9

99'990

F i g u r e 1D.

--

/',

f I l 5I0 ~ ~0 17 0 8 0I9 0I 100 I

}I0 2 0 5 0 4 0

)ol

C u m u l a t i v e m a r g i n a l d i s t r i b u t i o n s o f ATz(. ), N~( x ), N 3 ( A ) o n n o r m a l p r o b a b i l i t y paper.

N = 100. (f) (ii): Q1 -- Q2 - - Qs -- 0; F t = 0.240, F= = 0.333, F 3 = 0.427.

620

J O H N P. MULLOOLY

LITERATURE Bartlett, M. S., J. Brennan

and J. I~ Pollack. 1971. "Stochastic Analysis of Some

E x p e r i m e n t s on t h e Mating of Blowflies." Biometrics, 27, 725-730. Chiang, C. L. 1968. Introduction to Stochastic Processes in Biostatistics. N e w Y o r k : Wiley. McQuarrie, D. A. 1967. " S t o c h a s t i c A p p r o a c h to Chemical K i n e t i c s . " J. Appl. Prob., 4, 413-478. Mullooly, J . P., 1975. "Stochastic K i n e t i c Models F o r Small Reversible Systems: The E q u i l i b r i u m D i s t r i b u t i o n . " Bull. Math~ Biol., 37, 127-138.

RECErVED 1-9-76

APPENDIX

Verification of the closed-form solution ( 1 4 ) f o r three-component order-symmetriv systems. The verification t h a t (14) is the solution of (9) will be d e m o n s t r a t e d b y p r o v i n g t h a t an unnormalized version of (14) satisfies P - A s ( N ) = 0. Recalling the general inductive definition of states, we see t h a t for the c = 3 c o m p o n e n t system the kth-section consists of the following states: State index

Section ]z

N1

Lk Lk+ i :

]c--I k -- I : ....

Uk--1

]c--1

U~

k-1

N2

Ns

0 I :

N--}+l N-k :

N--k

1

.N-k+

1

(A1)

0

A form of (14) which is convenient for the d e v e l o p m e n t t h a t follows is

P(il, i2, is)

oc

(N,) 71., Q1'fl!," (N!) Q~.f212. (N~,.) Q3.fa,3.

(A2)

F r o m (A1) and (A2) it is seen t h a t the unnormalized probability for the i t h - s t a t e within t h e kth-section of states is given b y

pki=

[ N~]QI.fl~-I[

f o r k = 1,2 . . . . . N + I ; i

N~IQ2

t 1[

"f~-

N,

(N-k'-+~-i)!

]Q3.faN-k+~-'

---- 1, 2 . . . . . N - k + 2 .

Multiplication of the row vector P by the j~b-column within the ]r of As(N) gives

p(#-l)ffal(N - / c + 2 - j + --P~J"I l l s ( I t - 1 ) + f 1 3 ( k -

1) +P(~-1)(/+l)'f21(j)

of columns

+Pk(l-1)'fa2(N-- ]r~- 1 --j + 2)

1)+f21(j-- 1)+fa3(j-- 1)+f31(N--]c + 1 - - j + 1)

+f32(N -- k + 1 - j + 1 )] +Pkij+l) : r2a(j) + forj = 2,3 ..... N-k+l,k

(A3)

= 2 , 3 . . . . . N.

P(k+x)U-1).fl2(k) + p(~+l)/.fla(k)

= 0, (A4)

STOCHASTIC EQUILIBRIA FOR MULTIVARIATE REACTION SYSTEMS

621

The first column of the kth-section times the row vector P gives the equation P(k-1)I "f31( N -- ]e ~- 2) -}-P(k-1)2 "f21(1) --Pkl" [f12(/C-- 1) +f13(k -- 1) +fzx(N - k + 1) + f s u ( N - k + 1)] +P~2 "f23(1) +P(k+l)l "f13(k) = 0, while from the ( N - k + 2)nd-colurnn we obtain

(A5)

P(k-1)(N-k+2) "f31(1) +P(k-1)(N-k+8)'f21(N - k + 2) +Pk(N-k+l) "f32(1) -- Pk(N-k, 2)" [f12(/c -- 1) +f13(k -- 1 ) +f21(N - k + 1 ) +f28(N - k + 1)] +P(~+l)(~v-~+1)'fls(k) = 0.

(A6)

We now substitute the values of Pk~ given b y (A3) into (A4), (A5) a n d (A6), and p u t fts(l) = ftj .lO-~ for symmetric systems. E q u a t i o n (A5) becomes QI[N !]

(N -'-k--I- 1 ) !

9{f31 'Sl k-8 "fz N-k+2 -i-f21 "fl k-2 "f2 "f3 N-k+l

[ N, ]Q1 - (f12 -{-f13) "fl k-1 "f3 N-k+l ) -t- ~ - - ~ l

F N,

qQ3

[ N ! ] Q 2 | L ( N - ' " k) !J|

.{S2s.flk-l.f2.f8

N-~:

-{-f13 "fl k "f8 N - k -- (f31 "{-/32)"fl k - I "f3 N-k+l } = O. This equation is satisfied for a r b i t r a r y Qx, Q2, Q3 if and only if both of the curly-bracketed terms are zero, which is equivalent to requiring t h a t

A l . f ~ + fzl..f3 = (]12+As)'A, f13 "fl -t-]23 'f2 = (f31 -i-f32)"]3.

(A7)

B u t these relationships are consequences of the definitions of the composite rate constants fl,,f2, f3. Hence we have verified (A5). (A6) becomes (k-2)!]

(N-k+

1){

[N!]Q3"{f31"flk-2"f2N-k+I"f3

+f21 'fl k-2 "]2 N-'c+~ -- (f12 +f13) "]1k-1 ']2 N-~+I } 9

+A2 "A k ']2 x - ~ - (]21 +A3)"A k-i ']2 ~v-~§ } = 0. which is satisfied if a n d only if

]21 "]2+f81 .f8

= (A~+]13) "A,

f12 'fl +f32 "f3 = (fzl +f28) "f2-

(A8)

T h a t (A8) are identities is easily seen from the definitions off1 ,fu ,f3. Hence (A6) is verified. (A4) is seen to be [___~_r[ ] Q I [

(k= j

N!

]Q2

rLiN-N!2-j),J] Q3 "{f31"flk-2"f2i-1"]3N-Ir

"~'f21 "]1 k-2 "]2, "]3 N-k+2-j -- (f12 "q-]13)"fl k-1 "f2j-1 .]3N-k+2-j }

+

'(N - k + 2 - j ) !

"{f32 "]1 k-1 "f2t-2 "f3 N-k+3-j

]o,

+A~-"A ~ "]d -2 "]3 ~-~+2-j - (f21 +]23) "]1 x-x -]2 j-1 'f3 N-k+z-~}

I- ,lo,

+ t(~- l)!j

L ( J - 1)!J

(N-k-+ 1- j ) !

~

+f18 "]l k "]2s-1 "Is N-tc+l-j -- (f31 -I-]82)"fl k-1 "]2~-1 "]3 N-k+2-1 ) : O.

622

JOHN P. MULLOOLY

The three curly-bracketed terms are zero for all Q1, Q2, Q3 if a n d only if

f~ ./2+Ia~ "f3 = (/12+I~8)"/i, /12 "fl -t-/32 "f3 ---- (/21 -[-f23) "f2, f13 "fl +f2s "f2 = (f31 +fs~) "f3. (A9) Since these identities hold, we have verified (A4). This completes the proof for equations corresponding to the kth-section of columns of As(N), for k -- 2 . . . . . _hr. The verification of the equations derived from the first a n d last ( ( N + 1)st) sections is achieved in a similar m a n n e r .