Brownian dynamics simulations of diffusion-influenced reactions

Chemistry B 1997 101 (25), 5026-5030. Abstract | Full Text HTML | PDF ... Bruce H. Robinson , Colin Mailer , Gary Drobny. Annual Review of Biophys...
1 downloads 0 Views 534KB Size
J . Phys. Chem. 1990, 94, 7 133-7 136

7133

Brownian Dynamics Simulations of Diffusion-Influenced Reactions. Inclusion of Intrinsic Reactivity and Gating Stuart A. Allison,* Department of Chemistry, Georgia State University, Atlanta, Georgia 30303

J . Andrew McCammon, and Jacqueline J. Sines Department of Chemistry, University of Houston, Houston, Texas 77204-5641 (Received: January 24, 1990; In Final Form: April 2, 1990)

The Brownian dynamics algorithm for calculating steady-state bimolecular rate constants of diffusion-influenced reactions is applied to a number of simple model systems that are not fully diffusion controlled because of finite intrinsic reactivity and/or gating of the active site. Finite reactivity is dealt with by using procedures originally developed by Lamm and Schulten ( J . Chem. Phys. 1983, 78,2713) and extended by Northrup et al. ( J . Chem. Phys. 1986,84, 2196). According to Northrup et al., rate constants are derived from survival probabilities averaged over many independent Brownian dynamics trajectories. For the systems studied in this work, simulated rate constants are in very good agreement with analytic values. This is true even for extremely low intrinsic reactivities, where the reactions are far from diffusion controlled.

Introduction A Brownian dynamics approach first developed by Northrup et a1.l has been very useful in determining bimolecular rate constants for diffusion-controlled reactions. The original approach has undergone several refinements,2d and it has been applied to a number of complex bimolecular systems including superoxide dismutase (SOD)-superoxide,S-8 cytochrome ?cytochrome c p e r o ~ i d a s ecarbonic ,~ anhydrase-bicarbonate,lo and triose phosphate isomerase (TIM)-D-glyceraldehyde phosphate." Quite realistic modeling is now possible which includes taking account of the detailed topography of the enzymes (derived from crystal structures) as well as complex features of electrostatic interactions. In some but not all of these systems (Fe- and MnSOD, carbonic anhydrase, and TIM), simulated rate constants are substantially larger than experimental values, implying that the reactions are not strictly diffusion controlled. A number of factors might be responsible for this. One possibility is that fluctuations in structure are necessary to expose active sites. Such "gating" of active-site accessibility could be occurring in the FeSOD, MnSOD, and TIM systems mentioned above, as well as in oxygen binding to hemJ ~ substrate binding to lysozyme,14 to oglobin or m y ~ g l o b i n ' ~and name only a few. Alternatively, an activation barrier of some kind might have to be overcome in order for a reaction to take place. ( I ) Northrup, S. H.; Allison, S. A.; McCammon, J. A. J . Chem. Phys. 1984. 80. 1517. (2) Ailison, S. A.; Northrup. S. H.; McCammon, J. A. J . Chem. Phys. 1985, 83, 2894. (3) Ganti, G.; McCammon, J. A.; Allison, S. A. J . Phys. Chem. 1985,89,

7x99 --_,.

(4) Northrup, S. H.; Curvin, M. S.; Allison, S. A.; McCammon, J. A. J . Chem. Phys. 1986, 84, 2196. (5) Sharp, K.; Fine, R.; Schulten, K.; Honig, B. H. J . Phys. Chem. 1987, 91. 3624. (6) McCammon, J. A,; Bacquet, R. J.; Allison, S. A.; Northrup, S. H. Faraday Discuss. Chem. Soc. 1987, 83, 213. (7) Allison, S. A.: Bacquet, R. J.: McCammon, J. A. Biopolymers 1988, 27, 25 1. (8) Sines, J.; Allison, S. A.; Wicrtbicki, A.; McCammon, .I.A. J . Phys. Chem. 1990, 94, 959. (9) Northrup, S. H.; Boles, J. 0.; Reynolds, J. C. L. J . Phys. Chem. 1987, 91. 5991.

