- Email: [email protected]
On the stochastic theory of polymerization reactions
Physica 1llA (1982) 139-160 North-Holland
Publishing Co.
ON THE STOCHASTIC THEORY OF POLYMERIZATION REACTIONS M. MOREAU and L. VICENTE* Laboratoire
de Physique Thiorique des Liquides**, Unioersit& Pierre et Marie Curie, 4, Place Jussieu, Paris, France
Received 3 July 1981
We present a stochastic study of a polymerization reaction model. Using the WKB approximation and the method of factorial cumulants it is shown that the relaxation time of the second factorial cumulant is half the macroscopic relaxation time. Furthermore the system exhibits a globally non-Poissonian and locally Poissonian behavior in agreement with other results. The application of the model is illustrated for the particular case of dimer and trimer formation.
1. Introduction
The mechanism of polymer formation has been studied for a long tirnele3) and applied to a variety of problems such as nucleation in physics and chemistry’) and prebiotic formation in biologyMSu). In the first case the ideas are applied to a class of transformations that occurs by a growth mechanism, examples of which are boiling and freezing of a liquid and condensation of a vapor. So, in the homogeneous gas phase, nucleation of the system is described by a set of coupled differential equations (of the type of later eq. (7)) which generally are not easy to solve. Here it is necessary to know the thermodynamic properties of clusters, a serious difficulty being the accurate estimation of the free energy barrier’). The analysis of these problems is given by a deterministic causal description, provided by the equations of chemical kinetics. Our purpose here is to study a simple model of polymer chain formation according to the methods of stochastic chemical kinetics developed in recent years ‘-14) and extensively applied to study the fluctuations of the particle numbers and the transition to stable steady states far from equilibrium. We only consider a simple model of polymerization reactions in a closed system: the formation of a p-mer, where a p-mer is a polymer composed of p *Fellow of CONACYT and Facultad de Ciencias, UNAM, Mexico. **Equipe AssociCe au Centre National de la Recherche Scientifique.
037%4371/82/OOOO-OOOO/$O2.75 @ 1982 North-Holland
140
M. MOREAU
monomer
units,
number We
results
of monomer begin
studying
by
from
units
linking
results
which
the
macroscopic
evolution
permits
and is necessary
a monomer
conserved.
the case of an inhomogeneous formalism
L. VICENTE
being
analysing
the linearized
AND
equation system
behaviour
of the
for an homogeneous
is treated
a straightforward for the stochastic
The stochastic study is based first applied to the birth-and-death
to a (p-l)-mer,
the total system
system.
in the so-called
generalization
and Then
multivariate
of the
previous
description.
on the master description
equation formalism, which is of polymerisation in a homo-
geneous system. Using the generating function, the equations for the first two cumulants are derived and the asymptotic behaviour of the correlations between the number of molecules is deduced. But it must be noted15) that the fact of considering the system as homogeneous whole neglects all local effects and thus selects the limited class of big fluctuations at the level of the entire system. Therefore we undertake then the most complete treatment, considering local fluctuations by means of the multivariate master equation’.“) (MME) but making use of the foregoing case. Taking the continuous limit for the correlation we conclude that the system is dominated within each cell by Poissonian-type of fluctuations’).
fluctuations,
2. Macroscopic
in agreement
with the phase-space
description
theory
2.1. Homogeneous
closed system
2.1.1. Generalities It will component components
be
advantageous
for
our
purposes
fluid in which chemical reactions X,
to consider
a general
take place which according to
transform
multithe
m
c V”dXP+, p=l
iqx,,
a=1 ,...,
r.
(1)
The symbol X, is used both for the pth component and the corresponding particle number. The stoichiometric coefficients v;(F;) of X, in reaction CY correspond to reactants (products), and K”, Ka are the rate constants. With the forward and backward rates for the r reactions R,=K”Xf;q..
. X3,
l?,=K”X’;q..
.X3,
CW
STOCHASTIC
THEORY
OF POLYMERIZATION
REACTIONS
141
we write the standard equations of reaction kinetics’*ti)*
y$=z,(i$+$)(R,-R,).
(3)
We consider now the model of a chemical system in which a molecule of kind X1 (monomer) reacts with a polymer X, (chain of size p) producing a polymer X,,, (size p + l), p running from 1 to m always. Allowing also for the reverse process we have thus the scheme: X1+X1 e
x2, R’ KP
x*+x,
e
(4)
xp+,, Iv Km-1
x,+x,-,
e p-l
x,,
and this will be subject to the condition N = 2 pX, = total number of monomer units. p=l
(5)
From eqs. (2) and (4) we see that R,=K”X,X,,
m-l;
a=l,...,
~;=fi~,o+8~,,, F;=ap,a+,,
l?a=l-?‘X~+,, p=l,.
. . , m;
(6)
where S,, is the Kronecker symbol. Substituting this into eq. (3) yields d$=-2(K1X;-~1X3-~;
R,-,) -
$+=(R,-,-
(R,-I?,), (R,
-I?,).
In general a complete set of kinetic equations contains coupling of chemistry and hydrodynamics. Thus the hydrodynamic fluxes must be sufficiently small not to affect the chemistry appreciably. These conditions depend very much *Usually the laws of chemical kinetics then the forward rate is written:
r, = R,/fl = k” so that the constant
are written with the aid of the concentrations
v (xiJuq
0)
K” is related to the usual constant
K”=k”fl’-‘,‘P with analogous
relations
xi=X,/n;
k” by (2c)
for the backward
rate.
142
M. MOREAU
on the spatial return
and temporal
to this point
2.1.2. Stationary
when
scales
we consider
AND L. VICENTE
of the observed the effect
phenomena”).
We shall
of diffusion.
state
The stationary
state of the system
is given
by
!!&-J dt
and from eq. (7) X0,= ap(X$P,
(9)
where al= 1.
KP-'
K’ ap=-‘.K’
. -KP-”
p=2,.
. ,171.
(10)
which, together with eq. (5), defines the stationary detailed balance is then realized, this stationary confirmed later. 2.1.3. Perturbation
state uniquely. state is stable
Since the as will be
theory
In general the exact kinetic equations (7) cannot be solved exactly, excepted in some particular cases’). Thus it is interesting and necessary for our further purposes, to linearize these equations near the stationary state. By putting X, = Xi+ Yp we obtain the linear chemical system (II) with
-I, 7;=i+v;=
if p=as2, if p=o1+122, otherwise.
1, 0,
(12)
and 7:=2;
77=--l,
if ~~32;
CO,=K’X~X:=K’X:~,.
It is shown in appendix A that matrix R=(Rz) simple eigenvalues and may be diagonalized. The explicit solution of the linearized equation Y,(t)
=C c,U;I evf, llf0
has only
real,
non-positive,
is (13)
STOCHASTIC
THEORY
OF POLYMERIZATION
REACTIONS
where the sum is taken over all the negative eigenvalues right eigenvector corresponding to q, and
143
7 of R, (UT) is the
4
4
These results hold for any stationary state with detailed balance. In the case of the polymerization reactions an exact computation of the eigenvalues is generally not possible. But bounds may be found, and analytical or numerical approximations may be used*. Instead, we shall now particularize our system to the case m=3. 2.1.4. Particular case This situation has been studied both experimentally and theoretically by Turq16) from the macroscopic point of view. In many practical cases the forward rate constants are nearly the same for both reactions, and similarly for the backward rate constants, so that we can reasonably take K := K2= K, I?‘=Z?*=K. Introducing the scaled variables
x;=p’,
k'X2
k'X3
x;=j-gx.
(15)
The new variables obey equations similar to (7), with unit rate constants: K;=I?b= 1. Dropping the primes in the variables we obtain the following relation at the steady state: X*=(X9’,
x3= (X7)3.
In the neighborhood by setting Xp=X;+Yp,
(16)
of the steady state the evolution equations are linearized
p=l,2,3.
(17)
The linearized system reads %=c
4
R;Y,,
p=l,2,3,
(184
*It may be shown, for instance, that every vanishing eigenvalue is intermediary between two successive eigenvalues of the chemical system obtained when the concentration of monomers is kept constant”).
144
M. MOREAU
where
the matrix
R=(Rz)
is I
2-x’;
-x:(4+x:)
AND L. VICENTE
(1%) i (R;j) are
r),=o,
_ ?j? ,=S[-s+(&4?Y)“‘] _
(19)
where 6=x’j’+5x0+2
I
3 S’=
9x:‘+4x:‘+
1.
If (CJa)= U” is the eigenvector corresponding to qK, it is shown in appendix that (Y,) has no component on U’; thus the solution of Eq. (18a) is
It is interesting to consider the limit cases the first case X:+1; thus it is found that which implies a relaxation time of order I In the second case 7)?=-9 and v~-X$l, relaxation time of the component on U3 component equilibrium trimerc.
2.2. Di$usion
on U’=(O, 3, -2). value much more
A
of low and high concentrations. In qZ and q7 are nearly equal to -I, for both components on U’ and U’. which implies a very short and a relaxation time = l/9 of the
Thus the monomer concentration reach its rapidly than the concentrations of dimers and
in (1 non -homogeneous
closed
system
2.2.1. Cell description The foregoing analysis is a global one. which does not include the effect of diffusion. We proceed now to consider this effect in the deterministic equations by dividing the system into cells, within which the chemical reaction\ take place and between which matter is exchanged whereas homogeneity is maintained within each cell. The conditions are those of the MME and have been discussed elsewhere’). They imply that the length of each cell is at least of the order of the mean free path of molecules and far smaller than the mean free path between two reactive collisions. We assume these conditions to be realized. The interest of this description will be clear in section 3 when we treat the stochastic approach where the knowledge of the deterministic solution in every cell is needed.
STOCHASTIC
THEORY
OF POLYMERIZATION
145
REACTIONS
Then, we symbolise diffusion by a “reaction”. yp,i
SDP Yp,i',
(21)
where Yp,i is the p-mer number within cell i, and i’ denotes all cells adjacent to i; Dp is the diffusional constant which is related to Fick’s constant d, by the relation’)
4=+,
(22)
where 1 is the length of every cell. For the sake of simplicity we shall restrict ourselves to one-dimensional diffusion. The index i runs from 1 to n, 1 and n corresponding to the boundaries of the system and any interior cell i has two adjacent cells, i - 1 and i + 1. Throughout the analysis we will deal with bounded media and the boundary condition considered is zero fluxes at extremities. To this end we introduce two supplementary cells, 0 and n + 1 and the condition reads Yp.0= Yp.I,
Yp.n+l = Yp,n.
(22’)
2.2.2. Perturbation theory around the stationary state The stationary state of the system obviously corresponds to uniform chemical equilibrium and it is stable for the same reasons as previously. The evolution equation linearized near the stationary state is
where K%i = R;6{ + D$;Ai P.1
I.
R = (R;) is the “chemical matrix”, Si is the Kronecker symbol (also written 6ii according to convenience) and (A{) is the “diffusion matrix” whose elements are (A{) = 6i+l,i + 6i-l.j - 26ij = 6(+1+ s{_, - 26{, (A’,) = (A”,) = - 1. In terms of w, and U”, eigenvalues and eigenvectors of (K$), analogous to the solution of the homogeneous case*: Yp.i
=
C
Yp.i
(24) is exactly
(25)
C, ewKtU”.
I(
*But K now represents cell (i).
a pair of indices (p. i) corresponding
to the chemical
species (p) and to the
M.MOREAU
146
It has been vectors,
shown
chemical
ANDLVICENTE
by Lemarchand”) and diffusive:
that
U is the tensorial
CT= VW. In terms
product
(26)
(U,,) = (V,W,). The eigenvalue
equation
is then
C Kfl:{LJqj = 2 [(RzV,)W,fi{ q.i 4.1 = C RzV,Wi q Dividing
of two
of components
+ D,GzV,(A:.
Wi)]
+ D,V, 2 A:W, = wV,,W,. i
(27)
eq. (27) by U,W, we see that it implies
2 A(W, = AW,. 2 (R; + hD&)V,
(29)
= WV,.
4 By eq. (28) h is an eigenvalue of the diffusion matrix (A’,) and W the corresponding eigenvector. The well-known solutions of this eigenvalue problem, labelled by the index k, running from 0 to n ~ I are: Ak = *(cos+i)
Wt = \/~{expJ~+exp-J(i-,I)k?i]
(30)
where
Then eq. (29) allows us to find the set of eigenvalues {wf} and eigenvectors {Vk}. If the diffusional constants are the same for all species, that is D,, = D. we get from eq. (29) kWT - 77+ AkD.
(31)
where {q} are the eigenvalues identical to the eigenvectors (25) may then be written
given by eq. (19) and the eigenvectors V; of the chemical matrix (Rz). The
explicitly
V;l are solution
as
where
(WY]
c and G, is the concentration coefficient Coo vanishes.
at equilibrium.
As in the homogeneous
(33) case the
STOCHASTIC THEORY OF POLYMERIZATION
147
REACTIONS
2.2.3. Continuous limit If L is the total length of the system, the size of each cell is I = L/n, and to consider the continuous system we make si = iL/N -+ s and
Yp,i must be replaced by Yp(s) and the limn+,, X,/n. Thus, instead of eq. (33) we get
Wk(sW&,0) ds/
local
concentration
is G, =
xp qjoL(Wk(s))2ds,
(34)
where
Wk(s) = cos
F.
From eq. (32) we obtain lim A(k)D= _ !?$Z “+m and therefore,
instead of eq. (32) we have
Yp(s, t) = c VJ g c,,,k e(‘-kznzd’L*)tWk(S)+ z. V;jc,,Oe”‘W”(s). k#O ‘1
(35)
If we have the initial condition Y,(s, 0) = Y i(s) (deviation from equilibrium in the first cell), eq. (35) becomes Y,
=cr) V:[c
(VzY,/G,)/(z
+.X0v:[$ and approximating
Yp(s,t) =
(Va)‘/G,)]
CVYC,)/(~ (v:)~/G,)] e” P
(36)
the sum over k by an integral we obtain
C.Vz[T (V~Y,IG,)/(~ (v~)‘/G,)] :zo
e?’ ~oe-(kZ~ZdiLz)t cos?
P
P
V:[iC
P
(VZ/G,)/(C
P
evt{ Jge-S2/4dt}
(w’/G,)]
The first term is the typical result of diffusion and the to chemical contribution. Certainly this result may tinuous description of the spatial inhomogeneity but will be useful for the next section. A more interesting to assume the chemical equilibrium to be nearly
eq’.
(37)
second one corresponds be obtained by a conthe previous discussion approximation would be realized in every cell,
M. MOREAU
148
whereas
inhomogeneity
the present
subsists
approximation
AND L. VICENTE
between
will be sufficient
cells
(see
Appendix
B). However
for what follows.
3. Stochastic theory 3.1. introduction In this section
we want
to take
into
account
the
random
character
of
reactions, since the individual reaction between molecules is a random event and thus the number of molecule (or polymers) is a random variable, and we want to determine the probability distribution by means of the moment equations. First we consider the process in a global way and assume that reactions are spatially homogeneous and simulated by a stochastic process of birth-and-death type13). In fact this is only a first approximation because it neglects all local effects and thus selects the limited class of big fluctuations at the level of the entire system. The main advantage of such a description lies in its simplicity, which permits analytical approaches. The most complete treatment is made by considering local fluctuations, and this is done by means of the multivariate master equation (MME). In the treatment of the MME the system is divided into cells within which the chemical reactions take place (described by a birth-and-death process as before) and between which matter is exchanged by diffusion (simulated by a random walk). This representation’) supposes the characteristic dimension of cells to be of the order of the mean free path. It is considered that reactive collisions lead to a chemical process whereas elastic ones contribute to diffusion. Furthermore it is supposed that elastic collisions are much more frequent than reactive ones, which implies local equilibrium Finally let us note that the time scale of the diffusion process and the time scale of chemical reactions
are well separated.
3.2. Birth-death
theory
applied
to a homogeneous
system
3.2.1. The generating function method Let X, be the number of the p-mer molecules (p= 1,2,. . . , m) and . ,X,,,: t)= P{X,}) be the probability of having X, monomers, P(X,, . .
let Xz
dimers etc. at time t. We shall study the stochastic behavior of the system by the method of the generating function’,“), using the general notation of subsection 2.1.1 before particularizing to polymerization reactions. In the birth-and-death formalism, the evolution equation of any arbitrary smooth function Q({X,}) is’)
STOCHASTIC THEORY OF POLYMERIZATION REACTIONS
$(QWpIN=&I~
149
[Q({~,+~“,-v”,})-Q({~,})IW~((W,})P({X,}> P O1
+ x ~ u [Q({x,+u”,-V”,})-Q({Xp})lWOI((XP})P({XP}), FX P where the W’s are the transition probabilities
(38)
and are given by (39)
(40) According to the remark leading to eq. (2~) K” is used when the kinetic laws are written in terms of molecule numbers and k” when they are expressed in terms of concentrations. If we take Q({X,})=ll, sp, q=(Q) is the generating function*: (P({sp)r t,=(rJ
(41)
SP)
and we obtain, after manipulating eqs. (38) and (41)
(43) Let us write s=(s,, . . . ) s,) (with analogous notation for {X,}, {v;}, . . . ),
(4) Then eq. (43) may be written v=;
K”(s’=-- s+)&cp + 2 I-&s+ - s”“)&. cp. D
(45)
*cpis related to the average values of various powers of X, by acp ( as, >,r,+, = (xlJ1 with 8X, = x,
-
(X,).
(gk?,,,.,,
= @X,SX,) + W,)(X,) - W&xl,
(42)
150
M.MOREAUANDL.VICENTE
This
linear
cannot
partial
be solved
differential exactly.
so-called WKB approximation by Lemarchand”), writing: cp(s, t)=exp[Ns,
equation
However,
with
variable
for large
coefficients
values
in general
of R, we may
used by many authors’.“)
and developed
use the recently
(46)
t)l.
Substituting this approximation into eq. (45) and retaining dominant terms in 0, one converts this equation into a non-linear differential equation of the first order,
which
is similar
to a Hamilton-Jacobi
and analogous
equation:
for &a.
We may write eq. (47) in the form +)=K(s,
(48)
$“),
where K(s,
&)=c
Switching
a
k”(s”“-s”“)&+~
,$“(s”“-s~“)$QP.
(49)
to the variables (50)
&=ssp-I $ and K in powers
and expanding equations
of the distribution
of &, we obtain
factorial
from eq. (48) the evolution
cumulants:
dG,-K dt@$=T where
(R4,‘G,,,+R4,‘G,,,)+K,,,
p,q=
1,. . . , m,
(51b)
we have used the definitions
K=!E P
with
(Sla)
”
(86,>&=O’
t&g,,,
&q=(&,,, )=o) RZ= p q B
p
q
)=O’
P
(52)
STOCHASTIC THEORY OF POLYMERIZATION
151
REACTIONS
The factorial cumulants of order one and two are
G=?!k P
(atp)
t5p)=o
_ 0 ’
(534 (53b)
,=,=h [Wp~&>-Wp>4ql. Gpq=(&),, p 4 P
Thus G, is just the concentration and G,, is related directly to the deviation of the distribution from the Poissonian form’). From eq. (49) we obtain
Kp=C
(Wa-Ga)T;,
a
K,,=C (cd”-
W”)( q+“g- v;v”s-Ta~pq),
LI
R;=x
T;(w=v:-G%;)/G~,
(54)
OL
where T; is given in (12), and wa=ku n(GP)“B, P
O”=l;”
n(Gp)'?
(55)
P
Since G, = X,/n (eq. 53a), wa and W” are identical to the deterministic forward and backward rates r, and fa and (51a) indeed is the deterministic equation. 3.2.2. Evolution of the second cumulant If we assume that the macroscopic concentrations their equilibrium stationary values Gp, we have
have already reached
(56)
01a = w”(&‘w”((;‘),
&{Gi}.
Then the evolution equations (51b) are replaced by the homogeneous
(57) system
where &;q=ti;G;+&i3;. The properties of the asymptotic matrix a have been studied previously the forma1 solution of eq. (58) is G(t)=eR”‘G(0),
(59) and (60)
M. MOREAU
152
where e
W = (Z?yq). Using
the tensorial
AND L. VICENTE
product
we may write
RI = efif @ ,It.
(61)
If (7) are the distinct the eigenspace
eigenvalues
of 7 it is known
of k and P, is the projection
operator
on
that
(62) Then
the solution G,,,(t)=2
of eq. (58) may be written
(63)
c,, 1,e(“+““UaUi,
c~,=[z‘;i)O:G,.,(O),G~G~,,~
(U;U:)‘lG,G,]
P4
where
U”.”
as
=
~‘,
(64)
Pq
(UziJ
=
are the corresponding
are the eigenvectors eigenvalues”. Ix).
(UZU:)
of matrix
R and A,, ,. = rl+u
3.2.3. Stationary state If the evolution matrix 6 is regular at the stationary state, all its eigenvalues are strictly negative, and G,,(t) tends to 0 as t + m. Then the distribution becomes Poissonian and uncorrelated in agreement with the thermodynamical equilibrium of a grand ensemble. However, if some quantity is conserved in all reactions (and this is indeed the case for the number of polymer units in our model), i is singular and GP4(t) tends G&t)+
to CO,OU~U~ according Co,opq~&
to (63). Here
(t +m).
(65a)
with
(65b) which generally does not vanish. Thus the stationary state is generally non-Poissonian, unless of course the initial state is Poissonian. It may be observed that the variance of the total number of monomer units, N=C,pX,, is
W’J)2)=~q Pq@WXJ=fJ~
P4
pq(G,,+G,6,,),
(66)
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
153
so that by (65)
co.o=-cd’Gp+~(;~~p)l P
P
The asymptotic
distribution
((W’)=fl
c p2G,=C P
is only Poissonian
(67)
if
P’(&>.
(68)
P
In particular if N is known exactly, ((SN)2) vanishes and
Co,o=-(7 p2Gp)-'CO. Such deviations from Poissonian equilibrium have been predicted by Van Kampen20) and studied recently by numerical computation2’). They are clearly due to the incompatibility of the initial condition and the Poissonian distribution, which cannot disappear because of the conservation relation. 3.2.4. Asymptotic evolution On the other hand, if the macroscopic evolution matrix k of a chemical system is regular at the stationary state, the smallest (in absolute value) eigenvalue of the evolution matrix R for the second moments is twice the smallest eigenvalue of R, so that the relaxation time of (GP4) is half the macroscopic relaxation time, in agreement with other results’?. If k is singular, which is the case here and more generally when there are conservation relations the smallest non-vanishing eigenvalues of k and W are identical; however it is shown in appendix A that the previous result still holds. Thus the relaxation of the second moments is faster than the macroscopic relaxation to the stationary concentrations. However it should be remembered that the system is supposed here to be in macroscopic equilibrium. Otherwise the complete non-homogeneous eq. (51b) replaces eq. (58) and the evolution of (G,,) exhibits terms induced by the macroscopic evolution. 3.3. Multivariate
formalism
of a non-homogeneous
system
3.3.1. E#ect of diffusion Without writing the MME we can write immediately the equations for the first two moments, adding the effect of diffusion (random walk) and remembering that chemical reactions are described be the same birth-and-death process as before, but considering now the equations for every cell. We have
M. MOREAU
154
AND L. VICENTE
~=D,(G,i+i+G,,i~,-2G,i)+~~,i where
p denotes
p=l,m:i=l,..., and i the cell number.
the species
again the macroscopic 23). In the multivariate
solution
which
formalism,
p,p’,q
(70)
This equation
we have examined
in section
represents 2. (see eq.
eq. (58) becomes
. . ,111, i,j=l,...,
=I,.
n,
(71)
n,
where p, q and p’ denote species, i and j cells; RP,’ is again the matrix of linearized macroscopic equations. In the case of our system (cf. eq. (4)) and for equal
coefficients
GP,i,q.j(f) = 3
of diffusion
C,,kGv.h e
the solution
(v+u-(hk+Ah)D)f
v;l
of eq. (71) is
w;yy;, (72)
,a~h
with c7.k;v.h
=
z
Es
p,q,i,j
where
c
W;WfGp~i;q,j(0)
Gp
/[
Gq
q and v are eigenvalues
of the chemical
z
2
-&Wk)2
G, 1 [ i
p
matrix
’
’ 1 ’
and Ak, Ak eigenvalues
(73) of
diffusion; V;l, Vi are the corresponding chemical eigenvectors and W!, W: the corresponding diffusion eigenvectors. From eq. (72) we see that Gp,i.q.i state given by evolves toward a stationary
(74)
Generally this expression sonian. But if we consider treatment gpq(h
in section s2;
t)
=
z
2.2.3) we obtain v;v;
2
L
c p.9
does not vanish unless the initial condition is Poisthe limit of the continuous system (in analogy to the
v; v; G G
e
(t)+~~(k2+h*)nZdlL2)1~k(Si)~h(S2),
(75)
L
II 0
%k;+h
wk(sd
w”(sdgpqh,
~2; 0)
dsl
dsz
0
(76)
&k:v,h= [T
~]2[/oL(Wk(sd2
&I2
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
155
where g&l,
~2;O)=lim GP~i;~(“)=(SXp(S,)Sxq(S*))-(Xp)S~qS(S2-S,) n-o
(77)
and Sx, =x,-(x,). Thus g,, tends to the uniform value L &&It
s2;
aJ)
=
9
c P4
P4
II
L
gpqh,
~2;
0)
dsi
ds2.
(78) 2
As for homogeneous systems, the right-hand side of eq. (71) may be related to the variance ((SN)*) of the total number N of monomer units contained in the system, since L
((6N2)
=
c
pq
L
1
/
0
0
(S(x,(s,)Sxq(sz)) ds, dsz.
(79)
P
Thus, eq. (77) and eq. (78) yield
If N is exactly known, ((SN)2) vanishes, and the asymptotic value of g,, is given by the last term of eq. (80) which shows that the concentrations of the different chemical species are correlated and non-Poissonian. However for large systems (L + w), this term tends to 0, and the fluctuations become locally Poissonian, although the system is globally non-Poissonian, since LL
II 0
g&l,
~2;
~1
dsl
ds2
=
@~p~~q)-
Wp)&q
0
= _
PdWW,)
T
P~WP>
(81)
.
More generally if x,(sJ and xq(s2) are S-correlated, or at least if their correlation has a very short range, it is seen from eq. (77) and eq. (78) that g,,(sr, s2; 00)is of order I_-‘, and thus vanishes for large systems. Equivalently, this is true if the variance ((SN)2) of the total number of monomer units of the system is of order I., which corresponds to ordinary conditions.
M. MOREAU
156
AND L. VICENTE
Then the behaviour of the system is locally Poissonian, and globally non-Poissonian, in agreement with other results obtained analytically’5~2”) or by numerical This
simulation*‘).
conclusion
correlation
does
between
((6N)*) is of order a non-Poissonian
not hold
some
x,(s,)
L’; th en g,,(s,,
near
an instability,
and
xq(sZ) is of order
s2, m) remains
where
the range
of the
L, or equivalently,
finite, and the system
exhibits
local behaviour”).
4. Discussion The major difficulty in the study of chemical correlations appears generally in the approximations used to truncate systems of equations involving arbitrarily high moments. Indeed, as was pointed out by Nicolis’) the higher-order moments become increasingly important in the absence of asymptotically stable microscopic states. However this is not the case for the model used in this paper. Thus the use of a WKB approximation is justified according to the discussion of Kubo et al.‘*) in the limit where the volume 0 is considered to be macroscopic. The present study of a chain of polymerization reactions may be generalized in various directions. We have determined the behaviour of G,, supposing that the macroscopic concentrations have already reached their equilibrium stationary state. It is possible to apply the method to systems outside macroscopic equilibrium, and even to systems far from macroscopic equilibrium if the macroscopic evolution is known. Here it may be done -at least formally - near macroscopic equilibrium. In this case the inhomogeneous terms K,, of eq. (51.b) for the second cumulants G,, does not vanish; furthermore the matrix (Rz) is no longer constant. Then G,, indeed tends to the stationary value computed in section 3, but the inhomogeneous term K,, induces a coupling between the proper relaxation of the second cumulants and the macroscopic
relaxation,
as will be shown
more
generally
elsewhere.
Another possibility of extending the present method is outlined in appendix B: the system is supposed to be near local chemical equilibrium, but far from the stationary state of diffusion, that is, far from homogeneity. In this case the system is easily linearized since the diffusion-terms remain linear in the domain of validity of the multivariate master equation, and the previous results are partly generalized. It may be noticed that the correlations should be taken into account in the global equations of evolution for the average values, although they disappear in the local equations for concentrations. It would be interesting to introduce the effect of an electric field in our
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
157
model, namely an internal field due to charged particles, in order to study chemical reactions in electrolytes, conformily with the recent works of other authors16). In this case the correlations may have no serious effects. Research in this direction is in progress.
Appendix
A
The matrix R = (Rz) has only real, non-positive, simple eigenvalues and may be diagonalized. Indeed, its eigenvalues are also the eigenvalues of the symmetric matrix S = (Sz)
(A. 1) Furthermore,
the right eigenvectors
(Up) and (u,) of R and S are related by
Up = (x;)“*u,
(A.2)
and the left eigenvectors
(Up) and (up) by
UP =(X3_ ‘/*up = (jQ’/‘zup,
64.3)
so that UP = (x;))‘u, The eigenvalues c
s;u,u,
=
64.4) of R or S are non-positive, -2
0, n
PI4
I
x
7~u,(X~-"*
=C
1
7;u,(x;)-“2
co,
a
the zero eigenvalue corresponding
T
since
2
Up7”, = 0,
to a =
1,. . . , m - 1.
(A.@
P
Thus a left null vector Uop of R is also a left null vector of (7;) which implies, by (A.@, uop~p The corresponding
(A.7) right null vector of R being
u; cc pxo,. But clearly the perturbation
(A.@ (Y,) has no components
on (U3 since by eq. (5) (A.9)
M. MOREAU
158
In the same Indeed,
way the component
this component
AND L. VICENTE
of (G,,)
is given
on (CJ~U~) vanishes
for any
K#
0.
by
GO,", 3: 2 U",U"PWK~U~~'G~~~~.
(A. 10)
P’.4’ By eq. (53.b)
c
U”p’G,,,, = - u”q’&
(A.1 1)
P’
Using matrix
(A.ll) and relation (A.4) between the right R the right-hand side of (A.10) becomes
and
left eigenvectors
of
which vanishes if K f 0, because of the orthogonality of the right and left eigenvectors corresponding to different eigenvalues. Thus the eigenvalues 0 + K = K of matrix R = R @ I + 1 @R do not contribute to the relaxation of (G,,) and the “active” eigenvalue of R is twice the smallest
Appendix System
smallest non-zero
(in absolute value)eigenvalue of R.
B near local equilibrium
In the previous study of a non-homogeneous system, it has been assumed that the globally stationary state was nearly realized, which only permits to consider small deviations from chemical equilibrium and spatial homogeneity. However, the macroscopic chemical relaxation is faster than the diffusion relaxation. At intermediate times we may assume that the chemical equilibrium is approximately realized locally, whereas non-homogeneity subsists. If G = (G,,) represents macroscopic $G
the mean
evolution
equation
density
species
p in every
cell i, the
= R(G) + D(G),
(B. 1)
where R and D represent the chemical of monomer units in cell i is
Gi =
of every
is
and diffusion
operators.
The density
C PGp,i.
Gi
(B.2)
P
Since
the
number
of monomer
units
is conserved
in every
reaction,
Gi
STOCHASTIC
THEORY
OF POLYMERIZATION
159
REACTIONS
satisfies, by (B.l) and (B.2), %
=
C pDp(Gp,i+l+Gp,i-I- 2Gp.i).
(B.3)
P
If the diffusion coefficient is the same for all species, we obtain the discrete equivalent of the diffusion equation: (B.4).
& Gi = D(Gi+l + Gi-l- 2Gi). Let Go(t) be the state of local equilibrium, R( Go) = 0,
(B.9
satisfying the previous diffusion equation and its initial conditions. assuming that the actual state G is near Go, we write
Then,
G = Go+ G’ and linearize eq. (B.l) in the neighborhood
of Go, obtaining
where R” is the linearized chemical operator calculated at Go. The linear equation (B.6) may in principle be solved and yields a consistent result if the inhomogeneous terms remains small. We shall not discuss this equation here, but only consider the evolution of the second moments (Gpi,d) in the lowest approximation, where G is replaced by Go. Then the nonhomogeneous term Kpi,d of eq. (51b) vanishes: indeed the linear “reactions” representing diffusion do not contribute to Kpi,+ and neither do the local chemical reactions, since they are stationary for G = Go. For these reasons the linearized evolution matrix has the same properties as in the global stationary case. It now depends on time through Go, but changes slowly with respect to the chemical characteristic times since the evolution of Go depends on diffusion. Thus the evolution of the second moments is dominated locally by the chemical relaxations and in a first approximation it has the properties studied in section 3.
References 1) J.J. Burton, in Statistical Mechanics, part A, B.J. Berne, ed. (Plenum, 2) R.M. Ziff, J. Stat. Phys. 23 (1980). 241; and references therein. 3) F.C. Goodrich, Proc. Rov. Sot. A 277 (1964) 155. 4) A. Goldbeter and G. Nicolis, Biophysik 8 (1972) 212.
New York,
1977).
160
M. MOREAU
AND
L. VICENTE
5) A. Goldbeter and G. Nicolis. Progr. Theor. Biol.. vol. 4 (Academic Press. New York. 1976) p. 65. 6) A. Babloyantz, Biopolymers 11 (1972) 2349. 7) G. Nicolis and I. Prigogine. Self-Organization in Non-Equilibrium Systems (Wiley-lnterscience, New York, 1977). 8) H. Haken. Synergetics, an introduction, second ed. (Springer, Berlin. 1978). 9) H. Haken. Z. Phys. B20 (1975) 413. 10) C.W. Cardiner, K.J. McNeil, D.F. Walls and 1.S. Matheson, J. Stat. Phys. 14 (1976) 309. 11) H. Richter, I. Procaccia and J. Ross, Adv. Chem. Phys. 43 (1980) 217. 12) P. Glansdorff and I. Prigogine. Structure, Stabilite et Fluctuations (Masson. Paris, 1971). 13) D. McQuarrie. in Suppl. Rev. Her. Appl. Prob. (Methuen. London, 1967). 14) G. Nicolis and 1. Prigogine, Proc. Nat. Acad. Sci. (U.S.A.) 68 (1971) 2102. 15) C. Van den Broeck, J. Houard and M. Malek-Mansour, Physica 1OlA (1980) 167. 16) P. Turq, L. Orcil, J. Barthel and M. Chemla, Berichte der Bunzengesellschaft. Phys. Chem. 85 (1981) 535. 17) H. Lemarchand. Doctorat es-Sciences. Universite de Paris VI (1977). 18) R. Kubo. K. Matsuo and K. Kitahara. J. Stat. Phys. 9 (197.3) 51. 19) S. Dambrine and M. Moreau. Physica 106A (1981) 559-574. 20) N.G. Van Kampen, Physics Letters 59A (1976) 333. 21) W.H. Zurek and W.C. Shieve, J. Stat. Phys. 22 (1980) 289. 22) M. Moreau, unpublished. 23) A. Ikegami, Adv. Chem. Phys. 46 (1981) 363.
Publishing Co.
ON THE STOCHASTIC THEORY OF POLYMERIZATION REACTIONS M. MOREAU and L. VICENTE* Laboratoire
de Physique Thiorique des Liquides**, Unioersit& Pierre et Marie Curie, 4, Place Jussieu, Paris, France
Received 3 July 1981
We present a stochastic study of a polymerization reaction model. Using the WKB approximation and the method of factorial cumulants it is shown that the relaxation time of the second factorial cumulant is half the macroscopic relaxation time. Furthermore the system exhibits a globally non-Poissonian and locally Poissonian behavior in agreement with other results. The application of the model is illustrated for the particular case of dimer and trimer formation.
1. Introduction
The mechanism of polymer formation has been studied for a long tirnele3) and applied to a variety of problems such as nucleation in physics and chemistry’) and prebiotic formation in biologyMSu). In the first case the ideas are applied to a class of transformations that occurs by a growth mechanism, examples of which are boiling and freezing of a liquid and condensation of a vapor. So, in the homogeneous gas phase, nucleation of the system is described by a set of coupled differential equations (of the type of later eq. (7)) which generally are not easy to solve. Here it is necessary to know the thermodynamic properties of clusters, a serious difficulty being the accurate estimation of the free energy barrier’). The analysis of these problems is given by a deterministic causal description, provided by the equations of chemical kinetics. Our purpose here is to study a simple model of polymer chain formation according to the methods of stochastic chemical kinetics developed in recent years ‘-14) and extensively applied to study the fluctuations of the particle numbers and the transition to stable steady states far from equilibrium. We only consider a simple model of polymerization reactions in a closed system: the formation of a p-mer, where a p-mer is a polymer composed of p *Fellow of CONACYT and Facultad de Ciencias, UNAM, Mexico. **Equipe AssociCe au Centre National de la Recherche Scientifique.
037%4371/82/OOOO-OOOO/$O2.75 @ 1982 North-Holland
140
M. MOREAU
monomer
units,
number We
results
of monomer begin
studying
by
from
units
linking
results
which
the
macroscopic
evolution
permits
and is necessary
a monomer
conserved.
the case of an inhomogeneous formalism
L. VICENTE
being
analysing
the linearized
AND
equation system
behaviour
of the
for an homogeneous
is treated
a straightforward for the stochastic
The stochastic study is based first applied to the birth-and-death
to a (p-l)-mer,
the total system
system.
in the so-called
generalization
and Then
multivariate
of the
previous
description.
on the master description
equation formalism, which is of polymerisation in a homo-
geneous system. Using the generating function, the equations for the first two cumulants are derived and the asymptotic behaviour of the correlations between the number of molecules is deduced. But it must be noted15) that the fact of considering the system as homogeneous whole neglects all local effects and thus selects the limited class of big fluctuations at the level of the entire system. Therefore we undertake then the most complete treatment, considering local fluctuations by means of the multivariate master equation’.“) (MME) but making use of the foregoing case. Taking the continuous limit for the correlation we conclude that the system is dominated within each cell by Poissonian-type of fluctuations’).
fluctuations,
2. Macroscopic
in agreement
with the phase-space
description
theory
2.1. Homogeneous
closed system
2.1.1. Generalities It will component components
be
advantageous
for
our
purposes
fluid in which chemical reactions X,
to consider
a general
take place which according to
transform
multithe
m
c V”dXP+, p=l
iqx,,
a=1 ,...,
r.
(1)
The symbol X, is used both for the pth component and the corresponding particle number. The stoichiometric coefficients v;(F;) of X, in reaction CY correspond to reactants (products), and K”, Ka are the rate constants. With the forward and backward rates for the r reactions R,=K”Xf;q..
. X3,
l?,=K”X’;q..
.X3,
CW
STOCHASTIC
THEORY
OF POLYMERIZATION
REACTIONS
141
we write the standard equations of reaction kinetics’*ti)*
y$=z,(i$+$)(R,-R,).
(3)
We consider now the model of a chemical system in which a molecule of kind X1 (monomer) reacts with a polymer X, (chain of size p) producing a polymer X,,, (size p + l), p running from 1 to m always. Allowing also for the reverse process we have thus the scheme: X1+X1 e
x2, R’ KP
x*+x,
e
(4)
xp+,, Iv Km-1
x,+x,-,
e p-l
x,,
and this will be subject to the condition N = 2 pX, = total number of monomer units. p=l
(5)
From eqs. (2) and (4) we see that R,=K”X,X,,
m-l;
a=l,...,
~;=fi~,o+8~,,, F;=ap,a+,,
l?a=l-?‘X~+,, p=l,.
. . , m;
(6)
where S,, is the Kronecker symbol. Substituting this into eq. (3) yields d$=-2(K1X;-~1X3-~;
R,-,) -
$+=(R,-,-
(R,-I?,), (R,
-I?,).
In general a complete set of kinetic equations contains coupling of chemistry and hydrodynamics. Thus the hydrodynamic fluxes must be sufficiently small not to affect the chemistry appreciably. These conditions depend very much *Usually the laws of chemical kinetics then the forward rate is written:
r, = R,/fl = k” so that the constant
are written with the aid of the concentrations
v (xiJuq
0)
K” is related to the usual constant
K”=k”fl’-‘,‘P with analogous
relations
xi=X,/n;
k” by (2c)
for the backward
rate.
142
M. MOREAU
on the spatial return
and temporal
to this point
2.1.2. Stationary
when
scales
we consider
AND L. VICENTE
of the observed the effect
phenomena”).
We shall
of diffusion.
state
The stationary
state of the system
is given
by
!!&-J dt
and from eq. (7) X0,= ap(X$P,
(9)
where al= 1.
KP-'
K’ ap=-‘.K’
. -KP-”
p=2,.
. ,171.
(10)
which, together with eq. (5), defines the stationary detailed balance is then realized, this stationary confirmed later. 2.1.3. Perturbation
state uniquely. state is stable
Since the as will be
theory
In general the exact kinetic equations (7) cannot be solved exactly, excepted in some particular cases’). Thus it is interesting and necessary for our further purposes, to linearize these equations near the stationary state. By putting X, = Xi+ Yp we obtain the linear chemical system (II) with
-I, 7;=i+v;=
if p=as2, if p=o1+122, otherwise.
1, 0,
(12)
and 7:=2;
77=--l,
if ~~32;
CO,=K’X~X:=K’X:~,.
It is shown in appendix A that matrix R=(Rz) simple eigenvalues and may be diagonalized. The explicit solution of the linearized equation Y,(t)
=C c,U;I evf, llf0
has only
real,
non-positive,
is (13)
STOCHASTIC
THEORY
OF POLYMERIZATION
REACTIONS
where the sum is taken over all the negative eigenvalues right eigenvector corresponding to q, and
143
7 of R, (UT) is the
4
4
These results hold for any stationary state with detailed balance. In the case of the polymerization reactions an exact computation of the eigenvalues is generally not possible. But bounds may be found, and analytical or numerical approximations may be used*. Instead, we shall now particularize our system to the case m=3. 2.1.4. Particular case This situation has been studied both experimentally and theoretically by Turq16) from the macroscopic point of view. In many practical cases the forward rate constants are nearly the same for both reactions, and similarly for the backward rate constants, so that we can reasonably take K := K2= K, I?‘=Z?*=K. Introducing the scaled variables
x;=p’,
k'X2
k'X3
x;=j-gx.
(15)
The new variables obey equations similar to (7), with unit rate constants: K;=I?b= 1. Dropping the primes in the variables we obtain the following relation at the steady state: X*=(X9’,
x3= (X7)3.
In the neighborhood by setting Xp=X;+Yp,
(16)
of the steady state the evolution equations are linearized
p=l,2,3.
(17)
The linearized system reads %=c
4
R;Y,,
p=l,2,3,
(184
*It may be shown, for instance, that every vanishing eigenvalue is intermediary between two successive eigenvalues of the chemical system obtained when the concentration of monomers is kept constant”).
144
M. MOREAU
where
the matrix
R=(Rz)
is I
2-x’;
-x:(4+x:)
AND L. VICENTE
(1%) i (R;j) are
r),=o,
_ ?j? ,=S[-s+(&4?Y)“‘] _
(19)
where 6=x’j’+5x0+2
I
3 S’=
9x:‘+4x:‘+
1.
If (CJa)= U” is the eigenvector corresponding to qK, it is shown in appendix that (Y,) has no component on U’; thus the solution of Eq. (18a) is
It is interesting to consider the limit cases the first case X:+1; thus it is found that which implies a relaxation time of order I In the second case 7)?=-9 and v~-X$l, relaxation time of the component on U3 component equilibrium trimerc.
2.2. Di$usion
on U’=(O, 3, -2). value much more
A
of low and high concentrations. In qZ and q7 are nearly equal to -I, for both components on U’ and U’. which implies a very short and a relaxation time = l/9 of the
Thus the monomer concentration reach its rapidly than the concentrations of dimers and
in (1 non -homogeneous
closed
system
2.2.1. Cell description The foregoing analysis is a global one. which does not include the effect of diffusion. We proceed now to consider this effect in the deterministic equations by dividing the system into cells, within which the chemical reaction\ take place and between which matter is exchanged whereas homogeneity is maintained within each cell. The conditions are those of the MME and have been discussed elsewhere’). They imply that the length of each cell is at least of the order of the mean free path of molecules and far smaller than the mean free path between two reactive collisions. We assume these conditions to be realized. The interest of this description will be clear in section 3 when we treat the stochastic approach where the knowledge of the deterministic solution in every cell is needed.
STOCHASTIC
THEORY
OF POLYMERIZATION
145
REACTIONS
Then, we symbolise diffusion by a “reaction”. yp,i
SDP Yp,i',
(21)
where Yp,i is the p-mer number within cell i, and i’ denotes all cells adjacent to i; Dp is the diffusional constant which is related to Fick’s constant d, by the relation’)
4=+,
(22)
where 1 is the length of every cell. For the sake of simplicity we shall restrict ourselves to one-dimensional diffusion. The index i runs from 1 to n, 1 and n corresponding to the boundaries of the system and any interior cell i has two adjacent cells, i - 1 and i + 1. Throughout the analysis we will deal with bounded media and the boundary condition considered is zero fluxes at extremities. To this end we introduce two supplementary cells, 0 and n + 1 and the condition reads Yp.0= Yp.I,
Yp.n+l = Yp,n.
(22’)
2.2.2. Perturbation theory around the stationary state The stationary state of the system obviously corresponds to uniform chemical equilibrium and it is stable for the same reasons as previously. The evolution equation linearized near the stationary state is
where K%i = R;6{ + D$;Ai P.1
I.
R = (R;) is the “chemical matrix”, Si is the Kronecker symbol (also written 6ii according to convenience) and (A{) is the “diffusion matrix” whose elements are (A{) = 6i+l,i + 6i-l.j - 26ij = 6(+1+ s{_, - 26{, (A’,) = (A”,) = - 1. In terms of w, and U”, eigenvalues and eigenvectors of (K$), analogous to the solution of the homogeneous case*: Yp.i
=
C
Yp.i
(24) is exactly
(25)
C, ewKtU”.
I(
*But K now represents cell (i).
a pair of indices (p. i) corresponding
to the chemical
species (p) and to the
M.MOREAU
146
It has been vectors,
shown
chemical
ANDLVICENTE
by Lemarchand”) and diffusive:
that
U is the tensorial
CT= VW. In terms
product
(26)
(U,,) = (V,W,). The eigenvalue
equation
is then
C Kfl:{LJqj = 2 [(RzV,)W,fi{ q.i 4.1 = C RzV,Wi q Dividing
of two
of components
+ D,GzV,(A:.
Wi)]
+ D,V, 2 A:W, = wV,,W,. i
(27)
eq. (27) by U,W, we see that it implies
2 A(W, = AW,. 2 (R; + hD&)V,
(29)
= WV,.
4 By eq. (28) h is an eigenvalue of the diffusion matrix (A’,) and W the corresponding eigenvector. The well-known solutions of this eigenvalue problem, labelled by the index k, running from 0 to n ~ I are: Ak = *(cos+i)
Wt = \/~{expJ~+exp-J(i-,I)k?i]
(30)
where
Then eq. (29) allows us to find the set of eigenvalues {wf} and eigenvectors {Vk}. If the diffusional constants are the same for all species, that is D,, = D. we get from eq. (29) kWT - 77+ AkD.
(31)
where {q} are the eigenvalues identical to the eigenvectors (25) may then be written
given by eq. (19) and the eigenvectors V; of the chemical matrix (Rz). The
explicitly
V;l are solution
as
where
(WY]
c and G, is the concentration coefficient Coo vanishes.
at equilibrium.
As in the homogeneous
(33) case the
STOCHASTIC THEORY OF POLYMERIZATION
147
REACTIONS
2.2.3. Continuous limit If L is the total length of the system, the size of each cell is I = L/n, and to consider the continuous system we make si = iL/N -+ s and
Yp,i must be replaced by Yp(s) and the limn+,, X,/n. Thus, instead of eq. (33) we get
Wk(sW&,0) ds/
local
concentration
is G, =
xp qjoL(Wk(s))2ds,
(34)
where
Wk(s) = cos
F.
From eq. (32) we obtain lim A(k)D= _ !?$Z “+m and therefore,
instead of eq. (32) we have
Yp(s, t) = c VJ g c,,,k e(‘-kznzd’L*)tWk(S)+ z. V;jc,,Oe”‘W”(s). k#O ‘1
(35)
If we have the initial condition Y,(s, 0) = Y i(s) (deviation from equilibrium in the first cell), eq. (35) becomes Y,
=cr) V:[c
(VzY,/G,)/(z
+.X0v:[$ and approximating
Yp(s,t) =
(Va)‘/G,)]
CVYC,)/(~ (v:)~/G,)] e” P
(36)
the sum over k by an integral we obtain
C.Vz[T (V~Y,IG,)/(~ (v~)‘/G,)] :zo
e?’ ~oe-(kZ~ZdiLz)t cos?
P
P
V:[iC
P
(VZ/G,)/(C
P
evt{ Jge-S2/4dt}
(w’/G,)]
The first term is the typical result of diffusion and the to chemical contribution. Certainly this result may tinuous description of the spatial inhomogeneity but will be useful for the next section. A more interesting to assume the chemical equilibrium to be nearly
eq’.
(37)
second one corresponds be obtained by a conthe previous discussion approximation would be realized in every cell,
M. MOREAU
148
whereas
inhomogeneity
the present
subsists
approximation
AND L. VICENTE
between
will be sufficient
cells
(see
Appendix
B). However
for what follows.
3. Stochastic theory 3.1. introduction In this section
we want
to take
into
account
the
random
character
of
reactions, since the individual reaction between molecules is a random event and thus the number of molecule (or polymers) is a random variable, and we want to determine the probability distribution by means of the moment equations. First we consider the process in a global way and assume that reactions are spatially homogeneous and simulated by a stochastic process of birth-and-death type13). In fact this is only a first approximation because it neglects all local effects and thus selects the limited class of big fluctuations at the level of the entire system. The main advantage of such a description lies in its simplicity, which permits analytical approaches. The most complete treatment is made by considering local fluctuations, and this is done by means of the multivariate master equation (MME). In the treatment of the MME the system is divided into cells within which the chemical reactions take place (described by a birth-and-death process as before) and between which matter is exchanged by diffusion (simulated by a random walk). This representation’) supposes the characteristic dimension of cells to be of the order of the mean free path. It is considered that reactive collisions lead to a chemical process whereas elastic ones contribute to diffusion. Furthermore it is supposed that elastic collisions are much more frequent than reactive ones, which implies local equilibrium Finally let us note that the time scale of the diffusion process and the time scale of chemical reactions
are well separated.
3.2. Birth-death
theory
applied
to a homogeneous
system
3.2.1. The generating function method Let X, be the number of the p-mer molecules (p= 1,2,. . . , m) and . ,X,,,: t)= P{X,}) be the probability of having X, monomers, P(X,, . .
let Xz
dimers etc. at time t. We shall study the stochastic behavior of the system by the method of the generating function’,“), using the general notation of subsection 2.1.1 before particularizing to polymerization reactions. In the birth-and-death formalism, the evolution equation of any arbitrary smooth function Q({X,}) is’)
STOCHASTIC THEORY OF POLYMERIZATION REACTIONS
$(QWpIN=&I~
149
[Q({~,+~“,-v”,})-Q({~,})IW~((W,})P({X,}> P O1
+ x ~ u [Q({x,+u”,-V”,})-Q({Xp})lWOI((XP})P({XP}), FX P where the W’s are the transition probabilities
(38)
and are given by (39)
(40) According to the remark leading to eq. (2~) K” is used when the kinetic laws are written in terms of molecule numbers and k” when they are expressed in terms of concentrations. If we take Q({X,})=ll, sp, q=(Q) is the generating function*: (P({sp)r t,=(rJ
(41)
SP)
and we obtain, after manipulating eqs. (38) and (41)
(43) Let us write s=(s,, . . . ) s,) (with analogous notation for {X,}, {v;}, . . . ),
(4) Then eq. (43) may be written v=;
K”(s’=-- s+)&cp + 2 I-&s+ - s”“)&. cp. D
(45)
*cpis related to the average values of various powers of X, by acp ( as, >,r,+, = (xlJ1 with 8X, = x,
-
(X,).
(gk?,,,.,,
= @X,SX,) + W,)(X,) - W&xl,
(42)
150
M.MOREAUANDL.VICENTE
This
linear
cannot
partial
be solved
differential exactly.
so-called WKB approximation by Lemarchand”), writing: cp(s, t)=exp[Ns,
equation
However,
with
variable
for large
coefficients
values
in general
of R, we may
used by many authors’.“)
and developed
use the recently
(46)
t)l.
Substituting this approximation into eq. (45) and retaining dominant terms in 0, one converts this equation into a non-linear differential equation of the first order,
which
is similar
to a Hamilton-Jacobi
and analogous
equation:
for &a.
We may write eq. (47) in the form +)=K(s,
(48)
$“),
where K(s,
&)=c
Switching
a
k”(s”“-s”“)&+~
,$“(s”“-s~“)$QP.
(49)
to the variables (50)
&=ssp-I $ and K in powers
and expanding equations
of the distribution
of &, we obtain
factorial
from eq. (48) the evolution
cumulants:
dG,-K dt@$=T where
(R4,‘G,,,+R4,‘G,,,)+K,,,
p,q=
1,. . . , m,
(51b)
we have used the definitions
K=!E P
with
(Sla)
”
(86,>&=O’
t&g,,,
&q=(&,,, )=o) RZ= p q B
p
q
)=O’
P
(52)
STOCHASTIC THEORY OF POLYMERIZATION
151
REACTIONS
The factorial cumulants of order one and two are
G=?!k P
(atp)
t5p)=o
_
(534 (53b)
,=,=h [Wp~&>-Wp>4ql. Gpq=(&),, p 4 P
Thus G, is just the concentration and G,, is related directly to the deviation of the distribution from the Poissonian form’). From eq. (49) we obtain
Kp=C
(Wa-Ga)T;,
a
K,,=C (cd”-
W”)( q+“g- v;v”s-Ta~pq),
LI
R;=x
T;(w=v:-G%;)/G~,
(54)
OL
where T; is given in (12), and wa=ku n(GP)“B, P
O”=l;”
n(Gp)'?
(55)
P
Since G, = X,/n (eq. 53a), wa and W” are identical to the deterministic forward and backward rates r, and fa and (51a) indeed is the deterministic equation. 3.2.2. Evolution of the second cumulant If we assume that the macroscopic concentrations their equilibrium stationary values Gp, we have
have already reached
(56)
01a = w”(&‘w”((;‘),
&{Gi}.
Then the evolution equations (51b) are replaced by the homogeneous
(57) system
where &;q=ti;G;+&i3;. The properties of the asymptotic matrix a have been studied previously the forma1 solution of eq. (58) is G(t)=eR”‘G(0),
(59) and (60)
M. MOREAU
152
where e
W = (Z?yq). Using
the tensorial
AND L. VICENTE
product
we may write
RI = efif @ ,It.
(61)
If (7) are the distinct the eigenspace
eigenvalues
of 7 it is known
of k and P, is the projection
operator
on
that
(62) Then
the solution G,,,(t)=2
of eq. (58) may be written
(63)
c,, 1,e(“+““UaUi,
c~,=[z‘;i)O:G,.,(O),G~G~,,~
(U;U:)‘lG,G,]
P4
where
U”.”
as
=
~‘,
(64)
Pq
(UziJ
=
are the corresponding
are the eigenvectors eigenvalues”. Ix).
(UZU:)
of matrix
R and A,, ,. = rl+u
3.2.3. Stationary state If the evolution matrix 6 is regular at the stationary state, all its eigenvalues are strictly negative, and G,,(t) tends to 0 as t + m. Then the distribution becomes Poissonian and uncorrelated in agreement with the thermodynamical equilibrium of a grand ensemble. However, if some quantity is conserved in all reactions (and this is indeed the case for the number of polymer units in our model), i is singular and GP4(t) tends G&t)+
to CO,OU~U~ according Co,opq~&
to (63). Here
(t +m).
(65a)
with
(65b) which generally does not vanish. Thus the stationary state is generally non-Poissonian, unless of course the initial state is Poissonian. It may be observed that the variance of the total number of monomer units, N=C,pX,, is
W’J)2)=~q Pq@WXJ=fJ~
P4
pq(G,,+G,6,,),
(66)
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
153
so that by (65)
co.o=-cd’Gp+~(;~~p)l P
P
The asymptotic
distribution
((W’)=fl
c p2G,=C P
is only Poissonian
(67)
if
P’(&>.
(68)
P
In particular if N is known exactly, ((SN)2) vanishes and
Co,o=-(7 p2Gp)-'CO. Such deviations from Poissonian equilibrium have been predicted by Van Kampen20) and studied recently by numerical computation2’). They are clearly due to the incompatibility of the initial condition and the Poissonian distribution, which cannot disappear because of the conservation relation. 3.2.4. Asymptotic evolution On the other hand, if the macroscopic evolution matrix k of a chemical system is regular at the stationary state, the smallest (in absolute value) eigenvalue of the evolution matrix R for the second moments is twice the smallest eigenvalue of R, so that the relaxation time of (GP4) is half the macroscopic relaxation time, in agreement with other results’?. If k is singular, which is the case here and more generally when there are conservation relations the smallest non-vanishing eigenvalues of k and W are identical; however it is shown in appendix A that the previous result still holds. Thus the relaxation of the second moments is faster than the macroscopic relaxation to the stationary concentrations. However it should be remembered that the system is supposed here to be in macroscopic equilibrium. Otherwise the complete non-homogeneous eq. (51b) replaces eq. (58) and the evolution of (G,,) exhibits terms induced by the macroscopic evolution. 3.3. Multivariate
formalism
of a non-homogeneous
system
3.3.1. E#ect of diffusion Without writing the MME we can write immediately the equations for the first two moments, adding the effect of diffusion (random walk) and remembering that chemical reactions are described be the same birth-and-death process as before, but considering now the equations for every cell. We have
M. MOREAU
154
AND L. VICENTE
~=D,(G,i+i+G,,i~,-2G,i)+~~,i where
p denotes
p=l,m:i=l,..., and i the cell number.
the species
again the macroscopic 23). In the multivariate
solution
which
formalism,
p,p’,q
(70)
This equation
we have examined
in section
represents 2. (see eq.
eq. (58) becomes
. . ,111, i,j=l,...,
=I,.
n,
(71)
n,
where p, q and p’ denote species, i and j cells; RP,’ is again the matrix of linearized macroscopic equations. In the case of our system (cf. eq. (4)) and for equal
coefficients
GP,i,q.j(f) = 3
of diffusion
C,,kGv.h e
the solution
(v+u-(hk+Ah)D)f
v;l
of eq. (71) is
w;yy;, (72)
,a~h
with c7.k;v.h
=
z
Es
p,q,i,j
where
c
W;WfGp~i;q,j(0)
Gp
/[
Gq
q and v are eigenvalues
of the chemical
2
-&Wk)2
G, 1 [ i
p
matrix
’
’ 1 ’
and Ak, Ak eigenvalues
(73) of
diffusion; V;l, Vi are the corresponding chemical eigenvectors and W!, W: the corresponding diffusion eigenvectors. From eq. (72) we see that Gp,i.q.i state given by evolves toward a stationary
(74)
Generally this expression sonian. But if we consider treatment gpq(h
in section s2;
t)
=
z
2.2.3) we obtain v;v;
2
L
c p.9
does not vanish unless the initial condition is Poisthe limit of the continuous system (in analogy to the
v; v; G G
e
(t)+~~(k2+h*)nZdlL2)1~k(Si)~h(S2),
(75)
L
II 0
%k;+h
wk(sd
w”(sdgpqh,
~2; 0)
dsl
dsz
0
(76)
&k:v,h= [T
~]2[/oL(Wk(sd2
&I2
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
155
where g&l,
~2;O)=lim GP~i;~(“)=(SXp(S,)Sxq(S*))-(Xp)S~qS(S2-S,) n-o
(77)
and Sx, =x,-(x,). Thus g,, tends to the uniform value L &&It
s2;
aJ)
=
9
c P4
P4
II
L
gpqh,
~2;
0)
dsi
ds2.
(78) 2
As for homogeneous systems, the right-hand side of eq. (71) may be related to the variance ((SN)*) of the total number N of monomer units contained in the system, since L
((6N2)
=
c
pq
L
1
/
0
0
(S(x,(s,)Sxq(sz)) ds, dsz.
(79)
P
Thus, eq. (77) and eq. (78) yield
If N is exactly known, ((SN)2) vanishes, and the asymptotic value of g,, is given by the last term of eq. (80) which shows that the concentrations of the different chemical species are correlated and non-Poissonian. However for large systems (L + w), this term tends to 0, and the fluctuations become locally Poissonian, although the system is globally non-Poissonian, since LL
II 0
g&l,
~2;
~1
dsl
ds2
=
@~p~~q)-
Wp)&q
0
= _
PdWW,)
T
P~WP>
(81)
.
More generally if x,(sJ and xq(s2) are S-correlated, or at least if their correlation has a very short range, it is seen from eq. (77) and eq. (78) that g,,(sr, s2; 00)is of order I_-‘, and thus vanishes for large systems. Equivalently, this is true if the variance ((SN)2) of the total number of monomer units of the system is of order I., which corresponds to ordinary conditions.
M. MOREAU
156
AND L. VICENTE
Then the behaviour of the system is locally Poissonian, and globally non-Poissonian, in agreement with other results obtained analytically’5~2”) or by numerical This
simulation*‘).
conclusion
correlation
does
between
((6N)*) is of order a non-Poissonian
not hold
some
x,(s,)
L’; th en g,,(s,,
near
an instability,
and
xq(sZ) is of order
s2, m) remains
where
the range
of the
L, or equivalently,
finite, and the system
exhibits
local behaviour”).
4. Discussion The major difficulty in the study of chemical correlations appears generally in the approximations used to truncate systems of equations involving arbitrarily high moments. Indeed, as was pointed out by Nicolis’) the higher-order moments become increasingly important in the absence of asymptotically stable microscopic states. However this is not the case for the model used in this paper. Thus the use of a WKB approximation is justified according to the discussion of Kubo et al.‘*) in the limit where the volume 0 is considered to be macroscopic. The present study of a chain of polymerization reactions may be generalized in various directions. We have determined the behaviour of G,, supposing that the macroscopic concentrations have already reached their equilibrium stationary state. It is possible to apply the method to systems outside macroscopic equilibrium, and even to systems far from macroscopic equilibrium if the macroscopic evolution is known. Here it may be done -at least formally - near macroscopic equilibrium. In this case the inhomogeneous terms K,, of eq. (51.b) for the second cumulants G,, does not vanish; furthermore the matrix (Rz) is no longer constant. Then G,, indeed tends to the stationary value computed in section 3, but the inhomogeneous term K,, induces a coupling between the proper relaxation of the second cumulants and the macroscopic
relaxation,
as will be shown
more
generally
elsewhere.
Another possibility of extending the present method is outlined in appendix B: the system is supposed to be near local chemical equilibrium, but far from the stationary state of diffusion, that is, far from homogeneity. In this case the system is easily linearized since the diffusion-terms remain linear in the domain of validity of the multivariate master equation, and the previous results are partly generalized. It may be noticed that the correlations should be taken into account in the global equations of evolution for the average values, although they disappear in the local equations for concentrations. It would be interesting to introduce the effect of an electric field in our
STOCHASTIC THEORY OF POLYMERIZATION
REACTIONS
157
model, namely an internal field due to charged particles, in order to study chemical reactions in electrolytes, conformily with the recent works of other authors16). In this case the correlations may have no serious effects. Research in this direction is in progress.
Appendix
A
The matrix R = (Rz) has only real, non-positive, simple eigenvalues and may be diagonalized. Indeed, its eigenvalues are also the eigenvalues of the symmetric matrix S = (Sz)
(A. 1) Furthermore,
the right eigenvectors
(Up) and (u,) of R and S are related by
Up = (x;)“*u,
(A.2)
and the left eigenvectors
(Up) and (up) by
UP =(X3_ ‘/*up = (jQ’/‘zup,
64.3)
so that UP = (x;))‘u, The eigenvalues c
s;u,u,
=
64.4) of R or S are non-positive, -2
0, n
PI4
I
x
7~u,(X~-"*
=C
1
7;u,(x;)-“2
co,
a
the zero eigenvalue corresponding
T
since
2
Up7”, = 0,
to a =
1,. . . , m - 1.
(A.@
P
Thus a left null vector Uop of R is also a left null vector of (7;) which implies, by (A.@, uop~p The corresponding
(A.7) right null vector of R being
u; cc pxo,. But clearly the perturbation
(A.@ (Y,) has no components
on (U3 since by eq. (5) (A.9)
M. MOREAU
158
In the same Indeed,
way the component
this component
AND L. VICENTE
of (G,,)
is given
on (CJ~U~) vanishes
for any
K#
0.
by
GO,", 3: 2 U",U"PWK~U~~'G~~~~.
(A. 10)
P’.4’ By eq. (53.b)
c
U”p’G,,,, = - u”q’&
(A.1 1)
P’
Using matrix
(A.ll) and relation (A.4) between the right R the right-hand side of (A.10) becomes
and
left eigenvectors
of
which vanishes if K f 0, because of the orthogonality of the right and left eigenvectors corresponding to different eigenvalues. Thus the eigenvalues 0 + K = K of matrix R = R @ I + 1 @R do not contribute to the relaxation of (G,,) and the “active” eigenvalue of R is twice the smallest
Appendix System
smallest non-zero
(in absolute value)eigenvalue of R.
B near local equilibrium
In the previous study of a non-homogeneous system, it has been assumed that the globally stationary state was nearly realized, which only permits to consider small deviations from chemical equilibrium and spatial homogeneity. However, the macroscopic chemical relaxation is faster than the diffusion relaxation. At intermediate times we may assume that the chemical equilibrium is approximately realized locally, whereas non-homogeneity subsists. If G = (G,,) represents macroscopic $G
the mean
evolution
equation
density
species
p in every
cell i, the
= R(G) + D(G),
(B. 1)
where R and D represent the chemical of monomer units in cell i is
Gi =
of every
is
and diffusion
operators.
The density
C PGp,i.
Gi
(B.2)
P
Since
the
number
of monomer
units
is conserved
in every
reaction,
Gi
STOCHASTIC
THEORY
OF POLYMERIZATION
159
REACTIONS
satisfies, by (B.l) and (B.2), %
=
C pDp(Gp,i+l+Gp,i-I- 2Gp.i).
(B.3)
P
If the diffusion coefficient is the same for all species, we obtain the discrete equivalent of the diffusion equation: (B.4).
& Gi = D(Gi+l + Gi-l- 2Gi). Let Go(t) be the state of local equilibrium, R( Go) = 0,
(B.9
satisfying the previous diffusion equation and its initial conditions. assuming that the actual state G is near Go, we write
Then,
G = Go+ G’ and linearize eq. (B.l) in the neighborhood
of Go, obtaining
where R” is the linearized chemical operator calculated at Go. The linear equation (B.6) may in principle be solved and yields a consistent result if the inhomogeneous terms remains small. We shall not discuss this equation here, but only consider the evolution of the second moments (Gpi,d) in the lowest approximation, where G is replaced by Go. Then the nonhomogeneous term Kpi,d of eq. (51b) vanishes: indeed the linear “reactions” representing diffusion do not contribute to Kpi,+ and neither do the local chemical reactions, since they are stationary for G = Go. For these reasons the linearized evolution matrix has the same properties as in the global stationary case. It now depends on time through Go, but changes slowly with respect to the chemical characteristic times since the evolution of Go depends on diffusion. Thus the evolution of the second moments is dominated locally by the chemical relaxations and in a first approximation it has the properties studied in section 3.
References 1) J.J. Burton, in Statistical Mechanics, part A, B.J. Berne, ed. (Plenum, 2) R.M. Ziff, J. Stat. Phys. 23 (1980). 241; and references therein. 3) F.C. Goodrich, Proc. Rov. Sot. A 277 (1964) 155. 4) A. Goldbeter and G. Nicolis, Biophysik 8 (1972) 212.
New York,
1977).
160
M. MOREAU
AND
L. VICENTE
5) A. Goldbeter and G. Nicolis. Progr. Theor. Biol.. vol. 4 (Academic Press. New York. 1976) p. 65. 6) A. Babloyantz, Biopolymers 11 (1972) 2349. 7) G. Nicolis and I. Prigogine. Self-Organization in Non-Equilibrium Systems (Wiley-lnterscience, New York, 1977). 8) H. Haken. Synergetics, an introduction, second ed. (Springer, Berlin. 1978). 9) H. Haken. Z. Phys. B20 (1975) 413. 10) C.W. Cardiner, K.J. McNeil, D.F. Walls and 1.S. Matheson, J. Stat. Phys. 14 (1976) 309. 11) H. Richter, I. Procaccia and J. Ross, Adv. Chem. Phys. 43 (1980) 217. 12) P. Glansdorff and I. Prigogine. Structure, Stabilite et Fluctuations (Masson. Paris, 1971). 13) D. McQuarrie. in Suppl. Rev. Her. Appl. Prob. (Methuen. London, 1967). 14) G. Nicolis and 1. Prigogine, Proc. Nat. Acad. Sci. (U.S.A.) 68 (1971) 2102. 15) C. Van den Broeck, J. Houard and M. Malek-Mansour, Physica 1OlA (1980) 167. 16) P. Turq, L. Orcil, J. Barthel and M. Chemla, Berichte der Bunzengesellschaft. Phys. Chem. 85 (1981) 535. 17) H. Lemarchand. Doctorat es-Sciences. Universite de Paris VI (1977). 18) R. Kubo. K. Matsuo and K. Kitahara. J. Stat. Phys. 9 (197.3) 51. 19) S. Dambrine and M. Moreau. Physica 106A (1981) 559-574. 20) N.G. Van Kampen, Physics Letters 59A (1976) 333. 21) W.H. Zurek and W.C. Shieve, J. Stat. Phys. 22 (1980) 289. 22) M. Moreau, unpublished. 23) A. Ikegami, Adv. Chem. Phys. 46 (1981) 363.