Mic: _. . .
DEPARTQ‘HQNT (3F C‘,H&-_‘:\»,¢:_) z .4.
um mmsms. wow-em

ABSTRACT

ON THE APPLICATION OF STOCHASTIC PROCESSES
TO REACTION KINETICS

by Casimir J. Jachimowski

The theory of stochastic processes has been applied to the
problem of determining the inherent random fluctuations in reaction
systems. Such fluctuations, similar to the. fluctuations at equilibrium
treated by statistical thermodynamics, are inherent to some extent
in all chemical processes. The stochastic approach provides a well-
defined mathematical basis for analyzing the extent of inherent
fluctuations by assigning a time dependent probability function to the
composition variables.

The formulation of stochastic. reaction models was approached
in two ways. The first approach =assu1ned that the course of a chemi-
cal reaction could be treated as a Markov process. The basic notion
of a Markov process is the assignment of probabilities to the possible
transitions between different states of the system; the transitions
between different states being due to the occurrence of the various
elementary events that comprise the reaction scheme. From the
transition probabilities stochastic differential equations are obtained
which define the stochastic model. These stochastic equations parallel
the system of ordinary differential equations of the deterministic
approach.

Exact solutions for Markov models of the bimolecular reactions

2 A -—->- B and A + B -->- C were obtained and compared with their

Casimir J. Jachimowski

respective deterministic models to determine the extent to which the
deterministic models are applicable to small systems. The appli-
cation of the Markov approach to complex reaction systems was also
considered using approximate techniques which, in most cases,
easily generated the first and second moments of the concentration
variables. The reaction systems treated approximately included
radical chain polymerization, condensation polymerization, random
polymer degradation, chain reactions with and without branching, and
various enzyme-substrate reactions.

The second approach used to construct stochastic reaction
models was based on the premise that the possible states or compo-
sitions of a reaction system form a so-called "composition space. "
From the concepts of set theory and a frequency interpretation of
probability a probability function was formulated. The set theoretic
approach was applied to the first order processes A—->- B and
A: B and the second order process 2 A—->- B. The results were
compared with those obtained from the Markov approach.

Further extensions of stochastic methods to other types of
reaction processes and physico-chemical systems are proposed and

discus sed.

ON THE APPLICATION OF STOCHASTIC PROCESSES
TO REACTION KINETICS

BY

Casimir Joseph Jachimow ski

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Chemistry

1964

To Mary

ii

AC KNOW LED GMEN TS

The author expresses his appreciation to
Drs. M. E. Russell and D. A. McQuarrie for their
guidance and encouragement during the course of
this investigation.

Appreciation is extended to the National Insti-
tutes of Health for financial support in the form of a
research as sistantship.

>i< >I< >§<>i< >3 31‘ 31‘ *>€<* ** **

iii

TABLE OF CONT ENTS

Page
I. INTRODUCTION 0 O O O O O O I O O O O O O O O O O O . O O O 1

II. THEORY OF MARKOV PROCESSES. . . . . . . . . . . . . 6

A. IntrOduCtion O O O O O O O O O O O O O O O O O O O O O O 6
B. General Formalism for Reaction Systems. . . . . . 8

III. APPLICATION OF THE MARKOV TREATMENT. . . . . . 11

A.ExactTreatment.... 11
1.Reaction:ZA——->-B.............. ll
2.Reaction:A+B—-9-C............. 21

B. Approximate Stochastic Approach . . . . . . . . . . 26
1.Introduction................... 26
2. Application of Method One . . . . . . . . . . 30

a. Reaction: Radical Chain Polymerization 30
b. Reaction: Condensation Polymerization. 38
c. Reaction: Random Polymer Degradation 42
d. Chain Reaction in General. . . . . . . . 46
3. Application of Method Two . . . . . . . . . . . 51
a. Reaction: Autocatalytic A-—>— B. . . . . 51
b.Reaction:A—>—B—>C........ 52
c. Reaction: 2A-->- B—->- C. . . . . . . . 53
d.Reaction:ZA\:-_:B.......... 53
e.Reaction:A:B............ 54
f. Enzyme-Substrate Reactions. . . . . . . 55

IV. SET-THEORETIC APPROACH TO REACTION KINETICS . 63

A.Introduction......................63
B.Applications......................67
l. Irreversible First Order Process . . . . . . . 67
Z. Reversible First Order Process . . . . . . . . 69
3. Consecutive First Order Process . . . . . . . 70
4. ASecond Order Process . . . . . . . . . . . . 74

iv

TABLE OF CONTENTS - Continued

Page
V. DISCUSSIONANDSUMMARY............. 75
LIST OF REFERENCES. . . . . . . . ......... . . . . 79

APPENDICES 0000000 O O O O O O O O O O O O O O O O O O O 82

TABLE

II.

III.

IV.

LIST OF TABLES

Page
Fraction of Reactant for ZA-—>- B . . . . . . . . . 19
Fraction of Reactant for A+ B —->- C, A0 = Bo . . 25
Fraction of Reactant and Coefficient of Variation
for the Exact and Approximate Methods, 2A-—>- B 31

Fraction of Reactant and Coefficient of Variation
for the Exact and Approximate Methods, A + B —>— C,

A0 = BO 0 O O O O O O O O O O O O O O O O O O 0 O O O 32

vi

FIGURE

1.

LIST OF FIGURES

Page

Fraction of reactant when x0 = 10 for the reaction
2A-->-Basafunctionoftime............ 17

Fraction of reactant when x0 = 50 for the reaction
ZA-é-Basafunctionoftime............ 18

Coefficient of variation for the reaction 2A——>— B
as a function Of time. O O O O O O O O O O O O C O O O 20

Coefficient of variation for the reaction A + B —->- C,

A=B asafunctionoftime............ 24

O 0’

vii

I. INTRODUCTION

Chemical kinetics can be defined as the study of those systems
whose properties, or more specifically whose composition, is time
dependent. The most obvious goal of a kinetic study of a particular
system is the creation of a model that agrees to a reasonable degree
with the real system. When formulating a mathematical model, we
can select one of two approaches for the study of the reaction system
concerned. These two approaches, which are termed deterministic
and stochastic (or probabilistic), reflect the nature of the postulated
model which is usually expressed in mathematical form.

The mathematical theory of classical chemical kinetics is con-
cerned with the prediction of the concentration of the reacting Species
as a function of time. The concentration is treated as alcontinuous,
real-valued function of time and the mechanism of the reaction is
expressed by a set of differential equations which are obtained by
invoking the law of mass action. The integral solutions of these
equations, when they can be obtained, lead to ordinary analytic
functions of time. ‘On the other hand, the basis of the stochastic ap-
proach is to treat the concentration variable as a time varying,
discrete, random variable, the values of which at different times
being influenced by probability effects. To illustrate these two ap-
proaches consider the irreversible first order process, A-—>- B .

Let the function x (t), real valued and continuous, denote the
number of A molecules in the system at time t. According to the
classical approach if it assumed that at time t there are x molecules of
A in the system then the change in A during the time interval (t, t + A t)

is proportional to the number of A molecules present at time t.

Therefore we have the relation

Ax(t)=-kx(t)At

which leads to the differential equation

dx (t) _
dt - -kx(t) (1.1)
If we assume that x (0) = x 0' then the solution of (1. 1) is
-kt
x (t) = x e (1.2)

0

Such a mathematical approach results in a "deterministic" model

since once the initial value x is known expression (1.2) determines

o
the value for x (t) at subsequent times. This also implies that

whenever the initial value x is the same, the value of x (t) will

0
always be the same for a given time t.