(IO) Reynolds, J . C. L.;Cooke, K. F.; Northrup, S. H. J . Phys. Chem. 1990. 94. 985. (1'1) Madura, J. D.; McCammon, J . A. J . Phys. Chem. 1989, 93, 7285. (12) Perutz, M. F.; Mathews, F. S. J . Mol. Biol. 1966, 21, 199. ( 1 3 ) McCammon, J. A.; Harvey, S. C. I n Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, U.K., 1987. (14) McCammon, J. A.; Gelin, B. R.;Karplus, M.; Wolynes, P. G. Nature (London) 1976, 262, 325.

0022-3654/90/2094-7 133$02.50/0

Desolvation of the active-site or inner-sphere rearrangement of the kind seen in many electron-transfer processesI5 represent possible barriers. In any event, it will be necessary to generalize the simulation approach in order to account for these factors. Fluctuations in the active-site regions of enzymes have been studied directly by molecular dynamics simulation on the picosecond time ~ c a l e . ' ~ ,However, '~ because of the very long trajectory times of substrates diffusing in the neighborhood of the enzyme and the necessity of carrying out many trajectories in order to determine a rate constant, it is not generally feasible to monitor directly enzyme fluctuations in real time during a rate constant simulation. In the present work, Brownian dynamics simulation is used to study several simple model systems with two features that make them slower than purely diffusion controlled. One feature is simple stochastic gating to mimic the effect of fluctuations that alternately expose and cover the active The second feature is finite intrinsic reactivity of the active sites, so that not every encounter complex formed results in a reaction.20 Northrup et aL4 describe a Brownian dynamics rate constant algorithm with finite intrinsic reactivity, but to the best of our knowledge, this algorithm has not been applied previously to any system. It is worth noting that the reactivity is modulated by these features even though the enzyme and substrate topographies remain constant. This is important from a computational standpoint, since much longer dynamics time steps can be taken under these conditions. Several Brownian dynamics simulation studies have been reported in which a spatially varying rate constant is introduced to mimic Forster energy transfer2' and exchange In these studies however, the encounter complex, which can be viewed as a relative configuration or series of configurations the reactive species must pass through in order for a reaction to occur, is not defined. In the present work, intrinsic reactivity refers to the reactivity of the encounter complex and is defined in terms of boundary gating rates and/or an intrinsic rate constant, which enters into the solution of the diffusion (15) Newton, F. D.; Sutin, N . Annu. Rev. Phys. Chem. 1984, 35, 437. (16) Shen, J.; Subramaniam, S.; Wong, C. F.; McCammon, J. A. Biopolymers 1989, 28, 2085. (17) McCammon, J . A.; Northrup, S. H. Nature (London) 1981,293,316. (18) Szabo, A.; Shoup, D.; Northrup, S. H.; McCammon, J. A. J . Chem. Phys. 1982, 77, 4484. (19) Northrup, S. H.; Zarrin, F.; McCammon, J. A. J . Phys. Chem. 1982, 86, 2314. (20) Collins, F. C.; Kimball, G. E. J . Colloid Sci. 1949, 4 , 425. (21) Allison, S. A . J . Chem. Phys. 1986, 90, 4020. (22) Herbert, R. G.; Northrup, S. H. J . Mol. Liq. 1989, 41, 207.

0 1990 American Chemical Society

7134

Allison et al.

The Journal of Physical Chemistry, Vol. 94, No. 18, 1990

equation subject to reactive boundary conditions. The principal objective of this work is to show that the Brownian dynamics simulation technique can be used to determine bimolecular rate constants for model systems that are not strictly diffusion controlled because of finite intrinsic reactivity and/or gating of the active site. The algorithm is described in the next section with emphasis on new features such as the way in which diffusion near a reflecting surface is handled. The accuracy of the simulation technique is demonstrated next by comparing simulated rate constants with analytic values for several simple model systems. The final section summarizes the main findings and conclusions of this work.

Methods Rate constants are calculated by following the approach of Northrup et a1.l The diffusion-controlled rate constant k is given by k = kD(b) p ( b ) (1) where p ( b ) is the probability that reactant pairs, initially at interparticle separation b but otherwise randomly distributed, “react” rather than diffuse to infinite separation, and kD(b)is the diffusion-controlled bimolecular rate constant at which the reactive particles initially achieve separation b. In order for eq 1 to be strictly valid, the interparticle potential of mean force at separation r = b must be isotropic.23 If this is not the case, generalizations of eq 1 should be e m p l ~ y e d . ~For . ~ the cases considered in this work, where interparticle forces are centrosymmetric, eq 1 is exact. For centrosymmetric interparticle forces, k D ( b )is given byZ4

[

kD(b) = L - d r exp[V(r)/kB~/4.r2D(r)]’

(2)

k b is Boltzmann’s constant, Tis absolute temperature, V(r) is the interparticle potential of mean force, and D(r) is the relative translational diffusion coefficient. In the absence of hydrodynamic interaction, D(r) is simply the sum of the diffusion constants of the two interacting particles, which shall be denoted D. (D = Dl D2 where Di = kB7‘/6st)ai,t) is the solvent viscosity, and aiis the hydrodynamic radius of subunit i.) In the special case V(r) = 0, the Smoluchowski rate, k,(b), is obtained:

+

ks(b) = 4xDb

e k

(open)

(closed)

(5)

where k , and k , are rate constants for closing and opening. In this example, 9 can simply be chosen large enough that the time required to diffuse from separation 9 to R is long compared to the lifetimes of both open and closed states, which are proportional to 1 / k , and 1 / k c , respectively. Now p is calculated from the survival probability of a single trajectory, wi ( i denotes trajectory number), averaged over many (typically several thousand) independent Brownian dynamics trajectories (23) Zhou, H. X . J . Chem. Phys., in press. (24) Northrup, S. H.; Hynes. J . T. J . Chem. P h y s . 1979, 7 1 , 871.

r’(t

(6)

i

+ At)

= r(t)

+ (kB7‘)-lD f(t) + s

(7)

where r(t) is the position vector of “substrate” at time f in a body-fixed frame of reference with origin at the center of the “enzyme”, f(t) is the direct force on substrate at time t , and s is a vector of independent Gaussian random numbers of zero mean = 2DAt (j denotes x, y , or z component). with with variance (s:) As discussed elsewhere,“ eq 7 is strictly valid only in the absence of absorbing or reflecting boundaries. This limitation has been overcome in the past by taking short dynamics time steps near such boundaries. Northrup et aL4 have recently addressed this problem by adapting one-dimensional Brownian dynamics algorithms for reaction and reflection, developed by Lamm and Schulten,26to the three-dimensional problem of calculating bimolecular rate constants. Cf. also related work by Green.27 A similar approach is followed in this work although a different algorithm is used to generate displacements near a reflecting surface. Consider the special case of a reflecting infinite p z plane located at x = 0. Let d be the initial position of the particle (along the position x axis), x its position after a dynamics step of duration At, and p,,,(x,Af)d)dx the conditional probability that the particle initially at d has moved to the interval {x,x + dx] after time Ar. Then4-26 p,,dx,Atld) = pf(x,Atld) + egdpr(x,Atl-d) + O.Sge-g* erfc ((x + d - ~ ) / ( 4 D h t ) ’ / ~(8) ) where g = - - ( d ) / k B Ta, = gDAt, erfc denotes the complementary error function, and pris the corresponding probability distribution in the absence of a reflecting boundary: pf(x,Atld) = ( 4 ~ D A t ) - l /exp(-(x ~

-d

+ u ) ~ / ~ D A (9) ~)

As discussed by Northrup et al.,“ the third term on the right-hand side of eq 8 is small and can be ignored to a good approximation. In this work, the approximation is made

~,,n(x,Afld)= p&X,Atld) + cpf(x,Atl-d)

(4)

In order for eq 4 to be valid, q must be chosen sufficiently large that all trajectories reaching separation q would have an equal probability of ultimately reacting if they were continued. Because of this restriction, some care must be used in choosing q. A case in point is the gated reaction of two spheres.Is The two spheres react upon encounter (reaching separation R ) only if the gate is in the “open” position when an encounter occurs. The dynamics of the gate can be described by the simple kinetic scheme

1 - (Cwi)/N

where N is the number of trajectories in the simulation. In order to get an idea of the uncertainty in P and k , the simulation is usually broken down into four or five subsimulations from which standard deviations can be obtained. Each trajectory, in turn, consists of many individual dynamics steps. If hydrodynamic interaction is ignored, the trajectory propagation formula to generate unconstrained diplacements is given by the ErmakMcCammon algorithmZS

(3)

The direct result of Brownian dynamics simulation is a truncated recombination probability, @, that the reacting pair, initially at separation b, react before reaching some larger separation q. The quantity p ( b ) is related to @ by

p ( b ) = @ / ( I - ( 1 - @)kD(b)/kD(q))

P=

(10)

where c is a constant that can be determined by the normalization of prenon the interval x = 0 to infinity. In the limit of no direct forces, eq 10 is exact. It is straightforward to show that c

= [2 - erfc

((a-

d)(4DAt)-’/2)]/erfc ( ( a + d)(4DAt)-Il2) (1 1 )

