552
A. K. Thakur, A. Rescigno, and C. DeLisi
The Journal of Physical Chemistry, Vol. 82,No. 5, 1978
(3) H. Farber and S. Petrucci, J . Phys. Chem., 79, 1221 (1975). (4) W. E. Vaughan In "Dielectric Propertiesand Molecular Behaviour", Van Nostrand, Reinhold, New York, N.Y., 1969, p 146. (5) A. R. Davis, D. E. Irish, R. B. Roden, and A. J. Weerheim, Appl. Spectrosc., 28, 384 (1972). (6) D. E. Irish, S. Y. Tang, H. Talts, and S. Petrucci, to be published. (7) R. M. Fuoss, Chem. Rev., 17, 27 (1935); J. Am. Chem. Soc., 57, 2604 (1935). ( 8 ) C. J. F. Bottcher, "Theory of Electrical Polarization", Elsevier,
Amsterdam, 1952. (9) E. A. S. Cavell, P. C. Knight, and M. A. Sheik, Trans. Faraday Soc., 67, 2225 (1971). (10) M. Davies in "Dielectric Propertiesand Molecular Behavlour", Van Nostrand, Reinhold, 1969, p 298. (11) P. Debye, "Polar Molecules", Chemical Catalog Co., New York, N.Y., 1929. (12) F. Perrin, J. Phys. Radium Paris, 5 , 497 (1934); E. Fischer, Physica, 7, 60, 645 (1939).
Stochastic Theory of Second-Order Chemical Reactionst Ajlt K. Thakur," Aldo Resclgno, and Charles DeLisl Laboratory of Theoretical Biology, National Cancer Institute, National Institutes of Health, Bethesda, Maryland 200 14 (Received April 25, 1977; Revised Manuscript Received September 21, 1977) Publication costs assisted by the Natlonai Cancer Institute
A stochastic theory of the reversible bimolecular reaction 2A C is formulated and its thermodynamic and kinetic properties investigated in detail. The distribution function was obtained exactly by an analytical method with the system in equilibrium, and by numerical solution of the Kolmogorov equation for the more general time-dependent situation. The solutions provide a standard against which various approximation procedures of general interest are evaluated; in particular, calculation of the moments and inference of the distribution by an expansion over an appropriate orthonormal basis set. The determinants of an efficient approximation procedure in terms of constraints on the moments are delineated in a general way and then illustrated by application to other commonly studied bimolecular systems.
I. Introduction The stochastic theory of compartmental analysis in general and first-order chemical reactions in particular is well developed and has made possible a clear understanding of a number of simple chemical reactions both near and far from eq~ilibrium.l-~ The stochastic theory of second- and higher-order chemical reactions, on the other hand, is considerably less developed and consequently such processes are not as fully understood. Exact expressions for the stochastic mean and the generating function at equilibrium have been found for a number of reversible reaction ~ c h e m e s ~and , ~ *expressions ~ for time evolution of the stochastic mean and variance as well as the entire distribution have been obtained for some irreversible r e a ~ t i o n s . ~ J -Darvey l~ and Ninham'l have investigated some aspects of the reversible kinetic problem A B F? C D based on different constraints and firstand second-order perturbation methods. Besides that, little seems to have been done in the way of a complete analysis of both the equilibrium and kinetic stochastic properties of even a relatively simple reversible secondorder chemical reaction. The purpose of this paper is twofold. First, we present a detailed investigation of the reaction
+
+
2A* C
as a function of time, system size, and kinetic parameters. We obtain the exact distribution function at equilibrium by a combinatorial method which for this particular problem appears to be simpler and more direct than a +Partsof this work were presented at the 21st Annual Meeting of the Biophysical Society held at New Orleans, La., 1977 and appeared as an abstract. * Present Address: Endocrinology and Reproduction Research Branch, National Institute of Child Health and Human Development, NIH, Bethesda, Md. 20014. This article not subject to US.Copyright.
generating function approach. Far from equilibrium we have solved the master equation numerically. The other aspect of the problem is the use of these exact solutions as a standard against which the utility of approximate computational procedures of general interest are evaluated in particular, the calculation of moments and the inference of the distribution from an expansion over a Gaussian basis set.
11. Formulation of t h e Equations Let N 2C + A be the total number of units in the system. We take N to be time independent so that it is necessary to keep track of only one type which we choose as A. Let P,(t) be the probability that the system is in state i at time t ; i.e., Pi(t)is the probability of finding i units of type A in the system a t time t. Then if aAt and @At are, respectively, the transition probabilities for the formation and dissociation of C in the small time interval At
P i ( t + A t ) = a('
+
') Pi+ ,(t)A t (/3/2)(N- i 2)Pi.z(t)At + [l - a ( $ ) A t - ( / 3 / 2 ) ( N i ) A t ] P j ( t )+ O ( A t )
...
+ (1)
The first term on the right is the probability that a system in state i 2 will undergo a transition to a system in state i as the result of a pairwise interaction between two A's to form a C; the second the probability of a transition from state i - 2 to i as the result of dissociation of a C; and the third the probability that a system in state i remains in that state during At. Other terms, which are not written explicitly, represent multiple associations and/or dissociations and are proportional to terms in At higher than the first order. Rearranging eq 1 and letting At go to zero gives the forward Kolmogorov differential-difference equation for the system:
+
Published 1978 by the American Chemical Society
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978 553
Stochastic Theory of Second-Order Chemical Reactions
+
dPi(t)/dt = ( a / 2 ) ( i l ) ( i + 2)Pi+ z ( t ) + (0/2)(N - i + 2)P,.,(t) - [ (a/2)i(i - 1) (P/2)(N -
+
ill Pi(t)
1.5
t
(2)
with Pi(0) = 8 i ~ where , 6 is the Kronecker delta. Introducing the probability generating function G(s, t ) =
0
I=
Pi(t)si
where Is1 5 1, eq 2 becomes
aG/a t = ( a / 2 ) ( 1 - s2)(a2G/as2)+ ~ ( 0 / 2 ) ( 1s')(aG/as) ? ( 0 / 2 ) N ( 1 - s2)G
which is to be solved subject to the boundary condition that
The generating function for the equilibrium problem, i.e., the problem with aG/at = 0, has been found previously,l0 and from it an expression for the stochastic mean was obtained. In theory, the generating function uniquely defines the distribution; in practice, however, the inversion required to obtain it is often a formidable task. On the other hand, the distribution can be found exactly and directly by a statistical thermodynamic approach. We suppose that there is an ensemble of systems, each having exactly N units, and define the free energy zero of a system as the state i = 0; Le., the state in which all units are paired. Then a system with i = 2n unpaired units will have a statistical weight Knrelative to the reference state of zero free energy where = p/a =
exp(- AGIRT)
(4)
AG being the free energy change accompanying bond cleavage. The degeneracy factor Q ( N ,n) associated with this state is the ratio of the number of ways n C type molecules can dissociate to the number of ways in which the resulting 2n A type molecules can reassociate. Then the partition function X is of the form [NlZl
n=1+
x
n= 1
(5 1
n(N,n)Kn
and the probability that a randomly picked system will be in state N - 2n is Pzn
(Go)
= s-2 (N, n ) K n / n
(6)
The combinatorial term can be found as follows. Since there are N / 2 pairs in the reference state, there will be N / 2 ways in which one pair can dissociate, ( N / 2 ) ( N / 2- 1)ways in which two can dissociate, and ( N / 2 ) ( N / 2- 1)...[N / 2 (n- l)]ways in which n can dissociate. For reassociation in a system with 2n free units the first pair can be formed in ways, the second in ways, etc. Therefore
a ( N ,n )=
( N / 2 ) ( N / 2 - 1)...[N / 2 - (n - l)]
n!(2n- l)!!
(7)
The distribution function at equilibrium is therefore Pzn(-) =
,
I
103
200
300
I
4w
K
Figure 1. The difference between the stochastic mean (SEQ) and the deterministic value (DEQ) of the number of particles at equilibrium as a function of dissociation constant (K)for different values of Nfor the C. reaction 2A
*
G(s, 0) = s N
K
01 10
(3)
(N/2)(N/2 - 1)...[N / 2 - IZ n!(2n - l)!!Z
+ 1]P
(8)
and the mth moment is
(9)
Writing eq 9 with m = 1 in terms of confluent hypergeometric series we have M i = N K i F i ( - N / 2 + 1;3 / 2 ; -K/2)/1Fi(-N/2; 1/2;-K/2) (10) which is the same as the expression obtained by Darvey et al.1° using the probability generating function technique.
111. Moment Approximations for the Transient Reaction The equations for the moments are obtained by multiplying both sides of eq 2 by "i and summing over i. Thus dMm/dt = ( a / 2 ) [ Z i ( i - l ) ( i - 2)mPi(t) i
m + y i - l)Pi(t)] i
+ (@/2)[Z.(i+ 2 ) m ( N i
i)Pi(t) - Z : P ( N - i)Pi(t)] i
m = 1,2,3, ... (11)
Explicit equations for the first six moments are given in Appendix A. It is evident that although the system of differential equations is linear, the equation for the mth moment contains the (m list as a variable, so that some approximate closure condition must be introduced. Two different truncation methods were used. One was obtained by setting
+
Mm = C(N, IZ, K)MlM,-, in the ( m - l)stequation where C is to be evaluated by comparison. The other involves writing eq 11 in terms of central moments p,,, defined as p m = ( ( N - 2n - Ml)m)
(12) In this case with m even, the system is truncated by setting the succeeding odd central moment equal to zero (Appendix A).
IV. Results According to the results shown in Figure 1, for the equilibrium case, the stochastic mean number of free molecules always exceeds the deterministic mean with the difference rising to a maximum as K increases and then decreasing toward zero. This behavior can be understood by recognizing that for K infinite (no interaction) and K zero (infinite strength interaction) the deterministic and stochastic means must be the same when the reaction has gone to completion (for N even); in particular, there will be no free molecules for the infinite strength interaction and no bound molecules for zero strength interaction. In between fluctuations favor dissociation so the stochastic
554
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978
A. K. Thakur, A. Rescigno, and C. DeLisi
TABLE I: Comparison of Expectations with 40 U n i t s I n i t i a l l y in State A (Fkee) a n d K = 20
Time, s
Stochastic exact
Deterministic exact
0 0.2 0.4 0.8 1 1.2 1.6 2
40 29.701 24.912 21.426 20.812 20.498 20.223 20.144
40 29.537 24.883 21.391 20.755 20.412 20.124 20.037
Truncation
Truncation
Truncation
I
I1
I11
Truncation
IV
Truncation
M, = M,M,
M , = M,M,
M , = M,M,
Cl,=o
Clr=
40 29.719 25.140 21.789 21.183 20.844 20.531 20.425
40 29.703 25.014 21.307 20.596 20.216 19.917 19.842
40 29.783 25.099 21.773 21.191 20.856 20.516 20.378
40 29.708 25.065 21.543 20.895 20.543 20.243 20.151
40 29.627 24.774 21.058 20.405 20.069 19.812 19.749
V
0
0.12-
0.10 -
0.08 -
c d 5
P
0.06-
0.WC I
1
80
160
I
I
240
320
I
A
I 4w
K
Figure 2. The equilibrium variance as a function of dissociation constant for different values of N for the same reaction.
NUMBER OF UNITS IN THE FREE STATE (AI
Figure 4. Distribution functions for (from left to right) equilibrium ( t = a)and t = 0.2 s for N = 40, a = 0.05,and = 1 for the same reaction. See Figure 3 for explanations of symbols. Initial conditions are set as in Figure 3: p, = 0.004, p2 = 2.93 (equilibrium) and PI = 0.02, p2 = 2.92 ( t = 0.2 s).
NUMBER OF UNITS IN THE FREE STATE IAI
Flgure 3. Distribution functions for (from left to right) equilibrium ( t = m), t = 1 s, and t = 0.1 s for N = 20, a = 0.05, and p = 1. The solid lines are results obtained by exact numerical solution of eq 8 for the equilibrium case. The triangles (A)are points generated by reconstructing the distribution using a Gaussian basis using four central moments and the solid circles ( 0 )are points generated by assuming a Gaussian distribution. The initial conditions are set such that all units were in state A at t = 0 for the above reaction: (3, = 0.002, p2 2.85 (equilibrium); 0, = 0.005, p2 = 2.84 ( t = 1 s);and (3, = 0.7, p2 = 3.34 ( t = 0.1 s).
mean must increase relative to the deterministic value as K increases from zero, and then decrease as K continues to increase. Figure 1also shows that, as one would expect,
the magnitude of the difference decreases and its dependence on K becomes less pronounced, as N increases. The physical basis of these behavioral patterns suggests that the qualitative features of the results shown in Figure 1 may have general validity for reversible bimolecular processes. The equilibrium variance, as demonstrated in Figure 2, reaches a maximum for a certain value of K (depending on N) and slowly goes to a stable minimum for very high values of K , indicating the dispersion behavior of the system. Of course the dispersion is always larger for larger systems as is evident from the figure. The approach of the distribution toward its equilibrium form is shown in Figures 3 and 4. The details of this approach will depend on the initial condition of the system (in this case all molecules are of t y p e A at t = 0), although the equilibrium distribution should not. It is evident that for N = 20 the equilibrium distribution approaches a Gaussian whereas the nonequilibrium distribution does not, This observation has a bearing on the extent to which the distribution can be accurately inferred from a given number of moments. There are two aspects to the accuracy of the inference: the accuracy of the moments if an approximate method is being used to obtain them and the rate at which the expansion of the distribution over an orthonormal basis set (in this case Tchebycheff-Hermite polynomials, Appendix B) converges. With regard to the former, there appears to be little difference in the results obtained by
Stochastic Theory of Second-Order Chemical Reactions
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978 555
TABLE 11: Comparison of Expectations with 100 Units Initially in State A (Free) and K = 20 Truncation Truncation Truncation Truncation Stochastic DeterminI I1 I11 IV Time, s exact istic exact M, = M,M, M , = MlM4 M, "M,M, P'J= 0 100 100 100 100 100 0 100 54.086 54.266 54.268 53.923 53.875 0.2 54.058 42.439 40.975 40.084 41.488 42.272 0.4 41.043 36.901 36.065 36.681 35.273 36.798 0.8 36.417 36.305 36.613 36.939 34.810 36.212 1 36.074 36.068 36.114 36.618 34.718 35.98 1.2 35.964 35.935 36.110 36.475 34.814 35.85 1.6 35.917 35.914 36.110 36.463 34.916 35.83 2 35.912 TABLE 111: Comparison of Second Moments with 40 Units Initially in State A (Free) and K = 20 Truncation Truncation Truncation Truncation Stochastic I I1 I11 IV Time, s exact M, = M,M, M , = MlM4 M, = MIMs w 3 = 0 1600 1600 1600 1600 0 1600 894.30 893.44 896.13 0.2 893.87 890.43 634.43 641.22 649.74 631.17 0.4 633.59 477.33 464.03 490.15 466.83 0.8 472.28 449.90 440.33 460.00 441.87 1 446.41 435.36 429.64 442.92 429.35 1.2 433.46 423.18 420.72 418.70 428.62 1.6 422.34 419.49 416.92 424.62 414.8 2 419.17 TABLE IV: Comparison of Second Moments with 100 Units Initially in State A (Free) and K = 20 Truncation Truncation Truncation Truncation IV I1 I11 I Time, s Stochastic exact P3=0 M, =MIM, M4= M,M, M , = M,M4 104 104 104 104 104 0 2.9565 X l o 3 2.8848 x l o 3 2.9109 x l o 3 3.0115 x l o 3 0.2 2.9525 X l o 3 1.8300 x l o 3 1.6402 X l o 3 1.7799 x l o 3 1.9352 x l o 3 0.4 1.713 X i o 3 1.3899 x lo3 1.4246 x l o 3 1.4122 x l o 3 0.8 1.3542 x l o 3 1.3897 x l o 3 1.3462 X l o 3 1.3901 x l o 3 1.3593 x l o 3 1 1.3293 x l o 3 1.3466 x l o 3 1.3291 X l o 3 1.3511 X l o 3 1.3233 x l o 3 1.3411 x l o 3 1.2 1.3214 X l o 3 1.3195 X l o 3 1.3074 X l o 3 1.3082 x l o 3 1.3323 x l o 3 1.6 1.3181 x l o 3 1.3180 X l o 3 1.3042 X l o 3 1.33 X l o 3 2 1.3177 x l o 3 1.3077 X l o 3
the two truncation methods used here (section 111); truncation at different levels giving moments with high degree of accuracy in both cases (Tables I-IV). Higher moments can be obtained with comparable accuracy by retaining a correspondingly larger number of equations in the system. The closure approximation was also studied by looking for a coefficient C(N, n, K ) such that M,, = C(N, n, K ) . M,&f1, with the moments obtained from eq 8 when the system is at equilibrium, and by numerical solution of the master equation with the system away from equilibrium. The results (Tables V and VI) indicate that C is always close to one, and is relatively insensitive to time, equilibrium constant, order of truncation, and system size. The question which remains is how quickly the expansion of the distribution function converges. The rate of convergence is clearly going to depend on the set of orthonormal polynomials chosen as the basis set. Evidently, with the Gaussian basis used here, the procedure works well, even when there are as few as 20 molecules in the system, when the first four moments are used. This is perhaps not surprising for the equilibrium distribution since even a t N = 20, it is similar to a Gaussian. What is surprising is the accuracy of the expansion at small times when the distribution is highly asymmetric (Figures 3 and 4). This is in marked contrast to results obtained in studies of polymer chain distance distribution functions which indicate poor convergence for distributions that are not similar to a Gaussian.16 In addition, the convergence here appears to be asymptotic, with the results becoming more accurate as the number of moments increases to four and becoming worse thereafter. This is not
Truncation
V Ps=O
100 53.532 41.473 36.048 35.551 35.382 35.309 35.303
Truncation V PIS= 0 1600 914.91 665.96 489.1 458.39 442.3 429.5 426.16
Truncation V Ps= 0 104 3.041 X l o 3 1.863 X l o 3 1.395 x l o 3 1.352 X lo' 1.337 X IO3 1.33 X l o 3 1.329 X l o 3
TABLE V: Values of the C's at Different Times and also at Equilibrium with 20 Units Initially in State A (Free) for Two Different Values of K K = 20 K = 200
c,
Cl
0.2 0.6 1.2
1.03 1.06 1.13 1.17
1.03 1.05
Equilibrium
1.18 1.15
t, s 0.1
c 3 Cl Cl c 3 1.02 1.02 1.01 1.01
1.04
1.03
1.02
1.02
1.10 1.08 1.03 1.03 1.02 1.13 1.11 1.03 1.03 1.02 1.11 1.03 1.03 1.02
TABLE VI: Values of the C's at Different Times also at Equilibrium with 100 Units Initially in State A (Free) for Two Different Values of K
K = 200
K = 20
Cl
c,
0.1 0.2 0.6 1.2
1.03 1.04 1.09 1.10
1.02 1.04 1.07 1.08
Equilibrium
1.10 1.08
t, s
Cl Cl c 3 1.01 1.02 1.01 1.01
c 3
1.03 1.05 1.06
1.03 1.02 1.03 1.02 1.03 1.02
1.01
1.06
1.03 1.02 1.01
1.01 1.01
found in polymer problems in which the inference becomes increasingly more accurate as the number of moments increases.16 The origin of the apparently different behavior is not clear, but may reside in the use of both even and odd moments in the present problem, as opposed to only even moments in the polymer problem. In fact, the asymptotic convergence found with a Gaussian basis when both even
556
A. K. Thakur, A. Rescigno, and C. DeLisi
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978 1.2
-
1.0
-
.2 -
0.8 .15
G 0.6 0.4
-
0.2
-
,"
3.0
34
38
4.2
4.6
5.0
54
58
62
6.6
-
.l-
70
P .05
Flgure 5. PI and & phase-plane regions showing the positive definiteness (I) and unimodality (11) of distributions reconstructed by Gram-Charlier type A series method (reproduced from ref 20).
-
-
.w.a-
.02-
and odd moments are used has been noted previou~ly~'-~~ 30 O10 2 4 6 8 10 12 14 16 18 20 and is probably a property of the set (or subset) or orNO OF A PARTICLES thornormal polynomials used, rather than the physical Flgure 6. Distribution functions for the reaction A 4- B * C at t = problem. The efficiency of the inference has a direct 1 s with a = 0.1 and fl = 1. Initially there were 40 and 60 units in bearing on a number of biological problems including states A and B, respectively, and none in state C. pi = 0.153 and calculation of the free energy of RNA structures,15 Pa = 2.2 for this case. See Figure 3 for explanation of symbols. characterization of the affinity heterogeneity of various biologically active plasma membrane bound receptors,lg and models of immunological feedback in carcinogenesis.20 The criteria for accurate determination of the distribution from its moments, within the context of certain approximation procedures, has received considerable attention in the statistical literature. Barton and Dennis2I investigated the regions of positive definiteness and unimodality of Gram-Charlier type A series reconstruction (Appendix B) of distribution functions from their first four moments. Their results are shown in Figure 5 . The consequence of their finding is that before making conclusions about distributions obtained by this method the Pearson moment coefficients PI = pu32/var3 and P2 = p4/varz should be carefully examined. Unless their values lie within the prescribed ranges shown in the figure, any interpretation based on uni- or multimodality of such reconstructed distributions may be erroneous. In this respect the works of Fisherz2 and ShentonZ8are also noteworthy. They derived ranges of values of PI and P2 NO. OF A PARTICLES within which such reconstructions are 80% efficient. Figure 7. Distribution functions for the case as in Figure 6 at t = 0.5 Unfortunately there is no such general theory available at s: I, exact distribution; 11, distribution by Gram-Charlier method; and present for expansions based on higher order moments or 111, Gaussian approximation of the distribution. pi = 0.27, p2 = 1.4. on different basis sets. This problem is addressed in a different context el~ewhere.'~ The limitations are illustrated by application of the techniques to the reactions A B i=+C and A B + C D (Appendix C). For the first of these reactions near equilibrium, P1 and P2 deviate only slightly from the required range indicated in Figure 5 and the inference is reasonably accurate (Figure 6). Far from equilibrium, however, deviation is large and the inference fails (Figure 7). For the second reaction, the distribution cannot be ," .1 accurately approximated by a four-moment expansion, even at equilibrium as is evident in Figure 8. Thus, as one might expect, the accuracy of the technique is sensitive to the physical chemical nature of the reaction, through the parameters PI and &.
+
+
r
+
1
Appendix A Putting m = 1, 2, ..., 6 in eq 11 and remembering that the mth moment is defined as = Z
i= 0
impi
the following six linear ordinary differential equations can be written for the first six moments:
' 20
i
40 NUMBER OF A PARTICLES
ca
M,
,\-,
, 10
30
+
Figure 8. Equilibriumdistribution functions for the reaction A B C D with a = 0.1, p = 1, a = 40, b = 60, and c = d = 0. Explanation of symbols same as in Figure 6. p, = 18.19 and & = 44.9. I is obtained by solving eq D.2 numerically.
+
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978 557
Stochastic Theory of Second-Order Chemical Reactions
The Ai are determined from the orthonormality condition of Hj and are found to be
Where the Tchebycheff-Hermite polynomials are given by H~ (y = (- 1yeY I 2
0j
(e-Y2/2 )
(B.4)
The series is generally truncated at j = 4 to avoid problems involving negative frequencies. Explicitly the first few Hi and A j are given by
and
Ho(Y) = 1,Hl(Y) = Y HZ(Y) = Y - 1 y 3 - 3y, H4(y) = y 4 - 6yz + 3 9
9
1=
H3(Y
(B.5)
where p3 and p4 are the third and fourth moments around the mean (central moments). Substituting these into eq B.2 gives
1 P ( X ) = ~ e x p 42na A3H3
?+)+
Ad4
r+)]
(B.6)
Details of this method are given in ref 16-19.
Appendix C (i) T h e Reversible Reaction A + B + C. Following the same line of arguments presented in section 11, the forward Kolmogorov equation for the reversible second-order reaction A t B z C
P
can be written as dP,(t)/dt=a[(i+ l ) ( b - a
a
(c
if it is known that it is not far away from a Gaussian. P2,(X)is a polynomial of degree 2n. In this particular type of reconstruction procedure Pznis expressed as a sum of Tchebycheff-Hermite polynomials (Hj) and the distribution is written as P(X) = 1 exp
42nU
[-;
+ l)Pi+l(t)-i(b-
+ i)Pj(t)] + P [ ( c + a - i + l)P,-,(t)+ a - i)Pj(t)]
(C.1)
with initial condition Pi(0)= 6i,, where a and 0 are the forward and backward rate constants as before, Pi@)is the probability of finding i number of particles in state A at time t and a , b, and c are the numbers of particles present initially in states A, B, and C, respectively. (ii) T h e Reversible Reaction A B + C D. In a similar fashion the forward Kolmogorov equation for the reversible second-order reaction
+
+
A + B S C t D P
is
2iljHjlx j= 0
+
d P j ( t ) / d t = a [ ( i l ) ( b - a + i + l ) P j t , ( t )i(b - a i)Pj(t)] p [ ( c a - i l ) ( d a - i + 1 ) P i _ l ( t ) - (c a - l ) ( d a i)Pi(t)l
+
+
+
+
+ +
+
(C.2)
558
The Journal of Physical Chemistry, Vol. 82, No. 5, 1978
J. A. Amelse, J. 5.Butt, and L. H. Schwartz
with initial condition and other symbols same as before and a! is the number of particles present initially in state D.
(10) I. G. Darvey, B. W. Ninham, and P. J. Staff, J. Cbem. Pbys., 45, 2145 (1966). (11) I. G. Darvey and B. W. Ninham, J . Chem. Phys., 46, 1626 (1967). (12) A. Renyi, Magy. Tud. Akad. Alkalm. Mat. Int. Kozl., 2, 93 (1954). P. J. Staff, J. Chem. Phys., 46, 2209 (1967). C. J. Jachimowski, D. A. McQuarrie, and M. E. Russel, Blochemistry, 3, 1732 (1964). C. DeLisi, Biopo/ymers, 11, 2251 (1972). M. G.Kendaii, “The Advanced Theory of Statistics”, Voi. I, Griffin, London, 1943. H. Cramer, “Mathematical Methods of Statistics”, Princeton University Press, Princeton, N.J., 1946. W. P. Eiderton and N. L. Johnson, “Systems of Frequency Curves”, Cambridge University Press, London, 1969. A. K. Thakur and C. DeLisi, Biopo/ymers, in press. N. Dubin, “A Stochastic Model for Immunological Feedback in Carcinogenesis” in “Lecture Notes in Biomathematics”, Vol. 9, Springer-Veriag, Berlin, 1976, D. E. Barton and K. E. Dennis, Biometrika, 39, 425 (1952). R. A. Fisher, Phil. Trans. R. SOC.London, Ser. A, 222, 309 (1922). L. R. Shenton, Biometrlka, 38, 58 (1951).
References and Notes A. F. Barthoiomay, Bull. Math. Blopbys., 20, 175 (1958). A. T. Bharucha-Reid, “Elements of the Theory of Markov Processes and Their Applications”, McGraw-Hili, New York, N.Y., 1960. D. A. McQuarrie, “Methuen’s Monographs on Applied Probability and Statistics”, Voi. 8, Methuen, London, 1967. A. K. Thakur, A. Rescigno, and D. E. Schafer, Bull. Math. Biol., 34, 53 (1972). A. K. Thakur, A. Rescigno, and D. E. Schafer, Bull. Math. Biol., 35, 263 (1973). A. K. Thakur, Ph.D. Thesis, University of Minnesota, Minneapolis, Minn., 1975. D. A. McQuarrie, Adv. Cbem. Phys., XV, 149 (1969). D. A. McQuarrie, C. J. Jachimowski, and M. E. Russel, J . Cbem. Phys., 40, 2914 (1964). K. Ishida, J . Chem. Phys., 41, 2472 (1964).
Carburization of Supported Iron Synthesis Catalysts J. A. Amelse,+ J. 8. Butt,* and L. H. Schwartd Department of Chemlcal Englneerlng, and Department of Material Science and Engineerlng, North western Unlversity,Evanston, Illinois 6020 1 (Received October 17, 1977) Publication costs assisted by the Petroleum Research Fund
The chemical state of iron in a 4.94 wt % Fe/SiOz catalyst was examined by Mossbauer spectroscopy at various stages of calcination and reduction, and after use as a synthesis catalyst. About 90% of the iron in the initial oxide (a-Fez03)was reduced to a-Fe metal during 24 h reduction in Hzat 425 “C. Under reaction conditions the catalyst carburized within 90 min to the extent that no metallic iron could be detected in the Mossbauer spectra. The carbide formed is the 6’ form, as opposed to the x form previously reported for unpromoted fused or precipitated catalysts or the F form reported for KzO or Cu promoted catalysts. Reaction rates were measured in a differential reactor system. Methane formation rates increased monotonically as the catalyst carburized to a steady value of turnover frequency of 0.0037 f 0.0003 s-l. Product selectivity shifted toward a higher alkane content at the expense of methane formation as the catalyst carburized.
Introduction In the recently reawakened interest in the catalysis of Fischer-Tropsch ahd related synthesis reactions much work has centered on catalysts dispersed on high area supports such as silica or alumina rather than the fused or precipitated types originally employed. Primary topics of interest have been the determination of specific rates or turnover frequencies per unit surface area of the metal, the mechanism(s) of the synthesis reactions, and the role that electronic and alloying effects play in determining activity and ~electivity.l-~ Some caution must be exercised in the interpretation of recent attempts to deal with these topics, since much of the characterization which has been presented does not deal with the fact that Fischer-Tropsch catalysts apparently change their state under reaction conditions. It has been known for some time that the three principle catalysts, iron, cobalt, and nickel, all form carbides, at least in the fused or precipitated state.5 Controversy continues as to whether the synthesis reactions involve these carbides or whether some form of oxygenated species is the active *To whom correspondence should be addressed. Department of Chemical Engineering, Ipatieff Catalytic Laboratory. ‘Department of Chemical Engineering. %Department of Material Science and Engineering, Ipatieff Catalytic Laboratory, Materials Research Center. 0022-3654/78/2082-0558$01 .OO/O
intermediate.2i6-8 It is not our present purpose to engage in details of the synthesis reaction mechanism, rather we wish to investigate the possible formation of carbide phases under reaction conditions for a typical metal synthesis catalyst dispersed on a high area support. No information exists in the literature concerning the rate of formation or the form of carbide produced in supported catalysts. In the following we report the results of a detailed investigation of an iron on silica catalyst. Silica was chosen as the support material as it is generally considered more inert than alumina or silica/alumina. Mossbauer spectroscopy was used for characterization of the chemical state of the iron after various oxidation-reduction cycles and after use as a synthesis catalyst. Chemical reactivity studies were also carried out to provide information on synthesis rates and the rate of carburization as a function of time of catalyst use in addition to data on product selectivity.
Experimental Section Catalyst. The support material used was 80-120 mesh Davison Grade 62 silica gel, BET area 340 m2/g, average pore diameter 14 nm. The gel was repeatedly washed with redistilled water and 0.1 N HN09 to remove alkali metal impurities, then dried and impregnated with an Fe(N03)3 solution (Mallinckrodt Reagent Grade, Lot AVL) via the incipient wetness technique. About 1.1cm3of solution per 0 1978 American Chemical Society