Now consider the stochastic model of the reaction A ——>— B.
The stochastic model of this reaction was first studied by Bartholomay
(1). Let the integral valued random variable X (t) represent the
number of A molecules at time t, and assume X(O) = x0. In the
stochastic approach instead of deriving a functional equation for X(t)
an attempt is made to find an expression for the probability that at
time t the number of A molecules is x'. That is, we want to find
P (x , t) = Prob { X(t) = x3 . Bartholomay showed that the prob-
ability function is

-k -kt x -x
Xo e x 1:(1 '8 ) o (1-3)

X

P(x.t)=

An examination of the two models shows that in the deterministic
case (1.2), for k and x0 fixed, for every value of t there is associated

a real number x (t). While in the stochastic case (1. 3), for k and x0

fixed and for every pair ( x, t), there exists a number P(x, t),
0: P(x, t) _<_ l, which is the probability that the random variable
will assume the value x at time t. In view of this difference it can
be seen that the deterministic approach (1.2) describes the average
number of A molecules while the stochastic approach (1. 3) takes into
consideration random fluctuations. Since the stochastic approach
takes into consideration random fluctuations, it can be a valuable
tool for studying the kinetics of chemical reactions. Moreover, for
many reaction systems the deterministic average is not the same as
the stochastic mean, and therefore the stochastic approach provides
an excellent method for examining the extent to which the determin-
istic approach is applicable to small systems, i.e. , systems contain—
ing 100 or less molecules.

Fluctuations in the number of reactant species are inherent to
some extent in all chemical reactions, but because of the operation
of the law of large numbers such fluctuations are microscopic and
not noticeable. However, there are some chemical processes in
which one may expect macroscopic fluctuations to occur. Such fluc-
tuations in the number of reactants can be the cause fof irreproducible
reactions (2). For example, a reaction will be inherently irreproducible
if a measurable change of concentration of reactant depends on the
intervention of a process which is so slow that the number of ele-
mentary events in a macroscopic time interval is subject to large
relative fluctuations. , This is so in radioactive decomposition where
the events of a very slow process can be observed. Also a chemical
reaction may be such that the formation of one or a few reactive mole-
cules is followed by a diverging chain reaction, which leads to the
rapid (possibly explosive) completion of the reaction: if the formation
of the reactive molecule which initiates the chain is a ”rare event, "
i.e. , one preceded by a long, variable period of time, the course of

the reaction will be irreproducible (2).

The application of the stochastic method to reaction kinetics
has been considered by several workers. The first application
appears to be that of Delbruck (3). Delbruck studied a stochastic
model of the autocatalytic reaction A—->- B. He found that large
statistical fluctuations in the amount of reaction taking place in a
given time can be encountered during the initial stages of the reaction
and approach a constant value during the latter stages. Delbruck
showed that during the initial stages the fluctuations in the number
of reactant species was of the order of the square root of the mean
number of reactant molecules. In 1953 Singer (2) discussed the
application of the stochastic method to the study of irreproducible
chemical reactions. The reactions studied by Singer were a chain
reaction without branching and a chain reaction with branching.
Singer reported that a number of chemical reactions, for instance,
the oxidation of formic acid by potassium nitrate (4), the slow or
explosive decomposition of some solids, and the initial stages of
certain polymerization reactions (5) have been reported to show
fluctuations or “irreproducibilities. ” Singer concluded that in small
systems fluctuations in the number of reactant species could be
responsible for irreproducible reactions, together of course with
those due to the presence of impurities. Some time later Renyi (6)
treated the bimolecular reaction A + B —->- C and showed that the law
of mass action, which is a basic axiom in the deterministic theory,
is only approximately valid and does not hold for small systems
(see Appendix C). Recent applications have been done by Bartholomay
who discussed the irreversible unimolecular reaction (1), A—>— B
and also a stochastic model for the Michaelis-Menten formulation
of an enzyme—substrate reaction (7-9). From the unimolecular
stochastic model Bartholomay showed that the stochastic mean and the

deterministic average concentration were identical; however,

fluctuations were still present. From the stochastic Miclaelis-Menten
model he shOwed that the deterministic rate equations may be obtained
from the stochastic rate equations by a process of "averaging, " but
failed to calculate the extent of random fluctuations about the mean
number of enzyme and substrate molecules. Ishida (10) extended
Bartholomay's' method of treating unimolecular reactions to a more
general stochastic process having a time dependent rate constant.
Using this generalized stochastic approach she studied the bimolecular
process 2A—->- B (see Appendix B) and a multidimensional stochastic
model of the unimolecular decomposition theory of Kassel. The most
recent study was 'that of McQuarrie (11) in which he discussed the
application of stochastic methods to the study of the kinetics of small
systems. McQuarrie calculated the extent of fluctuations for the
reversible first order process A: B and for the first order
parallel reaction system A <: .

The previous investigations of stochastic models of chemical
reactions were based entirely on the assumption that the course of a
chemical reaction can be treated as a Markov process. This approach
has been successfully applied to various first order processes and to
some complex systems. However, the application of Markov theory to
second order processes and to complex reactions has received little
attention because of the inability to obtain explicit solutions for the
proposed stochastic reaction models. The present investigation was
undertaken first, to determine the feasibility of finding workable solu-
tions for stochastic Markov models of certain second order processes
using the probability generating function technique; second, to develop
well-defined mathematical techniques for easily obtaining approximate
solutions of stochastic models of complex reaction systems; and third,
to consider the feasibility of applying a stochastic set-theoretic

approach to the study of reaction kinetics.

II. THEORY OF MARKOV PROCESSES

A. Introduction

 

Our consideration of stochastic processes will be restricted to
what are called Markov processes, that is to say, to processes
where all future probability relations are fully determined by the
present state, analogous to classical mechanics where the present
state uniquely determines the future development of a system.
Furthermore, we shall focus our attention on the so—called discon-
tinuous type of Markov processes where the changes occur in jumps,
that is, the system remains for a time unchanged, then undergoes a
sudden change into another state.

Suppose that a physical system possesses a finite or denumer—
ably infinite number of states E1, E2, . . . , Ej, . . . , and let

X(t) = j if the system is in the state Ej at the instant t. Now we define

the transition probabilities

Prob {X(t) =k| X(s)= i} = P(i—>k: s,t) where s <t

to be the probability that at time t the system is in state Ek if at time
s it had been in the state Ei' Consider a time u in the interval (s,t)
that is, s _<_ u _<_ t and let us now consider the transition from the state
Ei at time s to the state Ek at time t through all times 11 in the inter-
val (s, t). The passage from E1 at time s to Ek at time t must occur
by way of some state Ej at time u, and for a Markov process the
probability that the passage goes through aparticular state Ej is

P(i ->j < s, 11) P (j —>- k : u, t). It follows that from the total prob-

ability theorem (12) we must have

P(i->k:s,t)=?P(i—->j:s,u)P(j->-k: u,t) (2.1)

6

This is the fundamental equation which may be taken as the descrip-
tion of a simple discontinuous Markov process; it is generally known
as the Chapman-KolmOgorov equation (13, 14).

The basic notion in the application of a Markov process to
chemical reactions is that of a system with finite or denumerably
infinite states, a continuous time parameter, and transition-prob-
abilities. Since it has been acknowledged, in general, that under a
given temperature and external parameters the rate of a chemical
reaction is determined only by the composition of the reaction system
at the time and does not depend on its previous history, the course
of a chemical reaction is a Markov process. Thus a reaction system
can be described by the Chapman-Kolmogorov equation for discon-
tinuous Markov processes. Furthermore we may also assume that
the transition probabilities P(i ->- k : s, t) are stationary in time since
the probability of a transition from one state to another state of the
system in a reaction system depends on the length of the timé.’ interval
alone, and not on the particular instant at which the present state was
reached. . With this further restriction the Chapman Kolmogorov
equation may be written as

P(i—>- k: t-s) = z P(i—>j : u-s) P(j—>- k: t-u) (2.2)
J

Analytically, a Markov process is completely determined by its
transition probabilities. . In addition, if the random variable X(t) has a
given initial value we can define the absolute probability of an event,

say E and therefore the probability function

j:
P(j,t) = Prob g X(t) = 3'}
If the initial distribution P(i, s) is known we can write

P(k, t) = E P(i, s) P (i-—>- k: t—s) (2. 3)
1

or from equation (2. 2) we get

P(k,t)=2 (E P(i,s) P(i—>j : u-s) ) P(j—> k: t- u)
j 1

which can be reduced to

P(k, t) = z P(j,u) P(j->k : t-u) (2.4)
3

By letting t = t + At and u = t equation (2.4) becomes

P(k,t + At) = 2 P(j, t) P(j ->— k: At) (2.5)
J

Equation (2. 5) is the Chapman-Kolmogorov equation for changes in

a system in the time interval (t, t + A't).

B. General Formalism for Reaction Systems

As mentioned previously the basis of the stochastic approach
is to treat the composition (defined as the number of molecules per
constant volume of the reaction system) as a time-varying, discrete,
random variable the values at different times being influenced by
probability effects. If the random variables N1(t), Nz(t), . . . , Nr(t)
represent the number of molecules of species 1, 2, . . . , r
respectively, in a multidimensional reaction process, the probability
that the system is in a state N(t) = n (where N(t) denotes the variables
N;(t), Nz(t), . . . , Nr(t)) is given by the function P(n,t) .

According to the Chapman- Kolmogorov equation (2. 5) the prob-

ability that the system is in the state n at time (t + At) is
P(n,t+At)= E P(n',t)P(n'—> n;At) (2.6)
nl

where P(nfit) is the probability that the system is in the state n' at t,
and P (n' '9- n ; At) is the probability of the transition from the state
n' to the state n in the interval (t, t + At). The transitions between

different states in a reaction system are due to the occurrence of the

various elementary events which form part of the reaction scheme.

Equation (2. 6) can also be written in the form

P(n,t+At)= Z P(n',t) P(n'—>' n; At)

n'afn
(2.7)
+ P (n.t) P(n->- 11; At)

where the summation is over all states n' except where n' = n. The
transition probability P(n —> n : At) is the probability that no change

occurs in the time interval (t, t + At), and therefore
P(n—>n;At)=l- I; P(n—9 n':At) (2.8)
n 11

since 23, P(n-—>- n' :At) = 1. In a chemical reaction the probability
of the occurrence of an elementary event is considered to be of the

form
P(n' ->- n ; At) = F(n',n) At + 0 (At) (2.9)
where 0(At) is probability of the multiple occurrence of the event

n" ->' n and is defined so that 0(At)/At ->- O as At ->- 0. Substituting

equations (2. 8) and (2.9) into (2. 7) gives

P(n,t + At) = n, 5n P(n'.t) F(n',n)At (2.10)

+ P(n,t) i 1 - ,5 F(n,.n')At} + 0(At)
n n
Transposing the term P(n, t) and dividing by At gives

P(n, t + At) - P(n, t)
At

 

'2 &P(n', t) F(n',n)
n In (2.11)
0(At)

At

 

- P(n, t)F(n, n'fg +-

Taking the limit we get

10

Eggs—Pl— .—. n'gn {LP(n'.t) F(n'.n) - P(n,t) F(n.n')} (2-13)

which is the fundamental equation which may be taken as a description
of the stochastic reaction system; it is often called the stochastic

differential difference equation. Equation (2. 12) defines the rate of

change of the probability function with time.

III. APPLICATION OF THE MARKOV TREATMENT

A. Exact Tr eatment

 

1. Reaction: 2A->r- B

 

According to the deterministic theory of reaction kinetics the

rate of this bimolecular reaction is given by the differential equation

_ d [A.]_= k [A]2 (3.1)

(it

From (3. l) the number of reactant molecules in the system at time
t is given by

,. . [Me
[A]- 1+[A]okt (3'2)

where [A]o is the initial concentration of reactant. It should be noted,
however, that for small systems (with the concentration expressed

as the number of molecules per constant volume of the reaction mix-
ture) the rate, according to the law of mass action, is proportional to
the concentrations of the reacting species, in other words, to the
number of ways a pair of reactant molecules can be chosen from a
total of A molecules. Hence for the reaction 2A—>— B the rate is
Mam-Ll

proportional to or

.. EléJ.=k[A] [11-1] (3.3)

(it

From (3. 3) the number of reactant molecules at time t is given by

[A],

[A] = [A]o + F- [A]; exp (- kt) (3'4)

11

12

However, the modified deterministic model expressed by (3. 3) still
predicts a precise value for the number of reactant molecules at
time t and neglects any fluctuations about this value. In order to
consider fluctuations, a stochastic model must be used.

The basis of the stochastic approach, as outlined previously,
is to consider the reaction 2A -—->- B as being a pure death process
with a continuous time parameter and transition probabilities for
the elementary events that make up the reaction process. Letting
the random variable X(t) be the number of A molecules in the system
at time t the stochastic model is then completely defined by the follow-
ing assumptions:

(1) The probability of the transition x + 2 -> _x in the

interval (t,t + At) is %— k (x + 2) (x + 1) At + 0(At) where
k is a constant and 0(At) is defined so that 0(At)/At ->- 0

as At->- 0.

(2) The probability of the transition x + j ->- x , j > 2 in
the interval (t, t + At) is 0(At).

(3) The probability of the transition x - j ->- x , j > O,

in the interval (t, t '1' At) is. zero.

(4) The probability of the transition 1): ->- x , in the interval

(t, t + At) is [l - %- kx (x- 1)At] + 0(At),
In view of these assumptions the following relation is obtained.

_ 1
P}, (t+At)-Px +2(t)[2—1<(it+ 2)(x+l)At] (3 5)
+ Px (t) [ 1 - 7;— kx( 9: -1)At]+ 0(At)
where Px (t) =Prob { X(t)= x , X =O,2,4, . . . , 1:0} . By
transposing the term Px'(t)' dividing by At and passing to the limit

one easily obtains the differential difference equation

13

(1 PK

--——'- = %k(x+2)(x+1)Px

dt

1
+2'2'k X(X llpx

(3.6)

Differential difference equations of this type are usually dealt with
by transformation into a partial differential equation by making use
of the probability generating function
00
F(s.t)= Px (t) s", Is )5 l (3.7)
x=0 ‘

where s is an arbitrary variable with no physical significance, intro-
duced purely for mathematical convenience. By means of (3.7)

equation (3. 6) may be transformed into the partial differential

equation

2
af 52—; <3...

‘5‘ ‘ i k ‘1 ’ 52’ a
Equation (3. 8) may be solved by separation of variables to give

00
F(s,t) = 2 An Sn(s) Tn(t) (3.9)
n=0

where 1

2-
Sn(S) = Cn (S)
Tn(t) = exp i - )— kn(n-1)t}

Sn(s) is a solution of the Gegenbauer equation (15-16)

2V V
(1-52)}3—‘5-27351” (2v+l)s 91.33.811.131.
+n(2v+n)C;(s)=O (3.10)

whenv= ~2-.

14

The solutions of (3. 10) are Gegenbauer polynomials defined for

. . . n. .
integral values of n to be the coeff1c1ent of h 1n the expansmn of

(l - 2hs + hz) - v in powers of h, that is
- v 00 v n
(l-Zhs+hz) = z Cn(s) h (3.11)
n = o
For the case when v = - i‘: the first few polynomials are

1

c;7(s) = 1
_1
Clz—(s) = -s
l
c;7(s) = i- (1- 82) (3.12)

The coefficients, An, can be determined most easily from the

boundary condition

a]? x '1
s at t=0 (3.13)

= X.
65 °

Applying (3.13) to (3.9)

00
2

0

 

1
__ '2'
x 5‘0 1 = A dcn (S) (3.14)

0 n= 1 n ds

From Whittaker and Watson (17) we find that

 

1
dggfis) __ i— -
ds - - Cn-l (s) - Pn-l (s) (3.15)

where Pn-1(s) is a Legendre Polynomial (18). Equation (3. 14) can

be written as

An Pn—l (S) (3.16)
1

15

Multiplying both sides of (3.16) by Pm_1(s) and integrating between
-1 and +1 yields

+1 .xO-l 00 +1
{ xos Pm_1ds=- 21A,, [1 Pn_1Pm_1ds
' n" ' (3.17)

Using the orthogonality properties of Legendre polynomials (l9)
equation (3. 17) can be written as

+1
-1

__ l-2n x0
An - 2 f1 xos Pn_1ds (3.18)

 

From Irving and Mullineux (20) we have the relation

 

+ 1 ' m-R+l
I = f sm PR(s) ds = I: r( 2 :R-i-3
" 1 2 (m-R)‘.r('Ln—'Z—-")

with the conditions

(a) m/and R_>_ O

(b) m 2 R with (m - R) even or zero
(c) I: 0 for (m - R) odd

(d) I: o for m < R

In order that An 7! 0 equation (3. 18) must satisfy conditions (a) and

(b). These conditions are satisfied if n = 2, 4, . . . , x0 where

X0 is even. Equation (3.18) becomes

l-Zn r( 9‘0 +1) r( .xn "zn'l' l)

 

 

 

 

A = (—-n---) , r (3-19)
11 2 - x +n+1
fix. -n+1)[‘(~—1’a )
n= 2, 4, . . . , x
o
By using the easily proved relations
2.
<x>= <81") .<X(x-1)>=(Q—§-)
as $51 as 8:1

(3. 20)

16

where < x > is the mean or first moment, and < x (x - 1) > is defined

as the second factorial moment, it can be easily shown that

.xo
<x > = z; A'n T(t), A'n = -An (3.21)
n=2
30 _
<x(x -1) > = z. _9_(_%_1)_ A1,, Tn(t) (3.22)

The measure of the inherent fluctuations is given by the variance.

In terms of the probability generating function the variance is

(931? OF (31“
Var(;x)= (--2) +(-————) - (_)Z
(9 5 3:1 (38 3:1 (33 8:1
=<x2>-<x>z (3.23)

The coefficient of variation, CV (x, t), can be used to measure
the relative extent of fluctuations. CV (x,t ) yields the fluctuations

relative to the mean and is defined as

 

1
cv (x, t) = V2332) T (3.24)

An alternate approach to the bimolecular reaction 2A ——>-1B
given by Ishida (10) is outlined in the Appendix.

The fraction of reactant molecules have been calculated at
various times for (x0 = 10, 50, and 100 from the deterministic,
the modified deterministic, and stochastic expressions. The results
are shown in Figures 1 and 2 and Table I. The coefficient of variation
is plotted in Figure 3. The results indicate that the stochastic and
deterministic means approach each other quite rapidly, though fluc-
tuations about the mean still exist. Figure 3 shows that the relative

fluctuations decrease as xo becomes large.

l7

 

 

 

 

 

 

 

 

O. 9 _ Stochastic
_______ Modified
\
——-— — ---—- Deterministic
\
\

0.8 -— \

O. 7 h
< x >

X

o
0.6 _-
\
0. 5 .... \
\
\
\
\
\
N
x
‘\
1 l , X
0 0 0. 50 1. 00 1, 50
1x0 kt

Figure 1. Fraction of reactant when x0 = 10 for the reaction 2A->- B
as a function of time.

18

 

 

 

 

 

 

 

0.9 _ f 4 Stochastic-
)
——_...---- Deterministic
0.8 '—
0.7 “'-
< x ‘ \
2
x — \
.) \
\,\
0.6 "-
o 5 L-
0” * I
(l I? 0.50

 

x0 kt

Figure 2. Fraction of reactant when. x0 = 50 for the reaction 2A->- B
as a function of time.

19

 

 

 

 

 

ooe.o moa.o woe.o ee¢.o ~o¢.o eoe.o 2~¢.o om.s
oom.o mom.o mom.o amm.o mom.o eom.o Hmm.o oo.s
Nem.o sem.o eem.o eoe.o vem.o eem.o mam.o me.o
see.o oas.o Hee.o moe.o oee.o oes.o eme.o om.o
oom.o sem.o mom.o mam.o Hom.e mow.o mzm.o m~.o
2: u om om u 0.x S n o... 2: .1. o... om u 0.x S n ox E o .x

owumwcfibuouofl Uofifipoz owpmmnogm

L

 

m .A.|.<N you «cmuomom mo cofiomuh .H 3an

20

 

 

 

0.5“-
0.4....
xo=10
S:
o
‘3
"-1 0.3"—
H
«3
:>
1....4
o
E
93
'o
3":
o
o
0 0.2 t—
x0=50
x0 - 100
0.1 -
o I I I - l l l
.0 0.50 1.00 1.50
x0 kt

Figure 3. Coefficient of variation for the reaction 2A—->- B as a

function of time.

I
‘\

 

21

2.. Reaction: A + B——>- C

 

Let X(t), the discrete, time varying, random variable, be the
number of A molecules at time t and Y(t) = Z0 + X(t), the number of
B molecules where Zo = Y(O) - X(O) = Bo - x0 and Bo _>_ x0 . The
stochastic model is based on the following assumptions:

(1) The probability of the transition x+ 1 -->- [X in the interval

(t, t + At) is k( x +1)(Zo + x + 1) At + 0(At) where k is a
constant and 0(At) refers to the probability of more than

one such transition.

(2) The probability of the transition x - j ->- ,X , j > 0 in

the interval (t, t + At) is zero.

(3) The probability of the transition x —->— x is the interval
(t, t+ At) is 1 - {1...(20 + x) At +0(At).

In view of these assumptions we obtain the differential difference

equation

dP

x' =k +1
dt (x )(zo+x+1)1:>x

 

+1— kx (Zo -x)Px (3.25)

where Px(t) =Prob {X(t)=x,x =0, l, 2, . . . , x0}

As before, we put
(I)
F(s,t) = 2 PK (t) s"
X = O

and from (3. 25) we obtain the partial differential equation

(91“ ’-
5';— = kS(I-S) 3% + k(Zo +1) (1 " S)? (3026)

Application of the method of separation of variables yields

22
(X)

F(s,t) = 2: An Sn(s) Tn(t) (3.27)
n=o

where Sn(s) = Jn(Zo, Z0 +. l, s) are Jacobi polynomials and solu-

tions to the differential equation (21)

 

 

2
S(1-S) d Jnégigi S) + [q _ (P + 1) S] Citing); q, s)
+ n(n+p)Jn(p, q, s) = 0 (3. 28)

and Tn(t) = exp { -n(n+Zo)kt3

Using the boundary condition given by (3. 17) together with the

relation (22)

 

d Jn(P:q. 5) n(n + E)
= - + +1

d5 q Jn-1 (P 2: q 9 S)
and the orthogonality relation (23)

l
I s q'lu-mp'q Jn(p. q, s) Jn(p,q, s)ds

 

 

 

0
we get
__(_'1)n (Zn‘tZQ) I-(n'l'Zo) eri-l) rxn'l'Zo +1)
P(n-H) I"(zo o+1) 1"(x0 -n+1) 1"(xo+zo C,+n+1)
wheren=l, 2, 3, . . . , x

o O

The first moment and the second factorial moment are

'< x > = g0 (2n+ZO) (Xb'l'l) nxO+ZQ+1)T T41“)
n=l (sec—n+1) I"'(xo + Z0 + 11 +1)

0 (n-1)(n+ZQ+1)(2n+Zo) r(§Q+l) fix o+zo+1)
T(xo -n+l) I‘D:o + 2‘0 +n+ l)

(3.29)

<x (x -1)> = 2
n=2

Tnm

(3. 30)

23

An alternate, but not as convenient, solution given by Renyi (6)
is outlined in the Appendix.

For the bimolecular reaction A + B -—->- C if it is assumed that
Z0 is very large, that is, Bo > > A0 = x0, it can be easily shown
that the bimolecular process changes to a pseudo first order process.
When Z0 is sufficiently large only the first term in the expansion of

(3. 29) is needed. The first term (n = 1) of (3.29) is

Z +1
x0 (“32%) exp { -Zokt} (3.31)

For large Zo this reduces to

xO exp {-Zokt} = X0 exp {L -k't}

01'

< x > = x0 exp { -k't} (3. 32)
In a similar manner equation (3. 30) becomes

< x (x .1) > = x0 (x0 -1) exp{ -2k't} (3.33)

The variance is

-k't -k't
e

<x (x-1)> +<x’>-<xfz> = x0 (l-e ) (3.34)

which is identical with the variance found by Bartholomay (1) and
McQuarrie (11) for the first order reaction A-—>- B.
The fraction of reactant molecules and the coefficient of variation

for B0 = A0 = x0 = 10, 50, and 100 have been calculated. The results

are given in Figure 4 and Table II.

24

 

 

 

 

 

 

0.5)-
0 4r—
:3
.9.
E3
'2 0 3L- x0- 10
a:
>
9.4
o
E
Q)
'8
I“:
8
o
0 0.20—
___x0 =50
0.1r— __ xo=10(i
0 l I L l l l
0 0.50 1.00 1.50

xokt

Figure 4. Coefficient of variation for the reaction A + B -->- C,
A0 = Bo, as a function of time.

25

Table II. Fraction of Reactant for A + B ->- C, A0: Bo

 

 

 

m
--—--—< x > Stochastic < x > Determin-
x kt "0 T istic
0 x0: 10 xo=50 x0: 100 o

0.25 0.798 0.799 0.800 0.800

0.50 0.662 0.664 0.666 0.667

0.75 0.565 0.567 0.571 0.572

1.00 0.491 0.498 0.499 0.500

1.50 0.389 0.398 0.399 0.400

 

26

An examination of the results obtained from the stochastic and
deterministic models for bimolecular systems shows that there is
no correspondence between the deterministic average and the stochastic
mean in reaction systems containing a small number of interacting
species. Furthermore, fluctuations in such systems are large com—
pared to the number of species present. Such a lack of correspond-
ence between the deterministic equations and stochastic equations
will result in all bimolecular processes. This lack of correspondence
is the result of the assumption made in the deterministic model; that the
higher moments can be expressed as products of the-lowermoments,
e.g., <az>=< a>2 and <a b> = <a> <b> . Such anassumption is

not valid for smallsy-stems as shown by the results in Tables I and II.

B. Approximate Stochastic Approach

 

1. Introduction

 

The stochastic approach is applicable to all types of reaction
systems. However, in many cases the resulting partial differential
equation which defines the stochastic model cannot be solved exactly.
For this reason some approximate methods must be used to determine
the first moment and the variance from the stochastic equations.
Since the quantities of main interest in the stochastic approach are
only the first moment, < x >, and the second moment, < x?" >
much effort can be saved by applying methods which produce these
lower moments without having to solve for the probability generating
function.

A set of equations involving the time derivatives of the moments
of the concentration can be generated from the differential difference
equation for the process under consideration. Since the moments of

the concentration variable are defined by

27
00
< x > = 2 x P (t) (3.35)

a differential equation involving the moments can be obtained by
multiplying the differential difference equation by x n and summing
over the variable x. This technique is equivalent to the use of the
cumulant generating function which generates the curnulants of the
process (24, 25).

Consider the bimolecular reaction 2A—>- B. Multiplication of

2

equation (3.6) by x and x respectively, followed by summation

over x gives

d<x>

dt =-k<x7‘> +k<x> (3-36)

d< x2)

dt = -2k<x3> +4k<x.Z>-2k<x> (3.37)

In all second order processes the set of equations obtained in this
manner cannot be solved because the right side always contains

higher moments than the left side. Therefore some relationship between
the higher and lower moments must be assumed. Two approximate
methods of determining the first and second moments have been formu-

lated and are given below.

1. An assumption frequently made is to express the higher

11 ..

moments as a product of the lower moments, e. g. , < [X > -
n-' .

< X J><_x']>(2,26). Assumingthat<px3>=<xz><x>and

substituting this into (3. 37) and putting < x 2 > = < x > 2 in equation

(3. 36) we obtain the approximate set of equations

<
SL313 =-k<x>2+ k<x> (3.38)
dt
d<xz> 2 z
-———-—————=-2k<x> <x> +4k<x >-2k<x> (3.39)

dt

28

The solution of (3. 38) is

° (3.40)

 

Note that (3.40) and (3.4) are identical. The assumption < x 2 > =
< x > z is equivalent to reducing the stochastic model to the determin-
istic model, i. e. , no fluctuations are permitted. In a similar manner

from (3. 39) and (3.40) we get the second moment
<xz>=<x>2 z(—£‘?—'—-1—)(eZkt.e.kt)+1 (3.41)
, '3 x0

The variance of X(t) as a function of time can therefore be expressed

as

Var( x) = -§- (_Z‘_Q£'_1__) (e -e ) (3.42)

2. The second approximate method formulated for determining
the first and second moments is somewhat empirical. An examination
of Figures 3 and 4 suggests that the coefficient of variation increases

exponentially with time. Hence, if it is assumed that

[CV(x ,t)]" = exp { pt} -1 (3.43)

where p is a constant, then

 

<xz> = <x>7' exp {pt} - (3.44)
Since
‘ <x"‘>-<x>Z
[CV( X ,t)]z = < x > ZT

Substituting (3.44) into (3. 36) and solving for the first moment

gives

xo(p+k)
Xokexp{pt} -(k xO-p-k)exp{-kt}

 

<x> =

(3.45)

29

The constant p can be determined from the equation involving the time
. . n
derivative of the second moment if 1t 15 assumed that < x > = x0

at time t = 0. Differentiating (3.44) with respect to t gives

2 < >
£5.12; = 1,002 ept + 2<x>ept _d_.x__
dt dt
or
z
9.15.; .-. p<x>2.Pt- ak<x><xz>ept +2k<X>Zept
(3.46)
By equating (3. 37) and (3.46), setting t = 0, it can be shown that
-l
x
0
Therefore
<x2>=<x>‘z exp i2(-¥-§;—1-)kt} (3.47)
o

For large x0, that is, as xo-—>- oo

<xz> —->- <x>z

x - l
since—44L;— ——>- 1 and kt O(—-.-—) —-—>— 0

9‘0 x0

as expected.

Now consider the reaction A + B -—>- C, for the case A0 = B

0'
Multiplication of (3.25) by x and x 2 respectively, and summing
over x we obtain the set of equations
25%.. = _ k < x3 >
(3.48)
d < x z > 3

By expressing the higher moments in terms of the lower moments

(approximate method 1) the set of equations (3.48) have the solutions

30

 

(3.49)

When (3.44) is substituted into (3.48) (approximate method 2) the
first and second moments of X(t) are

X
0

(X > = 1+ :50 [exp(kt)-l] (3°50)

 

<x"->=<x>‘2 expikt} (3.51)

The mean and the coefficient of variation have been calculated
from the approximate methods for the bimolecular reactions 2A-->— B
and A + B -->- C (A() = B0). The results are given in Tables III and
IV together with the results calculated from the exact solutions. The
subscripts s, 1, and 2 refer to the results obtained from the exact
equations, and the results obtained from approximate methods 1 and

2, respectively.

2. Application of Method One

 

a. Reaction: Radical Chain Polymerization

 

Olefinic substances, such as ethylene and styrene, usually
polymerize by free radical mechanisms. The reactions sometimes
occur in the gas phase and sometimes in solution; the free radical
processes are similar in all cases. The nature of the initiation process
varies with conditions (27). Sometimes, as with ethylene, the initiation
may be of the thermal type, the molecule forming the diradical “CHZCHz—.
In the case of vinyl acetate, for example, light of certain frequencies

favors the production of radicals and causes the substance to polymerize

at ordinary temperatures. Polymerization can also be induced by the

31

 

 

 

Table III. Fraction of Reactant and Coefficient of Variation for the
Exact and Approximate Methods. 2A-—>- B
<xs> <x1> (XE)

X0 kt x0 x0 XO CV5 CV1 CV2
x0= 10

0.25 0.815 0.818 0.802 0.196 0.213 0.214

0.50 0.686 0.695 0.684 0.258 0.304 0.306

0.75 0.593 0.606 0.592 0.300 0.375 0.377

1.00 0.521 0.539 0.513 0.332 0.435 0.444

1.50 0.421 0.444 0.411 0.382 0.542 0.548
x0 = 50

0.25 0.803 0.803 0.801 0.091 0.099 0.101

0.50 0.670 0.671 0.670 0.119 0.141 0.141

0.75 0.576 0.577 0.575 0.137 0.172 0.172

1.00 0.504 0.508 0.502 0.152 0.199 0.200

1.50 0.404 0.408 0.401 0.176 0.244 0.246
x0 = 100

0.25 0.801 0.801 0.801 0.059 0.070 0.070

0.50 0.669 0.670 0.668 0.082 0.100 0.100

0.75 0.574 0.574 0.573 0.097 0.122 0.122

1.00 0.502 0.503 0.501 0.109 0.141 0.141

1.50 0.402 0.402 0.401 0.125 0.173 0.173

 

32

Table IV. Fraction of Reactant and Coefficient of Variation for the
Exact and Approximate Methods. A + B —>- C A0 B

 

 

0
(XS) < X!) <xz >

xO k t x0 x0 x0 CVS CV1 = CV2

x. = 10

o

0.25 0.798 0.800 0.799 0.121 0.158
0.50 0.662 0.667 0.661 0.192 0.226
0.75 0.565 0.572 0.562 0.226 0.278
1.00 0.491 0.500 0.487 0.250 0.324
1.50 0.389 0.400 0.382 0.293 0.390

xo==50
0.25 0.799 0.800 0.799 0.056 0.071
0.50 0.664 0.667 0.664 0.091 0.105
0.75 0.567 0.572 0.567 0.101 0.122
1.00 0.498 0.500 0.497 0.105 0.144
1.50 0.398 0.400 0.396 0.125 0.174

X0 = 100
0.25 0.800 0.800 0.800 0.037 0.050
0.50 0.666 0.667 0.667 0.057 0.071
0.75 0.571 0.572 0.571 0.070 0.086
1.00 0.499 0.500 0.498 0.076 0.100
1.50 0.399 0.400 0.398 0.088 0.122

 

33

introduction of free radicals into a monomer; hydrogen atoms and
methyl radicals, for example, bring about the polymerization of
vinyl acetate. In all the above cases the propagation and termination
steps are of the same character, and a general scheme of reactions
will be considered, without the nature of the initiating reaction being
specified.

The type of propagation reaction that is involved is the addition

of a radical to a double bond:

R + R' CH = CHz-—>- R' CH-CHz R

This produces another radical, which in turn can add on to a double
bond. This process can continue with the formation of large radicals,
which finally react with another radical to give a molecule. The
general scheme for this type of polymerization reaction can be repre-

sented as follows (28);

Initiation ——->- R1
M + R1 "_>- R2
Propagation M + Rz—>- R3

M + RIF-9. Rn+1

Term1nat1on Rn + Rm "9’ Pn +m

Here M represents a monomer molecule; R1 an initiating free radical,
and Rn a radical consisting of R1 added to a chain of n-l monomer
molecules. The rate of initiation of R, will be written as v, and the
rate constants for the propagation reactions, assumed to be the same
regardless of chain length, will be written as kp. The termination step

involves the mutual annihilation of two radicals with the rate constant kt°

34

Let P ( {r3 , m, Sij} , t) be the probability that the system

is in the state Ri(t) = r

i' M(t) = m, and Pj(t) = pj where i = 1, 2, 3, . . . ,

j: Z. 3. . - . ;and I'j. m, Pj = 0. 1. 3. ... . . The transitions

between different states in the reaction system are due to the occurrence
of the various elementary events (reactions) which form part of the
reaction scheme of the polymerization process. The elementary events,

Ek, in the radical chain polymerization process are:

E1: 1'1-1 -->-r1

E2: m+1,rn+l, r -l—->- m, rn, rn+1

n+1

E3: rn+ l, rm+ 1, pn+m - l—>~rn, rm, pn+m

E4: m, r pj-—>- m, ri, pi

i:
The transition probabilities for these events are

P(Ez) = kp(m+l)(rn+l)At + 0(At)

kp(rn+l)(rm+l)At + 0(At), when n 7! m

P(E3) :{%-kp(rn+2)(rn+1)At + 0(At). when n = m

P(E4) = 1 - vAt - kp(m)(rn)At - {kp(rn)(rm)At + 0(At)
i—kp(rn)(rn-1)At

From the above we get the relation

P( {r1}, m, {p93, t+ At): vAt P({ri} ' , r1 - l, m, {pig , t)
00

I
+ kp(m +1)At Z (rx+1) P( {rik , rx+1, rx+1-1, m+1, {pj} , t)

x=1
00

I I
+kt (r1+1)At 2 (rx+l) P( {r11 , r,+1, rx+l, m,{pj‘( , pr-l, t)
x=2

(I) I t
+ k?(rz+1)At g3 (rx+1) P( {r95 , rz+l, rx+1, m, {p3} , px+2 -1, t)

m I
+ . . . +i- kFAt E (rx+2)(rx+1) P( {r5} , rx+2, m, {pj},, pzx - l, t)
x=l oo oo 00
+ [ 1- 12412 - kp(m)At 2'3 (rx) - kt(r1)At )3 (rx) - kt(rz)At E (rx)
x=1 x=2 x=3

oo
—. . . - -},_- kyAt 421”") rx-1)] P((r.} . m. {P11 . t) + 0(At) (3.52)

35

Where P({fj} , m, {pj} , t) = P (r19 r2! 9 e o 9 m9 p29 p39 0 o o 9 t)

and P( {r25 , r1 - 1, m, {p}; , t) = P (r1 '19 r29 r39 0 0 0 9 m9 P29 p3
. . . , 1:)

By transposing the term P( {r1} , m, {p} , t), dividing by At and taking

the limit one gets the differential difference equation

dP({r3 , m, {p} , t)

 

; dt = VP(£riI}'9 1'1'1, m9 {Pj} 9 t)
00
+ kp(m+l) 12 (rx+1) P({ri’lfl rx+1, rx+1-1, m+1, ([le , t)
00
+ kt (r1+1) Z (rx+l) P({riifl r1+1, rx+1, m, {pfis ', px+1- l, t)
2
00

+ . . . +%- kt E (rx+2)(rx+l) P({ri} ', rx+2, m, {pfifi pZX-l, t)
1

oo

.. vP({ri} , m, {pj‘g , t) - kp(m) Z (rx)P({ri‘g , m, {pjl , t)
1

oo
- kt(r1) i (rx) P ({IBS , m, {pj} , t) .. , . . ,
oo
- i- kt >31 (rx) (rx+1) P ({ri} . m. {133-}. t) (3.53)

Multiplication of (3.53) by p2, p3, . . . , pn, . . . , respectively,

followed by summation gives the equations

d<P2> 1 2 1
T=fkt<rl>-§-kt<rl>
d<23>=

dt kt<r1rz>

< >
i525__:kt<rlr3>+%-kt<rz(rz-l)>

and generally for n even

n-1
if} >=l . , _l.
dtn— :- kt J§1<r1rn1> 2 kt<rn/z> (3.54)

and for n odd

36

n—l
d < Er: > = 1. < .
dt 2 kt jZil rJrn-j > (3. 55)

Multiplication of (3.53) by pg, pg , . . . , pnz, . . . , respectively,

followed by summation gives the equations

d < 2 >
+ = 2' kt <2 p2r1(r1-1)-r1(r1-1)>

< 7' >
.d_:1_tB3___ = kt <2p3rlr2+r1rz >

and generally, for 11 even

n-l

kt < (2pn+1) {g- s rjrmj} > (3.56)

d<pr7i> _
(it

and for n odd

n-l

d< z>
-——dz;n——— 3 kt<(2pn +1) i%‘ 21 rjrn_j " i“ rn/Z > (3957)

Putting < pn Z rj n-j > = < Pn > <2 rj ran >

and _
< Pnrn/z > " < Pn > < rn/z >

we get from (3.54), (3.55), (3.56) and (3.57) that

 

2
d<p > = {2<Pn> +1} d<d2tn>

 

dt
or
Z 2
d<2n> = d<p9> +d<2n> (3.58)
dt dt dt

By imposing the initial conditions, < pn > = 0 on (3.58) and solving

we obtain the variance of Pn(t)

Var (pn) = (pn> (3,59)

37

In a similar fashion multiplication of (3.53) by m and ma,

respectively, followed by summation gives the equations

d < m) 00
dt kp m Elrn (3.60)

In most radical chain polymerization reactions m > > rx so that (3.60)

and (3. 61) may be written as

d<m> ‘30

= - k < > < > 0
(it p m ? rn (3 62)
< 3> 00 oo
£1———r5‘——=-2k<riri”><>:i~n>+1< <m><2r> (3.63)
dt P 1 P l n
Multiplication of (3.53) by r1, r2, . . , rn, . . ., respectively

followed by summation yields the set of equations

d<r> 00
—a-t-l—-=v-kp<mrl>-kt<r1.i‘3 rn>+kt<r1>
d<r >

71L=kp<mrl>-kp<mrz>

oo
'kt<1'2 E rx>+kt<rz>
1

. . . (3.64)

Summation of these equations gives

(I)
d<Z rn>

1 00 2 0°
dt =v-kt<(?rn)>+kt<213rn> (3.65)

 

38

Note that the deterministic counterpart of (3. 65) is

 

CD
d[21rn] (I) 2
dt - v-kt[231rn] (3.66)

When the steady state condition is imposed on (3.65) then
00 oo
v=kt<(2rn)2>-kt<21rn> (3.67)
1
In order to evaluate Var (m) we must assume that

V: g‘2 = k <Zr >2 (3.68)

which is the deterministic steady state condition. From (3.62),

(3. 63) and (3.68)

d < m >
cit p g m (3 69)
d < m2 > _

If the system reaches the steady state rapidly, for all practical pur-
poses g can be considered to be constant at t: 0. This assumption
permits solution of equations (3. 69) and (3.70). The mean and

variance of the variable M(t) are

<m > = m0 e-gkpt (3-71)

Var (m) = moe”gkpt (1 - e-gkpt)

b. Reaction: Condensation Polymerization

Here we consider a mechanism involving successive addition of
monomer to the same growing chain. A reaction which is an example
of the type to be considered is the synthesis of polypeptides by poly-
merization of N-carboxyanhydrides of amino acids, the simplest

possible mechanism of which is (29)

39

o
//
RCH-C \ RCH-CO-R'
( O + R'H—-> \ +COz
NH-C / NHZ
\\
o
,o
RCH-C/\ RC‘SH-COR' RCH-CO-NH-RCH-CO-R'
_>_ +
i /O + NH; NH; C02
NH-C\\
0
etc.

We shall consider a stochastic model for a mechanism represented by

the following scheme:

A + M1-—->- Mz
A + Mz—"> M3
A + M3—>' M4

A + MX——>- Mx+1

Here let A(t) and Mx(t) represent, respectively, the number of A and

Mx (x = 1, 2, . . .) molecules at time t. The rate constants for each

step in the process are assumed to be the same regardless of chain
length. The stochastic model of this scheme is defined by the follow-

ing differential difference equation.

(I)

dP(a-9imj19t)= k (21+ 1) 21(mx +1)P(a+ l,{mj}‘, mx+1, t)
x:

(it

 

oo
- k (a) 21 (mx) P(a, imj} , t) (3.72)
x:

whereP(a, {In}, t)= P (a, m1, m2, . . . , mx, . . . , t)

anda=0,1,2,...,ao;mj=0,1,2,...,;j=1,2,3,....

40

Let C() be the initial number of M1 molecules. The mechanism
suggested by the condensation scheme conserves molecules of the
type Mx' that is, each molecule M. destroyed is replaced by a mole-

J
cule Mj 4; 1. Therefore Co represents the total number of Mx mole-

oo
cules, i. e. , E mx = C0. The rate of change of A is obtained by

multiplying equation (3. 72) by a and summing. We have that

. 00
< >
LL=_k<aZm >= -kC <8.) (3.73)
dt 1 X 0
with the solution
<a> = a0 exp {-kcot} (3.74)

Again, multiplication of (3.72) by a2 followed by summation over a
gives

d<az> _ 2
___Tfi__..-2k<a>Co+k<a>Co (3.75)

It can be shown that

-2
z e kCOt

-kCot -2kCot
0 e

<az>=a +a0e -a0

and therefore the variance of a is

e-kCot -kCot

Var (a) = a0 (l-e ) (3.76)

Multiplication of (3.72) by m1, m2, . . . , respectively and summing

gives the rate equations

d<m1>

dt = - k <am1 > (3.77)
and for x>1
d< >
___I§J{_=k<amx_l>-k<amx> (3.78)

dt

41

In most condensation reactions of this type the species Mx do
not reach a steady state, therefore, the steady state condition cannot
be used. Instead, each of the rate equations must be solved to give
each < mx > as a function of time. The standard method of treating
(3. 77) and (3. 78) is to transform these equations so that in! in =
a - <a>

—0—C—-——— +1, is the independent variable instead of t. in is the degree
0

of polymerization (30).
Then we have

- l
= .. —- < > = < > .
d Xn Co d a k a dt (3 79)

Combining (3.7.9) with (3.77) and (3. 78) we have

 

 

d<m1> _ (Ulla)
an .. <a> (3°80)
d < m2; > = <am1_1 > - < a mx > (3.81)
an <a>. <a>

No simple solution can be found for < m1 > and < mx > unless some
relation between <amj > and < a > < mj > is assumed. In most cases,
especially during the early stages of the reaction process, < a > > >

< mj > so that < am- > < a > < mj > . Equations (3.80) and (3.81)

J
become
d<m1> .
—————-——=—,< >
di’n “‘1
. 2

—————————d<mx>=< > < > (38)
(1)—(n mx_1 - Inx

With the initial condition in = 1 and < m1 > = C0 the equations (3. 82)
have the solution
-1 _
Co Yx e y
(x-l)'.

 

<1'nx> = (3.83)

where y: Xn - l.

42

Similarly it can be shown that

d < m} >
dt

~2k<am§> +k<am1> (3.84)

and forx >1

d<mxz>

dt

-2k<am§>+k<amx>

(3.85)
+ 2k< amx_1mX > + k<amx_1>

By assuming < amj > = < a > <mj > and <mjmj_1>’R-’v<mj> <mj_1>.

equations (3. 84) and (3. 85) can be easily solved. The variance of

M1(t) and Mx(t) are

Var (m1) = C0 e'y (1 - e'y)

_ x1 (3.86)
Co e yy

(x-l)‘.

'

-<mx>

 

Var (mx) =

c. Reaction: Random Polymer Degfadation

In its most general form the degradation process can be illus-

trated by the scheme (31)

M——>-MX+M

Y y-x

We assume that each link in any chain has the same rate constant.
Let Mx(t) and My(t) be the discrete time varying random variables
representing the number of Mx and My molecules, respectively, at

time t. The transitions between various states are given by the ele-

mentary events,

E1: mx-l, my +1—9- mx, my
E2: mx+1 -—>- mx

E3: mx, my->- m m

X’ y

43

and the transition probabilities by

P(El) = 2k(my + 1) At + 0(At)
P(Ez) = k(mx + 1)At + 0(At)
P(E3) = 1 - 2k(my)At - k(mx)At + 0(At)

The factor two in the transition probability P(El) is necessary because
there are two ways an x-mer can be formed from a y-mer (from
either end of the y-mer). The stochastic differential difference

equation for this degradation process is

 

dP(mx,m ,t)
‘f z - +
2ky__2_3x‘{_1(my-)-1)P(mx l, rny 1,t)

dt
x-l
+ k jEI (rnx +1) P (rnx +1, my, t) (3.87)
x-l
- 2k 23 + k2 P , ,t
{ y=x+1(my) j=1mx} (mx my )

where P (mx,my, t) = Prob { Mx(t) = mx, My(t) = my} .

Multiplication of (3.87) by m5: and mi, respectively, followed

by summation give 5

 

 

d<m >
X =2k2 <rn >-k(x-l)<m > (3.88)
dt y y x
< 2 h
d mX) =4k<m 2‘. m >+2k<2 m > (3.89)

- 2k(x-1) <mxz >+ k(x-1) <mx >

The solutions of (3.88) and (3. 89) depend on the initial conditions.
Here we shall assume an initially uniform polymer, 1. e. , we shall
suppose that at t = 0 all < mx >o = 0 unless x= n. < rnn > = m0 is
then the total concentration of polymer molecules at the beginning
of the degradation process. No molecules with x > n appear at any

time. Then when x = n equation (3.88) becomes

44

d<mn>
dt

 

= -k(n-1) <rnn >

and
< rnn > = Ino exp {-k (n-l)t}
Equation (3.88) for x = n-l is

d < mn-l >
dt

 

= 2k<mn>- k(n-2)<mn_1>

Substituting (3. 90) in (3. 91) and solving we get

-k -2 -kt
<mn_1> = 2 moe (n )t(1-e )

Continuing in this fashion we get that (except where x = n)

kt

e-k(x-1)t(1_e'kt) [z +(n-x-1)(1'e- )1

< mx > -_- InO
Similarly when xan equation (3. 89) becomes

d<mz> 2 ..
tn—+2k(n-1)<mn> -k(n-1) <mn>

 

Combining (3.91) and (3. 94) and solving we get

- ..1 -k -1 ~2k -1t
<m§>= mge 2k(n )t‘l'moe (n )t-moe (n )

and therefore

-k(n-l)t

Var(mn) = moe [ l-e-km-l)t

I

In order to determine Var (mx) we proceed as follows. From

equation (3. 88) and the relation

 

z
d<m-Z-2<m > d<dmt >
x

 

dt
wehave
d<m >2 . z _
dt -——+ 2k(x-l)<mx> - 4k<mx> 2:3,<my>

01'

(3.90)

(3.91)

(3.92)

(3.93)

(3. 94)

(3. 95)

45

d 2k(x-1)t _ 2k(x-1)t
Ei<mx>ze }-4k<mx>§<my>e (3.96)

Similarly (3. 89) can be written as

-—d-- <m >62k(x-1)t = 4k<m 23m >
(it X x y Y

2k(x- l)t

+2k<zm >+k(x-l)<mx>}e (3.97)

Y Y

If we assume < m E m > .'—'-' < m > < 2 m > and subtract equation
x Y Y x y Y

(3. 96) from equation (3.97) we get

i { Var (m ) e2k(x-l)t} 3 {2k 23 <m > + k(x-1)<1m >} eZMX-l)t
dt x y x

Y
(3.98)
When x = n-l integration of (3.98) gives
-k - -k -l
Var (mn_1) = rno {2e (11 2)t - 2e (11 )t} (3.99)
For x = n-2, n-3, n-4, and so on we get
Var(m ) ... m £3e-k(n-3)t -4e-k(n-2)t + e«km—1n}
n-2 0
Var(mn_3)= mo {4e'k(n'4)t ~6e 'kin'm + .2e'l'dn'a)t }
Var(mn_4) = m0 SL5e-k(n-5)t --8e-k(n-4)t +3e-k(n-3)t}
and generally (except where x = n)
Var(mx) :3 mo $‘(n-x+l)e"k(x-l)t -2(n-x)e—kXt
.. +1
+ (n-x-1)e k(x )t}
which reduces to
Var(mx) = < mx > (3.100)

Therefore the variation in mean is of the order of the square root of

the mean.

46

d. Chain Reactions in General

 

The essence of the chain mechanism (32) is the existence of
unstable intermediates or chain carriers that take part in the chain
propagation steps and are regenerated. If a single carrier is re-
generated in the propagation step the reaction is called a chain re-
action without branching. If more than one carrier is produced from
the original one, a branched chain reaction is taking place. In either
type of chain mechanism the rate of reaction especially during the
early stages depends on the concentration of chain carriers and their
rate of initiation. A striking feature of some chain reactions is the
occurrence of an induction period during which the concentration of
chain carriers is increasing. Any inherent fluctuations in the concen-
tration of the chain carriers may result in irreproducible reactions.

Instead of presenting a stochastic treatment of particular
examples of chain mechanisms, the present aim will be to treat
generalized schematic mechanisms of non-branching and branching
chain mechanisms. In the branching chain mechanism no distinction
is made between the various reactive chain carriers even though they

may be chemically different (33).

A Chain Reaction without Branching

 

A-l-e- R1
2.
Rl‘i" B'—'>" Ra‘i’P
3
Rz+A—-)‘ R1+I

Here A and B are the reactants, R1 and R2 are the reactive chain
carriers and P and I are the product molecules. Let A(t), B(t),
R; (t), Rz(t), P(t) and I(t) be the discrete random variables. The

stochastic differential equation for this process is

47

dP(a, r1, b, r2, p, i, t)
dt

 

= k1(a+l)P(a+1, r1- 1, b, r2, p, i, t)

+ kz(r1+l)(b+l)P(a, r1+1, b+l, rz-l, p- l, i, t)
+ k3(rz+1)(a+1)P(a+l, rl-l, b, rz+l, p, i-l, t)
- {k1(a)+kz(ri)(b)+k3(rz)(a)} P(a. r1. b. rz.p. 1.1:) (3. 101)

The set of differential equations obtained from (3. 101) for the first

and second moments of R1(t) and Rz(t) are

d <r1>

dt = k1< a > ' k2 <(r1)(b) > + k3 <(I‘2) (a) > (3.102)
2
ELSE—L)— = 2k1< (a)(r1)> + k, < a > - 2k; <(r1)z(b) >
+kz <(I‘1)(b) > + 2k3<(r1)(rz)(a) > + k3 <(rz)(a) > (3-103)
1.552. = k, <(r1)(b) > - k3 <(r2)(a) > (3.104)
d <r§ >
T = 2k2<(r1)(b)(rz)> + kz<(1'1)(b)>

- 2k3<(rz)z(a) > + k3<(rz)(a)> (30105)

During the initial stages of the reaction, when r1, r2, < < a0, be
one can use the approximation

a = a0 b = bo
With this approximation equations (3.102) - (3. 105) contain only r1
and r2 as variables. In order to find <r1> we substitute (3.104) in

(3.102) which gives

d <r1 > d<1' >
dt 8 mo - “"th (3.106)

Solution of (3.106) under the boundary conditions <r1> = 0, <r2> = 0 at
t = 0, is

<f1> +<rz> 3 k1 3.0 t (3.107)

48

Equation (3. 107) can be used to eliminate <rz> in (3.102) to give

d<r1>
dt

 

+ Xo<r1> = k1a0(1-k3aot)

(where X0 = kzbo + k3a0) which upon integration gives

 

<1~1> = k1a° (1 - e’XOt) + 51-53-39— (t+ e‘X9t - 1) (3.108)
x0 x0

The variance of r; can be readily calculated from equations

(3.102) and (3.103). Multiplication of (3.102) by 2 <r1> and since

 

d <r >2 d<r >
52:513.: +2kb<r>z =2ka<r>+2ka<r><r> (3109)
dt 2 o 1 1 0 1 3 0 1 2 '

Similarly, equation (3.103) can be written as

d<r 2)
“3%” + 2kzbo<rf> = 2k1a0<r1> + klao +kzbo<r1>
+ 2k3ao<(r1)(rz)> +k3ao<r2> (3.110)

Subtraction of (3. 109) from (3. 110), followed by multiplication with

. . 2 .
the 1ntegrat1ng factor e g1ves

k3b0t1= d_ [<rl>82kzbot

d 2
— [Var(r1)e dt

dt 1 (3.111)

2k
+ 21.3.40 [<(r1)(r3)>-<r1><rz>]e zbot

In order to find an easy solution to (3.111) we make the assumption

that <(r1)(rz)>= <r1><rz>. Then

2
Zkzbot] = E? [(1:96 kzbot

l

d
a? [Var(r1)e

and with the initial conditions

49

Var(r1)=0 <r1> =0
we get
Var(r1)= (1'1) (3.112)

In a similar fashion from (3.104) and (3.105) we have

d

2k a t 2k a t

a: [Var(rz)e 3 0)= [kzbo<r1> + k.a.<r.> 1 e 3 0 (31131

From (3.104) we see that
d< >

kzbo<r1> + k3ao<r2> s "EEL + 2k3ao<rz>
which further simplifies (3.113) into

(1 21(33. t _. (1 21(33 t

a? [Var(rz) e O] "‘ a? ["--<rz> e 0]
Therefore, the variance of Rz(t) is equal to the mean of Rz(t).

Var(rz) = <r2> (3.114)

A Chain Reaction with Branching
1
A—-—>— R
R + B-—2—>- P + a R

3
R -->- destruction

Here A and B are the reactants, R is the reactive chain carrier, P is
the product and a is the number of chain carriers formed from one

carrier in the chain propagating step. The stochastic differential

difference equation describing the reaction scheme is

50

dP(a, r, b,p, t)

dt = k1(a+l)P(a+1,r-1,b,p,t)

+k2(r- o.+1)(b+1) P(a,r- (1+1, b+ l, p- l, t)
+ k3(r + 1) P (a, r+1, b, p, t)

- {k1(a) + kz(r)(b) + k3(r)} P(a, r, b,p,t) (3.115)

Multiplication of (3. 115) by r and r3 respectively, followed by sum-

mation gives

 

 

‘12:) = k1<a> + k2( o. -1)<(r)(b)> - k3<r> (3.116)
2
d: > =2k1<(a)(r)> + k1<a> + 2( a-1)kz<(r)z(b) >
+ kz(o. -l)z<(r)(b)> - 2k3<(r)z> + k3<r> (3.117)

Putting a = a0 and b =bo integration of (3. 116) gives

<r> = 3:181:9— (l - e-Ct) (3.118)

where C = (k3 - kz(o. -1)bo). The variance, which is calculated from

(3.116) and (3.117), is

- t - z
Var(r) = e ZCt of [klaO +kz( o. -1)2bo<r> +k3<r> ] e Ct dt

ka -2Ct kao -2Ct
= —21-C-Q-[ I-e ]+[k3+kz( Q'l)2bo] {—217.52 (1- e )

_ k2? (e-Ct_e-2Ct,) }

or in terms of the mean <r>, the variance of R(t) is

1 -Ct z ' <r>?-
= — < > 1+ -
Var(r) 2 r ( e ) +Qk3 +k3( a 1) b0} 2k1ao

 

51

3. The Application of Approximate Method Two

In this section the application of approximate method two to
other reaction systems is illustrated. This particular method,

because of its simplicity, is easily applicable to any reaction scheme.

a. Reaction: Autocatalytic A—->- B

Let X(t) be the number of A molecules at time t and Y(t) = Y(o)

+ X(o) - X(t) = Z0 - X(t) the number of B molecules. The relation

Px (t + At) = k(x +1)(Zo - x - 1)At Px+1(t)
+ [1 - kx (Zo-x)At] Px(t) + 0(At)

yields the differential difference equation
...—x: k(x + 1)(Zo-x-1)PX+1 - kx(Zo--x)Px (3.119)

' Multiplication of (3.119) by x and x‘2 respectively, followed by sum—

mation over x gives the differential equations

d<x>

 

 

dt = k<xz> - Zok<x> ’ (3.120)
d<xz> 3 2
dt = 2k<x > + Zo<x> - (1 + ZZO)<x > (3.121)

Substituting <xz> = <x>2 exp {pt} into (3.130) and (3.121) yields

the first and second moments of X(t)

(XQP - Zak) (3.122)

< x > (p-20k + xok) exp {Zokt} ~xok exp { pt}

A

K

V
I

2 -
“<x>zexp R(E—QJT-ch) kt}
O

52

1 z
b. Reaction: A—->- B ——>- C

 

Let A(t), B(t), and C(t) be the number of A, B, and C molecules,
respectively, at time t. The differential difference equation for this

multidimensional stochastic process is

dP(a. 1). C. t)
dt

 

= k1(a+1)P(a+l,b-1,c,t)
+ kz(b+1)P(a, b+1,c-1,t)
- (k1a+k2b)P(a, b, c, t) (3.133)

where P(a,b,c,t) = Prob {A(t) ='- a, B(t) = b,C(t) = C } .

Multiplying (3.123) by a and a‘2 respectively, and summing over a

 

gives
< >
d d: s -k,< a > (3.124)
d< a2 > z
T 2 -2k1< a > + kl< a. > (30125)

Putting < 32 > = < a >3 exp { pt} , equation (3.125) is satisfied when
p = k1 /ao. The first and second moments of B(t) are found in a
similar manner.

-kzt e'k‘t -15"th

 

 

<b>=boe +klao[ k2_k1 ]
< b2 > = < b >Z eXp); 7391(1th kalt }
O

In a similar manner it can be shown that the second moment of C(t) is

<cz>=<c>z exp $33-3t

53

c. Reaction: 2A-‘-—>'- 133-)- c

 

Let A(t), B(t), and C(t) be the number of A, B, and C molecules
respectively at time t. The differential-difference equation for this

process is

 

dP(Zlb,C,t) :1 %_ k1(a+2)(a+1)P(a+2,b-1,C,t)

+ kz(b+1)P(a, b+l, c-l, t)
- [7;- k1 a(a-l) + kzb] P(a,b, c,t)

Proceeding as described above a set of equations for the first

moments is

< >
...—......dd: :-kl<az>+k1<a>

< >
51—5-3— =i—k,<a2>-%—kl<a>-kz<b>
d<c>
—————- = < >

dt k2 b

The second moments of A(t), B(t), and C(t) were found to be

<a2>= <a>z exp {Mag-Li) 1(1):}
0

2
< b2 > a < b >2 exp {klao - ;£a%+ 21(sz t}
O

 

<c‘7‘>=<c>Z exp{-li%PzQ- 1%
o

d. Reaction: 2A..<—_..a B

 

Let A(t) and B(t) be the random variables. The differential dif-

ference equation is

54

+ kz(b+1)P(a-2,b+1,t) (3.126)

- [i—k; a(a-1) + kzb] P(a,b, t)

The set of differential equations obtained from (3. 126) for the first

and second moments of A(t) are

< >
Lfihkfiau +kl<a>+2kz<b> (3.127)
< 2>

+ 4k2<ab>+ 4kz<b>

Putting < a2 > = < a >‘2 exp {pt} in (3. 128‘) and solving for p gives

2
<az> =<a>z exp {2k1ao “2.31.3:33 A4323; t}
o

_ 3....
e. Reaction: A72— B

Here we assume that at time t = 0, A(O) = a0 and B(0) 2 0.

The stochastic differential difference equation for this reversible first

order reaction is

dP(acll’tb’ Q" = k1(a+1)P(a+1. b-l. t) + kz(b+1)P(a-1, b+1.t)

- {k1a+ 1:26} P(a,b,t)

The differential equation for the second moment of A(t) is

d< 2>
—-—d%——— z—Zkl<az>+k1<a>+2kz<ab>+kz<b>

with the solution

k
<a2> = <a>‘Z exp{—1}-}
a0

55

Hence the variance of A is

 

Var (a) = <a>7‘ {expflilt )-1}
o

McQuarrie (11) has shown that the variance of A for this

reversible reaction is

Var (a) = <a>2{ kd1:(:+k)t

" (k1+kz)t

J

+k2 I

During the early stages of the reaction we can assume
e-(kI-i-k‘z)t - l - (k1+kz)t

McQuarrie's expression for the variance of A then becomes

k t
= < 2 1
Var (a) a > Sac“ _ kit);

which is in good agreement with the variance obtained by method two.

 

f. Enzyme-Substrate Re actions

 

Relatively few kinetic studies have been carried out on enzyme-
catalyzed reactions, and much remains to be done before the details
of the mechanisms will be fully known. However, several mechanisms
have been postulated and found to describe certain enzyme-substrate
reactions. It is for these mechanisms that stochastic models were
constructed. Bartholomay (8, 9) studied the Michaelis-Menten model

through a stochastic model but did not determine the fluctuations.

1. Here we consider the Michaelis-Menten scheme

.1...

s+E,,z_c

3
C——>E+P

56

where S and E are the substrate and enzyme, P the product, and

C the enzyme-substrate complex. Letting E(t), S(t), C(t), P(t),
represent the number of E, S, C and P molecules, respectively, at
time t the stochastic model is then defined by the differential difference
equation.

dP(s, e, c, p, t) =
dt

 

k1(s+1)(e+l)P(s+l, e+1, c- 1, p, t)

-k3(c+l)P(s-1, e-l, c+1, p,t) + k3(c+1)P(s, e-l, c+1,p-1, t)

- {k,(s)(e) + 1t,(e) + 143(6)} P(s, e, c,p,t) (3.119)

Multiplication of (3.129) by e and e2 respectively, followed by sum-

mation over e gives

d<c>

dt = -k1<(e)(s) > + (kz‘i'k3) < C > (3.130)
2
Eli-{:3- = -2k1<(e)z(s) > + k1 <(e)(s) > + 2 kz < (e)(c) >
+(kz+k3) < C > + 2 k3 < (6)(C) > (3.131)

If, as before, we let < e2 > = < e >2 exp I 0.1 t} differentiating

with respect to t gives

d<e2> z d<c>
Cit (11(6) expi It} +2<e>expi1t} (11'.
(3.132)

By equating (3.131) and (3.132), setting t = 0, and assuming
E(0) = e0, 8(0) = so, C(O) = 0, P(O) = 0, it can be easily shown

that

0.1 = (So/eo)k1

Multiplication of (3.129) by s and 57‘, respectively, followed by sum-

mation over 8 gives

57

d<s>

dt = R2 < c > - k1<(e)(s)> (3.133)
2
355%: = -21.. < (e)(s)z > + k. < (e)(s) >
+ 2k1<(c)(s)> +k2<c> (3.134)

Again if we let < s2 > = < s >2 exp {(13 t} , equation (3.134) is
satisfied if 0. 2 = (eO/so) k1

It is obvious that the calculation of the variance depends on the
knowledge of the first moment or mean. The rate equations (3.111)
and (3. 114) cannot be solved unless it is assumed that < (e)(s)> =
< e > 4 s >. For most enzyme substrate reactions, especially during
the early stages of the reaction when 8 > > 'e,, such an assumption is
essentially correct. Therefore (3.130) and (3.133) become identical
with the deterministic rate expressions

iLe—L =k1[e}[s]+(k3+k3)[c]

dt

iL94:-ktlc1-ktteltt1

(it
and for early stages in the reaction good approximations to < e > and
< s > can be found (34) (see Appendix A).

The enzyme-substrate systems usually studied are examples
of closed systems; matter neither leaves nor enters the system and
the total mass of the system remains constant while the composition
of the system changes as a result of the reactions taking place. This
conservation condition is usually implicit in the reaction scheme
and can be used to determine a relation between the variances of
the species taking part in the reaction process. For the Michaelis-
Menten mechanism we have the following conservation conditions:

P(t) = so - S(t) - C(t)

E(t) = e0 - C(t) (3.135)

58

From (3.135) we get relations

Var (p) = Var (s + c) (3.136)
Var (e) = Var (c)
Furthermore we know that (12)
Var (s + c) = 2 Cov (s,c) + Var (s) + Var (c) (3.137)
where Cov(s,c) = <(s)(c)> - <s ><c>
Since 8 > > c Cov (s, c) 0 , then

Var (p) = Var (s) + Var (c) = Var (s) + Var (e)

2. The Michaelis-Menten mechanism assumes that only one inter-
mediate is involved in the reaction process. However, in many
enzyme-substrate reactions more than one intermediate is found.
In this section we discuss a process in which n intermediates are

involved. Such a process can be represented by the scheme

I)?

k k
H$ Cn—-E>E+P

k k
E + i =2;- -
8 k1' C1 kz' C2 3 kn'

F

If we again let E(t), S(t), C1(t), . . . , Cn(t), P(t) be the number of

E, S, C1, . . . , C and P molecules, respectively, at time t, then

n

we get

dP(e9 S9C19 ° 0
dt

 

- » end“): k1(e+1)(s+1)P(e+1.3+1.Ci"11”-1Gn'P't)

+ k1'(c1+1)P(e-1, 5'19 Cl'l9 ° ° 0 9 Cn: pot)

+. k'(Cj-1+]-) P (e, S,C1,.o.,Cj_1+1, C

n

1332 1 --1.....cn.p.t)
n

E

J

+ k'j (Cj +1) P (e, 39 C19 9009Cj-1-19 Cj+11-°°'Cn'p't)

i=7-

+ kp(cn+1)P(e-1,s,c1,...,cn+ 1, p-l, t)
n n
- {k.(e)(s) + 14' (c1) + {:72 kj(c,-1) + £2 k1,- (c,)

+ kp {CID} P (C, 8, C19 0 0 0 9 C119 P: t) (30138)

59

Multiplication of (3.138) by e and e2 respectively, followed by sum-
mation over 8 gives

d<e>

dt ='k1<(e)(5)> +k1'<Ci>+lfi3<Cn>

d<ez>

dt = - 2k1<(e)(s)z > + k1 < (e)(S) > + Zki' (190(9) >

+ k'1<c1>+2k <(cn)(e)>+ kp<cn>

P

Again if we put < e2 > =< e >2 exp {(11 t} we get 0.1 = (so/eo)k1 and
similarly for the substrate we have < s2 > = < s >3 exp {(eo/soHfit} .
The results show that the relative fluctuations in the number of
enzyme and substrate molecules during the early stages of the reaction
do not depend on the number of intermediate complexes, but only on
the ratio of the initial enzyme and substrate concentrations. There-
fore even though a reaction process involves more than one inter-
mediate an analysis of the fluctuations in the system can be accomp-

lished through the Michaelis-Menten model.

3. Here we consider an enzyme-substrate reaction in the presence

of an inhibitor (35)

E+1—-> C;
<—

5
Let E(t), S(t), I(t), C1(t), Cz(t), P(t) be the number of E, S, I, Cl, C;
and P molecules in the system at time t. The stochastic differential

difference equation for this process is

60

C113 9 9.9 9 9 1t '
(e S ldtCL CZ E ) = k1(e+1) (8+1) P(e+1,S+1,1.C1’11021P1 t)

+ k2 (c1+l) P(e-l, s-1,i, m+1, cz,p,t)
+ k3 (c1+1) P(e-l, s, i, cl-l, c2, p-1,t)
+ k‘ (e+l) (i+1) P(e+1, s,i+1,c1,cz-1,p, t)
+ k5 (ca-+1) P (e-l, s, i-l, c1, cz+1,p, t)
-~((k.(e)(s) + k.(c.) + k.(c.) + k4(e)(i)
+ 1‘5 (C2)) P(e, 8: it C19 021 P. t) (3.139)

Proceeding as before it can be shown that
< e2 > = < e >’- exp [(11150 + know/e0]
< s2 > =< 3 >2 exp [ (eo/so)k1t]

Using the conservation conditions

e

. E10 + C10) + c.(t)

II

S

O S(t) + C1(t) + P(t)

1. = I(t) + cam

we get
Var (e) = Var (c; + c?)
Var (s) = Var (p + c1)

Var (i) = Var (c2)

4. Enzyme catalyzing two reactions simultaneously (36)

1 3
E+A-—>-C1—>E+PA

(.._..

2

4 6

——>-

E+B Cz—+E+PB

61

The stochastic differential difference equation for this process is

dP(e9 8'9 b9 C19 C29 pA9 PB: t)
dt

 

= k1(e+l)(a+l)P(e+l, a+1, b, cl- 1, c2, pA, PB: t)
+ kz(C1+1) P(e'19 a"'19 b9 Cl+19 C29 pA9 PB: t)

+ k3 (C1+1) P(e'19 a9 b9 CI+19 C29 pA'li pB! t)

+ k4(e+l)(b+1)P(e+1,a, b+1,C1,Cz'1,pA, pB, t)

+ k5 (cz+1) P(e-l, a, b, c1, cz+l, pA, p'B-l, t)
-{k. (e) (a) + k. (c.) + k. (ct) + k. (e) (b) + k. (...) + k. (4)}

P(e,a,b,c1,cz,pA,pB,t) (3.140)

If we again proceed as before the second moments of the variables

A(t), B(t) and E(t) are

< a2 > = < a >2 exp [(eo/ao) klt]
< bz‘> 3 < b >2 exp [ (eO/bo) lgt]
< e2 > = < e >Z exp [(klao ‘1' k‘bOM/eol

From the conservation conditions

a0 = A(t) + C1(t) + PA (t)

b0 : B(t) + Cz(t) + PB (t)

e0 2 E(t) + C1(t) + Cz(t)
We get

Var (pA) = Var (a + Cl)

Var (b + C2)

Var (pB)

Var (e)

Var (C1 + C2)

62

The stochastic models developed here for various types of
enzyrne-substrate systems illustrates again that the stochastic rate
equations are not consistent with the deterministic rate expressions.
Bartholomay (8) in his treatment of the Michaelis-Menten scheme
suggested that the stochastic rate equations agree "on the average"
with the deterministic rate equations. It should be pointed out,
however, that the rate equations agree only in the sense that they
parallel one another. Only under certain conditions, for instance,
in the Michaelis-Menten scheme when < s > > > < e >, will the
equations agree.

The particular interest that is attached to the probabalistic
approach to enzyme-substrate reactions is in the study of reactions
taking place in the cells of living organisms where low concentrations
of enzyme and substrate are found. Such systems have, in fact,
received special mathematical treatment (deterministic) by Strauss

(37) and Goldstein (38).

IV. SET-THEORETIC APPROACH TO REACTION KINETICS

A. Introduction

 

Since the rate of a chemical reaction at a time t is determined
only by the composition of the system at time t and does not depend
on its previous history it was shown that the reaction process could
be treated stochastically as a simple Markov process. The present
investigation of reactions is considered from the point of view of a
set theoretic description of the reaction process (39, 40). Whereas
the basis of the Markov approach depends on transition probabilities,
the set theory approach is based on the assumption that the molecules
in the reaction system can be treated as elements of what might be
called a molecular composition space. The concepts to be presented
in the remainder of this section provide a framework for the so-called
set theory approach to kinetics. We shall begin by defining those con-
cepts of set theory that will be used in this mathematical approach.

When formulating a probabilistic model of a random phenomenon,
such as the composition of a reaction system, the notion of the sample
description space of the random phenomenon must be defined. The
sample description space of a random phenomenon is the space of
descriptions of all possible outcomes of the phenomenon (41). It should
be noted that the word “space" is used to mean a set that is in some
sense complete. For example, the sample description space of the
reaction system A -->- B which initially contains 11 molecules of A,
contains descriptions which contain n elements. At time t = 0 all
the elements correspond to the A molecules only. At subsequent

times some of these elements correspond to A molecules and some

63

64

to the product or B molecules. Such a space is said to consist of
descriptions which are n-tuples. Each of these descriptions corres-
pond to the composition of the reaction system. For the case where
say 11 = 3 the sample description space (or more properly, the
composition space) R, consists of 3-typles of the form (x1,xz,x3),
where x1, x2, x3, are either the elements a or b. Such a space has

8 members

R: { (a.a.a). (a,a.b). (a,b,a), (b.a.a)
(a,b,b), (b.b.a). (b. a. b). um. 11)}

each member corresponding to a possible composition of the system.
At present we shall consider such descriptions as (a, a, b), (a, b, a),
(b, a, a) as being distinct even though they represent the same compo-
sition, namely, two molecules of A and one molecule of B.

Just as in the Markov approach, the aim of the set theoretic
approach is to determine the probability that the system has a certain
composition at time t. Since in the set theoretic description of
reactions several descriptions may represent the same composition,
we introduce the notion of'an event. We define an event as a set‘of
descriptions. In the example cited above the event two A molecules

and one B molecules contains three descriptions

E: i(a,a,b), (a,b, a), (b, a, 61)}

To say an event E has occurred is to say that the outcome of the random
situation under consideration has a description that is a member of E.
The notion of a sample description space and the notion of an
event provide a framework for the development of the notion of the
probability of a random event; in other words, set theory provides a
method for determining the probability that the reaction system has a

certain composition at a time t. The development of the probability

65

function will first be formulated from set theory for reaction systems
in general.

Consider a reaction system in which the total number of mole-
cules remains constant, say the number n. The composition space
of this system, denoted by R, therefore contains descriptions that are
n-tuples (x1, x2, . . . , xn) where the members xj represent a particu-
lar type of molecule in the reaction mixture. Further we define a
probability function on the space R that satisfies the following con-

ditions (42);
1. P(E) Z O for every event E
2. P(R) = 1 for the certain event R

3. P(E or F) = P(E) + P(F) for the probability of the occurrence

of the event E or F

In order to determine the exact form of the probability function P on R
we proceed in the following way.

Let X1, X2, X3, . . . , Xn be ndescription spaces (which may
be alike) on which are defined probability functions P1, P2, . . . , Pn'
We form the composition space R, which consists of n-tuples
(x1, x2, . . . , xn), by taking for the first component x1 any descrip-
tion in X1, by taking for the second component x2 any description in
X2, and so on. In this way all the descriptions in R can be formed
from the spaces X1, X2, . . . , Xn. We introduce a notation to express

these facts; we write

R=Xl® XZ®° 0 '® Xn

which is read "R is a combinatorial product of the spaces X1, X2, . . . ,
Xn.” In order to clarify the foregoing ideas consider again the reaction
A—->- B when n = 3. All the descriptions in the composition space R

can be formed by forming a combinatorial product with the spaces

66

X1 = X; = X3 = {(a), (b)} where (a), and (b) are the descriptions

in X1, X2 and X3. Using the combinatorial product notation we have

R: {(4), (15)} ® {6). (b1) 69 {(a). (b)}

Using the rules for the formation of any or all the descriptions in R
all descriptions in R are, of course, 3-tuples (x1, x3, x3). Since there
are two possible descriptions that can be assigned to x; ( (a) or (b) )
and similarly two for x; and x3 there are 2 . 2 - 2 = 23 or 8 descrip-
tions that can be formed.

In general then for a composition space consisting of n-tuples

we will arbitrarily list the descriptions in X1, X2, . . . , Xn as

_ J J J
Xj '- iDl , Dz , o o o , DNj }

Where j = 1, 2, . . . , n. A typical description in R can then be

written as

( 5,51), Dig) , . . ., 6.51) 1

where forj = l, 2, . . . , n Dijm represents a description in Xj

and ii is some integer l to Nj.
Since an event is defined as a set consisting of descriptions

found in the description space,we can define a combinatorial product

event of the events on the spaces X1, X2, . . . , X IfCl, C3, . . . ,

n.
Cn are events defined on the spaces X1, X2, . . . . Xn,then any event

on R can be written as

C 2 C1 ® C2 ® 0 o n ® Cn
Furthermore we choose the events C1, C2, . . . , Cn so that

P(C) = P(C1)P(Cz) . o o P(Cn)

67

When the event C on the composition space R contains only one
(1) 11“") , Dim)
n

12 000,

description say {D11

then

mph”), D-(Z), . . . . 69);)

12
P({D‘”})P({D1‘Z’})~ P({Di(2)} )

It should be noted that the function P ( 3.1311(1); Diiz), . . . ,

Di( n)

} ) refers to the probability of the occurrence of an event which

is defined by the particular sequence

(1) (2)
(Di-I ’Diz 900°9D'
For example, in the system A -—->- B, the function P ({a, b, c} )
refers the sequence (a, b, a) and not to (a, a, b) or (b, a, a).
Therefore to determine the probability that the system has two mole-
cules of A and one molecule of B, regardless of the sequence we use

the total probability

ProbiA=2, B: l} =P({a, b, a})+P({b, a, a})
+P({a, a, b})

B. Applications

 

1. Irreversible First Order Process

 

Consider the reaction system A-—>- B where at time t = 0,
A(0) == nA, B(O) = 0. Let R be the molecular composition space (sample
description space) where each description is an nA-tuple. Define R as
the combinatorial product of the spaces X1, X2, . . . , XPA' Further-

more, list the descriptions in X1, X2, . . . , X as

nA
: iDlJ 9 DZJ}

68

where j = 1, 2, . . . , nA. Since the reaction system contains only
two types of molecules put D1J = a and Dzj = b for all j. This means
A) in R each xj is
assigned the description a or b, and therefore for each single member

that in an nA-tuple description (x1, x2, . . . , xn

event in R

P({(x1, x2, . . . , XnA)}) =P({X1})P({xz} )° ° ' P({xng )

(4.1)
where P( {xj} ) = P (a) or P (b)

The probabilities P(a) and P(b) may be interpreted as the prob-
ability that there is an A molecule and a B molecule in the system,
respectively. From the deterministic model of the reaction A—59- B
the ratio of the number of A molecules at time t to the number of

A molecules initially is

exp ( - kt)

where k is the rate constant. By the frequency interpretation of

probability w e put

P(a) = exp (- kt) (4. 2)
and

P(b) = 1 - exp (- kt) (4.3)

The probability that the reaction system contains x molecules of A

and nA“ x molecules of B at a time t is the total probability or sum

of the single event probabilities of the form (4. 1) that contain x

elements assigned the description a and nA-x elements assigned the

description b. It can be easily shownthat there are ‘11:.) = nAi
(HA-X)! x!

such descriptions. Therefore, from (4.1), (4. 2) and (4. 3) H '

69

P(A(t) = x) = (21A) [18(4)]x (P(b))““‘*'x

or
- _ nA X _ nA-x
P(A(t) - x) - (x ) [P(a) 1 [1 P(a)] (4.4)
The total probability function given by (4. 4) is identical with the
binomial law. The mean and variance of x, the number of A mole-
cules, with respect to a binomial type probability function can be

shown to be
<x>=nAP(a) (4.5)

Var (x)::nA P(a)[l - P(a)] (4.6)

From (4.2), (4.5) and (4.6) we have

- kt

<x>=.nAe (4.7)

-k -k
Var (x): nAe t (l - e 1:) (4.8)
These results are identical with the mean and variance obtained

through the Markov approach (1, 11).

2. Reversible First Order Process

1
Consider the reaction A T2... B where at t = 0 A(0)_= nA,

 

B(O) =-0. Again we define the molecular composition space R, as a
combinatorial product of the spaces X1, X2, . . . , Xn with the

descriptions in Xj as

X1 ‘iDlJ' DZJ} = 15"")
The probability that there are x molecules of A at time t is
given by the binomial probability law. From the deterministic model

the ratio of A molecules at time t to the initial number of A mole-

cules is again set equal to P(a).

70

- (k1+k3)t
k1 + k2

kle + k2

P(a) =

 

(4.9)

The probability of finding a B molecule amon g the “A molecules

in the system is
nA - A
P(b) = ---—.-———- = 1 - P(a) (4.10)
nA
The mean and variance are

- (k1+1<2)t

 

 

 

< X > : nA (kle +k2) / (kl'i'kz) (4.11)
“(k1+k2)t

nA -(k1+k2)t } { kle +151)
-.- _ .12
Var (x) k1+kz {kle +k2 1 k1+kz 3(4 )

- +
Putting )1 = kl/kz, ct) = )te (k1 kz)t+ 1 in (4.12) gives
HALO (.0

Var (x)-m— (1- “A ) (4.13)

which is identical with the result obtained through the Markov approach
by McQuarrie (11).

3. Consecutive First Order Process

 

1
Consider the reaction A—->- B 39- C, where initially A(O) =
nA, B(O) = 0, C(0) = 0. Define the composition space as the combi-
natorial product of the spaces X1, X2, . . . , anyvith the descriptions in

X1

Xj z {131‘}, DZJ, D3J} = {3, b, C}

Therefore in an nA-tuple description (x1, x3, . . . , an) in R each xj
is assigned the description a, b, or c, and for each single member

event in R

P( ((4.14. ...,an)} )=P({x.1)P((x.} )-~P({xn}5)

(4.14)
where P( 5‘ij ) = P(a), P(b), or P(c).

71

Since the numbe r bf descriptions in the space Xj is greater than
two the probability that there are xA molecules of A, xB molecules of
B and xc = “A - XA — xB molecules of C is given by the multinomial

probability law

P(A(t) = xA, B(t) = x3, C(t) = xC ) =

x .

 

n A x
x I; x [P(a)] [P(b)] B [P(cnxc
A c r (4. 15)
11A _ nA'.
where XA XB xc _ xA! xB! xc'. is the

number of nA-tuple descriptions in which XA elements are assigned
the description a, xB elements the description b, and xC = nA - xA - XB
elements the description 0.

The calculation of the moments of the composition variables
for the process described by a multinomial probability law can be
accomplished through the use of the probability generating function (43).

Let the multinomial be represented in general by

n n n n
P(nl’nz’ °'°’nk) = n1n2n3...nk P11P22.0.Pkk

(4.16)

Multiplying both sides of (4. 16) by SlnlSznz° ° ' Sknk, where Sj is
an arbitrary parameter, and summing over all nj (j = 1, 2, . . . , k)

we define the probability generating function g as

n n n n
g(SI, Sz, 0 o o , Sk) = E P(nl, n2, . . . , nk) SI 152 z. o o Sk k

1(4. 17)
Equation (4. 17) may also be written as
n
_ n . n1 n2 nk
g(ShSZa . . . , sk) — z ulna. . . nk (P181) (P282) . - -(PkSk)

(4.18)

7?.

or

n
g(SI,Sz, o o o , Sk) ={P151'l' P252 + ' ° ° + Pksk} (4,19)

Differentiating (4. 18) with respect to Sj gives

n

—a—g— = ‘ n1. 0 0 0n. n°-l. o o nk

853' Z nJ n1n2° ° ° nk (P151) PJ Jsj J (Pksk)
Setting 5; = $2 = ° - . = 5k = 1 gives

fig. n

asj 1 = 2 nj P(n1,nz ~ . . nk) = <nj>
From (4.19)

5%; =nPj(P1+Pz+°°°+Pk)n=nPj

J
1

since P; + 3+- ‘ ° ' + Pk=1° Therefore

< nj > = n Pj (4. 20)
In a similar fashion it can be shown that

< an > = n(n - 1) sz + nPj (4.21)

< ninj > = n(n - l) 13in (4.22)

The mean and the variance of xA and xB, the numb er of A

molecules and B molecules respectively is from (4.20) and (4. 21)
< XA > = nA P(a) , < xB > =nAP(b) (4.23)

Var(xA) = nA P(a) [ 1 - P(a) ] , Var(xB) = nA P(b) [ l .. P(b)]
(4.24)
P(a), P(b) and P(c) can be determined very easily from the determin-

istic model.

73

The Markov approach to the consecutive first order reaction
requires the solution of a multidimensional stochastic process. The

reaction is defined by the differential difference equation

dP (xA, xB, xC, t)

dt -.- k1(xA + 1) P(xA+1,xB-1. xc.t)

+ kz(xB+l) P(xA,xB+l, xC-l, t)
- (kle+ksz) P (xA, xB,xC,t) (4. 25)

The rate equation containing the derivative of the second moment of

3:13 derived from equation (4. 25) is

d<fo>
dt

 

: < > <
Zkl XAXB + k1 XA>

- 2k < 7' > + < > 4. 26
z XB k2 XB ( )
Equation (4. 26) cannot be easily solved unless some relation between

A
< xAxB> = < XA> < xB > but from equation (4. 22) we have

< x x B> and < xA > < XB > is known. Previously, we assumed that

> (1-3:)

<xAxB> = < xA >< XB
Using this relationship and equations (4.21) and (.4. 23) it can be

easily shown that

> (1,232.)

< z>=< >3+ <

is a solution of (4. 26). _

The application of the set theoretic method to other first order
processes will yield exact solutions for the mean and variance as
long as the total number of molecules remains constant. The investi-
gation of simple second order reactions requires that certain assump-

tions be made about the set theoretic description of the composition space.

74

4. A Second Order Process

 

Consider the reaction ZA-—>- B where at time t = 0 there are

a0 = 2xo molecules of A. By arbitrarily forming x0 pairs of A mole-

cules and by considering these pairs as single units, say X, the re-

action can be treated as the first order process

X ——> B
From part 1 we know that
< x > = x0 P(x)

Var (x) = x0 P(x) [ 1 - P(x) ]

From the deterministic model for 2A——->- B we have

 

p43) - fiil. _ 2291 ._.___3;___._. pqx)
30 X0 1+aokt
then
ao
< > = 2< )3
a x 1+aokt

Var(a) 4 Var(x) = < a >2 i 2kt}

These results are identical with Ishida's (12) treatment of 2A-->- B
using a similar pairing of reactant molecules but using a Markov

treatment (Appendix B).

V. DISCUSSION AND SUMMARY

The theory of stochastic processes has been applied to the
problem of determining the inherent statistical fluctuations in
reaction systems. Such fluctuations, similar to the fluctuations at
equilibrium treated by statistical thermodynamics (44), are inherent
to some extent in all chemical processes. The main purpose for
doing a stochastic calculation was the fact that the stochastic approach
explicitly acknowledges the presence of fluctuations and provides a
method of analyzing these fluctuations by taking as its aim, not the
specification of the concentration variable as an ordinary function of
time, but the assignment of a time dependent probability function to
the concentration variable.

The problem of determining the probability function was
approached by assuming that the course of a reaction could be treated
as a Markov process and also through a set-theoretic formulation of
the composition of the reaction system. The basic notion of a Markov
reaction system was the assignment of probabilities to the transitions
between different states of the system; the transitions between dif-
ferent states being due to the occurrence of the various elementary
events that comprise the reaction scheme. From the transition prob-
abilities a stochastic differential equation was obtained which gave
the rate of change of the probability function. This differential equation
parallels the system of ordinary differential equations of the determin-
istic approach.

Exact solutions for Markov models of the bimolecular reactions
2A-—>- B and A + B —->- C were obtained and compared with the

respective deterministic models to determine the extent to which the

75

76

deterministic models were applicable to small systems. The appli-
cation of the Markov approach to complex reaction systems was also
considered. The extension to complex systems was accomplished by
using mathematical techniques which concentrated on producing
approximate expressions for the first and second moments of the
concentration variables. The first mathematical technique considered,
termed approximate method one, generated a set of equations involving
time derivatives of the first and second moments. These equations
were then solved by assuming that the higher moments could be ex-
pressed as products of the lower moments. The reaction systems
treated by this technique included radical chain polymerization, con-
densation polymerization, polymer degradation, and chain reactions.
The second technique considered, termed approximate method two,

was based on the assumption that the second moment could be ex-
pressed as a product of the square of the first moment and an exponential
function, i. e. , < x2 > = < x >‘2 exp {at} . The reaction systems studied
by this technique were: A: B, 2A—>- B, A+ B-->- C, the auto-
catalytic reaction A—9— B, the consecutive reactions A——>— B —->- C
and 2A—9- B -—>- C, the reaction 2A: B, and various enzyme-sub-
strate reactions.

The set-theoretic approach to the study of reaction systems was
based on the premise that the possible states or compositions of a
reaction system formed a so-called "composition space. " From the
concepts of set theory and a frequency interpretation of probability
a probability function was formulated. The set theoretic approach
was applied to the first order processes A—->- B and A If: B
and the second order process 2A——> B. The results were com-

pared with those obtained from the stochastic Markov approach.

77

There are several directions in which the stochastic study of
reaction kinetics could be extended. First, in the development of
the stochastic reaction models it was tacitly assumed that the initial
concentration of reactants could be determined exactly. This meant
that the variance in the initial composition was zero. However, due
to inaccuracies in weighing, etc. , one can only talk about an average
or mean initial concentration. Therefore, in order to include the
uncertainty in the initial composition one has to consider the prob-
ability that the system has a given initial composition. As an illus-
trative example consider the system A-->- B. The rate equations

for the first and second moments of A(t) are

d<a>
dt

d<az>
dt

= -k<a>
= -2k<az>+k<a>

'When exactly a0 molecules of A are present initially the mean and

variance are

-kt

<a> = aoe

-k -
aoe t(1- e kt)

Var (a; t)

However, when only the mean number of A molecules at t= 0 is
known, say < a0 >, the mean and the variance of A(t) are

-kt
<a> =<ao>e

-k -k
Var(a;t)=Var(a;0)+<ao>e t(l-e t)

where

Var(a;0)=<aoz> -<a0>3

78

Therefore, when fluctuations in the initial composition are present
the variance is larger by a factor of Var (a; 0). An investigation
of the initial random variation of reactant composition would be of
considerable interest and would provide a more complete and
realistic analysis of fluctuations, especially in small systems.

Another aspect of the stochastic treatment of reaction kinetics
that might be of theoretical interest and possibly of some practical
importance would be the construction of stochastic reaction models
in which the composition variable is treated as a continuous random
variable, while still retaining the notion of assigning a probability
function to the composition variables. Such an approach, termed a
Diffusion Markov Approach, has been used extensively to describe
empirical phenomenon that occur in Physics (45) and Biology (46),
but has not been used to describe chemical processes.

The stochastic models studied in this investigation have been
limited to certain reaction schemes. Not all reaction schemes were,
of course, studied, therefore, the methods developed here can be
applied to other types of reactions that might be of interest. The
extension of the stochastic methods presented here to other physico
chemical processes, such as chemisorption, and nucleation or crystal
growth is possible. Furthermore, the stochastic Markov approach
can be used to adequately describe many phenomena which are
executed in space rather than in time: in such cases the parameter
t can stand for volume, area, length, force, and the like, instead of
time. A stochastic treatment of chemical processes dependent on such
parameters apparently has not been attempted and therefore represents

another area that might be investigated.

10.

ll.

12.

l3.

14.

15.

16.

LIST OF REFERENCES

. A. F. Bartholomay, BullI Math, Biophys,, 22, 175 (1958).

. K. Singer, J, Roy. Statist. Soc., B15, 92 (1953).

M. Delbruck, J. Chem. Phys., g, 120 (1940).

A. Quartaroli, Gazz. chirn. ita1., _5__3_, 345 (1923).

. P. J. Flory, J. Am. Chem. Soc., 6_2_, 2261 (1940).

 

. A. Renyi, Magyar Tudomanyos Akad. Kem. Tudomanyok

Osztalyanak, Kozlemenyei _2_, 93 (1954).

 

 

A. F. Bartholomay, Bull. Math. Biophys., __Z_I_, 363 (1959).

. A. F. Bartholomay, Ann. N. Y. Acad. Sci., 26, 897 (1962).

 

. A. F. Bartholomay, Biochem., i, 223 (1962).

K. Ishida, Bull. Chem. Soc. Japan, _3_3_, 433 (1963).

 

 

D. A. McQuarrie, J. Chem. Phys., _3__8_. 433 (1963).

E. Parzen, .Modern Probability Theory and Its Applications,
John Wiley and Sons, Inc., New York, 1960, chap. 2.

 

A. T. Bharucha-Reid, Elements of the Theory of Markov Processes

 

and Their Applications, McGraw-Hill Book Company, Inc. , New

 

York, 1960, chap. 2.

W. Feller, An Introduction to Probability Theory and Its Appli-

 

cations, John Wiley and Sons, Inc., New York, 1960, chap. 17,

E. T. Whittaker and G. N. Watson, A Course of Modern Analysis,
Cambridge University Press, Cambridge, England, 1952, p. 329.

 

A. Erdelyi, Higher Transcendental Functions, McGraw—Hill Book
Company, Inc., New York, 1953, Vol. I, chap. 3.

 

79

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

28.

29.

30.

31.

32.

33.

34.

80

E. T. Whittaker and G. N. Watson, 22. git” p. 330.

J. Irving and N. Mullineux, Mathematics in Physics and Engineering,
Academic Press, New York, 1959, chap. 3.

 

P. Morse and H. Feshbach, Methods of Theoretical Physics,
McGraw-Hill Book Company, Inc., New York, 1953, Part I, p. 751.

 

J. Irving and N. Mullineux, 22. c_i__., p. 180.
P. Morse and H. Feshbach, 22. 93., p. 781.
Ibid., Part II, p. 1754.

Ibid., p. 1756.

A. T. Bharucha-Reid, _op. (33., p. 21.

A. C. Aitken, Statistical Mathematics, Interscience Publishers,
Inc., New York, 1957, p. 22.

 

E. R. Cohen and H. Reiss, J. Chem. Phys., _3_§_, 680 (1963).

 

P. J. Flory, Principles of Polymer Chemistry, Cornell
University Press, Ithaca, New York, 1953, chap. 4.

 

C. Tanford, Physical Chemistry of Macromolecules, John Wiley
and Sons, Inc., New York, 1961, p, 596.

 

Ibid., p. 608.
Ibid., p. 609.
E. Parzen, 22. it.” p. 150.

F. S. Dainton, Chain Reactions: An Introduction, John Wiley and
and Sons, Inc., New York, 1956, chap. 3.

 

J. Rose, Dynamic Physical Chemistry, Sir Isaac Pitman and Sons,
Ltd., London, 1961, chap. 18.

 

M. Dixon and E. C. Webb, Enzymes, Academic Press Inc.,
New York, 1958, chap. 4.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

81

Ibid., p. 173.
Ibid., p. 91.

O. H. Strauss and A. Goldstein, J. Gen. Physiol., _2_§_, 559 (1943).

 

A. Goldstein, J. Gen. Physiol., 21, 529 (1944).

 

P. R. Halmos, Measure Theory, D. Van Nostrand Company, Inc.,
Princeton, N. J., 1950, chap. 1.

 

E. Parzen, 22' git” chap. 3.
Ibid., p. 9.
Ibid., p. 18.

F. N. David and D. E. Barton, Combinatorial Chance, Hafner Pub-
lishing Company, New York, 1962, chap. 2.

 

T. L. Hill, An Introduction to Statistical Thermodynamics,
Addison-Wesley Publishing Company, Inc. , Reading, Mass. ,
1960, chap. 2.

 

S. Chandrasekhar, "Stochastic Problems in Physics and Astronomy, "
Reviews of Modern Physics, _1_5_, 1 (1943).

 

W. Feller, “Diffusion Processes in Genetics',‘ Proc. Second
Berkeley Symposium on Math. Statistics and Probability, 227 (1951).

APPENDICES

82

APPENDIX A

 

Here we illustrate how the Michaelis-Menten rate equations
can be treated to give good approximations for the enzyme and sub-
strate concentrations. The deterministic rate expressions for the

enzyme, the substrate, and the complex are

iii-L=”k1[e][81+(kz+k3)[C] (1)

m=kz[c]-ki[e][s] <2)

dt

9—L9—1—=k1[e1[s1-<k.+k.>[c1 <3)

dt

Application of the steady state treatment gives rise to

_ k1[e][sl
[CJ- k,+k, <4)

From the conservation condition [ eo ] == [ e ] + [ c ] equation (4) can

be written as

_ kle ][S]
”1‘ (kz+1k3;)+k1[51 (5)

Substituting (5) into (2) gives

d[s] g _ antenna]

dt (kz+k3)+kl[S]

 

which has the solution

(k‘l’k) S 1
’t—t‘klk. lni‘dT +1: ‘[81'[Sol>=-t

83

84

During the early stages of the reaction when [ so ] is not much larger

than [ s ], the logarithm term can be omitted, therefore

[S]=[So]'k3t (6)
The concentration of enzyme is

_ _, (k +k)[eol
[el—[eol-[cl- (kzj-ij-fklfs] (7)

APPENDIX B

 

Ishida's treatment of the system 2A—->- B is based on the
assumption that the stochastic differential difference equation for a
unimolecular process can be used to study the bimolecular system.
Ishida proceeds in the following manner (10).

The deterministic rate of the bimolecular reaction 2A —->- B
is expressed in the form

(1 2n)

d1: = - k2(2n)"- (1)

where (2n) denotes the number of reactant molecules at time t and k2
is the rate constant. Solving (1) under the condition that 2n = 2no at

t = 0, we obtain

 

_ 2no
2“ " 1 + kz(2n)t (2)

or if we set 2n = N, this becomes

No

N: l+kzNot

 

(3)

From equation (2) the probability that any one of n reactant molecules
undergoes a chemical transformation during the time interval (t, t + At)

is given by

 

 

I An I _ kzNo
n - 1 + kaNot At + 0(At) (4)
. kzNQ . . . . .
Ishida then sets k(t) = , which 18 substituted into the first

order stochastic differential difference equation

85

86
9%42 = - k(t)nP(n,t) + k(t) (n+1) p (n+1,t) (5)

where k(t) is defined as the "instantaneous" value of the rate constant.

We then have

dP(n t) k N
___...1_.. z _ _.___L_Q... _
dt 1+kzNot n P(n't)
+ ———l—9———k N (n+1) P (n+1 t) (6)
1+k2Not ’

From equation (7) we obtain in the usual manner the rate equations

d<n>_ kzNo

 

__ _ .. < >
dt 1+kzNot n
(7)
d < n2 > 2k N kzN‘J
+ Z O W < 3 > = < >
dt 1+ kzNot n 1+ kzNot n
whose solutions give
< N > = J9— (8)
1+kzNot

Var (N) = < N >2 (2kzt)

APPENDIX C

 

Here we outline the method used by Renyi (6) for treating the
reaction A + B -->- C.

Let the random variables X1(t), Xz(t), X3(t) be the number of
A, B, and C molecules at time t, respectively and let x1, x2, x,
(xi 2 0) denote the values that these random variables can assume.
If at time t = 0, X,(0) = x‘o and Xz(0) :3 x20 then X1(t) = x10 - X3(t) and
x,(t) = x20 - X3(t). Let Px3(t) = Prob {x,(t) = x,} where
x3: 0, 1, 2, . . . , x: and x: = min (x10, x30).

From the relation

Px3(t + At) = [ 1 - 1. h ( x3 )At ] Px3(t)

(1)
+ x h ( x3 — 1)At Px3_1(t) + 0(At)

where h( x3) 5‘— (x1o - x3) (xzo - x3), we get the differential difference

equations

tip 9
€519— : i [h(x3-l) px,-1<t) - tux.) Px3(t)] (2)

a):
X3=1, 2,...,X

185151;)— : " A X10 X20 P0(t) (3)

Renyi solves this system by using the Laplace transformation.

By applying the transformation to equation (2) we get

- .. hix '110 .-1(S)
QX3(S) “ L ) Pam} ‘ 33+ xh (i3) (4)

 

With the initial condition PO(0) = 1 from (3) we have

87

88

Po(t) = exp i- A x10 x20 t} (5)

The Laplace transform of Po(t) is

 

 

90(5) = s + x X“; x20 <6)
Qx3(s) can be obtained by iteration using (4) and (6); therefore
X3 X3
1 ,Tr h (i - 1)
i=1
QX3(S) = X3

1‘20 [s+).h(i)]

By applying the inverse transformation Renyi obtained

. x1
Px3(t) = [1: exp i - i h (i) t} ]

[2: J L: W3

Multiplication of (8) by x3 and x32, followed by summation over x3

(7)

 

gives the first and second moments. The mean < x3 > satisfies the

differential equation

< >
i—fi— = x< (x1) (x2) >
which differs from the rate equation based on the deterministic model

< >
51%: x<x1><xz>.

LHIWHSTH‘ Ublwh’l’

 

f)“; - N

"I7'11111:1111:E11117?