Equation 7 is used to generate a displacement according to p,(x’,Atld). This step is accepted if x’> 0. Otherwise, the area, CP, under the pr(x,Atld) distribution from x’to 0 is first computed 9= [erfc ((x’+ a - d)(4DAt)-’/2) - erfc ( ( a - d)(4DAf)-’l2)]/2 (12) Then the particle is reflected back into the x > 0 domain by finding that x, distributed according to cpt(x,Atl-d), which yields the same area according to

(25) Ermak, D. L.; McCammon, J . A. J . Chem. Phys. 1978, 69, 1352. (26) Lamm, G.; Schulten, K. J . Chem. Phys. 1983, 78, 2713. (27) Green, N . J. €3. Mol. Phys. 1988, 65, 1399.

The Journal of Physical Chemistry. Vol. 94, No. 18, 1990 7135

Brownian Dynamics Simulations Solving for x x=

-d - a

+ (4DAt)’/2erfc-l

[erfc ( ( d

+ a)(4DAt)-’l2) - 2@/c] (14)

where erfc-’ denotes the inverse of the complementary error function. (Lamm and Schulten26describe a simple algorithm to compute i = erfc-’ (p) where z is an algebraic function of p for 0 C p C I . In the event 1 < p C 2, one sets erfc-’ (p) = -z.) In this work, the reflecting surfaces are not planar but spherical and the problem involves three dimensions and not one. Since At is small during dynamics steps near reflecting boundaries, the reflecting surface can be assumed planar to a good appr~ximation.~ Also in this work, the spherical reflecting surfaces are centered at the origin of the body-fixed coordinate system, which makes the connection with three-dimensional diffusion very simple. At the beginning of a dynamics step, d = Ir(t)l - R where R is the radius of the reflecting surface. If r’ = Ir’(t + Ar)1 > R, the unconstrained step is accepted as the new position, r(r + Af). Otherwise, x is computed according to eq 14 and we set r(t

+ At) = r’(t + At) + (x + R - r’) n(t)

(15)

where n(t) = r(t)/lr(t)l is the initial outward normal from the reflecting surface. The survival probability of a single trajectory, wi, can be expressed as the product of the survival probabilities for every dynamics step of trajectory i

wj = nwjk k

(16)

where k denotes a particular dynamics step. For the sake of brevity, the ik subscript shall be omitted. As discussed by Lamm and Schulten26and Northrup et al.? one may generate dynamics steps using the reflecting algorithm discussed above but assign a statistical survival probability for each step given by w = Preac/Pren

(17)

In this work, eqs IO and 1 1 are used to calculate plefl. Appropriate equations for prcacare summarized in the Appendix. In the case of a partially reactive boundary, prcacis the one-dimensional Green’s function solution of the diffusion equation subject to a radiative boundary condition at the reactive surface

W a / a x + g)Preac = kbreac

(18)

where k’is a measure of boundary reactivity. In the limit k’m every encounter results in a reaction ( w can go to zero for steps taken near the boundary) whereas the limiting case k’= 0 corresponds to no reaction. Equation 18 can be generalized to “gate” the boundary by replacing k’on the right-hand side with k’h(t) where h(t) is a characteristic function equal to 1 or 0 depending upon whether the gate is open or closed. If the gate is open during the dynamics step, w is calculated in the same way as before, but if the gate is closed, then w = I . In this work, we shall restrict ourselves to stochastic gating. Suppose that the gate was open/closed during the previous dynamics step and that a dynamics step of duration At is now to be taken. The probability, p g , that the gate remains open/closed is

pe = (k,

+ k,

pg = (kc +

kO

exp(-(ko

+ kc)At))/(ko + k,)

(19a)

exp(-(kO + kC)Ar))/(kO + kc)

where 19a/ 19b correspond to open/closed initial states. Note that k, and k, refer to gate opening and closing rates following the kinetic scheme of eq 5. A uniform random number @ is generated on the interval (0,l). If p g exceeds @, the state of the gate remains unchanged; otherwise, it is switched. These are several factors in the simulation algorithm that are very important in connection with computational efficiency. Let d denote the distance of the substrate from the nearest reactive boundary. Since survival probabilities, w, only deviate from 1 when d is close to 0, it is unnecessary to go through the lengthy calculation of w much of time. At the start of the simulation, a cutoff

TABLE I: Comparison of Simulated and Analytic Rates for the Model of Szabo et aLI8 with k. >> k. Z krd(theory)O krd(sim)* Z kd(theory)“ k,d(sim)b 0.01 0.03 0.1 0.3

