J. Phys. Chem. 1996, 100, 2597-2604
2597
Theory and Simulation of Stochastically-Gated Diffusion-Influenced Reactions Huan-Xiang Zhou* Department of Biochemistry, Hong Kong UniVersity of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
Attila Szabo* Laboratory of Chemical Physics, National Institute of Diabetes and DigestiVe and Kidney Diseases, Room 138, Building 5, National Institutes of Health, Bethesda, Maryland 20892-0520 ReceiVed: August 15, 1995X
The kinetics of the irreversible diffusion-influenced reaction between a protein (P) and a ligand (L) is studied when [L] . [P] and the reactivity is stochastically gated due to conformational fluctuations of one of the species. If gating is due to the ligand, we show that the Smoluchowski rate equation, d[P(t)]/dt ) -k(t)[L][P(t)], can be generalized by simply using a stochastically-gated time-dependent rate coefficient, ksg(t). However, if gating is due to the protein, this is no longer true, except when the gating dynamics is sufficiently fast or the ligand concentration is very low. The dynamics of all the ligands around a protein become correlated even when they diffuse independently. An approximate theory for the kinetics of proteingated reactions that is exact in both the fast and slow gating limits is developed. In order to test this theory, a Brownian dynamics simulation algorithm based on a path-integral formulation is introduced to calculate both ksg(t) and the time dependence of the protein concentration. Illustrative simulations using a simple model are carried out for a variety of gating rates. The results are in good agreement with the approximate theory.
I. Introduction Conformational fluctuations of a macromolecule or a ligand can play a significant role in determining the kinetics of ligand binding. A classical example is the entrance of oxygen into the heme pocket of myoglobin. In the static X-ray structure of myoglobin, there is no “hole” for the ligand to enter, and it was suggested that the side chains blocking the entrance “might act as a gate which has to swing out to admit the ligand”.1 Another example is provided by a ligand (e.g., a peptide) that can interconvert among various conformations with differing reactivity for the binding site of a protein. Perhaps the simplest way to handle such conformational fluctuations is to model them as a stochastic gate. The influence of gating on the kinetics of diffusion-influenced reactions is the focus of this paper. In the simplest case, the gate has two states: an open state representing a reactive conformation and a closed state representing a nonreactive conformation. McCammon and Northrup2 analyzed the influence of such a gate when the transitions between the states are deterministic. For stochastic gating, Szabo et al.3 showed how one can generalize the Smoluchowski calculation of the steady-state bimolecular rate constant. This work was based on the assumption that it does not matter whether it is the protein or the ligand that interconverts between the two conformations. While such symmetry indeed exists for the reaction between a single protein and a single ligand (i.e., the isolated pair problem), we will show that it breaks down for the situation where a single protein is surrounded by many ligands. The basic idea of the Smoluchowski approach to the kinetics of diffusion-influenced bimolecular reactions is to reduce the many-particle problem to a pair problem. Since the pair problem has the above symmetry, one naively expects that this symmetry carries over to the many-particle problem. In other words, the kinetics of the diffusive reaction between a X
Abstract published in AdVance ACS Abstracts, January 15, 1996.
This article not subject to U.S. Copyright.
gated protein with a large number of ligands appears at first sight to be the same as the kinetics of the diffusive reaction between a protein with gated ligands. We will show that this is true only if the gating dynamics is sufficiently fast or the ligand concentration is very low and, in general, the naive generalization of the Smoluchowski theory for the kinetics is only useful when the ligands are gated. When the protein is gated, the dynamics of the ligands become correlated even if it is assumed that they diffuse independently. This correlation arises because the switching of the protein from a reactive conformation to a nonreactive one (or vice versa) is “felt” simultaneously by all the ligands. For the kinetics of diffusion-influenced reactions in which the ligands are gated, we will show that the familiar Smoluchowski theory can be generalized by simply using a stochastically-gated version of the time-dependent rate coefficient. It is the long-time limit of this quantity that was obtained by Szabo et al.3 To derive a theory for the kinetics of diffusion-influenced reactions in which the protein is gated, we use the approach of Waite,4 which was generalized to reversible reactions by Lee and Karplus.5 In this approach, one first obtains an infinite hierarchy of equations for the reactive reduced distribution functions and then truncates this hierarchy at the pair level using a superposition approximation. In order to test the above theory, we performed many-particle Brownian dynamics simulations of the reaction between a gated protein and a collection of diffusing ligands. These simulations are based on a path-integral formulation of the many-particle survival probability, in which the independence of the diffusive motion of the ligands is fully exploited. II. Theory As noted in the Introduction, the effects of gating on the kinetics of the diffusion-influenced reaction between a macromolecule and a collection of ligands are quite different,
Published 1996 by the American Chemical Society
2598 J. Phys. Chem., Vol. 100, No. 7, 1996
Zhou and Szabo
depending on whether the macromolecule or the ligands are gated. To introduce notation and for future reference, first we treat the reaction of an isolated pair. Then, we treat the two cases of gating in the many-particle system. 1. Isolated Gated Pair. The diffusion-influenced reaction between a single macromolecule and a single ligand when either is stochastically gated is considered. κ(r,g) denotes the reactivity for a given value r of the intermolecular displacement vector and a given gating state g. The probability p(r,g,t) of finding the pair at a displacement r and in gating state g at time t satisfies
∂p(r,g,t)/∂t ) [Lr + Lg - κ(r,g)]p(r,g,t)
(2.1)
In this equation,
Lr ) ∇De-βU(r)‚∇eβU(r)
(2.2)
where D is the sum of the diffusion coefficients of the macromolecule and the ligand, U(r) is the interaction potential between the pair, and Lg is an operator describing the fluctuation of the gating states. If the gating states are continuous, then
Lg )
∂ ∂ 1 D p (g) ∂g g eq ∂g peq(g)
(2.3)
where Dg is an effective diffusion coefficient and peq(g) is the equilibrium distribution of the gating states. Continuous gating states have been previously considered by Lee and Karplus.6 If the gating states are discrete, then Lg is a matrix. For example, when the gate fluctuates between an open and a closed state, one has
Lg )
[
-a +a
+b -b
]
a
c (closed state) o (open state) a b
(2.5)
Often the reactivity κ(r,g) can be factored into the form h(r) f(g). For the two-state gating model represented by eq 2.5, f(g) ) 1 when g ) o and 0 when g ) c. If reaction occurs only at contact, then the reactivity is a δ function
k0
δ(r - R)f(g) 4πR2
(2.6)
where R is the contact distance. When combined with a reflecting boundary condition at contact, eq 2.6 is equivalent to using the Collins and Kimball radiation boundary condition.7,8 We will develop our theory for a general gating operator Lg and a general reactivity κ(r,g), but particular attention will be paid to the two-state gating model and the above δ-function reactivity. To formally generalize the Smoluchowski approach to the calculation of the time-dependent rate coefficient, one solves eq 2.1 subject to the initial condition p(r,g,0) ) exp[-βU(r)]peq(g) and then defines the stochastic-gated version of the rate coefficient as
[
]
∂p(r,g,t) ksg(t) ) ∫∫d r dg ∂t 3
1 1 a ) + b(s + a + b)k ˆ (s + a + b) skˆ sg(s) skˆ (s)
(2.8)
which reduces to the result of Szabo et al. when s f 0. The results for ksg(t) in the limits that the gating dynamics is very fast and very slow will find use later in the section. In the fast gating limit, the distribution of gating states immediately relaxes to the equilibrium one and p(r,g,t) can be approximated as peq(g) pj(r,t). Using this in eq 2.1 and integrating both sides over g, one finds
∂pj(r,t)/∂t ) [Lr - κj(r)]pj(r,t)
(2.9)
where κj(r) ) ∫dg κ(r,g) peq(g). The result for the rate coefficient is
kfast (t) ) ∫d3r κj(r) pj(r,t)
(2.10)
In the slow gating limit, the gate is frozen in its initial state and thus the distribution of gating states is kept at the initial equilibrium one. For a given gating state g, the probability of finding the pair at a displacement r at time t, p˜ (r,t|g) ) p(r,g,t)/ peq(g), satisfies eq 2.1 with Lg set to zero, i.e.,
∂p˜ (r,t|g)/∂t ) [Lr - κ(r,g)]p˜ (r,t|g)
(2.11)
The rate coefficient for a given gating state g is
(2.4)
where a and b are the rate constants in the kinetic scheme
κ(r,g) )
involving the two-state gating model and the δ-function reactivity in the absence of an interaction potential. One can readily extend this result to obtain the Laplace transform of ksg(t) (i.e., kˆ sg(s) ) ∫∞0 dt exp(-st) ksg(t)). In terms of the Laplace transform of the rate coefficient k(t) for the ungated problem, one has
) ∫∫d3r dg κ(r,g) p(r,g,t)
(2.7)
Szabo et al.3 obtained the t f ∞ limit of ksg(t) for the case
kslow(t|g) ) ∫d3r κ(r,g) p˜ (r,t|g)
(2.12)
For an equilibrium distribution of gating states, the rate coefficient is
kslow(t) ) ∫dg peq(g)kslow(t|g)
(2.13)
Note that both pj(r,t) and p˜ (r,t|g) have the same initial value exp[-βU(r)]. One can also use the survival probability S(t|r,g) to describe the pair problem. This is the probability that the pair, initially at a displacement r and in gating state g, does not bind by time t. From the fact that the Green function of eq 2.1 satisfies the detailed balance condition, one can easily show that p(r,g,t) and S(t|r,g) are related by
p(r,g,t) ) e-βU(r) peq(g) S(t|r,g)
(2.14)
The rate coefficient defined in eq 2.7 becomes
ksg(t) ) ∫∫d3r dg κ(r,g)e-βU(r) peq(g) S(t|r,g)
(2.15)
In particular, for the two-state gating model, one has
ksg(t) ) ∫d3r h(r) e-βU(r) peq(o) S(t|r,o)
(2.16)
where peq(o) ) b/(a + b) and S(t|r,o) is the survival probability of the pair starting at a displacement r and in the open gating state. 2. Many-Particle System: Gated Ligands. We now consider the kinetics of the diffusion-influenced reaction between the macromolecule at a small concentration and the ligand at concentration C when each ligand undergoes confor-
Stochastically-Gated Diffusion-Influenced Reactions
J. Phys. Chem., Vol. 100, No. 7, 1996 2599
mational fluctuations or is gated. When the concentration of the macromolecule is sufficiently small, we can focus on the fate of a single macromolecule surrounded by N ligands in a volume V provided both N and V are large and N/V ) C. We assume that the ligands do not interact with each other but do interact with the macromolecule via a potential of mean force U(r) and that the gating dynamics of each ligand is described by the operator Lg. Within a Born-Oppenheimer-like approximation,8,9 we can fix the macromolecule at the origin and let each ligand diffuse around it with the diffusion coefficient D. After these simplifications, the probability P(r1,g1,...,rN,gN,t) of finding the ith ligand at a displacement ri from the macromolecule and in gating state gi at time t satisfies
which has the familiar Smoluchowski form but with the stochastically-gated time-dependent rate coefficient. 3. Many-Particle System: Gated Macromolecule. We turn to the case where it is the macromolecule that undergoes conformational fluctuations or is gated. In this case, we need to deal with the probability P(r1,...,rN,g,t) of finding the macromolecule in gating state g and the ith ligand at a displacement ri from it at time t. This satisfies
∂
N
∂t
P(r1,...,rN,g,t) ) {∑[Lri - κ(ri,g)] + Lg}P(r1,...,rN,g,t) i)1
(2.24)
and has the initial value
∂
P(r1,g1,...,rN,gN,t) ) N
[Lr + Lg - κ(ri,gi)]P(r1,g1,...,rN,gN,t) ∑ i)1 i
i
(2.17)
We are interested in the fate of the macromolecule after surrounding it with an equilibrium distribution of ligands at t ) 0, i.e.,
e-βU(ri)
N
P(r1,g1,...rN,gN,0) ) ∏
e-βU(ri)
i)1
V
Again the quantity of interest is the survival probability of the macromolecule, given by
SN(t) ) ∫d3r1 ... ∫d3rN∫dg P(r1,...,rN,g,t)
i
i)1
e-βU(ri) peq(gi) V i)1 N
N
) ∏[p(ri,gi,0)/V]
(2.18)
i)1
Clearly eq 2.17 when subject to the initial condition of eq 2.18 is separable, and the solution is N
P(r1,g1,...,rN,gN,t) ) ∏[p(ri,gi,t)/V]
(2.19)
i)1
The survival probability of the macromolecule in the presence of N ligands at time t is
SN(t) ) ∫∫d r1 dg1 ... ∫∫d rN dgNP(r1,g1,...,rN,gN,t) ) [∫∫d3r dgp(r,g,t)/V]N
(2.20)
Differentiating with respect to t, one finds
]
dSN(t) ∂p(r,g,t) N SN-1(t) (2.21) ) - ∫∫d3r dg dt V ∂t Using the definition of the stochastically-gated time-dependent rate coefficient in eq 2.7 and noting that SN-1(t) ≈ SN(t) when N is large, we have
dSN(t)/dt ) -Cksg(t) SN(t)
(2.22)
The final result for the survival probability of the macromolecule is
SN(t) ) exp[-C∫0 dt′ ksg(t′)] t
(2.23)
e-βU(ri) V (2.26b)
where S(t|r1,...,rN,g) is the survival probability starting from a particular configuration (r1,...,rN,g). Even though the initial value of P(r1,...,rN,g,t) is factorable, the coupling between g and ri in the reactivity κ(ri,g) prevents eq 2.24 from being separable. This is directly related to the fact that each conformational change of the macromolecule influences the dynamics of all the ligands simultaneously. We have not been able to find the exact solution of eq 2.24, except in the limits that the gating dynamics is very fast and very slow. In the fast gating limit, one can approximate h (r1,...,rN,t). Using this in eq 2.24 and P(r1,...rN,g,t) as peq(g) P intergrating both sides over g, one finds that P h (r1,...rN,t) satisfies a separable equation involving the averaged reactivity κj(ri) ) ∫dg κ(ri,g) peq(g) for the ith ligand. In analogy to the derivation of eq 2.23, one can easily show that the survival probability of the macromolecule in this limit is
SN(t) ) exp[-C∫0 dt′ kfast(t′)] t
3
[
(2.26a) N
≈∏
3
(2.25)
) ∫d3r1 ... ∫d3rN∫dg S(t|r1,...,rN,g) peq(g)∏
peq(gi)
∫d3ri e-βU(r )
i)1
N
P(r1,...rN,g,0) ) peq(g)∏
∂t
(2.27)
Thus, when the gating dynamics is very fast, eqs 2.17 and 2.24 predict the same result for the survival probability of the macromolecule and it does not matter whether the macromolecule or the ligands are gated. In the slow gating limit, the operator Lg in eq 2.24 can be set to zero, leading again to a separable equation for each gating state g. The resulting survival probability of the macromolecule is
SN(t) ) ∫dg peq(g) exp[-C∫0 dt′ kslow(t′|g)] (2.28) t
When the exponential functions of eqs 2.23 and 2.28 are expanded in power series, agreement is found only to the first order in ligand concentration. Thus, except at low ligand concentrations and short times, the survival probability is different depending on whether the macromolecule or the ligands are gated if the gating dynamics is slow. As an illustration, consider the long-time behavior of the survival probability when
2600 J. Phys. Chem., Vol. 100, No. 7, 1996
Zhou and Szabo
the gating model is the two-state one represented by eq 2.5. The long-time limits of kslow(t|o) and kslow(t) ) peq(o) kslow(t|o) are in general nonzero, and thus, the long-time limit of SN(t) predicted by eq 2.23 is zero. On the other hand, as kslow(t|c) ) 0, eq 2.28 predicts the long-time limit of SN(t) to be peq(c). To develop a tractable theory for arbitrary gating dynamics, one needs to make approximations. We will employ the standard strategy of converting eq 2.24 into a hierarchy of equations for the reduced distribution functions and then truncating the hierarchy. The pair and triplet distribution functions are10
F (r,g,t) ) (2)
F (r,r′,g,t) ) (3)
1
N
∫d r1 ... ∫d rN P(r1,...,rN,g,t)∑δ(ri-r) 3
3
C 1
i)1
(2.29)
∫d r1 ... ∫d rN × 3
3
C2
N
N
P(r1,...,rN,g,t) ∑∑δ(ri-r) δ(rj-r) (2.30) i)1 j)1 j*i
∫d3r1
∫d3rN
Let Q(g,t) ) ... P(r1,...,rN,g,t), which gives the survival probability of the maromolecule through
SN(t) ) ∫dg Q(g,t)
(2.31)
One can derive an equation for Q(g,t) by integrating both sides of eq 2.24 over r1, ..., rN. After writing κ(ri,g) as ∫d r κ(r,g) δ(ri-r), one finds
∂Q(g,t)/∂t ) LgQ(g,t) - C∫d3r κ(r,g) F(2)(r,g,t)
(2.32)
N By multiplying both sides of eq 2.24 by (1/C)∑j)1 δ(rj - r) and then integrating over r1, ..., rN, one finds the following equation for F(2)(r,g,t):
∂F(2)(r,g,t)/∂t ) [Lr + Lg - κ(r,g)]F(2)(r,g,t) C∫d3r′ κ(r′,g) F(3)(r,r′,g,t) (2.33) Similarly, one can derive an equation for F(3) involving F(4) and so on. We want to truncate this hierarchy by constructing an approximate relation between F(3) and F(2). As usual, we assume that the triplet distribution function can be written as a product of pair distribution functions so that F(3)(r,r′,g,t) is proportional to F(2)(r,g,t) F(2)(r′,g,t). To find the appropriate proportionality factor, we note that
Q(g,t) ) ∫d3r1 ... ∫d3rN P(r1,...,rN,g,t) 1 3 (2) ∫d r F (r,g,t) V 1 ) 2 ∫d3r ∫d3r′ F(3)(r,r′,g,t) V (1 - 1/N) 1 ≈ 2∫d3r ∫d3r′ F(3)(r,r′,g,t) (2.34a) V )
Consequently
∫d3r ∫d3r′ F(3)(r,r′,g,t) Q(g,t) ) ∫d3r F(2)(r,g,t)∫d3r′ F(2)(r′,g,t)
(2.34b)
Thus, we make the superposition approximation
F(3)(r,r′,g,t) ) F(2)(r,g,t) F(2)(r′,g,t)/Q(g,t)
(2.35)
Use of eq 2.35 in eq 2.33 leads to
∂F(2)(r,g,t) ) [Lr + Lg - κ(r,g)]F(2)(r,g,t) ∂t C ∫d3r′ κ(r′,g) F(2)(r′,g,t) F(2)(r,g,t) (2.36) Q(g,t) As a consistency check, we note that eq 2.36 reduces to eq 2.32 when both sides are integrated over r. Defining W(r,g,t) ) F(2)(r,g,t)/Q(g,t) and rearranging eqs 2.32 and 2.36, one finds
∂Q(g,t)/∂t ) [Lg - C∫d3r (κ(r,g) W(r,g,t)]Q(g,t)]
(2.37)
∂W(r,g,t)/∂t ) [Lr - κ(r,g)]W(r,g,t) + HgW(r,g,t)
(2.38)
where
HgW(r,g,t) )
1 [L Q(g,t) W(r,g,t) - W(r,g,t)LgQ(g,t)] Q(g,t) g (2.39)
The nonlinear coupling between W(r,g,t) and Q(g,t) in the last term of eq 2.38 prevents the equation from being easily solvable. We thus have tried to approximate this term so that a simple solution is possible. To this end, we note that, in the fast gating limit, Q(g,t) ) peq(g) SN(t) and thus
HgW(r,g,t) )
1 L p (g) W(r,g,t) peq(g) g eq
(2.40)
Now in the slow gating limit, Lg can be set to zero anyway, so eq 2.40 is also correct in this limit. With this approximation for the last term of eq 2.38, one can readily check that the solution of this equation is just W(r,g,t) ) p(r,g,t)/peq(g) ) exp[-βU(r)]S(t|r,g), where S(t|r,g) is the survival probability of an isolated gated pair. Using this in eq 2.37, one finds
∂Q(g,t)/∂t ) [Lg - C∫d3r κ(r,g) e-βU(r) S(t|r,g)]Q(g,t) (2.41) This equation provides the simplest possible description of the kinetics in the case where the maromolecule is gated.11 As expected, it correctly predicts the survival probability SN(t) in both the fast and the slow gating limits. We now specialize to the two-state gating model, for which κ(r,o) ) h(r) and κ(r,c) ) 0. Using eq 2.4 for Lg and eq 2.16 to express the integral over r in the open gating state in terms of ksg(t), one has
[ ] [
-a - Cksg(t)/peq(o) d Q(o,t) ) dt Q(c,t) +a
+b -b
][ ]
Q(o,t) Q(c,t) (2.42)
The initial conditions are Q(o,0) ) peq(o) ) b/(a + b) and Q(c,0) ) peq(c) ) a/(a + b). The sum of Q(o,t) and Q(c,t) gives the survival probability of the macromolecule. In order to test the utility of the theory developed here for the kinetics of the diffusion-influenced reaction between a gated macromolecule and a collection of ligands, we need to compare its predictions with accurate results for the many-particle system.
Stochastically-Gated Diffusion-Influenced Reactions
J. Phys. Chem., Vol. 100, No. 7, 1996 2601
In the next section we discuss a Brownian dynamics simulation algorithm for obtaining such results. III. Simulation Algorithm for Calculating ksg(t) and SN(t) Both ksg(t) and SN(t) are given by survival probabilities starting from particular initial distributions (see eqs 2.15 and 2.26b). Now we introduce a Brownian dynamics simulation algorithm for calculating survival probabilities. This algorithm, based on the work of Lamm and Schulten,12 uses the path-integral formulation of Brownian motion and requires the generation of nonreactiVe trajectories for its implementation. We consider a system obeying a general reaction-diffusion equation of the form
∂P(x,t)/∂t ) [Lx - κ(x)]P(x,t)
(3.1)
where x represents the position, gating state, orientation, etc., of one or more particles in the system. Note that eqs 2.1, 2.17, and 2.24 all have this form. In terms of the Green function of eq 3.1, G(x,t|x0,0), the survival probabilty of the system starting at x0 is
S(t|x0) ) ∫dx G(x,t|x0,0)
(3.2)
By breaking the time interval [0,t] into many small segments and repeatedly using the Chapman-Kolmogorov equation,13 one finds M
S(t|x0) ) ∏∫dxm G(xm,tm|xm-1,tm-1)
function (1 before the termination and 0 afterward) if the termination occurs before t or just 1 if termination has not occurred by time t. The survival probability S(t|x0) is obtained by generating a sufficient number of trajectories from x0 and averaging the results from individual trajectories. The pathintegral procedure has certain advantages over the react-andterminate one, which will become clear when we discuss specific systems. It should be emphasized that the two procedures described above are completely general and should prove useful for calculating any quantity that can be expressed in terms of the survival probability of a system starting from a particular initial distribution. For example, within the Smoluchowski approach, a time-dependent rate coefficient can be defined in terms of the dynamics of a pair of particles as
k(t) ) ∫dx0 peq(x0) κ(x0) S(t|x0)
where peq(x) is the equilibrium distribution (i.e., Lxpeq(x) ) 0) that satisfies the normalization condition limr0f∞ ∫dx δ(r-r0) peq(x) ) 1, and x represents all relevant variables that describe the dynamics of the pair, e.g., their displacement r, gating states, and orientations. Dividing both sides of eq 3.6 by
k(0) ) ∫dx0 peq(x0) κ(x0) k(t)/k(0) ) ∫dx0 F(x0) S(t|x0)
(3.3)
G(xm,tm|xm-1,tm-1) ) e-∆t[κ(xm)+κ(xm-1)]/2 G0(xm,tm|xm-1,tm-1) (3.4) Substituting this into eq 3.3 and taking the ∆t f 0 limit, one finds
S(t|x0) ) 〈exp{-∫0 dt′ κ[x(t′)]}〉x0 t
(3.5)
where 〈...〉x0 represents an average over nonreactive trajectories started from x0. In practice, eq 3.5 can be implemented as follows. (1) A set of nonreactive trajectories is propagated starting from x0. (2) For each trajectory, κ[x(t′)], is calculated numerically and integrated up to time t using the trapezoidal rule, and the exponential of the negative of the integral is evaluated. (3) The result from each trajectory is added up, and the sum is divided by the number of trajectories. When the same time steps are used for all the components of x, instead of saving trajectories and processing them later, one can simply accumulate ∆t[κ(xm) + κ(xm-1)]/2 for each trajectory and use the sum for calculating the survival probability. The above path-integral procedure can be compared with the react-and-terminate procedure developed previously by Zhou.14,15 For the system described by eq 3.1, the react-andterminate procedure can be recast as follows. One starts a nonreactive trajectory from x0 and, after each step, compares exp{-∆t[κ(xm) + κ(xm-1)]/2} with a random number that is uniformly distributed between 0 and 1. If the random number is larger, the trajectory is terminated at this step; otherwise the trajectory is propagated. Between 0 and t, the survival probability calculated from each trajectory is either a step
(3.7a)
one obtains
) ∫dx0 F(x0)〈exp{-∫0 dt′ κ[x(t′)]}〉x0 t
m)1
where (xM,tM) ) (x,t) and t0 ) 0. If the increment ∆t ) tm tm-1 is small, G(xm,tm|xm-1,tm-1) is related to the Green function G0(xm,tm|xm-1,tm-1) of eq 3.1 with κ(x) ) 0 through12
(3.6)
(3.7b)
where
F(x0) ) peq(x0) κ(x0)/∫dx0 peq(x0) κ(x0)
(3.7c)
is a normalized distribution. Thus, the time-dependent rate coefficient can be obtained by choosing initial configurations of the pair from F(x0) and then calculating its survival probability by either the path-integral procedure or the react-and-terminate procedure. In most applications, κ(x) is localized in a small region of configurational space and thus one needs to initiate trajectories only from this reactive region. We note in passing that the react-and-terminate procedure shares some features with the algorithm of Allison et al.16 for calculating the long-time limit of k(t). We now consider in detail two systems on which we have carried out Brownian dynamics simulations in order to test our theory for the kinetics of the diffusion-influenced reaction between a gated macromolecule and a collection of ligands. The first, described by eq 2.1, consists of a gated spherical macromolecule and a ligand.17 The second, described by eq 2.24, consists of a gated spherical macromolecule and N noninteracting ligands. The pair system is simulated because its dynamics determines the time-dependent rate coefficient ksg(t), which appears in our theory. The gating model is the two-state one represented by eq 2.5. For this model, a trajectory of the gating dynamics consists of a series of intervals spent in the open state interrupted by intervals spent in the closed state. The intervals spent in the open state have the distribution exp(-ato); thus if one starts the trajectory from the open state, the time interval before switching to the closed state can be generated by
to ) -a-1 ln rand
(3.8a)
where rand is a random number that is uniformly distributed
2602 J. Phys. Chem., Vol. 100, No. 7, 1996
Zhou and Szabo
between 0 and 1. Similarly, the time interval spent in the closed state can be generated by
tc ) -b-1 ln rand
(3.8b)
The reactivity is κ(r,c) ) 0, and
κ(r,o) ) k0/4πR2, R < r < R1 ) 0,
r > R1
(3.9)
When R1 - R ≡ f 0, this reduces to the δ-function reactivity given in eq 2.6. The trajectory of a ligand’s diffusive motion can be generated by18
r ) r0 - β∇r0U(r0)D∆t + s(2D∆t)1/2
(3.10)
where ∆t is the time step and s is a vector of Gaussian random numbers. For the two systems simulated, U(r) was set to zero. While more sophisticated procedures are available,12 we treated each reflecting boundary by taking small time steps while near the boundary and putting a ligand back to its previous position if a step brought it across the boundary. 1. Pair System and Calculation of ksg(t). To calculate ksg(t) using the path-integral procedure, one propagates gating trajectories from the open state and ligand trajectories from initial positions distributed uniformly in the reactive region R < r < R1 (see eq 3.7c).19 In propagation of ligand trajectories, in general variable time steps are used, but numerical operations are simplified enormously if a fixed time step (say, ∆treac) is used in the reactive region. Then, for the jth gating trajectory one records a function fj(t) that is 1 in the open state and 0 in the closed state, and for the jth ligand trajectory one records the time for each step at which the ligand is inside the reactive region (τ(l) j is the time for the lth such step). If J gating trajectories and J ligand trajectories have been generated, then the time-dependent rate coefficient calculated from them is20 J
[
ksg(t) ) ksg(0)J-1∑exp j)1
k0∆treac 2
]
∑ fj(τj(l))
4πR l,τj(l)