0.0102 (0.0004) 0.0300 (0.0010) 0.0935 (0.0028) 0.236 (0.006)

0.0099 0.0291 0.0909 0.231

Llk,hsory = Z/(l

+ Z).

1 3 10 m

0.500 0.750 0.909 1.000

0.507 (0.008) 0.759 (0.008) 0.921 (0.01I ) 1.020 (0.014)

bStandard deviations are in parentheses.

distance, d*, is determined such that (1 - w ( d * ) ) / (1 - w ( 0 ) ) I Then during any dynamics trajectory w is only calculated for d < d*; otherwise, it is set equal to 1. When several rates are being computed in a simulation, a separate d* is computed for each case. Another important factor is the dynamics time step, A!. Even in the first simulation studies,’ it was recognized that a variable At can drastically influence computational efficiency. Basically, short time steps should be taken where there are absorbing/reflecting boundaries nearby or there are steep gradients in interparticle force fields, but substantially longer time steps can be taken otherwise. Longer time steps can be taken when the survival probability approach is used4 than in the simpler “hit or miss” approach used initially.’ In this work, we set At =

+

r

r(d + R)2/D q2/(6000D) lo

d < d, d , < d < d2 d, < d < 4 - R

(20)

+

where d , = d* 2(2Dt0)1/2and d2 = q / ( 1 3 ( 6 ~ ) l / -~ )R. For a range of different rate calculations, d , is computed by using the largest d*. The adjustable parameters to and T determine At. Results and Discussion Consider first the case of two identical uniformly reactive spheres of radius a without direct forces or hydrodynamic interaction. The reaction is stochastically gated with first-order opening and closing rate constants k, and kc, respectively (eq 5). The reaction radius is taken to be at R = 2a, where the radiative boundary condition (eq 18), with k’replaced by k’h(t)) is assumed to hold. This special case was solved analytically by Szabo et al.’* and will serve as one of the benchmarks for simulations. To simplify the notation, let kd = k/ks(R) where k is the bimolecular rate constant and ks is the Smoluchowski rate (eq 3). Also let k’ = ks(R) 2/4aR2 where now the intrinsic reactivity is contained in the dimensionless quantity 2 . In this notation kred = Z(koz/kc)/[(z

+ 2 ) + ( 1 + z)(koz/kc)l

(21)

where z= 1

+ ((k, + kc)R2/D)’I2

(22)

In the limit ko >> k,, which corresponds to the case of an always open gate, krd = Z/( 1 E), which is a result derived by Collins and Kimbal120 and Noyes.** A comparison of simulated and analytical rates for the limiting case k, >> k, is presented in Table I. In this example, T = 293 K, 7 = 1 cP, a = 1 A, R = 2 A, 6 = 5 A, q = 8 A, T = 0.001, and to = 0.001 ps. The CPU time for the 8000 trajectory simulation was 26 minutes on an Amdahl 5880 mainframe computer. The simulations are able to accurately reproduce the theoretical reduced rates, verifying the soundness of the basic algorithm. Agreement can be slightly improved by decreasing T and/or to, but this increases CPU time proportionately. A comparison of simulated and theoretical rates for a more general example is given in Table 11. In this case, k, = 0.5D/R2, k, = 1.5D/R2, q = 25 A, T = 0.003, to = 0.003 ps, and other parameters are the same as above. Agreement between theory and simulations is excellent. These results confirm the validity of the simulation method with regard to gating and finite reactivity in the absence of direct forces. A simple model system with direct forces consists of two uniformly reactive spheres interacting via a centrosymmetric potential of mean force, U ( r ) ,with intrinsic reactivity k ’ = ks(R)

+

(28) Noyes, R. M . Prog. React. Kinet. 1961, 1, 129.

7136

The Journal of Physical Chemistrv, Vol. 94, No. 18, 1990

TABLE 11: Comparison of Simulated and Analytic Rates for the Model of Szabo et al." with k . = 0.5 and k . = 1.5 (in units of D / R * ) Z krd(theory)O k,,(sim) Z krd(theory)" krd(sim) 0,160 0.159 (0.007) 0.01 0.0025 0.00244 (0.00012) I 0.280 0.279 (0.012) 0.03 0.0074 0.00725 (0.00036) 3 IO 0.378 0.379 (0.015) 0.1 0.0237 0.0234 (0.0011) 0.446 (0.016) 0.0635 (0.0029) m 0446 0.3 0.0642

From eqs 21 and 22 in text

TABLE 111: Comparison of Simulated and Analytic Rates for Two Charged Spheres (+e. -e)

10-8

I 0-7 10-6 I 0-5 I 0-4 0.001 0.003 0.0 I 0.03 0. I 0.3 1

3 10 m I,

3.86 x 3.86 X 3.86 x 3.86 x 3.85 X 3.82 X 0.112 0.350 0.884 I .90 2.83 3.42 3.63 3.7 1 3.75

10-7

IO" 10-5 10-4

IO-'

4.04 (0.15) X 4.05 (0.15) X 4.06 (0.17) X 4.06 (0.17) X 4.06 (0.16) X 4.02 (0.15) X 0.1 I8 (0.005) 0.367 (0.010) 0.920 (0.035) 1.95 (0.03) 2.87 (0.04) 3.44 (0.04) 3.65 (0.04) 3.73 (0.04) 3.76 (0.04)

Allison et al. k'(0r E), or Some combination of the two factors is a matter of convenience. Summary

The results summarized in Tables 1-111 show that the Brownian dynamics simulation approach can yield accurate rate constants for models where the active-site reactivity is modulated by finite intrinsic reactivity and/or by gating. Furthermore, the algorithm works well even for extremely low intrinsic reactivities. This accuracy is made possible by using the survival probability app r ~ a c hinstead ~ , ~ ~ of the simpler "hit or miss" strategy employed in early studies. In future work, we plan to apply this methodology to complex enzyme substrate systems such as the reaction of superoxide with iron and manganese superoxide dismutases.8 Acknowledgment. S.A.A. acknowledges the NSF for support through a Presidential Young Investigator Award. Work at the University of Houston was supported by grants from the National Institutes of Health, the Robert A. Welch Foundation, and the Texas Advanced Research Program. J.A.M. is the recipient of the 1987 George H . Hitchings Award from the Burroughs Wellcome Fund.

IOd IO-' IO4

I OT2

Appendix I n case of finite k".25 pr,,,(k',d + Ax,Atld) = ~ r ( d+ Ax,At(d) + egd p f ( d + Ax,Atl-d) + p2(kr,d + Ax,At(d) ( A l ) where pz(k',d

From eqs 23 and 24.

Z/4nR2 as in the previous example. In the ungated case (the limit k, >> kc), it can be shown that kred = x / ( l

+ x exp(U(R)/kBr)Z)

(23)

where

Consider two identical spheres of radius a containing charges q1 and q2 at their respective centers and diffusing in a solution of dielectric constant t (assumed here to equal 78, corresponding to water) with no added salt. In this case, the potential of mean force is given by Coulomb's law, U(r) = q1q2/tr. Under these conditions, eq 23 reduces to the Debye result29in the limit of Z approaching infinity. Table 111 compares simulated rates with those predicted from eqs 23 and 24 for the case a = 1 A, R = 2 A, b = 5 A, q = I O A, 7 = 0.003, to = 0.003 ps, q1= +e, and q2 = -e. The CPU time for 5000 trajectories was 58 min. It should be emphasized that the rates for the 16 different L: values in the table were all derived from the same simulation. During each trajectory, recombination probabilities are computed for each Z in the manner described under Methods. The agreement between analytic and simulated rate constants, which is as good as the no-force cases considered previously, shows that the approximate reflective and reactive probability distributions yield accurate rate constants. It is also evident from Table 111 that the algorithm works well even when the intrinsic reactivity is very low. For L: less than 0.001, the asymptotic form of p,(k? given by eq A3 must be used because of approximations in the algorithms used to compute error functions. It is not always possible to distinguish between ungated models with finite Z (type I) and gated models with infinite L: (type 11). A case in point is the stochastically gated force-free model described by Szabo et a1.,I8 where type I and I 1 cases that satisfy the condition 2; =

k,z/kc

(25)

have the same rates and are indistinguishable. I n such cases, whether one modulates active-site reactivity by gating, varying (29) Debye. P. Trans. Electrochem. Soc. 1942, 82. 265

+ Ax,Atld) = 0.5(g - 2kr/D) exp[-g(Ax + d) +

(k'/D)(Ax + 2d + (k'/D - g)DAt)] X erfc [ ( A x 2d (2kr/D - g)DAt)/(4DAt)1/2] (A2)

+

+

and pf is given by eq 9 in the text. Simple algebraic expressions can be used to evaluate e r f ~ . In~ the ~ limit of low reactivity (Rk'/D